text:
using Test

@class Parent begin
    pfield::Int
    pfun(self) = self.pfield
    pfunAdd(self,x) = self.pfield + x
end

@class Parent2 begin
    pfield2::Int
    pfun(self) = self.pfield2
    pfunAdd(self,x) = self.pfield2 + x
end

@class Child(Parent, Parent2) begin
    cfield::Int
    pfun(self) = self.cfield
    pfunAdd(self,x) = self.cfield + x
    cfunSuper(self) = @super Parent pfun(self)
    cfunAddSuper(self, x) = @super Parent pfunAdd(self, x)
    cfunSuper2(self) = @super Parent2 pfun(self)
    cfunAddSuper2(self, x) = @super Parent2 pfunAdd(self, x)
end

p = Parent(0)
c = Child(0,1,2)
@test c.pfield == 0
@test c.pfield2 == 1
@test c.cfield == 2
pvalue = p.pfield
pvalue2 = c.pfield2
cvalue = c.cfield
@test pfun(p) == pvalue
@test pfun(c) == cvalue
@test pfunAdd(p,1) == pvalue + 1
@test pfunAdd(c,1) == cvalue + 1
@test_throws(MethodError, pfunAdd(c,"a"))

@test cfunSuper(c) == pvalue
@test cfunAddSuper(c,1) == pvalue+1

@test cfunSuper2(c) == pvalue2
@test cfunAddSuper2(c,1) == pvalue2+1


text:
module Coriolis

export
    FPlane, BetaPlane,
    x_f_cross_U, y_f_cross_U, z_f_cross_U

# Physical constants for FPlane and BetaPlane constructors.
const Ω_Earth = 7.292115e-5 # [s⁻¹] https://en.wikipedia.org/wiki/Earth%27s_rotation#Angular_speed
const R_Earth = 6371.0e3    # Mean radius of the Earth [m] https://en.wikipedia.org/wiki/Earth

"""
    AbstractRotation

Abstract supertype for parameters related to background rotation rates.
"""
abstract type AbstractRotation end

include("no_rotation.jl")
include("f_plane.jl")
include("beta_plane.jl")

end

ext: jl


text:
module ImplicitPlots

export implicit_plot, implicit_plot!, implicit_curve!, ImplicitFunction

import Contour
using RecipesBase
import MultivariatePolynomials, StaticPolynomials
import StaticArrays: SVector
const MP = MultivariatePolynomials
const SP = StaticPolynomials

using Requires: @require
function __init__()
    @require HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327" include("hc.jl")
end

struct ImplicitFunction{N,F}
    f::F
end
ImplicitFunction{N}(f) where {N} = ImplicitFunction{N,typeof(f)}(f)
(IF::ImplicitFunction{2})(x, y) = IF.f(x, y)
(IF::ImplicitFunction{3})(x, y, z) = IF.f(x, y)

ImplicitFunction(f::MP.AbstractPolynomialLike) = ImplicitFunction(SP.Polynomial(f))
function ImplicitFunction(f::SP.Polynomial)
    if SP.nvariables(f) == 2
        ImplicitFunction{2}((x, y) -> f(SVector(x, y)))
    elseif SP.nvariables(f) == 3
        ImplicitFunction{3}((x, y, z) -> f(SVector(x, y, z)))
    else
        throw(ArgumentError("Given polynomial doesn't have 2 or 3 variables."))
    end
end
function ImplicitFunction(f)
    nargs = 0
    try
        f(1.0, 1.0)
        nargs = 2
    catch e
        try
            f(1.0, 1.0, 1.0)
            nargs = 3
        catch e
            throw(ArgumentError("Provided function does not accept 2 or 3 arguments."))
        end
    end
    ImplicitFunction{nargs}(f)
end

@recipe function implicit(f::ImplicitFunction{2}; aspect_ratio = :equal, resolution = 400)
    xlims --> (-5.0, 5.0)
    xlims = plotattributes[:xlims]
    ylims --> xlims
    ylims = plotattributes[:ylims]

    linewidth --> 1
    grid --> true
    aspect_ratio := aspect_ratio

    rx = range(xlims[1]; stop = xlims[2], length = resolution)
    ry = range(ylims[1]; stop = ylims[2], length = resolution)
    z = [f(x, y) for x in rx, y in ry]

    nplot = plotattributes[:plot_object].n
    lvl = Contour.contour(collect(rx), collect(ry), z, 0.0)
    lines = Contour.lines(lvl)
    !isempty(lines) || return p

    clr = get(plotattributes, :linecolor, :dodgerblue)
    for (k, line) in enumerate(lines)
        xs, ys = Contour.coordinates(line)
        @series begin
            seriestype := :path
            linecolor --> clr
            if k > 1
                label := ""
            end
            xs, ys
        end
    end
end
implicit_plot(f; kwargs...) = RecipesBase.plot(ImplicitFunction(f); kwargs...)
implicit_plot!(f; kwargs...) = RecipesBase.plot!(ImplicitFunction(f); kwargs...)
implicit_plot!(p::RecipesBase.AbstractPlot, f; kwargs...) =
    RecipesBase.plot!(p, ImplicitFunction(f); kwargs...)

end # module

max_issues_repo_name: saschatimme/ImplicitPlots

max_forks_repo_name: saschatimme/ImplicitPlots


text:
module TestIVFADC

using Test
using IVFADC
using Distances

include("index.jl")
include("utils.jl")
include("search.jl")
include("persistency.jl")

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "SCIP_PaPILO"

version = v"0.1.0"

sources = [
    ArchiveSource("https://scipopt.org/download/release/scipoptsuite-8.0.0.tgz", "9b85283db0ac939b2d8eb3475067c8e1164b239e0c78e68f55dcc55859b78b2d"),
]

# Bash recipe for building across all platforms
script = raw"""
cd scipoptsuite*
mkdir build
cd build/
cmake -DCMAKE_INSTALL_PREFIX=$prefix\
  -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN}\
  -DCMAKE_BUILD_TYPE=Release\
  -DZIMPL=OFF\
  -DUG=0\
  -DAMPL=0\
  -DSYM=bliss\
  -DIPOPT_DIR=${prefix} -DIPOPT_LIBRARIES=${libdir} ..
make -j${nproc} scip
make -j${nproc} gcg
make papilo-executable

make install
cp bin/papilo "${bindir}/papilo${exeext}"

mkdir -p ${prefix}/share/licenses/SCIP_PaPILO
for dir in papilo scip soplex; do
    cp $WORKSPACE/srcdir/scipoptsuite*/${dir}/COPYING ${prefix}/share/licenses/SCIP_PaPILO/LICENSE_${dir}
done
cp $WORKSPACE/srcdir/scipoptsuite*/gcg/LICENSE ${prefix}/share/licenses/SCIP_PaPILO/LICENSE_gcg
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_gfortran_versions(expand_cxxstring_abis(supported_platforms(; experimental=true)))

filter!(platforms) do p
    arch(p) ∉ ("armv6l", "armv7l") && !Sys.iswindows(p) && libgfortran_version(p) >= v"4" && libc(p) != "musl"
end

# The products that we will ensure are always built
products = [
    ExecutableProduct("papilo", :papilo),
    ExecutableProduct("scip", :scip),
    LibraryProduct("libscip", :libscip),
    LibraryProduct("libgcg", :libgcg),
]

dependencies = [
    Dependency(PackageSpec(name="bliss_jll", uuid="508c9074-7a14-5c94-9582-3d4bc1871065")),
    Dependency(PackageSpec(name="boost_jll", uuid="28df3c45-c428-5900-9ff8-a3135698ca75"); compat="=1.76.0"),
    Dependency(PackageSpec(name="Bzip2_jll", uuid="6e34b625-4abd-537c-b88f-471c36dfa7a0"); compat="1.0.8"),
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
    Dependency(PackageSpec(name="GMP_jll", uuid="781609d7-10c4-51f6-84f2-b8444358ff6d"), v"6.2.0"),
    Dependency(PackageSpec(name="Ipopt_jll", uuid="9cc047cb-c261-5740-88fc-0cf96f7bdcc7")),
    Dependency(PackageSpec(name="oneTBB_jll", uuid="1317d2d5-d96f-522e-a858-c73665f53c3e"); compat="2021.4.1"),
    Dependency(PackageSpec(name="Readline_jll", uuid="05236dd9-4125-5232-aa7c-9ec0c9b2c25a")),
    Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(
    ARGS,
    name,
    version,
    sources,
    script,
    platforms,
    products,
    dependencies;
    preferred_gcc_version=v"7",
    julia_compat="1.6",
)

max_stars_repo_name: cometscome/Yggdrasil

max_issues_repo_name: cometscome/Yggdrasil

max_forks_repo_name: cometscome/Yggdrasil


text:
export NormalOp, normalOperator

struct NormalOp{S,D}
  parent::S
  weights::D
end

function Base.copy(S::NormalOp)
  return NormalOp(copy(S.parent), S.weights)
end

function normalOperator(S, W=I)
  return NormalOp(S,W)
end

# Generic fallback -> TODO avoid allocations
function Base.:*(N::NormalOp, x::AbstractVector)

  #@info size(x)  size(N.parent)   size(N.weights)

  d = adjoint(N.parent)*(N.weights*(N.parent*x))
  return d
end


text:
# Base model struct
Base.@kwdef mutable struct graded{T <: AbstractFloat} <: IRTmodel
    a::AbstractVector{T} = [0.0]
    d::AbstractVector{T} = [0.0]
    group::Vector{Int64} = [1]
    fixed::Bool = false
    # graded(;a = [0.0], d = [0.0], group = [1], fixed = true) = new(a, d, group, fixed)
end

function _distribute!(new, old::graded)
    old.a = new[begin:length(old.a)]
    old.d = new[length(old.a)+1:end]
end

function _distribute!(new, old::graded, q)
    old.a = new[begin:length(old.a)][q]
    old.d = new[length(old.a)+1:end]
end

function _distribute(new, old::graded)
    out = copy(old)
    out.a = new[begin:length(out.a)]
    out.d = new[length(out.a)+1:end]
    return out
end

function _distribute(new, old::graded, q)
    out = copy(old)
    out.a = new[begin:length(out.a)][q]
    out.d = new[length(out.a)+1:end]
    return out
end

function _take(par::graded)
    return [par.a; par.d]
end

function _take(par::graded, q)
    return [par.a[q]; par.d]
end

function _copy(x::graded)
    graded(a = x.a, d = x.d, group = x.group, fixed = x.fixed)
end

function _rescale!(par::graded, μ, Σ, q)
    if length(q) == 1
        # This item follows unidimensional model.
        A = 1/ sqrt(Σ[1]); K = -A*μ[1]
        par.a = par.a / A
        par.d = @. par.d - par.a*K/A
    else
        # This item seems to follow the multidimensional model.
        L = cholesky(Hermitian(Σ)).L
        tmp = par.a' * L
        par.a =  tmp'# Dimension mismatch ?
        par.d = par.d .+ par.a' * μ
    end
end

function _irf(p::graded, θ, u, q)::Float64
    if u == 0
        return 1.0 - logistic(p.a[q]'θ[q] + p.d[u+1])
    elseif u == length(p.d)
        return logistic(p.a[q]'θ[q] + p.d[u])
    else
        return logistic(p.a[q]'θ[q] + p.d[u]) - logistic(p.a[q]'θ[q] + p.d[u+1])
    end
end

function _irf(p, m::graded, θ, u, q)::eltype(p)
    a = p[1:length(m.a[q])]
    d = p[length(m.a[q])+1:end]
    if u == 0
        return 1.0 - logistic(a'θ[q] + d[u+1])
    elseif u == length(d)
        return logistic(a'θ[q] + d[u])
    else
        return logistic(a'θ[q] + d[u]) - logistic(a'θ[q] + d[u+1])
    end
end


text:
const STDLIB_UUID_STRS = Dict(
    "Pkg"              => "44cfe95a-1eb2-52ea-b672-e2afdf69b78f",
    "Statistics"       => "10745b16-79ce-11e8-11f9-7d13ad32a3b2",
    "Test"             => "8dfed614-e22c-5e08-85e1-65c5234f0b40",
    "CRC32c"           => "8bf52ea8-c179-5cab-976a-9e18b702a9bc",
    "Random"           => "9a3f8284-a2c9-5f02-9a11-845980a1fd5c",
    "SuiteSparse"      => "4607b0f0-06f3-5cda-b6b1-a6196a1729e9",
    "Libdl"            => "8f399da3-3557-5675-b5ff-fb832c97cbdb",
    "UUIDs"            => "cf7118a7-6976-5b1a-9a39-7adc72f591a4",
    "Distributed"      => "8ba89e20-285c-5b6f-9357-94700520ee1b",
    "Serialization"    => "9e88b42a-f829-5b0c-bbe9-9e923198166b",
    "SHA"              => "ea8e919c-243c-51af-8825-aaa63cd721ce",
    "DelimitedFiles"   => "8bb1440f-4735-579b-a4ab-409b98df4dab",
    "LinearAlgebra"    => "37e2e46d-f89d-539d-b4ee-838fcccc9c8e",
    "REPL"             => "3fa0cd96-eef1-5676-8a61-b3b8758bbffb",
    "FileWatching"     => "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee",
    "SharedArrays"     => "1a1011a3-84de-559e-8e89-a11a2f7dc383",
    "LibGit2"          => "76f85450-5226-5b5a-8eaa-529ad045b433",
    "Base64"           => "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f",
    "SparseArrays"     => "2f01184e-e22b-5df5-ae63-d93ebab69eaf",
    "Mmap"             => "a63ad114-7e13-5084-954f-fe012c677804",
    "Profile"          => "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79",
    "Unicode"          => "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5",
    "Dates"            => "ade2ca70-3891-5945-98fb-dc099432e06a",
    "InteractiveUtils" => "b77e0a4c-d291-57a0-90e8-8db25a27a240",
    "Future"           => "9fa8497b-333b-5362-9e8d-4d0656e87820",
    "Sockets"          => "6462fe0b-24de-5631-8697-dd941f90decc",
    "Logging"          => "56ddb016-857b-54e1-b83d-db4d58db5568",
    "Printf"           => "de0858da-6303-5e67-8744-51eddeeeb8d7",
    "Markdown"         => "d6f4376e-aef5-505a-96c1-9c027394607a"
)

const STDLIBS = collect(keys(STDLIB_UUID_STRS))

function uses(repo::String, tree::String, lib::String)
    pattern = string(raw"\b(import|using)\s+((\w|\.)+\s*,\s*)*", lib, raw"\b")
    success(`git -C $repo grep -Eq $pattern $tree`)
end

function gitmeta(pkgs::Dict{String,Package})
    @assert length(STDLIBS) ≤ 64 # use 64 bits to encode usage
    fd = joinpath(@__DIR__, "sha1map.toml")
    fs = joinpath(@__DIR__, "stdlib.toml")
    d = ispath(fd) ? TOML.parsefile(fd) : Dict()
    s = Dict()
    if ispath(fs)
        s = TOML.parsefile(fs)
        get(s, "STDLIBS", nothing) == STDLIBS || empty!(s)
        mv(fs, "$fs.old", force=true)
    end
    s["STDLIBS"] = STDLIBS
    io = open(fs, "w")
    println(io, "STDLIBS = [")
    for lib in STDLIBS
        println(io, "    ", repr(lib), ",")
    end
    println(io, "]")
    println(io)
    for (pkg, p) in sort!(collect(pkgs), by=first)
        (pkg == "julia" || isempty(p.versions)) && continue
        uuid = string(p.uuid)
        @info "Package [$uuid] $pkg"
        haskey(d, uuid) || (d[uuid] = Dict())
        haskey(s, uuid) || (s[uuid] = Dict())
        updated = false
        repo_path = joinpath(homedir(), ".julia", "clones", uuid)
        repo = nothing
        for (ver, v) in p.versions
            haskey(d[uuid], v.sha1) &&
            if v"0.7" ∉ v.requires["julia"].versions ||
                haskey(s[uuid], v.sha1)
                continue
            end
            if repo == nothing
                repo = ispath(repo_path) ? LibGit2.GitRepo(repo_path) : begin
                    updated = true
                    @info "Cloning [$uuid] $pkg"
                    LibGit2.clone(p.url, repo_path, isbare=true)
                end
            end
            git_commit_hash = LibGit2.GitHash(v.sha1)
            if !updated
                try LibGit2.GitObject(repo, git_commit_hash)
                catch err
                    if !(err isa LibGit2.GitError && err.code == LibGit2.Error.ENOTFOUND)
                        rethrow(err)
                    end
                    @info "Updating $pkg from $(p.url)"
                    LibGit2.fetch(
                        repo,
                        remoteurl=p.url,
                        refspecs=["+refs/*:refs/remotes/cache/*"]
                    )
                end
            end
            failed = false
            git_commit = try LibGit2.GitObject(repo, git_commit_hash)
            catch err
                failed = true
                if !(err isa LibGit2.GitError && err.code == LibGit2.Error.ENOTFOUND)
                    rethrow(err)
                end
                @error("$pkg: git object $(v.sha1) could not be found")
            end
	    failed && continue
            if !(git_commit isa LibGit2.GitCommit || git_commit isa LibGit2.GitTag)
                error("$pkg: git object $(v.sha1) not a commit – $(typeof(git_commit))")
            end
            git_tree = LibGit2.peel(LibGit2.GitTree, git_commit)
            @assert git_tree isa LibGit2.GitTree
            git_tree_hash = string(LibGit2.GitHash(git_tree))
            d[uuid][v.sha1] = git_tree_hash
            # scan for stdlib dependencies
            v"0.7" in v.requires["julia"].versions || continue
            haskey(s[uuid], v.sha1) && continue
            libs = [uses(repo_path, git_tree_hash, lib) for lib in STDLIBS]
            s[uuid][v.sha1] = sum(d*2^(i-1) for (i,d) in enumerate(libs))
        end
        isempty(s[uuid]) && continue
        println(io, "[$uuid]")
        for (sha1, n) in sort!(collect(s[uuid]), by=string∘first)
            println(io, "$sha1 = $n")
        end
        println(io)
        flush(io)
    end
    open(fd, "w") do io
        TOML.print(io, d, sorted=true)
    end
    close(io)
    rm("$fs.old", force=true)
    return d, s
end


text:
# calculate the least squares line of the errors
function calc_line()
  dirs = Array(1:2)
  npts = length(dirs)
  println("maximum(dirs) = ", maximum(dirs))
  err_vals = zeros(maximum(dirs))
  h_vals = zeros(maximum(dirs))
  for i in dirs
    println("i = ", i)
    fname = "m$i/error_calc.dat"
#    fname2 = "m$i/counts_0.txt"
    val_i = readdlm(fname)
    err_vals[i] = val_i[1]

#    counts_i = readdlm(fname2)
#    h_vals[i] = 1/(sqrt(counts_i[3]/2))
    h_vals[i] = val_i[2]

  end

  h_vals = [0.5 0.2 0.1]

  min_idx = minimum(dirs)
  max_idx = maximum(dirs)
  err_vals = err_vals[min_idx:max_idx]
  h_vals = h_vals[min_idx:max_idx]
  npts = max_idx - min_idx + 1

  println("h values = ", h_vals)
  println("error values = ", err_vals)

  err_vals_log = zeros(npts)
  h_vals_log = zeros(npts)
  for i=1:npts
    err_vals_log[i] = log(err_vals[i])
    h_vals_log[i] = log(h_vals[i])
  end

  # now form system

  A = ones(npts, 2)
  A[:, 2] = h_vals_log

  # solve system
  AAt = A.'*A
  rhs = A.'*err_vals_log

  x = AAt\rhs

  println("x = ", x)

  for i=1:npts
    h_i = h_vals_log[i]
    err_exp = x[1] + h_i*x[2]
    err_act = err_vals_log[i]
    err_pct = 100*(err_exp - err_act)/err_act
    println("point $i percent error of best fit line = ", err_pct, "%")
  end

  # calculate slope from last two points
  slope2 = (err_vals_log[npts - 1] - err_vals_log[npts])/(h_vals_log[npts - 1] - h_vals_log[npts])

  println("2 point slope: ", slope2)
  println("least squares slope: ", x[2])

  # write data to file for plotting
  data_arr = hcat(h_vals, err_vals)
  writedlm("err_data.dat", data_arr)

  return x[2]  # return least squares slope
end


text:
import Pkg
Pkg.activate(".")

using PAM
print("Username: ")
username = readline()
password = read(Base.getpass("Password"), String)
println(authenticate(username, password))


text:
module NARemovers

using Random
using DataFrames: DataFrame, nrow

using ..AbsTypes
using ..BaseFilters
using ..Utils

import ..AbsTypes: fit, fit!, transform, transform!
export fit, fit!, transform, transform!
export NARemover

"""
    NARemover(
      Dict(
        :name => "nadetect",
        :acceptance => 0.10 # tolerable NAs percentage
      )
    )

Removes columns with NAs greater than acceptance rate.
This assumes that it processes columns of features.
The output column should not be part of input to avoid
it being excluded if it fails the acceptance critera.

Implements `fit!` and `transform!`.
"""
mutable struct NARemover <: Transformer
   name::String
   model::Dict{Symbol,Any}

   function NARemover(args::Dict = Dict())
      default_args = Dict{Symbol,Any}(
         :name => "nadetect",
         :acceptance => 0.10
      )
      cargs=nested_dict_merge(default_args,args)
      cargs[:name] = cargs[:name]*"_"*randstring(3)
      new(cargs[:name],cargs)
   end
end

"""
NARemover(acceptance::Float64)

Helper function for NARemover.
"""
NARemover(acceptance::Float64) = NARemover(Dict(:acceptance => acceptance))

"""
    fit!(nad::NARemover,features::DataFrame,labels::Vector=[])

Checks and exit of df is empty

# Arguments
- `nad::NARemover`: custom type
- `features::DataFrame`: input
- `labels::Vector=[]`:
"""
function fit!(nad::NARemover, features::DataFrame, labels::Vector=[])::Nothing
   return nothing
end

function fit(nad::NARemover, features::DataFrame, labels::Vector=[])::NARemover
   fit!(nad,features,labels)
   return deepcopy(nad)
end

"""
    transform!(nad::NARemover,nfeatures::DataFrame)

Removes columns with NAs greater than acceptance rate.

# Arguments
- `nad::NARemover`: custom type
- `nfeatures::DataFrame`: input
"""
function transform!(nad::NARemover, nfeatures::DataFrame)::DataFrame
   isempty(nfeatures) && return DataFrame()
   features = deepcopy(nfeatures)
   sz = nrow(features)
   tol = nad.model[:acceptance]
   colnames = []
   for (colname,dat) in collect(pairs(eachcol(features)))
      if sum(ismissing.(dat)) < tol*sz
         push!(colnames,colname)
      end
   end
   xtr =  features[:,colnames]
   return xtr
end

function transform(nad::NARemover, nfeatures::DataFrame)::DataFrame
   return transform!(nad,nfeatures)
end

end


text:
###############################
# functionalize
###############################
struct _functionalize{F} <: Function; f::F; end
@inline (f::_functionalize)(a...; ka...) = f.f(a...; ka...)

@inline functionalize(f) = _functionalize(f)
@inline functionalize(f::Function) = f
@inline functionalize(s::Symbol) = o -> getproperty(o, s)
@inline functionalize(s::Union{Tuple, AbstractArray, Base.Generator, Base.ValueIterator}) =
    (a...; ka...) -> map(f->functionalize(f)(a...; ka...), s)
@inline functionalize(s::Union{NamedTuple, AbstractDict}) = functionalize(values(s))

###############################
# arg
###############################

struct _Args{A,KA} <: Function
    a::A
    ka::KA
end
@inline (a::_Args)(f) = functionalize(f)(a.a...; a.ka...)

@inline arg(a...; ka...) = _Args(a, ka)

###############################
# grb, mth, asn, grbs, asns
###############################

#=
x |> grb.a = x.a
x |> grb[10] = x[10]
x |> grb.a[10] = x.a[10]
x |> grb(:a) = x.a

x |> asn.a(3) ===> (x.a = 3; x)
x |> asn.a[10](3) ===> (x.a[10] = 3; x)
x |> asn[10](3) ===> (x[10] = 3; x)
x |> asn(:a)(3) ===> (x.a = 3; x)

x |> mth.a(10) == x.a(10)
x |> mth[10](3) == x[10](3)
x |> mth(:a)(3) == x.a(3)

x |> grbs[5, 3] ===> [x[5], x[3]]
x |> grbs(:a, :b) ===> [x.a, x.b]

x |> asns[5, 3](10, 20) ===> (x[5] = 10; x[3] = 20; x)
x |> asns(:a, :b)(10, 20) ===> (x.a = 10; x.b = 20; x)
=#

abstract type _AbstGet <: Function end

Base.getproperty(g::_AbstGet, a::Symbol) =
    begin
        itms = Base.getfield(g, :__itms);
        a == :__itms ? itms : (append!(itms, [a]); g)
    end
Base.getindex(g::_AbstGet, a...) = (append!(g.__itms, a); g)

_grb_prp_idx(o, a::Symbol) = Base.getproperty(o, a)
_grb_prp_idx(o, a) = Base.getindex(o, a...)
_asn_prp_idx(o, a::Symbol, x) = Base.setproperty!(o, a, x)
_asn_prp_idx(o, a, x) = Base.setindex!(o, x, a...)

########### grb
struct _GrbSingleton end
const grb = _GrbSingleton()
(s::_GrbSingleton)(atr::Symbol) = Base.getproperty(s, atr)
Base.getproperty(g::_GrbSingleton, a::Symbol) = _Grb([a])
Base.getindex(g::_GrbSingleton, a...) = _Grb([a])
mutable struct _Grb <: _AbstGet __itms::Vector{Any} end
(g::_Grb)(obj::Any) = reduce(_grb_prp_idx, g.__itms; init=obj)

########### asn
struct _AsnSingleton end
const asn = _AsnSingleton()
(::_AsnSingleton)(atr::Symbol) = x -> (o -> Base.setproperty!(o, atr, x))
Base.getproperty(g::_AsnSingleton, a::Symbol) = _Asn([a])
Base.getindex(g::_AsnSingleton, a...) = _Asn([a])
mutable struct _Asn <: _AbstGet __itms::Vector{Any} end
(g::_Asn)(x) =
    obj -> (y = reduce(_grb_prp_idx, g.__itms[1:end-1]; init=obj);
            _asn_prp_idx(y, g.__itms[end], x);
            obj)

########### grbs
struct _GrbsSingleton end
const grbs = _GrbsSingleton()
(::_GrbsSingleton)(a...) = o -> map(x -> Base.getproperty(o, x), a)
Base.getindex(g::_GrbsSingleton, a...) = o -> map(x -> Base.getindex(o, x), a)

########### asns
struct _AsnsSingleton end
const asns = _AsnsSingleton()
(::_AsnsSingleton)(a...) =
    (b...) -> (o -> (map((atr, x) -> Base.setproperty!(o, atr, x), zip(a, b)); o))
Base.getindex(g::_AsnsSingleton, a...) =
    (b...) -> (o -> (map((atr, x) -> Base.setindex!(o, atr, x), zip(a, b)); o))

########### mth
struct _MthSingleton end
const mth = _MthSingleton()
(s::_MthSingleton)(atr::Symbol) =
    (a...; ka...) -> (o -> Base.getproperty(o, atr)(a...; ka...))
Base.getproperty(g::_MthSingleton, a::Symbol) = _Mth([a])
Base.getindex(g::_MthSingleton, a...) = _Mth([a])
mutable struct _Mth <: _AbstGet __itms::Vector{Any} end
(g::_Mth)(a...; ka...) =
    obj -> reduce(_grb_prp_idx, g.__itms; init=obj)(a...;ka...)

###########
# cry, wc, wd
###########

cry(f, n::Integer=1) =
    (x...; xkwa...) ->
        ((y...; ykwa...) ->
             f(x[1:(n==0 ? 0 : (n>0 ? n-1 : end+n+1))]...,
               y...,
               x[(n==0 ? 1 : (n>0 ? n : end+n+2)):end]...;
               xkwa..., ykwa...))
               # merge(n>=0 ? xkwa : ykwa, n>=0 ? ykwa : xkwa)...))
    #=
    begin
        if n < 0
            (x...; xkwa...) -> ((y...; ykwa...) -> f(x[1:end+n+1]...,
                                                     y...,
                                                     x[end+n+2:end]...;
                                                     merge(xkwa, ykwa)...))
        elseif n > 0
            (x...; xkwa...) -> ((y...; ykwa...) ->f(x[1:n-1]...,
                                                    y...,
                                                    x[n:end]...;
                                                    merge(xkwa, ykwa)...))
        else # n == 0
            cry(f,1)
        end
    end
    =#

Base.getproperty(::typeof(cry), p::Symbol) = cry(eval(p))

wc(f) = cry(f)
Base.getproperty(::typeof(wc), p::Symbol) = wc(eval(p))

wd(f) = cry(f, -1)
Base.getproperty(::typeof(wd), p::Symbol) = wd(eval(p))

macro cry(f)
    esc( :($f(a...;ka...) = x -> $f(x, a...; ka...)))
end


text:
using Gadfly, RDatasets, Dates

set_default_plot_size(6inch, 3inch)

approval = dataset("Zelig", "approval")
dates = Date[Date(y, m)
             for (y, m) in zip(approval.Year, approval.Month)]
approval.Date = dates

p = plot(approval, x="Date", y="Approve", Geom.line)


text:
#=
Parallel direct estimation of the probability for the Schlögl chemical reaction network
with X(0)= x₀ to reach X(T) ∈ A, the target set.

This reaction network corresponds to:
A + 2S → 3S
3S → A + 2S
B → S
S → B

Parameters were chosen so that the system is bistable, as in:
https://doi.org/10.1063/1.5017955
=#
using Distributed
# nw = 4; # number of workers
# addprocs(nw);

using Statistics
using HypothesisTests
using Printf

@everywhere include("schlogl_setup.jl")

Random.seed!(100)
@everywhere n_trials = 10^6;
n_in_set = @sync @distributed (+) for j in 1:n_trials
    X₀ = copy(x₀);
    f(integrate_schlogl(X₀, Tmax));
end
CI = confint(BinomialTest(n_in_set, n_trials));
@printf("P Estimate = %g\n",n_in_set/n_trials);
@printf("95%% CI = (%g,%g)\n", CI[1], CI[2])


text:
using Documenter
using Polynomials

DocMeta.setdocmeta!(Polynomials, :DocTestSetup, :(using Polynomials); recursive = true)

makedocs(
    modules = [Polynomials],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "Polynomials.jl",
    authors = "Jameson Nash, Keno Fischer, and other contributors",
    pages = [
        "Home" => "index.md",
        "Reference/API" => "reference.md",
        "Polynomial Types" => [
            "Polynomial" => "polynomials/polynomial.md",
            "Chebyshev" => "polynomials/chebyshev.md",
        ],
        "Extending" => "extending.md",
    ],
)

deploydocs(repo = "github.com/JuliaMath/Polynomials.jl.git")


text:
function tableaus_vprk_glrk()
    (
        ( TableauVPGLRK(1),          "vprk_glrk1" ),
        ( TableauVPGLRK(2),          "vprk_glrk2" ),
        ( TableauVPGLRK(3),          "vprk_glrk3" ),
        ( TableauVPGLRK(4),          "vprk_glrk4" ),
        ( TableauVPGLRK(5),          "vprk_glrk5" ),
        ( TableauVPGLRK(6),          "vprk_glrk6" ),
    )
end

function tableaus_vprk_lobatto_IIIA_IIIB()
    (
        ( TableauVPLobattoIIIA(2),    "vprk_lobatto_IIIA_IIIB2" ),
        ( TableauVPLobattoIIIA(3),    "vprk_lobatto_IIIA_IIIB3" ),
        ( TableauVPLobattoIIIA(4),    "vprk_lobatto_IIIA_IIIB4" ),
        ( TableauVPLobattoIIIA(5),    "vprk_lobatto_IIIA_IIIB5" ),
    )
end

function tableaus_vprk_lobatto_IIIB_IIIA()
    (
        ( TableauVPLobattoIIIB(2),    "vprk_lobatto_IIIB_IIIA2" ),
        ( TableauVPLobattoIIIB(3),    "vprk_lobatto_IIIB_IIIA3" ),
        ( TableauVPLobattoIIIB(4),    "vprk_lobatto_IIIB_IIIA4" ),
        ( TableauVPLobattoIIIB(5),    "vprk_lobatto_IIIB_IIIA5" ),
    )
end

function tableaus_vprk_radau()
    (
        ( TableauVPRadauIIAIIA(2),     "vprk_radau_IIA2" ),
        ( TableauVPRadauIIAIIA(3),     "vprk_radau_IIA3" ),
        ( TableauVPRadauIIAIIA(4),     "vprk_radau_IIA4" ),
        ( TableauVPRadauIIAIIA(5),     "vprk_radau_IIA5" ),
    )
end

function tableaus_srk_glrk()
    (
        ( TableauGLRK(1),          "srk_glrk1",      IntegratorSRKimplicit ),
        ( TableauGLRK(2),          "srk_glrk2",      IntegratorSRKimplicit ),
        ( TableauGLRK(3),          "srk_glrk3",      IntegratorSRKimplicit ),
        ( TableauGLRK(4),          "srk_glrk4",      IntegratorSRKimplicit ),
        ( TableauGLRK(5),          "srk_glrk5",      IntegratorSRKimplicit ),
        ( TableauGLRK(6),          "srk_glrk6",      IntegratorSRKimplicit ),
    )
end

function tableaus_firk_glrk()
    (
        ( TableauGLRK(1),          "firk_glrk1",     IntegratorFIRKimplicit ),
        ( TableauGLRK(2),          "firk_glrk2",     IntegratorFIRKimplicit ),
        ( TableauGLRK(3),          "firk_glrk3",     IntegratorFIRKimplicit ),
        ( TableauGLRK(4),          "firk_glrk4",     IntegratorFIRKimplicit ),
        ( TableauGLRK(5),          "firk_glrk5",     IntegratorFIRKimplicit ),
        ( TableauGLRK(6),          "firk_glrk6",     IntegratorFIRKimplicit ),
    )
end


text:
using DynamicGrids, Dispersal, Test, Unitful, Dates
using DynamicGrids: applyrule, SimData, extent
using Unitful: d

@testset "loglogistic survival" begin
    init =  [10.0 40.0 70.0;
             20.0 50.0 80.0;
             30.0 60.0 90.0]

    test1 = [10.0 40.0 70.0;
             20.0 50.0 80.0;
             30.0 60.0 90.0]
    test2 = [ 3.0327   12.1308  21.2289;
              6.06539  15.1635  24.2616;
              9.09809  18.1962  27.2943]
    test3 = [ 0.919724  3.6789   6.43807;
              1.83945   4.59862  7.35779;
              2.75917   5.51835  8.27752]

    output = ArrayOutput(init; tspan=1:3)
    rule = Ruleset(SurvLogLogistic(exposure = 100, LC50=50, hillcoefficient = 1.2))
    sim!(output, rule)

    @test output[1] == test1
    @test output[2] ≈ test2 atol=1e-4
    @test output[3] ≈ test3 atol=1e-4

end

@testset "loglogistic survival with exposure from exposure map" begin
    init =  [10.0 40.0 70.0;
            20.0 50.0 80.0;
            30.0 60.0 90.0]

    suit = [100.0 100.0 100.0;
            50.0  50.0  50.0;
            20.0  20.0  20.0]

    output = ArrayOutput(init; tspan=1:3, aux=(suit=suit,))
    output.extent
    rule = Ruleset(SurvLogLogisticMap(layerkey=Val(:suit)))
    sim!(output, rule)

    @test output[1] == [10.0 40.0 70.0;
                        20.0 50.0 80.0;
                        30.0 60.0 90.0]
    @test output[2] ≈  [4.42688  17.7075  30.9882;
                        9.19701  22.9925  36.7881;
                        14.4803   28.9607  43.441] atol=1e-4
    @test output[3] ≈  [1.95973   7.83892  13.7181;
                        4.22925  10.5731   16.917;
                        6.98935  13.9787   20.968] atol=1e-4

    @test DynamicGrids.normalise.(output[3], 0.25, 4) ≈ [0.455928  1.0  1.0;
                                                          1.0       1.0  1.0;
                                                          1.0       1.0  1.0] atol=1e-4

    output = ArrayOutput(init; tspan=1:3, aux=(suit=suit,))
    rule = Ruleset(SurvLogLogisticMap(layerkey=Val(:suit), LC50=10))
    sim!(output, rule)

    @test output[1] == [10.0 40.0 70.0;
                        20.0 50.0 80.0;
                        30.0 60.0 90.0]
    @test output[2] ≈ [4.42688  17.7075  30.9882;
                       9.19701  22.9925  36.7881;
                      14.4803   28.9607  43.441] atol=1e-4
    @test output[3] ≈ [1.95973   7.83892  13.7181
                       4.22925  10.5731   16.917
                       6.98935  13.9787   20.968] atol=1e-4
end

@testset "loglogistic survival with ArrayOutput layer" begin
    init =  [10.0 40.0 70.0;
            20.0 50.0 80.0;
            30.0 60.0 90.0]

    suit = [100.0 100.0 100.0;
            50.0  50.0  50.0;
            20.0  20.0  20.0]

    # suit = output
    # outputStack = ArrayOutput(init; tspan=1:3, aux=(suit=suit,index=1:3))
    # rule = Ruleset(SurvLogLogisticMap(layerkey=:suit, timeindex=:index))
    # sim!(outputStack, rule)

end

@testset "survival map masking" begin
    init = [1.0 1.0 1.0;
        1.0 1.0 1.0;
        1.0 1.0 1.0]

    suit =  [1.0 1.0 2.0;
        2.0 2.0 0.5;
        1.0 1.0 0.5]

    test1 = [1.0 1.0 1.0;
        1.0 1.0 1.0;
        1.0 1.0 1.0]

    test2 = [0.0 0.0 1.0;
            1.0 1.0 0.0;
            0.0 0.0 0.0]

    output = ArrayOutput(init; tspan=1:3, aux=(suit=suit,))
    maskrule = MaskSurvMap(layerkey=Val(:suit), threshold=1.1)
    ruleset = Ruleset(maskrule)
    data = SimData(extent(output), ruleset)

    @test applyrule(data, maskrule, 5, (1, 1)) === 0
    @test applyrule(data, maskrule, 5, (2, 2)) === 5

    @test applyrule(data, maskrule, 5.0, (1, 1)) === 0.0
    @test applyrule(data, maskrule, 5.0, (2, 2)) === 5.0

    sim!(output, Ruleset(maskrule))

    @test output[1] == test1
    @test output[2] == test2
end

@testset "Test double layers ExactLogisticGrowthMap2" begin
    popSizeInit = [ 1.0 4.0 7.0;
                    2.0 5.0 8.0;
                    3.0 6.0 9.0]

    Exposure = repeat([0.0  0.0  0.0;
                        20.0 20.0 20.0;
                        50.0 50.0 50.0],
                    inner=(1, 1,3))

    LC50 = cat( [0.0 0.0 0.0;
                0.0 0.0 0.0;
                0.0 0.0 0.0],
                [ 20.0 20.0 20.0;
                20.0 20.0 20.0;
                20.0 20.0 20.0],
                [ 10.0 10.0 10.0;
                10.0 10.0 10.0;
                10.0 10.0 10.0]; dims=3)

    survParameter = cat(Exposure, LC50; dims = 4)

    popSizeGrids = ArrayOutput(popSizeInit; tspan=1:6, aux=(survParameter=survParameter,));
    survRule = Ruleset(SurvLogLogisticMap2(layerkey=:survParameter));
    sim!(popSizeGrids, survRule);
end

@testset "Test double grids DiscreteGrowth2" begin
    init = (pop1 = [.5 0. 0.;
                0. 0. 1.;
                0. 0. 0.],
            pop2 = [.5 0. 0.;
                0. 1. 0.;
                0. 0. 0.],
            pop3 = [.5 0. 0.;
                0. 1. 0.;
                0. 1. 0.],)

    exposure = repeat([0.0  0.0  0.0;
                20.0 20.0 20.0;
                50.0 50.0 50.0],
                inner=(1, 1,3))

    output = ArrayOutput(init; tspan=1:3, aux=(exposure=exposure,))

    ruleGrowthSurv =  GrowthSurvLogLogisticMap3{Tuple{:pop1,:pop2,:pop3},Tuple{:pop1,:pop2,:pop3}}(
        layerkey=:exposure,
        intrinsicrate=1.5,
        hillcoefficient=2.5,
        LC501=50,
        LC502=40,
        LC503=30,
        carrycap = 5.)

    sim!(output,ruleGrowthSurv);

    @test output[1] == init

    @test output[2].pop1 ≈ [.525 0. 0.;
                        0. 0. 1.08973;
                        0. 0. 0.] atol=1e-4
    @test output[2].pop2 ≈ [.525 0. 0.;
                        0. 0.764801  0.;
                        0. 0. 0.] atol=1e-4
    @test output[2].pop3 ≈ [.525 0. 0.;
                        0. 0.660363 0.;
                        0. 0.26166 0.] atol=1e-4
end


text:
"""Provides an abstract type for spectroscopic graphs,
allowing for other representations of spectroscopic graphs
to be developed.
"""
abstract type SpectroscopicGraph end

struct LevelsGraph <: SpectroscopicGraph
	levels
	graph
end

"""A projected spectroscopic monograph, where nodes represent
energy levels, and edges between nodes represent transitions
that connect them.
"""
struct ProjectedSG <: SpectroscopicGraph
  levels
  transitions
  graph
end

"""A bipartite spectroscopic graph, with two classes of nodes:
energy levels and transitions. Transition nodes are connected to
two energy levels nodes.
"""
struct BipartiteSG <: SpectroscopicGraph
  levels
  transitions
  graph
end

Base.length(sg::SpectroscopicGraph) = length(sg.levels) + length(sg.transitions)
offset(sg::SpectroscopicGraph) = length(sg.levels)
indices(sg::SpectroscopicGraph) = 1:length(sg)

function ProjectedSG(ls, ts)
	sg = MetaGraph(length(ls))
	for t ∈ ts
        properties = Dict(:I => t.I, :ν => t.ν)
		lq, uq = lower_state_qnos(t), upper_state_qnos(t)
		l_index, u_index = match_encoding(ls, lq), match_encoding(ls, uq)
        MetaGraphs.add_edge!(sg, l_index, u_index, properties)
	end
	return ProjectedSG(ls, ts, sg)
end

function BipartiteSG(ls, ts)
	sg = MetaGraph(length(ls) + length(ts))
	for (t_index, t) ∈ enumerate(ts)
		t_offset = length(ls) + t_index + 1
        # get the transition strength and frequency, potentially for plotting
        properties = Dict(:I => t.I, :ν => t.ν)
		lq, uq = lower_state_qnos(t), upper_state_qnos(t)
		l_index, u_index = match_encoding(ls, lq), match_encoding(ls, uq)
		for index in [l_index, u_index]
          MetaGraphs.add_edge!(sg, t_offset, index, properties)
		end
	end
	return BipartiteSG(ls, ts, sg)
end

function LevelsGraph(ls)
	sg = MetaGraph(length(ls))
	return LevelsGraph(ls, sg)
end


text:
mutable struct DataTable

    dims        :: Dict
    grid_dims   :: Dict
    grid_dims2_str  :: Dict
    data_units  :: Dict
    missing_idx1 :: Dict   # with shape (z, x, y)
    missing_idx2 :: Dict   # with shape (x, y, z)
    flags       :: Dict    # What are flags for?

    function DataTable(;
        Nx,
        Ny,
        Nz,
        mask_sT    :: Union{Nothing, AbstractArray{Float64, 3}} = nothing,
        mask_T     :: Union{Nothing, AbstractArray{Float64, 3}} = nothing,
    )

        N1 = 1
        Nyp1 = Ny+1
        Nzp1 = Nz+1

        dims = Dict(
            "N1" => N1,
            "Nx" => Nx,
            "Ny" => Ny,
            "Nz" => Nz,
            "Nyp1" => Nyp1,
            "Nzp1" => Nzp1,
        )

        grid_dims = Dict(
            :T  => (Nz,   Nx  , Ny  ),
            :U  => (Nz,   Nx  , Ny  ),
            :V  => (Nz,   Nx  , Nyp1),
            :W  => (Nzp1, Nx  , Ny  ),
            :sT => (N1, Nx, Ny  ),
            :sU => (N1, Nx, Ny  ),
            :sV => (N1, Nx, Nyp1),
            :SCALAR => (N1, N1, N1),
            :cW => (Nzp1, N1, N1),
            :cT => (Nz, N1, N1),
        )

        # This is used for RecordTool output
        # Notice that this is oriented dimension so z is the last one
        grid_dims2_str = Dict(
            :T  => ("Nx", "Ny", "Nz"),
            :U  => ("Nx", "Ny", "Nz"),
            :V  => ("Nx", "Nyp1", "Nz"),
            :W  => ("Nx", "Ny", "Nzp1"),
            :sT => ("Nx", "Ny", "N1"),
            :sU => ("Nx", "Ny", "N1"),
            :sV => ("Nx", "Nyp1", "N1"),
            :SCALAR => ("N1", "N1", "N1"),
            :cW => ("N1", "N1", "Nzp1",),
            :cT => ("N1", "N1", "Nz",),
        )

        data_units = Dict()
        flags      = Dict()

        missing_idx1 = Dict(
            :sT => (mask_sT == nothing) ? nothing : mask_sT .== 0,
            :T  => (mask_T == nothing) ? nothing : mask_T .== 0,
            :U => nothing,
            :V => nothing,
            :W => nothing,
            :UV => nothing,
            nothing => nothing,
        )

        missing_idx2 = Dict()
        for (k,v) in missing_idx1
            missing_idx2[k] = (v != nothing) ? permutedims(v, [2, 3, 1]) : nothing
        end

        return new(
            dims,
            grid_dims,
            grid_dims2_str,
            data_units,
            missing_idx1,
            missing_idx2,
            flags,
        )

    end

end

function regVariable!(
    dt       :: DataTable,
    id       :: Union{Symbol, String},
    grid     :: Symbol,
    mask     :: Union{Symbol, Nothing},
    data     :: AbstractArray{T},
) where T

    N1 = dt.dims["N1"]
    Nx = dt.dims["Nx"]
    Ny = dt.dims["Ny"]
    Nz = dt.dims["Nz"]
    Nyp1 = dt.dims["Nyp1"]
    Nzp1 = dt.dims["Nzp1"]

    if haskey(dt.data_units, id)
        throw(ErrorException("Error: variable id " * String(id) *  " already exists."))
    end

    dim = dt.grid_dims[grid]

    dtype = eltype(data)
    if ! (dtype in (Float64, Int64))
        throw(ErrorException("Invalid data type. Only Float64 and Int64 are accepted"))
    end

    if Tuple(dim) != size(data)
        println("Expect ", dim)
        println("Get ", size(data))
        throw(ErrorException("Provided data does not have correct dimension: " * string(id)))
    end

    if dtype != T
        throw(ErrorException("dtype and provided data does not match."))
    end

    dt.data_units[id] = DataUnit(
        dt,
        id,
        grid,
        mask,
        data,
    )

    dt.flags[id] = 0
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
"""
    is_graphical(degs::Vector{Int})

Check whether the degree sequence `degs` is graphical, according to
[Erdös-Gallai condition](http://mathworld.wolfram.com/GraphicSequence.html).

Time complexity: O(length(degs)^2)
"""
function is_graphical(degs::Vector{Int})
    iseven(sum(degs)) || return false
    n = length(degs)
    cumdeg = cumsum(degs)
    for r=1:n-1
        cond = cumdeg[r] <= r*(r-1) + sum(i->min(r,degs[i]), r+1:n)
        cond || return false
    end
    return true
end

"""
    has_self_loops(g)

Returns true if `g` has any self loops.
"""
has_self_loops(g::AGraphOrDiGraph) = any(v->has_edge(g, v, v), vertices(g))

"""
    num_self_loops(g)

Returns the number of self loops in `g`.
"""
num_self_loops(g::AGraphOrDiGraph) = count(v->has_edge(g, v, v), vertices(g))


text:
const COVERAGE_GITIGNORE = ["*.jl.cov", "*.jl.*.cov", "*.jl.mem"]

"""
    Codecov(; file=nothing)

Sets up code coverage submission from CI to [Codecov](https://codecov.io).

## Keyword Arguments
- `file::Union{AbstractString, Nothing}`: Template file for `.codecov.yml`,
  or `nothing` to create no file.
"""
@plugin struct Codecov <: FilePlugin
    file::Union{String, Nothing} = nothing
end

source(p::Codecov) = p.file
destination(::Codecov) = ".codecov.yml"

badges(::Codecov) = Badge(
    "Coverage",
    "https://codecov.io/gh/{{{USER}}}/{{{PKG}}}.jl/branch/master/graph/badge.svg",
    "https://codecov.io/gh/{{{USER}}}/{{{PKG}}}.jl",
)

"""
    Coveralls(; file=nothing)

Sets up code coverage submission from CI to [Coveralls](https://coveralls.io).

## Keyword Arguments
- `file::Union{AbstractString, Nothing}`: Template file for `.coveralls.yml`,
  or `nothing` to create no file.
"""
@plugin struct Coveralls <: FilePlugin
    file::Union{String, Nothing} = nothing
end

source(p::Coveralls) = p.file
destination(::Coveralls) = ".coveralls.yml"

badges(::Coveralls) = Badge(
    "Coverage",
    "https://coveralls.io/repos/github/{{{USER}}}/{{{PKG}}}.jl/badge.svg?branch=master",
    "https://coveralls.io/github/{{{USER}}}/{{{PKG}}}.jl?branch=master",
)

gitignore(::Union{Codecov, Coveralls}) = COVERAGE_GITIGNORE
view(::Union{Codecov, Coveralls}, t::Template, pkg::AbstractString) = Dict(
    "PKG" => pkg,
    "USER" => t.user,
)

"""
    is_coverage(::Plugin) -> Bool

Determine whether or not a plugin is a coverage plugin.
If you are adding a coverage plugin, you should implement this function and return `true`.
"""
is_coverage(::Plugin) = false
is_coverage(::Union{Codecov, Coveralls}) = true

needs_username(::Union{Codecov, Coveralls}) = true


text:
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = page_model()
    include("../src/components/RCPSSPScenario.jl")
    include("../src/components/LGemissions.jl")

    rcpsspscenario = addrcpsspscenario(m, scenario)
    lgemit = add_comp!(m, LGemissions)

    lgemit[:er_LGemissionsgrowth] = rcpsspscenario[:er_LGemissionsgrowth]
    set_param!(m, :LGemissions, :e0_baselineLGemissions, readpagedata(m,"data/e0_baselineLGemissions.csv"))

    # run Model
    run(m)

    emissions= m[:LGemissions,  :e_globalLGemissions]
    emissions_compare=readpagedata(m, "test/validationdata/$valdir/e_globalLGemissions.csv")

    @test emissions ≈ emissions_compare rtol=1e-3
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
module FastaTests

using FastaIO
using GZip
using Test

const fastadata_ascii = Any[
    ("A0ADS9_STRAM/3-104",
     "-------------------------------------------STVELTKEN-F--D-Q-" *
     "-T-V----T-D---------N-----------E--F-V--LI--------D-----F--W" *
     "----A----S-W---------C---G------P----C-----R----Q---F------A" *
     "------P---V----Y-------E---K---A---A--E---A---NP------------" *
     "----------------DL---V--F--G----K-------V-------D-----------" *
     "T-----------E--------A------Q-----------------P------E------" *
     "------L----A--------Q-----A------F------G--------I------Q---" *
     "-----S------I------P-T---L--M------I---V-------R---D------Q-" *
     "---V----A-VF-----------AQP-GA----LP-EE-ALTDVI---GQA---------" *
     "------------------------------------------------------------" *
     "---------"),
    ("A0AHX4_LISW6/1-102",
     "-------------------------------------------MVKEITDAT-F--E-Q-" *
     "-E-T----S-E---------------------G--L-V--LT--------D-----F--W" *
     "----A----T-W---------C---G------P----C-----R----M---V------A" *
     "------P---V----L-------E---E---I---Q--E---E---RGE-----------" *
     "----------------AL---K--I--V----K-------M-------D-----------" *
     "V-----------D--------E------N-----------------P------E------" *
     "------T----P--------G-----S------F------G--------V------M---" *
     "-----S------I------P-T---L--L------I---K-------K---D------G-" *
     "---E----V-VE-----------TII-GY----HP-KE-ELDEVINK-Y-----------" *
     "------------------------------------------------------------" *
     "---------"),
    ("A0AJ61_LISW6/3-103",
     "-----------------------------------NLESVEQFD----------------" *
     "---V----I-K-------S-E-----------G--K-S--VF--------M-----F--S" *
     "----A----D-W---------C---G------D----C-----K----Y---I------E" *
     "------P---V----M-------P---E---I---E--A---E---NE------------" *
     "----------------DF---T--F--Y----H-------V-------D-----------" *
     "R-----------D--------E------F-----------------I------D------" *
     "------L----C--------A-----D------L------A--------I------F---" *
     "-----G------I------P-S---F--L------V---F-------E---D------G-" *
     "---E----E-VG-----------RFV-SKD--RKT-KE-EINDFLA--AI----------" *
     "------------------------------------------------------------" *
     "---------"),
    ("A0AK08_LISW6/41-151",
     "-----------------------------------------FLN--TISTKD-F--K-Q-" *
     "-Q-M----A-D---------K-----------T--T-G--FV--------Y-----V--G" *
     "----R----P-T---------C---E------D----C-----Q----A---F------Q" *
     "------P---I----L-------K---K---E---L--K---E---RKLN----------" *
     "---------------QNM---N--Y--Y----N-------T-------D-----------" *
     "K-----------A--------S------E-----------------K------SRDD---" *
     "---MIAL----L--------K-----K------M------D--------I------D---" *
     "-----S------V------P-T---M--V------Y---L-------K---D------G-" *
     "---K----V-AS-----------TYA-AT----DE-PE-KLTHWMNK-V-----------" *
     "------------------------------------------------------------" *
     "---------")]

const fastadata_uint8 = map(x->(x[1],Vector{UInt8}(codeunits(x[2]))), fastadata_ascii)
const fastadata_char = map(x->(x[1],Vector{Char}(x[2])), fastadata_ascii)

function test_fastaread(T::Type, infile, fastadata)
    FastaReader(infile, T) do fr
    @testset "collect(fr)" begin
        @test collect(fr) == fastadata
    end

        rewind(fr)

    @testset "iterate(fr)" begin
        for (desc, seq) in fr
            @test desc == fastadata[fr.num_parsed][1]
            @test seq == fastadata[fr.num_parsed][2]
        end
    end

        rewind(fr)
    @testset "readentry(fr)" begin
        while !eof(fr)
            desc, seq = readentry(fr)
            @test desc == fastadata[fr.num_parsed][1]
            @test seq == fastadata[fr.num_parsed][2]
        end
    end
    end
end

function test_fastawrite(infile, outfile, fastadata)
    # char-by-char mode
    FastaWriter(outfile) do fw
        gzopen(infile) do f
            while !eof(f)
                write(fw, read(f, UInt8))
            end
        end
    end
    @test readfasta(outfile) == fastadata

    # writeentry mode
    FastaWriter(outfile) do fw
        for (desc, seq) in fastadata
            writeentry(fw, desc, seq)
        end
    end
    @test readfasta(outfile) == fastadata

    # one element at a time mode
    FastaWriter(outfile) do fw
        fd_plain = Any[]
        for (desc, seq) in fastadata
            push!(fd_plain, ">" * desc)
            i = 0
            while i < length(seq)
                j = rand(i+1:length(seq))
                push!(fd_plain, seq[i+1:j])
                i = j
            end
            #push!(fd_plain, seq)
        end
        write(fw, fd_plain)
    end
    @test readfasta(outfile) == fastadata

    # writefasta (vector)
    writefasta(outfile, fastadata)
    @test readfasta(outfile) == fastadata

    # writefasta (dict)
    ## how to do this in a backward-compatible way?
    ##   fd_dict = Dict(desc=>seq for (desc, seq) in fastadata)
    fd_dict = Dict(); for (desc, seq) in fastadata; fd_dict[desc]=seq; end
    writefasta(outfile, fd_dict)
    @test sort!(readfasta(outfile)) == sort(fastadata)
    return
end

@testset "reading test.fasta$(suffix) as $T" for
    suffix in ["", ".gz", ".win", ".win.gz", ".no_eof", ".no_eof.gz", ".blanklines", ".blanklines.gz"],
    (T, fastadata) in [(Vector{UInt8}, fastadata_uint8),
                       (Vector{Char}, fastadata_char),
                       (String, fastadata_ascii)]

    infile = joinpath(dirname(Base.source_path()), "test.fasta" * suffix)
    outfile = joinpath(dirname(Base.source_path()), "test_out.fasta" * suffix)

    test_fastaread(T, infile, fastadata)

    try
        test_fastawrite(infile, outfile, fastadata_ascii)
    finally
        isfile(outfile) && rm(outfile)
    end
end

@testset "reading invalid fasta #$i" for i in 1:4
    infile = joinpath(dirname(Base.source_path()), "invalid_test_$i.fasta.gz")

    @test_throws Exception readfasta(infile)
    @test_throws Exception FastaReader(infile) do fr
        while !eof(fr)
            desc, seq = readentry(fr)
        end
    end
end

outfile = joinpath(dirname(Base.source_path()), "invalid_test_out.fasta.gz")

@testset "writing invalid fasta #$i" for i in 2:6
    infile = joinpath(dirname(Base.source_path()), "invalid_test_$i.fasta.gz")

    try
        @test_throws ErrorException FastaWriter(outfile) do fw
            gzopen(infile) do f
                while !eof(f)
                    write(fw, read(f, UInt8))
                end
            end
        end
    finally
        isfile(outfile) && rm(outfile)
    end
end

const invalid_fastadata_entries = [("DE\nSC", "DATA" ),
                                   (""      , "DATA" ),
                                   ("DESC"  , ""     ),
                                   ("ΔΕΣΞ"  , "DATA" ),
                                   ("DESC"  , "ΔΑΤΑ" ),
                                   ("DESC"  , "DA>TA")]

for (desc, data) in invalid_fastadata_entries
    try
        @test_throws Exception FastaWriter(outfile) do fw
            writeentry(fw, desc, data)
        end
    finally
        isfile(outfile) && rm(outfile)
    end
end

for invalid_fastadata in [[ife] for ife in invalid_fastadata_entries]
    try
        @test_throws Exception writefasta(outfile, invalid_fastadata)
    finally
        isfile(outfile) && rm(outfile)
    end
end

end


text:
mutable struct Node
	x :: Int
	y :: Int
	value :: Int
	bestdist :: Int
end

struct Edge
	sideA :: Node
	sideB :: Node
end

other(edge :: Edge, node :: Node) = node == edge.sideA ? edge.sideB : node == edge.sideB ? edge.sideA : error("node not on either edge")
getedges(n :: Node, edges :: Vector{Edge}) = [edge for edge in edges if edge.sideA == n || edge.sideB == n]
makenode(value :: Int, x :: Int, y :: Int) = Node(x, y, value, typemax(Int))

function create_edges(grid_to_node :: Dict{Tuple{Int, Int}, Node})
	nodes = values(grid_to_node)
	edges_from_node = Dict{Node, Vector{Edge}}(n => Edge[] for n in nodes)
	for node in nodes
		if haskey(grid_to_node, (node.x + 1, node.y))
			target_node = grid_to_node[(node.x + 1, node.y)]
			new_edge = Edge(node, target_node)
			push!(edges_from_node[node], new_edge)
			push!(edges_from_node[target_node], new_edge)
		end

		if haskey(grid_to_node, (node.x, node.y + 1))
			target_node = grid_to_node[(node.x, node.y + 1)]
			new_edge = Edge(node, target_node)
			push!(edges_from_node[node], new_edge)
			push!(edges_from_node[target_node], new_edge)
		end
	end

	return edges_from_node
end

function dijkstra!(grid_to_node :: Dict{Tuple{Int, Int}, Node}, edges_from_node :: Dict{Node, Vector{Edge}})
	grid_to_node[(1, 1)].bestdist = 0
	eligibles = Set{Node}([grid_to_node[(1, 1)]])
	counter = 0

	while counter < length(grid_to_node)
		current = argmin(n -> n.bestdist, eligibles)
		node_edges = edges_from_node[current]
		for other_edge in node_edges
			target = other(other_edge, current)
			if current.bestdist + target.value < target.bestdist
				target.bestdist = current.bestdist + target.value
				push!(eligibles, target)
			end
		end

		pop!(eligibles, current)
		counter += 1
	end
	# modifies in-place, return not needed
end

function part1()
	lines = readlines("inputD15.txt")

	grid_to_node = Dict{Tuple{Int, Int}, Node}()
	for (i, line) in enumerate(lines)
		for (j, value) in enumerate(line)
			grid_to_node[(i, j)] = makenode(parse(Int, value), i, j)
		end
	end

	edges_from_node = create_edges(grid_to_node)
	dijkstra!(grid_to_node, edges_from_node)

	return grid_to_node[(100, 100)].bestdist
end

function add(n :: Int, amount :: Int)
	for _ in 1:amount
		n == 9 ? n = 1 : n += 1
	end

	return n
end

function part2()
	lines = readlines("inputD15.txt")

	grid_to_node = Dict{Tuple{Int, Int}, Node}()
	for (i, line) in enumerate(lines)
		for (j, value) in enumerate(line)
			value = parse(Int, value)
			for m in 0:4
				for n in 0:4
					grid_to_node[(i+n*100, j+m*100)] = makenode(add(value, n+m), i+n*100, j+m*100)
				end
			end
		end
	end

	edges_from_node = create_edges(grid_to_node)
	dijkstra!(grid_to_node, edges_from_node)

	return grid_to_node[(500, 500)].bestdist
end

println(part1())
println(part2())


text:
module Mock2

# add dependencies here
using DelimitedFiles, ForwardDiff, TickTock

# functions to export
export addone, addtwo, my_f, my_f_der

function addone(x)
    return x+1
end

function addtwo(x)
    return x+2
end

my_f(x)=x^2-5x+1;
my_f_der(x)=ForwardDiff.derivative(my_f,x)

end


text:
export Astrobee2D

mutable struct Astrobee2D{T<:AbstractFloat} <: Robot
  mass::T
  J::T
  Jinv::T
  n_thrusters::Int
  s::T
  r::T
  hard_limit_vel::T
  hard_limit_accel::T
  hard_limit_ω::T
  hard_limit_α::T
  btCollisionObject

  xb::Vector{T}
  Jcollision
end
function Astrobee2D{T}() where T
  n_thrusters = 12
  s = 0.5*0.305   # each side of cube is 30.5cm
  r = sqrt(2)*s   # inflate to sphere

  # ground robot param: freeflyer/astrobee/config/worlds/granite.config
  mass = 14.4
  hard_limit_vel = 0.20
  hard_limit_accel = 0.02
  hard_limit_ω = 10*π/180
  hard_limit_α = 10*π/180

  J = 0.1083
  Jinv = inv(J)

  xb = [0.; 0.15; 0.]
  Jcollision = []

  # btCollisionObject = BulletCollision.convex_hull_cylinder(SVector{3}(zeros(3)), SVector{3}(0.,0.,1.), r)
  btCollisionObjects = BulletCollision.BulletCollisionObjectPtr[]
  push!(btCollisionObjects, BulletCollision.convex_hull_cylinder(SVector{3}(zeros(3)), SVector{3}(0.,0.,1.), r))
  push!(btCollisionObjects, BulletCollision.convex_hull_cylinder(SVector{3}(xb...), SVector{3}(0.,0.,1.), r))

  btCollisionObject = BulletCollision.compound_collision_object(btCollisionObjects)

  # new astrobee instance
  return Astrobee2D{T}(mass, J, Jinv, n_thrusters, s, r, hard_limit_vel, hard_limit_accel, hard_limit_ω, hard_limit_α, btCollisionObject, xb, Jcollision)
end
Astrobee2D(::Type{T} = Float64; kwargs...) where {T} = Astrobee2D{T}(; kwargs...)


text:
f = open("jawiki-country.txt", "r")

for line in readlines(f)
  if ismatch(r"\[\[Category:", line)
    print(line)
  end
end


text:
using Test

import Functional
using Functional: or, default, when
const F = Functional

struct TestStruct
    x::Int
    s::String
end

mutable struct TestMutableStruct
    x::Int
    s::String
end

join_test() = (["a", "b", "c"] |> F.join(",")) == "a,b,c"
sum_test() = ([1,2,3,4] |> F.sum) == 10
enumerate_test() = (["a", "b", "c"] |> F.enumerate |> collect) == [(1, "a"), (2, "b"), (3, "c")]

take_test() = ([1,2,3] |> F.take(2) |> collect) == [1,2]
takewhile_test() = ([1,2,3,4] |> F.takewhile(x -> x < 3) |> collect) == [1,2]

skip_test() = ([1,2,3] |> F.skip(2) |> collect) == [3]
skipwhile_test() = ([1,2,3,4] |> F.skipwhile(x -> x < 3) |> collect) == [3,4]

map_test() = ([1,2,3] |> F.map(x -> 2 * x) |> collect) == [2,4,6]
zip_test() = (([1,2], ["a", "b"]) |> x -> F.zip(x) |> collect) == [(1, "a"), (2, "b")]
reverse_test() = ([1,2,3,4] |> F.reverse |> collect) == [4,3,2,1]

all_test() = ([1,2,3,4] |> F.all(x -> x != 3)) == false
any_test() = ([1,2,3,4] |> F.any(x -> x == 3)) == true

groupby_test() = begin
    a = TestStruct(10, "a")
    b = TestStruct(10, "b")
    c = TestStruct(11, "c")
    d = TestStruct(11, "d")
    ([a,b,c,d] |> F.groupby(e -> e.x) |> collect) == [[a,b], [c,d]]
end

apply_test() = ([TestMutableStruct(1, "ciao"),
                TestMutableStruct(2, "abaco")] |> F.apply(x -> x.x = 2)) === nothing

@testset "FunctionTests" begin
    @test join_test()
    @test sum_test()
    @test enumerate_test()

    @test take_test()
    @test takewhile_test()

    @test skip_test()
    @test skipwhile_test()

    @test map_test()
    @test zip_test()
    @test reverse_test()

    @test all_test()
    @test any_test()

    @test groupby_test()

    @test apply_test()
end

matching_test() = 50 |> F.match(
    10 => 10,
    20 => 20,
    50 => 50
) == 50

matching_or_test() = 50 |> F.match(
    or(10, 20) => 10,
    or(20, 50) => 20
) == 20

matching_when_test() = 50 |> F.match(
    (50 |> when(x -> x < 50)) => 10,
    (50 |> when(x -> x >= 50)) => 20
) == 20

matching_with_default() = 50 |> F.match(
    20 => 20,
    default() => 10
)  == 10

@testset "MatchingTests" begin
    @test matching_test()
    @test matching_or_test()
    @test matching_when_test()
    @test matching_with_default()
end


text:
export HCAS_MDP

# Define MDP
mutable struct HCAS_MDP <: MDP{stateType, actType}
    ranges::Array{Float64,1}
    thetas::Array{Float64,1}
    psis::Array{Float64,1}
    vowns::Array{Float64,1}
    vints::Array{Float64,1}
    pras::Array{Int64,1}
    discount_factor::Float64
    turns::Dict{Int,Tuple{Vector{Float64},Vector{Float64}}}
    currentTau::Float64
end
HCAS_MDP() = HCAS_MDP(RANGES,THETAS,PSIS,OWNSPEEDS,INTRSPEEDS,ACTIONS,discount_f,turns,0.0)
#HCAS_MDP(tau::Float64) = HCAS_MDP(RANGES,THETAS,PSIS,OWNSPEES,INTRSPEEDS,discount_f,turns,tau)

# Define necessary functions for HCAS MDP
POMDPs.actionindex(::HCAS_MDP, a::actType) = a + 1
POMDPs.actions(mdp::HCAS_MDP)  = ACTIONS
POMDPs.discount(mdp::HCAS_MDP) = mdp.discount_factor
POMDPs.n_actions(::HCAS_MDP)   = length(ACTIONS)

function POMDPs.convert_s(::Type{V} where V <: AbstractVector{Float64}, s::stateType, mdp::HCAS_MDP)
    v = [s[1],s[2],s[3],s[4],s[5],convert(Float64,s[6])]
    return v
end

function POMDPs.convert_s(::Type{stateType}, v::AbstractVector{Float64}, mdp::HCAS_MDP)
    s = (v[1],v[2],v[3],v[4],v[5],convert(Int,v[6]))
    return s
end

max_stars_repo_name: cherrywoods/HorizontalCAS

max_issues_repo_name: ahmed-irfan/HorizontalCAS

max_forks_repo_name: ahmed-irfan/HorizontalCAS


text:
using RootHitFiles, BenchmarkTools

const eventpath = realpath(joinpath(dirname(pathof(RootHitFiles)), "..", "test", "test.root.hits"))

SUITE = BenchmarkGroup()

SUITE["construct"] = @benchmarkable RootHitFile($eventpath)

SUITE["read"] = @benchmarkable read(RootHitFile($eventpath))

SUITE["collect"] = @benchmarkable collect(RootHitFile($eventpath))


text:
using Stheno, Plots, Random
plotly();

# Define model.
σ_noise, ω, T = 1e0, 1.0, 25.0;
gpc = GPC();
f₁ = GP(periodic(EQ(), ω), gpc);
f₂ = GP(scale(EQ(), 0.1), gpc);
f₃ = f₁ + f₂;
y = f₃ + GP(Noise(σ_noise^2), gpc);

# Sample from marginal process to generate toy data.
rng, S = MersenneTwister(123456), 25;
X, Xp = linspace(0.0, T, 100), linspace(-2.5, T + 2.5, 500);
ŷ = rand(rng, y(X));

# Plots for stuff.
f₁′_plot, f₂′_plot, f₃′_plot = plot(), plot(), plot();

posterior_plot = plot();
approx_in_marginal_posterior_plot = plot();
approx_in_latents_posterior_plot = plot();

##################################### Exact Inference ######################################

# Compute the posterior processes, sample from them, and compute marginals.
@show logpdf(y(X), ŷ);
f₁′, f₂′, f₃′ = (f₁, f₂, f₃) | (y(X) ← ŷ);
f₁′Xp, f₂′Xp, f₃′Xp = rand(rng, [f₁′(Xp), f₂′(Xp), f₃′(Xp)], S);
μ₁′, σ₁′ = marginals(f₁′(Xp));
μ₂′, σ₂′ = marginals(f₂′(Xp));
μ₃′, σ₃′ = marginals(f₃′(Xp));

items = [
    (μ₁′, σ₁′, f₁′Xp, :red, "exact", f₁′_plot),
    (μ₂′, σ₂′, f₂′Xp, :red, "", f₂′_plot),
    (μ₃′, σ₃′, f₃′Xp, :red, "", f₃′_plot),
];

# Posterior marginal variance.
for (μ, σ, _, colour, name, plt) in items
    plot!(plt, Xp, [μ, μ];
        linewidth=0.0,
        fillrange=[μ .- 3 .* σ, μ .+ 3 * σ],
        fillalpha=0.3,
        fillcolor=colour,
        label="");
end

# Posterior samples.
for (μ, σ, f, colour, name, plt) in items
    plot!(plt, Xp, f,
        linecolor=colour,
        linealpha=0.2,
        label="");
end

# Posterior mean.
for (μ, σ, f, colour, name, plt) in items
    plot!(plt, Xp, μ;
        linecolor=colour,
        linewidth=2.0,
        label=name);
end

# Plot observations.
scatter!(posterior_plot, X, ŷ;
    markercolor=:blue,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.7,
    label="");

##################### Approximate inference with pseudo-data in f₃ #########################

# Compute approximate posterior process.
M₃ = 15;
Z₃ = linspace(0.0, T, M₃);
Z₃ = 0:T
Z₃ = 0:0.5:T
# Z₃ = 0:0.25:
u₃ = f₃(Z₃);
μ′u₃, Σ′u₃ = Stheno.optimal_q(f₃(X), ŷ, u₃, σ_noise);
q_u₃ = Stheno.Titsias(u₃, μ′u₃, Σ′u₃, gpc);
f₁′u₃, f₂′u₃, f₃′u₃ = (f₁, f₂, f₃) | q_u₃;
f₁′u₃Xp, f₂′u₃Xp, f₃′u₃Xp = rand(rng, [f₁′u₃(Xp), f₂′u₃(Xp), f₃′u₃(Xp)], S);
μ₁′u₃, σ₁′u₃ = marginals(f₁′u₃(Xp));
μ₂′u₃, σ₂′u₃ = marginals(f₂′u₃(Xp));
μ₃′u₃, σ₃′u₃ = marginals(f₃′u₃(Xp));

@show elbo(f₃(X), ŷ, u₃, σ_noise);

items = [
    (μ₁′u₃, σ₁′u₃, f₁′u₃Xp, :green, "Z in f3", f₁′_plot),
    (μ₂′u₃, σ₂′u₃, f₂′u₃Xp, :green, "", f₂′_plot),
    (μ₃′u₃, σ₃′u₃, f₃′u₃Xp, :green, "", f₃′_plot),
];

# Posterior marginal variance.
for (μ, σ, _, colour, name, plt) in items
    plot!(plt, Xp, [μ, μ];
        linewidth=0.0,
        fillrange=[μ .- 3 .* σ, μ .+ 3 * σ],
        fillalpha=0.3,
        fillcolor=colour,
        label="");
end

# # Posterior samples.
# for (μ, σ, f, colour, name, plt) in items
#     plot!(plt, Xp, f,
#         linecolor=colour,
#         linealpha=0.2,
#         label="");
# end

# Posterior mean.
for (μ, σ, f, colour, name, plt) in items
    plot!(plt, Xp, μ;
        linecolor=colour,
        linewidth=2.0,
        label=name);
end

# Plot observations and pseudo-input locations.
scatter!(approx_in_marginal_posterior_plot, X, ŷ;
    markercolor=:blue,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.7,
    label="");
scatter!(approx_in_marginal_posterior_plot, Z₃, zeros(M₃);
    markercolor=:black,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.8,
    label="Z");

############# Perform approximate inference by placing pseudo-data in f₁ and f₂ ############

# Compute approximate posterior process.
M₁, M₂ = 10, 5;
Z₁, Z₂ = linspace(0.0, 1 / ω, M₁), linspace(0.0, T, M₂);
u₁₂ = BlockGP([f₁(Z₁), f₂(Z₂)]);
μ′u, Σ′u = Stheno.optimal_q(f₃(X), ŷ, u₁₂, σ_noise);
conditioner = Stheno.Titsias(u₁₂, μ′u, Σ′u, gpc);
f₁′u₁₂, f₂′u₁₂, f₃′u₁₂ = (f₁, f₂, f₃) | conditioner;

f₁′u₁₂Xp, f₂′u₁₂Xp, f₃′u₁₂Xp = rand(rng, [f₁′u₁₂(Xp), f₂′u₁₂(Xp), f₃′u₁₂(Xp)], S);
μ₁′u₁₂, σ₁′u₁₂ = marginals(f₁′u₁₂(Xp));
μ₂′u₁₂, σ₂′u₁₂ = marginals(f₂′u₁₂(Xp));
μ₃′u₁₂, σ₃′u₁₂ = marginals(f₃′u₁₂(Xp));

@show elbo(f₃(X), ŷ, u₁₂, σ_noise);

items = [
    (μ₁′u₁₂, σ₁′u₁₂, f₁′u₁₂Xp, :blue, "Z in f1 and f2", f₁′_plot),
    (μ₂′u₁₂, σ₂′u₁₂, f₂′u₁₂Xp, :blue, "", f₂′_plot),
    (μ₃′u₁₂, σ₃′u₁₂, f₃′u₁₂Xp, :blue, "", f₃′_plot),
];

# Posterior marginal variance.
for (μ, σ, _, colour, name, plt) in items
    plot!(plt, Xp, [μ, μ];
        linewidth=0.0,
        fillrange=[μ .- 3 .* σ, μ .+ 3 * σ],
        fillalpha=0.3,
        fillcolor=colour,
        label="");
end

# # Posterior samples.
# for (μ, σ, f, colour, name, plt) in items
#     plot!(plt, Xp, f,
#         linecolor=colour,
#         linealpha=0.2,
#         label="");
# end

# Posterior mean.
for (μ, σ, f, colour, name, plt) in items
    plot!(plt, Xp, μ;
        linecolor=colour,
        linewidth=2.0,
        label=name);
end

# Plot observations and pseudo-input locations.
scatter!(approx_in_latents_posterior_plot, X, ŷ;
    markercolor=:blue,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.7,
    label="");
scatter!(approx_in_latents_posterior_plot, Z₁, zeros(M₁);
    markercolor=:black,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.8,
    label="Z₁");
scatter!(approx_in_latents_posterior_plot, Z₂, zeros(M₂);
    markercolor=:magenta,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.8,
    label="Z₂");

# display(posterior_plot);
# display(approx_in_marginal_posterior_plot);
# display(approx_in_latents_posterior_plot);

plot(f₁′_plot, f₂′_plot, f₃′_plot; layout=(3, 1))


text:
module ClassicalLMG

export hamiltonian,ClassicalLMGSystem,P_of_ϵ
    import ..ClassicalSystems
    using ..PhaseSpaces
    using ParameterizedFunctions

    LMG_step = @ode_def begin
        dQ = integate_backwards*P*(Ω - ξ*Q^2/2) # =dH/dP
        dP = integate_backwards*(Q*(ξ*P^2/2 - (2*ξ + Ω)) + ξ*Q^3) # =-dH/dQ
    end Ω ξ integate_backwards
    function out_of_bounds_LMG(u,_,t)
        return 4-u[2]^2-u[1]^2<=0
    end
    varnames=[:Q,:P]
    """
    ```julia
    struct ClassicalLMGSystem <: ClassicalSystems.ClassicalSystem
    ```
    Subtype of [`ClassicalSystems.ClassicalSystem`](@ref) which represents the
    classical LMG model with the given parameters ``Ω`` and ``ξ``. See Eq. (2) of Ref. [Pilatowsky2020](@cite).
    To generate this struct, use the constructor
    ```julia
        ClassicalLMGSystem(;Ω::Real,ξ::Real)
    ```
    For example, `system = ClassicalLMGSystem(Ω=1, ξ=1)`.

    This struct may be passed to all functions in this module that require an instance of `ClassicalLMG.ClassicalLMGSystem`,
    as well as functions in other modules that require the abstract [`ClassicalSystems.ClassicalSystem`](@ref),
    such as  [`ClassicalSystems.integrate`](@ref ClassicalSystems.integrate(::ClassicalSystems.ClassicalSystem,::AbstractVector{<:Real},::Real)).
    """
    struct ClassicalLMGSystem <: ClassicalSystems.ClassicalSystem
        parameters::Vector{Float64}
        ClassicalLMGSystem(;Ω::Real,ξ::Real)= new([Float64(Ω),Float64(ξ)])

    end
    function ClassicalSystems.parameters(system::ClassicalLMGSystem)
        return system.parameters
    end
    function ClassicalSystems.step(system::ClassicalLMGSystem)
        return LMG_step
    end
    function ClassicalSystems.out_of_bounds(system::ClassicalLMGSystem)
        return out_of_bounds_LMG
    end
    function ClassicalSystems.varnames(system::ClassicalLMGSystem)
        return varnames
    end
    function Base.show(io::IO, system::ClassicalLMGSystem)
        Ω,ξ=ClassicalSystems.parameters(system)
        print(io,"ClassicalLMGSystem(Ω = $Ω, ξ = $ξ)")
    end
    """
    ```julia
    function hamiltonian(system::ClassicalLMGSystem)
    ```
    Returns a classical Hamiltonian function `h(x)` where `x=[Q,P]`, which is given by
    Eq. (2) of Ref. [Pilatowsky2020](@cite).

    # Arguments
    - `system` should be generated with [`ClassicalLMG.ClassicalLMGSystem`](@ref).
    """
    function hamiltonian(system::ClassicalLMGSystem)
        Ω,ξ=ClassicalSystems.parameters(system)
        function H(u)
            Q,P=u
            Ω*(Q^2 + P^2)/2 + ξ*(Q^2 - P^2*Q^2/4 - Q^4/4) - Ω
        end
        return H
    end
    """
    ```julia
    function discriminant_of_P_solution(system::ClassicalLMGSystem,Q::Real,ϵ::Real)
    ```
    Returns the discriminant of the second degree equation in ``P`` given by
    ```math
        h_\\text{cl}(Q,P)=\\epsilon,
    ```
    where ``h_\\text{cl}`` is given by Eq. (2) of Ref. [Pilatowsky2020](@cite).

    # Arguments
    - `system` should be generated with [`ClassicalLMGSystem`](@ref).
    - `Q` and `ϵ` are the values of ``Q`` and ``\\epsilon``, respectively.
    """
    function discriminant_of_P_solution(system::ClassicalLMGSystem,Q::Real,ϵ::Real)
        Ω,ξ=ClassicalSystems.parameters(system)
        return (-4*ϵ + 4*Q^2*ξ - Q^4*ξ - 4*Ω + 2*Q^2*Ω)/(Q^2*ξ - 2*Ω)
    end
    """
    ```julia
    function P_of_ϵ(system::ClassicalLMGSystem;
        Q::Real,
        ϵ::Real,
        sgn::Union{typeof(-),typeof(+)}=+,
        returnNaNonError::Bool=true)
    ```
    Returns the solutions ``P_\\pm`` of the second degree equation in ``P`` given by
    ```math
        h_\\text{cl}(Q,P)=\\epsilon,
    ```
    where ``h_\\text{cl}`` is given by Eq. (2) of Ref. [Pilatowsky2020](@cite).

    # Arguments
    - `system` should be generated with [`ClassicalLMGSystem`](@ref).
    # Keyword arguments
    - `Q` and `ϵ` are values of ``Q`` and ``\\epsilon``, respectively.
    - `sgn` is `+` for ``P_+`` and `-` for ``P_-``
    - If `returnNaNonError` is `true`, then `NaN` is returned if there are no solutions. If it is `false`, and error is raised.
    """
    function P_of_ϵ(system::ClassicalLMGSystem;
        Q::Real,
        ϵ::Real,
        sgn::Union{typeof(-),typeof(+)}=+,
        returnNaNonError::Bool=true)
        Δ=discriminant_of_P_solution(system,Q,ϵ)
        if Δ<0
            if returnNaNonError
                return NaN
            else
                ϵmin= minimum_ϵ_for(system;Q=Q)
                error("The minimum energy for the given parameters is $ϵmin. There is no P such that (Q=$Q,P) has ϵ=$ϵ.")
            end
        end
        return sgn(0.0,sqrt(Δ))
    end
    """
    ```julia
    function minimum_ϵ_for(system::ClassicalLMGSystem;
        Q::Union{Real,Nothing}=nothing,
        P::Union{Real,Nothing}=nothing)
    ```
    Returns the minimum energy ``\\epsilon`` when constraining the system to
    one fixed value of the coordinates ``Q`` or ``P``.

    # Arguments
    - `system` should be generated with [`ClassicalLMGSystem`](@ref).
    # Keyword arguments
    - You may pass either ``P`` or ``Q``.
    """
    function minimum_ϵ_for(system::ClassicalLMGSystem;
        Q::Union{Real,Nothing}=nothing,
        P::Union{Real,Nothing}=nothing)
        if count(x->x==nothing,[Q,P])!=1
            error("You have to pass one of Q or P")
        end
        Ω,ξ=ClassicalSystems.parameters(system)
        if P==nothing
            return (2*Ω*Q^2 + 4*Q^2*ξ - Q^4*ξ)/4 - Ω
        else
            return Ω*P^2/2 - Ω
        end
    end

    """
    ```julia
    function Point(;Q::Real,P::Real)
    ```
    Returns the list `[Q,P]`
    """
    Point(;Q::Real,P::Real)=Float64[Q,P]
    """
    ```julia
    function Pointθϕ(;θ::Real,ϕ::Real)
    ```
    Returns a list `[Q,P]`, where `Q` and `P` are calculated from `θ` and `ϕ` using [`PhaseSpaces.Q_of_θϕ`](@ref) and  [`PhaseSpaces.P_of_θϕ`](@ref).
    """
    Pointθϕ(;θ::Real,ϕ::Real)=Point(Q=Q_of_θφ(θ,ϕ),P=P_of_θφ(θ,ϕ))
    Pointθφ(;θ::Real,φ::Real)=Pointθϕ(θ=θ,ϕ=φ)
    """
    ```julia
    function Point(system::ClassicalLMGSystem;
        Q::Real,
        ϵ::Real,
        sgn::Union{typeof(-),typeof(+)} = +)
    ```
    Returns a list `[Q,P]`, where `P` is calculated with [`P_of_ϵ`](@ref).
    If there are no solutions for ``P``, an error is raised.
    """
    Point(system::ClassicalLMGSystem;
        Q::Real,
        ϵ::Real,
        sgn::Union{typeof(-),typeof(+)}=+)=Point(Q=Q,P=P_of_ϵ(system;Q=Q,ϵ=ϵ,sgn=sgn,returnNaNonError=false))
end


text:
using Mamba

## Model and User-Defined Sampler Specifications

model = Model(

  y = Stochastic(1,
    (mu, s2) ->  MvNormal(mu, sqrt(s2)),
    false
  ),

  mu = Logical(1,
    (xmat, beta) -> xmat * beta,
    false
  ),

  beta = Stochastic(1,
    () -> MvNormal(2, sqrt(1000))
  ),

  s2 = Stochastic(
    () -> InverseGamma(0.001, 0.001)
  )

)

Gibbs_beta = Sampler([:beta],
  (beta, s2, xmat, y) ->
    begin
      beta_mean = mean(beta.distr)
      beta_invcov = invcov(beta.distr)
      Sigma = inv(Symmetric(xmat' * xmat / s2 + beta_invcov))
      mu = Sigma * (xmat' * y / s2 + beta_invcov * beta_mean)
      rand(MvNormal(mu, Sigma))
    end
)

Gibbs_s2 = Sampler([:s2],
  (mu, s2, y) ->
    begin
      a = length(y) / 2.0 + shape(s2.distr)
      b = sumabs2(y - mu) / 2.0 + scale(s2.distr)
      rand(InverseGamma(a, b))
    end
)

## Hybrid No-U-Turn and Slice Sampling Scheme
scheme1 = [NUTS(:beta),
           Slice(:s2, 3.0)]

## No-U-Turn Sampling Scheme
scheme2 = [NUTS([:beta, :s2])]

## User-Defined Sampling Scheme
scheme3 = [Gibbs_beta, Gibbs_s2]

## Sampling Scheme Assignment
setsamplers!(model, scheme1)

## Graph Representation of Nodes
draw(model)
draw(model, filename="lineDAG.dot")

## Data
line = Dict{Symbol, Any}(
  :x => [1, 2, 3, 4, 5],
  :y => [1, 3, 3, 3, 5]
)
line[:xmat] = [ones(5) line[:x]]

## Set Random Number Generator Seed
srand(123)

## Initial Values
inits = [
  Dict{Symbol, Any}(
    :y => line[:y],
    :beta => rand(Normal(0, 1), 2),
    :s2 => rand(Gamma(1, 1))
  )
  for i in 1:3
]

## MCMC Simulations

setsamplers!(model, scheme1)
sim1 = mcmc(model, line, inits, 10000, burnin=250, thin=2, chains=3)

setsamplers!(model, scheme2)
sim2 = mcmc(model, line, inits, 10000, burnin=250, thin=2, chains=3)

setsamplers!(model, scheme3)
sim3 = mcmc(model, line, inits, 10000, burnin=250, thin=2, chains=3)

## Gelman, Rubin, and Brooks Convergence Diagnostic
gelmandiag(sim1, mpsrf=true, transform=true) |> showall

## Geweke Convergence Diagnostic
gewekediag(sim1) |> showall

## Heidelberger-Welch Diagnostic
heideldiag(sim1) |> showall

## Raftery-Lewis Convergence Diagnostic
rafterydiag(sim1) |> showall

## Summary Statistics
describe(sim1)

## Highest Posterior Density Intervals
hpd(sim1) |> show

## Cross-Correlations
cor(sim1) |> show

## Lag-Autocorrelations
autocor(sim1) |> show

## State Space Change Rate (per Iteration)
changerate(sim1) |> show

## Deviance Information Criterion
dic(sim1) |> show

## Subset Sampler Output
sim = sim1[1000:5000, ["beta[1]", "beta[2]"], :]
describe(sim)

## Write to and Read from an External File
write("sim1.jls", sim1)
sim1 = read("sim1.jls", ModelChains)

## Restart the Sampler
sim = mcmc(sim1, 5000)
describe(sim)

## Plotting

## Default summary plot (trace and density plots)
p = plot(sim1)

## Write plot to file
draw(p, filename="summaryplot.svg")
#draw(p, filename="summaryplot.pdf", fmt=:pdf)

## Autocorrelation and running mean plots
p = plot(sim1, [:autocor, :mean], legend=true)
draw(p, nrow=3, ncol=2, filename="autocormeanplot.svg")
#draw(p, nrow=3, ncol=2, filename="autocormeanplot.pdf", fmt=:pdf)

## Pairwise contour plots
p = plot(sim1, :contour)
draw(p, nrow=2, ncol=2, filename="contourplot.svg")
#draw(p, nrow=2, ncol=2, filename="contourplot.pdf", fmt=:pdf)

## Development and Testing

setinputs!(model, line)             # Set input node values
setinits!(model, inits[1])          # Set initial values
setsamplers!(model, scheme1)        # Set sampling scheme

showall(model)                      # Show detailed node information

logpdf(model, 1)                    # Log-density sum for block 1
logpdf(model, 2)                    # Block 2
logpdf(model)                       # All blocks

sample!(model, 1)                  # Sample values for block 1
sample!(model, 2)                  # Block 2
sample!(model)                     # All blocks


text:
using Base.Test

using Graft

import Graft: eltypes, randindxs, completeindxs, LabelMap, IdentityLM, DictLM,
bfs, bfs_list, bfs_tree, bfs_subgraph

###
# TEST FILES
###
include("util.jl")

include("sparsematrix.jl")

include("edgeiter.jl")

include("labelmap.jl")

include("graph.jl")

include("generator.jl")

include("combinatorial.jl")

include("vdata.jl")

include("edata.jl")

include("mutation.jl")

include("display.jl")

include("graphio.jl")

include("subgraph.jl")

include("algorithms.jl")

include("operations.jl")

include("query.jl")


text:
# develop cascade tree init

# using Caesar, RoME, Images

function csmHexDelayTest(injectDelayBefore)

  fg = generateCanonicalFG_Hexagonal(graphinit=false)

  getSolverParams(fg).treeinit = true
  getSolverParams(fg).graphinit = false
  getSolverParams(fg).limititers = 100
  getSolverParams(fg).async = true

  # limitcliqs = [:x0=>8;:x4=>12;:l1=>21;:x1=>21;:x5=>50;:x3=>50] # breaks
  # limitcliqs = [:x0=>8;:x4=>13;:l1=>21;:x1=>21;:x5=>60;:x3=>60] # 50 # doesnt break, blocks

  # injectDelayBefore=[5=>(canCliqMargRecycle_StateMachine=>5); ]
  # injectDelayBefore = nothing

  mkpath(getLogPath(fg))

  fid = open(joinLogPath(fg,"delays.dat"), "w")
  println(fid, injectDelayBefore)
  close(fid)

  verbosefid = open(joinLogPath(fg, "csmVerbose.log"),"w")
  # verbosefid = stdout
  tree = solveTree!(fg, recordcliqs=ls(fg), timeout=40, verbose=true, verbosefid=verbosefid, injectDelayBefore=injectDelayBefore); #, limititercliqs=limitcliqs);

  ## wait for async to complete

  println("Going to sleep an long period so that solve can finish or fail")
  sleep(180)

  flush(verbosefid)
  close(verbosefid)

  didFail = 0 < ((smt .|> x->x.state == :failed) |> sum)

  fid = open(joinLogPath(fg,"didFail.log"), "w")
  println(fid, "did the solve fail, $didFail")
  close(fid)

  # async case
  fetchCliqHistoryAll!(smt, hists)

  fid = open(joinLogPath(fg, "csmSequ.log"),"w")
  printCliqHistorySequential(hists, nothing, fid)
  close(fid)

  fid = open(joinLogPath(fg, "csmLogi.log"),"w")
  printCSMHistoryLogical(hists, fid)
  close(fid)

  return fg, didFail
end

#


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _geo2mag{T<:AbstractFloat}(lat::T, long::T, pole_lat::T, pole_long::T)
    r       = 1.0 # Distance from planet center.  Value unimportant -- just need
                  # a length for conversion to rectangular coordinates
    lat  = deg2rad(lat)
    long = deg2rad(long)
    alt  = lat*0.0 + r
    x = alt*cos(lat)*cos(long)
    y = alt*cos(lat)*sin(long)
    z = alt*sin(lat)

    # Compute first rotation matrix: rotation around plane of the equator, from
    # the Greenwich meridian to the meridian containing the magnetic dipole
    # pole.
    geolong2maglong = Array(T, 3, 3)
    geolong2maglong[:,1] = [cos(pole_long), -sin(pole_long), 0.0]
    geolong2maglong[:,2] = [sin(pole_long),  cos(pole_long), 0.0]
    geolong2maglong[:,3] = [           0.0,             0.0, 1.0]
    out = geolong2maglong * [x, y, z]

    # Second rotation: in the plane of the current meridian from geographic pole
    # to magnetic dipole pole.
    tomaglat = Array(T, 3, 3)
    tomaglat[:,1] = [ cos(pi/2 - pole_lat), 0.0, sin(pi/2 - pole_lat)]
    tomaglat[:,2] = [                  0.0, 1.0,                  0.0]
    tomaglat[:,3] = [-sin(pi/2 - pole_lat), 0.0, cos(pi/2 - pole_lat)]
    out = tomaglat * out

    maglat  = rad2deg(atan2(out[3], hypot(out[1], out[2])))
    maglong = rad2deg(atan2(out[2], out[1]))
    # I don't care about that one...just put it there for completeness' sake
    # magalt  = vecnorm(out) - r
    return maglat, maglong
end

"""
    geo2mag(latitude, longitude[, year]) -> geomagnetic_latitude, geomagnetic_longitude

### Purpose ###

Convert from geographic to geomagnetic coordinates.

### Explanation ###

Converts from geographic (latitude, longitude) to geomagnetic (latitude,
longitude).  Altitude is not involved in this function.

### Arguments ###

* `latitude`: geographic latitude (North), in degrees.
* `longitude`: geographic longitude (East), in degrees.
* `year` (optional numerical argument): the year in which to perform conversion.
  If omitted, defaults to current year.

The coordinates can be passed as arrays of the same length.

### Output ###

The 2-tuple of magnetic (latitude, longitude) coordinates, in degrees.

If geographical coordinates are given as arrays, a 2-tuple of arrays of the same
length is returned.

### Example ###

Kyoto has geographic coordinates 35° 00' 42'' N, 135° 46' 06'' E, find its
geomagnetic coordinates in 2016:

``` julia
geo2mag(ten(35,0,42), ten(135,46,6), 2016)
# => (36.86579228937769,-60.184060536651614)
```

### Notes ###

This function uses list of North Magnetic Pole positions provided by World
Magnetic Model (https://www.ngdc.noaa.gov/geomag/data/poles/NP.xy).

`mag2geo` converts geomagnetical coordinates to geographic coordinates.

Code of this function is based on IDL Astronomy User's Library.
"""
geo2mag(lat::Real, long::Real, year::Real=Dates.year(Dates.now())) =
    _geo2mag(promote(float(lat), float(long),
                     deg2rad(POLELATLONG[year][1]::AbstractFloat),
                     deg2rad(POLELATLONG[year][2]::AbstractFloat))...)

function geo2mag{LA<:Real, LO<:Real}(lat::AbstractArray{LA},
                                     long::AbstractArray{LO},
                                     year::Real=Dates.year(Dates.now()))
    @assert length(lat) == length(long)
    typela   = typeof(float(one(LA)))
    maglat   = similar(lat, typela)
    maglong  = similar(lat, typela)
    polelat  = deg2rad(POLELATLONG[year][1])
    polelong = deg2rad(POLELATLONG[year][2])
    for i in eachindex(lat)
        maglat[i], maglong[i] =
            _geo2mag(promote(float(lat[i]), float(long[i]),
                             polelat, polelong)...)
    end
    return maglat, maglong
end


text:
module TestDensity
using Test
using LinearAlgebra
using GaussianDistributions
using GaussianDistributions: pushforward, nlml
using ForwardDiff
using CovarianceFunctions

@testset "density" begin
    @testset "univariate" begin
        μ, σ = randn(), exp(randn())
        G = Gaussian(μ, σ^2)
        @test nld(G, μ) ≈ log(sqrt(2π)*σ)
        @test exp(-nld(G, μ)) ≈ pdf(G, μ)
        x = randn()
        @test exp(-nld(G, x)) ≈ pdf(G, x)
        @test pdf(G)(x) ≈ pdf(G, x)
        @test nld(G)(x) ≈ nld(G, x)
        @test cdf(G)(μ) ≈ 1/2
        @test cdf(G)(Inf) == 1
        @test cdf(G)(-Inf) == 0
    end

    @testset "multivariate" begin
        d = 3
        μ = randn(d)
        Σ = randn(d, d)
        Σ = Σ'Σ
        G = Gaussian(μ, Σ)
        @test nld(G, μ) ≈ logdet(2π*Σ) / 2
        @test exp(-nld(G, μ)) ≈ pdf(G, μ)
        x = randn(d)
        @test exp(-nld(G, x)) ≈ pdf(G, x)
    end

    @testset "process" begin
        μ(x) = zero(x)
        Σ(θ) = (x, y) -> exp(-(x-y)^2/2exp(2θ))
        G = Gaussian(μ, Σ(0.))
        σ = .1
        n = 3
        x = randn(n)
        y = randn(n)
        @test nlml(G, x, y, σ^2) isa Real
        K = cov(marginal(G, x)) + σ^2*I
        @test nlml(G, x, zero(y), σ^2) ≈ (logdet(K) + length(x) * log(2π)) / 2

        # test gradient
        dΣ(θ) = (x, y) -> ForwardDiff.derivative(z->Σ(z)(x, y), θ)
        θ = 1.
        dk = dΣ(θ)
        dK = CovarianceFunctions.gramian(dk, x)
        val, push = pushforward(nlml, K, y)
        @test val ≈ nlml(K, y)
        @test push(dK) isa Real
    end
end

end


text:
"""
Extension of product_distribution from Distributions.jl to allow for taking products between arbitrary distributions, not just univariate.
"""
struct GeneralProduct <: ContinuousMultivariateDistribution
    v::AbstractVector{<:Union{ContinuousMultivariateDistribution, ContinuousUnivariateDistribution}}
    function GeneralProduct(v::AbstractVector{<:Union{ContinuousMultivariateDistribution, ContinuousUnivariateDistribution}})
        if all(x->x isa AbstractMvNormal, v)
            MvNormal(mapreduce(mean, vcat, v), HealthyCovariance(reduce(BlockMatrix, map(HealthyCovariance∘cov, v))))
        else
            new(v)
        end
    end
end
GeneralProduct(DM::Union{AbstractDataModel,AbstractDataSet}) = GeneralProduct([xdist(DM),ydist(DM)])

import Base.length
length(P::GeneralProduct) = sum(length(P.v[i]) for i in 1:length(P.v))

import Distributions: insupport, mean, cov, invcov, pdf, logpdf, gradlogpdf, product_distribution
# insupport(P::GeneralProduct, X::AbstractVector)::Bool = all([insupport(P.v[i], X[P.lengths[i],]) for i in 1:length(P.lengths)])
# sum(!insupport(P.v[i],X[i]) for i in 1:length(P.v)) == 0
mean(P::GeneralProduct) = reduce(vcat, map(GetMean, P.v))
cov(P::GeneralProduct) = BlockMatrix(map(Sigma, P.v)...)

for F = (:logpdf, :gradlogpdf)
    eval(quote
        # Base.$F(a::MyNumber) = MyNumber($F(a.x))
        $F(P::GeneralProduct, X::AbstractVector) = sum($F(P.v[i], X[i]) for i in 1:length(P.v))
        $F(P::GeneralProduct, X::AbstractVector{<:Real}) = $F(P, X, Val(length(P.v)))
        $F(P::GeneralProduct, X::AbstractVector{<:Real}, ::Val{1}) = $F(P.v[1], X)
        $F(P::GeneralProduct, X::AbstractVector{<:Real}, ::Val{2}) = $F(P.v[1], view(X,1:length(P.v[1]))) + $F(P.v[2], view(X,length(P.v[1])+1:lastindex(X)))
        # There has to be a more performant way than this!
        function $F(P::GeneralProduct, X::AbstractVector{<:Real}, ::Val)
            C = vcat([0],cumsum(length.(P.v)))
            sum($F(P.v[i], X[C[i]+1:C[i+1]]) for i in 1:length(P.v))
        end
    end)
end

pdf(P::GeneralProduct, X::AbstractVector) = exp(logpdf(P,X))
invcov(P::GeneralProduct) = BlockMatrix(map(InvCov,P.v)...)
product_distribution(X::AbstractVector{<:ContinuousMultivariateDistribution}) = GeneralProduct(X)

Sigma(P::GeneralProduct) = cov(P)

# LogLike(P::GeneralProduct, args...) = logpdf(P,[args...])

# Define as subtype of continuous distribution to get accepted by methods more seamlessly
# although it is actually discontinuous.
"""
pdf at μ takes value 1 and 0 everywhere else.
"""
struct Dirac <: ContinuousMultivariateDistribution
    μ::AbstractVector{<:Number}
    Dirac(μ) = new(floatify(Unwind(μ)))
end

import Base.length
length(d::InformationGeometry.Dirac) = length(d.μ)

import Distributions: insupport, mean, cov, invcov, pdf, logpdf
insupport(d::InformationGeometry.Dirac, x::AbstractVector) = length(d) == length(x) && all(isfinite, x)
mean(d::InformationGeometry.Dirac) = d.μ
cov(d::InformationGeometry.Dirac) = Diagonal(zeros(length(d)))
invcov(d::InformationGeometry.Dirac) = Diagonal([Inf for i in 1:length(d)])
pdf(d::InformationGeometry.Dirac, x::AbstractVector{<:Number}) = x == mean(d) ? 1.0 : 0.0
logpdf(d::InformationGeometry.Dirac, x::AbstractVector{<:Number}) = log(pdf(d, x))

# Fix gradlogpdf for Cauchy distribution and product distributions in general
import Distributions: gradlogpdf
gradlogpdf(P::Cauchy,x::Real) = gradlogpdf(TDist(1), (x - P.μ) / P.σ) / P.σ
gradlogpdf(P::Product,x::AbstractVector) = [gradlogpdf(P.v[i],x[i]) for i in 1:length(x)]


text:
@testset "rcParams" begin
    @testset "rcParams" begin
        @test rcParams isa ArviZ.RcParams
        @test pyisinstance(PyObject(rcParams), ArviZ.arviz.rcparams.RcParams)
        pyrcParams = ArviZ.arviz.rcParams
        @test rcParams == pyrcParams
        @test ArviZ.RcParams(pyrcParams) isa ArviZ.RcParams{Any,Any}
        @test isa(
            convert(ArviZ.RcParams{String,Union{Int64,String}}, pyrcParams),
            ArviZ.RcParams{String,Union{Int64,String}},
        )
        @test convert(ArviZ.RcParams, pyrcParams) isa ArviZ.RcParams
        @test haskey(rcParams, "plot.backend")
        def_backend = rcParams["plot.backend"]
        @test ("plot.backend" => def_backend) ∈ rcParams
        rcParams["plot.backend"] = "matplotlib"
        @test rcParams["plot.backend"] == "matplotlib"
        rcParams["plot.backend"] = "bokeh"
        @test rcParams["plot.backend"] == "bokeh"
        @test_throws KeyError rcParams["blah"]
        @test_throws KeyError rcParams["blah"] = 3
        @test_throws ErrorException rcParams["plot.backend"] = "blah"
        @test get(rcParams, "blah", "def") == "def"
        @test Dict(map(p -> Pair(p...), zip(keys(rcParams), values(rcParams)))) == rcParams
        rcParams["plot.backend"] = def_backend
    end

    @testset "defaults" begin
        @test rcParams["data.index_origin"] == 1
        @test rcParams["plot.bokeh.show"] == false
    end

    @testset "with_rc_context" begin
        def_backend = rcParams["plot.backend"]
        rcParams["plot.backend"] = "matplotlib"
        with_rc_context(; rc=Dict("plot.backend" => "bokeh")) do
            @test rcParams["plot.backend"] == "bokeh"
            return nothing
        end
        rcParams["plot.backend"] = def_backend
    end
end


text:
#### Number theory and combinatorics.

import Primes
import Combinatorics
import Polynomials

### BinarySearch

@mkapprule BinarySearch

@sjdoc BinarySearch """
    BinarySearch(list,item)

returns a list of the indices at which `item` appears in the sorted list `list`.
"""

@doap BinarySearch(a::ListT,x) = tolist(searchsorted(margs(a),x))

## TODO: We have a julia package that is much more efficient
## at functions related to primes and factoring. Find a way to make an optional Symata package.

### FactorInteger

@sjdoc FactorInteger """
    FactorInteger(n)

give a list of prime factors of `n` and their multiplicities.
"""
@mkapprule FactorInteger nargs => 1:2
@doap FactorInteger(x) = tolistoflistsfixed(Primes.factor(x))

### PrimeList

@sjdoc PrimeList """
    PrimeList(n)

return the list of all prime numbers `<= n`.
"""

## Convert `n` to Int to work around julia bug
@mkapprule PrimeList nargs => 1
@doap PrimeList(n::Integer) = tolistfixed(Primes.primes(convert(Int,n)))

### Prime

@sjdoc Prime """
    Prime(n)

returns the `n`th prime number.
"""
@mkapprule Prime :nargs => 1
@doap Prime(n::Integer) = sympy.prime(n)

### PrimePi

@mkapprule PrimePi nargs => 1  """
    PrimePi(n)

returns the number of primes less than or equal to `n`.
"""

@doap PrimePi(n::Integer) = sympy.primepi(n)

### Permutations

@mkapprule Permutations nargs => 1:2  """
    Permutations(expr)

give a list of all permutations of elements in `expr`.
"""

@doap Permutations(x::Mxpr) = tolistoflistsfixed(Combinatorics.permutations(margs(x)))

### IntegerPartitions

## The Mathematica "way" is not very convenient. It is clumsy to controll what you
## want (regardless of the implmentation here, it is clumsy for the user.) But, it
## is not clear how to avoid this kind of interface to julia `partitions`

## Setting "fixed" below increases speed by 3 or 5 times in some tests.

@mkapprule IntegerPartitions
@doap IntegerPartitions(n::Integer) = tolistoflistsfixed(Combinatorics.partitions(n))
@doap IntegerPartitions(n::Integer, m::Integer) = _integer_partitions_range(n,1,m)

function _integer_partitions_range(n,m1,m2)
    a = collect(Combinatorics.partitions(n,m1))
    for i in (m1+1):m2
        c = collect(Combinatorics.partitions(n,i))
        isempty(c) && continue
        append!(a,c)
    end
    tolistoflists(a)
end

@doap IntegerPartitions(n::Integer, ls::ListT) = listofintegersq(ls) ?  _integer_partitions_list(mx,n,margs(ls)...) : mx
_integer_partitions_list(mx,n,m) = tolistoflistsfixed(collect(Combinatorics.partitions(n,m)))
_integer_partitions_list(mx,n,m1,m2) = _integer_partitions_range(n,m1,m2)
_integer_partitions_list(mx,n,m,args...) = mx

@mkapprule NumberOfPartitions nargs => 1  """
    NumberOfPartitions(n)

gives the number of integer partitions of `n`.
"""

@doap NumberOfPartitions(n::Integer) = n <= 405 ? length(Combinatorics.partitions(n)) : bnpartitions(n)

## copied from partitions.jl. We just add the BigInt dict
## If n > 405, we need BigInt.
## We intentionally eschew type stability.
bnpartitions(n::Integer) = _bnpartitions_func(convert(Int,n))

let _bnpartitions = Dict{Int,BigInt}()
    global _bnpartitions_func
    function _bnpartitions_func(n::Int)
        if n < 0
            0
        elseif n < 2
            1
        elseif (np = get(_bnpartitions, n, 0)) > 0
            np
        else
            np = 0
            sgn = 1
            for k = 1:n
                np += sgn * (bnpartitions(n - (k * (3k - 1)) >> 1) + bnpartitions(n - (k * (3k + 1)) >> 1))
                sgn = -sgn
            end
            _bnpartitions[convert(Int,n)] = np
        end
    end
end

### Named polynomial sequences.

## @runtime_include is defined in util.jl

@runtime_include "PolynomialSequences.jl" _init_polys

### Fibonaccci

@mkapprule Fibbonaci nargs => 1:2

@doap function Fibbonaci(n::Integer)
    n < 0 && return mx
    Combinatorics.fibonaccinum(n)
end

# The choice of n <= 92 is arbitrary. This is where evaluation with x = 1 causes Int overflow
@doap function Fibbonaci(n::Integer,x)
    n < 0 && return mx
    _init_polys() || return mx
    p =  n <= 92 ? fibpoly(n) :
          fibpoly(big(n))
    fromjuliaPolynomial(p,x)
end

### LucasL

@mkapprule LucasL nargs => 1:2

@doap function LucasL(n::Integer)
    n < 0 && return mx
    Combinatorics.lucasnum(n)
end

@doap function LucasL(n::Integer,x)
    n < 0 && return mx
    _init_polys() || return mx
    p =  n <= 90 ? lucaspoly(n) :
          lucaspoly(big(n))
    fromjuliaPolynomial(p,x)
end

### Interface to Julia Polynomials.jl

## TODO: why do we have to do mergesyms here, but not when "var" is a Mxpr, as below
## ... a method is missing in evaluation.jl ?
function fromjuliaPolynomial(p::Polynomials.Poly, var::Symbol)
    res = _fromjuliaPolynomial(p,var)
    mergesyms(res,var)  # force evaluation when var changes.
    if isa(res,Mxpr)
        a = margs(res)
        foreach(setfixed, a)
    end
    setfixed(res)
end

## symbols are merged correctly automatically.
## TODO: It is probably faster to merge symbols here
function fromjuliaPolynomial(p::Polynomials.Poly, var)
    _fromjuliaPolynomial(p,var)
end

function _fromjuliaPolynomial(p::Polynomials.Poly, var)
    cs = Polynomials.coeffs(p)
    nc = length(cs)
    nargs = newargs()
    if nc > 0
        if cs[1] != 0
            push!(nargs,cs[1])
        end
    end
    if nc > 1
        if cs[2] > 1
            push!(nargs, mmul(cs[2], var))
        elseif cs[2] == 1
            push!(nargs,var)
        end
    end
    for i in 3:length(cs)
        c = cs[i]
        c == 0 && continue
        n = i - 1
        push!(nargs, _monomial(var, c, n))
    end
    mxpra(:Plus,nargs)
end

## assumes n != 0,  n != 1
function _monomial(var,c,n)
    if c == 1
        mpow(var,n)
    else
        mmul(c, mpow(var, n))
    end
end


text:
## Module associated with the MeshArrays package

module MeshArrays

import Pkg
thistoml=joinpath(dirname(pathof(MeshArrays)), "..", "Project.toml")
thisversion=Pkg.TOML.parsefile(thistoml)["version"]

using LazyArtifacts
p=dirname(pathof(MeshArrays))
artifact_toml = joinpath(p, "../Artifacts.toml")

GRID_LLC90_hash = artifact_hash("GRID_LLC90", artifact_toml)
GRID_LLC90 = joinpath(artifact_path(GRID_LLC90_hash)*"/","GRID_LLC90-1.1/")
GRID_LLC90_download() = artifact"GRID_LLC90"
GRID_LL360_hash = artifact_hash("GRID_LL360", artifact_toml)
GRID_LL360 = joinpath(artifact_path(GRID_LL360_hash)*"/","GRID_LL360-1.0/")
GRID_LL360_download() = artifact"GRID_LL360"
GRID_CS32_hash = artifact_hash("GRID_CS32", artifact_toml)
GRID_CS32 = joinpath(artifact_path(GRID_CS32_hash)*"/","GRID_CS32-1.1/")
GRID_CS32_download() = artifact"GRID_CS32"

include("Types.jl")
include("Grids.jl")
include("GridPaths.jl")
include("Operations.jl")
include("Exchanges.jl")
include("ReadWrite.jl")
include("Solvers.jl")
include("ReIndexing.jl")
include("Interpolation.jl")
include("VerticalDimension.jl")

export AbstractMeshArray, MeshArray, gcmgrid, varmeta
export GridSpec, GridLoad, GridLoadVar, UnitGrid, simple_periodic_domain
export exchange, Tiles, Interpolate, InterpolationFactors, knn
#The following exch_UV differs from normal exchange; incl. exch_UV_N
export exch_UV

#export InnerArray, OuterArray
#export smooth, mask
#export GridOfOnes
#export GridAddWS!

export ScalarPotential, VectorPotential
export UVtoUEVN, UVtoTransport, gradient, curl, convergence
export LatitudeCircles, Transect, ThroughFlow
export StereographicProjection, isosurface

#export location_is_out, NeighborTileIndices_dpdo, NeighborTileIndices_cs, RelocationFunctions_cs
#export update_location_cs!, update_location_llc!, update_location_dpdo!

end # module