text:
module Firestore

using Dates
using JSON3
using HTTP
using StructTypes: StructTypes

export fval, Fields

function __init__()
    if haskey(ENV, "FIRESTORE_PROJECT")
        set_project!(ENV["FIRESTORE_PROJECT"])
    else
        @warn "No FIRESTORE_PROJECT key found in ENV.  `set_project!` must be called before this package will work properly."
    end
    if haskey(ENV, "FIRESTORE_API_KEY")
        set_api_key!(ENV["FIRESTORE_API_KEY"])
    else
        @warn "No FIRESTORE_API_KEY key found in ENV.  `set_api_key!` must be called before authentication will work."
    end
end

project = Ref{String}()
api_key = Ref{String}()
token = Ref{String}()
db = Ref{String}("(default)")

set_project!(s::String) = (project[] = s)
set_api_key!(s::String) = (api_key[] = s)
set_token!(s::String) = (token[] = s)
set_db!(s::String) = (db[] = s)

abstract type FirestoreType end
StructTypes.StructType(::Type{<:FirestoreType}) = StructTypes.Struct()

fval(x) = error("No Julia-to-Firestore value mapping for $(typeof(x)).")

#-----------------------------------------------------------------------------# Fields/Values
struct Fields <: FirestoreType
    fields::Dict{String, FirestoreType}
    Fields(x::AbstractDict) = new(Dict(string(k) => fval(v) for (k,v) in pairs(x)))
end
Fields(; kw...) = Fields(Dict(kw))

struct Values <: FirestoreType
    values::Vector{FirestoreType}
end
Values(args...) = Values(collect(args))

#-----------------------------------------------------------------------------# ArrayValue
struct ArrayValue <: FirestoreType
    arrayValue::Values
end
fval(x::AbstractVector) = ArrayValue(Values(collect(fval.(x))))

#-----------------------------------------------------------------------------# BooleanValue
struct BooleanValue <: FirestoreType
    booleanValue::Bool
end
fval(x::Bool) = BooleanValue(x)

#-----------------------------------------------------------------------------# DoubleValue
struct DoubleValue <: FirestoreType
    doubleValue::Float64
end
fval(x::AbstractFloat) = DoubleValue(Float64(x))

#-----------------------------------------------------------------------------# IntegerValue
struct IntegerValue <: FirestoreType
    integerValue::String
end
fval(x::Int) = IntegerValue(string(x))

#-----------------------------------------------------------------------------# MapValue
struct MapValue <: FirestoreType
    mapValue::Fields
end
fval(x::AbstractDict) = MapValue(Fields(x))
fval(x::NamedTuple) = fval(Dict(pairs(x)))

#-----------------------------------------------------------------------------# NullValue
struct NullValue <: FirestoreType
    nullValue::Nothing
    NullValue() = new(nothing)
end
fval(x::Nothing) = NullValue()
fval(x::Missing) = NullValue()

#-----------------------------------------------------------------------------# StringValue
struct StringValue <: FirestoreType
    stringValue::String
end
fval(x::AbstractString) = StringValue(String(x))
fval(x::Symbol) = fval(string(x))

#-----------------------------------------------------------------------------# TimestampValue
"A timestamp with microsecond precision."
struct TimestampValue <: FirestoreType
    timestampValue::String
end
fval(x::Dates.TimeType) = TimestampValue(string(DateTime(x)) * "Z")

#-----------------------------------------------------------------------------# BytesValue
"A Base64-encoded string."
struct BytesValue <: FirestoreType
    bytesValue::String
end

#-----------------------------------------------------------------------------# ReferenceValue
"""
A Reference to another Firestore document. E.g.

    projects/{project_id}/databases/{databaseId}/documents/{document_path}
"""
struct ReferenceValue <: FirestoreType
    referenceValue::String
end

#-----------------------------------------------------------------------------# GeoPointValue
struct LatLng <: FirestoreType
    latitude::Float64
    longitude::Float64
end

"A geo point value representing a point on the surface of Earth."
struct GeoPointValue <: FirestoreType
    geoPointValue::LatLng
end

fval(x::LatLng) = GeoPointValue(x)

#-----------------------------------------------------------------------------# Auth
function get_token!()
    email = ENV["FIRESTORE_EMAIL"]
    password = ENV["FIRESTORE_PASSWORD"]
    url = "https://identitytoolkit.googleapis.com/v1/accounts:signInWithPassword?key=$(api_key[])"
	data = JSON3.write(Dict(:email => email, :password => password, :returnSecureToken => true))
	res = HTTP.post(url, ["Content-Type: application/json"], data)
	set_token!(JSON3.read(res.body).idToken)
end

#-----------------------------------------------------------------------------# API
# API methods (https://firebase.google.com/docs/firestore/reference/rest/)
function _url(path::String)
    "https://firestore.googleapis.com/v1beta1/projects/$(project[])/databases/$(db[])/documents/$path"
end
_data(x) = JSON3.write(Fields(x))

function createDocument(path::String, doc::AbstractDict)
    split_path = split(path, '/'; keepempty=false)
    documentId = split_path[end]
    h = ["Authorization" => "Bearer $(token[])", "Content-Type: application/json"]
    HTTP.post(
        _url(join(split_path[1:end-1], '/') * "?documentId=$documentId"),
        h,
        _data(doc)
    )
end

delete(path::String) = HTTP.delete(_url(path), ["Authorization" => "Bearer $(token[])"])

function get(path)
    res = HTTP.get(_url(path))
    fields = JSON3.read(res.body).fields
    Dict(k => f2j(v) for (k,v) in pairs(fields))
end

function patch(path::String, doc::AbstractDict)
    h = ["Authorization" => "Bearer $(token[])", "Content-Type: application/json"]
    HTTP.patch(_url(path), h, _data(doc))
end

#-----------------------------------------------------------------------------# Firestore-to-Julia
function f2j(x::JSON3.Object)
    prop = propertynames(x)[1]
    val = getproperty(x, prop)
    if prop === :arrayValue
        firestore2julia.(val.values)
    elseif prop === :mapValue
        Dict(field => firestore2julia(val.fields[field]) for field in propertynames(val.fields))
    elseif prop === :integerValue
        parse(Int, val)
    elseif prop === :nullValue
        nothing
    elseif prop === :timestampValue
        DateTime(val[1:end-1])
    else
        val
    end
end

end # module


text:
using ConfParser
using Test

conf = ConfParse(joinpath(@__DIR__, "confs", "config.simple"))
parse_conf!(conf)

@test retrieve(conf, "protocol") == "kreberos"
@test retrieve(conf, "port")     == "6643"
@test erase!(conf, "protocol")
@test haskey(conf, "port") == true
@test haskey(conf, "nothere") == false
save!(conf, joinpath(@__DIR__, "confs", "out.conf"))

conf = ConfParse(joinpath(@__DIR__, "confs", "out.conf"))
parse_conf!(conf)
@test retrieve(conf, "port")     == "6643"

ext: jl


text:
using Ipopt
using PowerModelsDistribution
const PMD = PowerModelsDistribution
# documentation https://lanl-ansi.github.io/PowerModelsDistribution.jl

eng = parse_file("./test/data/opendss/case3_unbalanced.dss")
math = parse_file("./test/data/opendss/case3_unbalanced.dss"; data_model=MATHEMATICAL)

ipopt_solver = Ipopt.Optimizer
result = run_mc_opf(eng, ACPPowerModel, ipopt_solver)

result = run_mc_opf(math, ACPPowerModel, ipopt_solver)


text:
# Checklists which prevent already-ran
# experiments from being re-ran. Useful when
# ComputeCanada experiments timeout

function mark_config(self::ConfigManager)
    chk_lst = get_checklist(self)
    chk_lst[self["param_setting"]] = 1
end

function check_config(self::ConfigManager)
    chk_lst = get_checklist(self)
    if chk_lst[self["param_setting"]] == 1
        println(string("Already ran param setting ", self["param_setting"], " ", self["run"]))
        exit(0)
    end
end

function get_checklist(self::ConfigManager)
    filename = join(split(self[ "save_path"],"/"), "_")
    path = joinpath(_checklists(self), string(filename ,".bin"))

    if !isfile(path)
        touch(path)
    end

    f = open(path,"r+")
    return Mmap.mmap(f, Vector{UInt8}, self.total_combinations, (self["run"] - 1) * self.total_combinations * sizeof(UInt8))
end

function unmarked_configs(self::ConfigManager)
    chkLst = get_checklist(self)
    indices = Int[]
    for i=1:total_combinations(self)
        if chkLst[i] != 1
            push!(indices, i)
        end
    end
    return indices
end


text:
struct ConvOp{T} <: AbstractArray{T,3}
    data::Array{T,3}
    k0::Array{Int,2}
    k1::Array{Int,2}
    tail::Array{T,2}
    length::Int
end

function Base.size(obj)
    return (size(obj.data)[2:3]...,obj.length)
end

function Base.getindex(obj::ConvOp, m::Int, n::Int, k::Int)
    if k < obj.k0[m,n]
        return zero(eltype(obj.data))
    end

    if k > obj.k1[m,n]
        return obj.tail[m,n]
    end

    return obj.data[k-obj.k0[m,n]+1,m,n]
end

function convolve!(y, Z::ConvOp, x, X, j, k_start, k_stop=size(Z,3))
    for n in axes(x,1)
        for m in axes(y,1)
            k0 = Z.k0[m,n]
            k1 = Z.k1[m,n]
            for k in max(k0,k_start):min(k1,k_stop)
                p = k - k0 + 1
                j-k < 0 && continue
                y[m] += Z.data[p,m,n] * x[n,j-k+1]
            end

            j-k1 > 0 || continue
            y[m] += Z.tail[m,n] * X[n,j-k1]
        end
    end
end

function polyeig(Z::ConvOp)
    kmax = maximum(Z.k1)
    M = size(Z,1)
    Q = zeros(eltype(Z), 2M, 2M, kmax+1)
    for k in 1:kmax
        Q[1:M,1:M,k] .= Z[:,:,k]
    end
    Id = Matrix{eltype(Z)}(LinearAlgebra.I,M,M)
    Q[M+1:2M,1:M,1] .= -Id
    Q[M+1:2M,M+1:2M,1] .= Id
    Q[M+1:2M,M+1:2M,2] .= -Id
    Q[1:M,M+1:2M,kmax+1] .= Z[:,:,kmax+1]
    return eigvals(companion(Q)), Q
    # return Q
end

function polyeig(Z)
    return eigvals(companion(Z))
end

function Base.:+(a::ConvOp, b::ConvOp)

    @assert size(a.data) == size(b.data)
    M,N = size(a.data)
    T = promote_type(eltype(a.data), eltype(b.data))

    k0 = min.(a.k0, b.k0)
    k1 = max.(a.k1, b.k1)

    bandwidth = maximum(k1 - k0) + 1
    data = zeros(T, bandwidth, M, N)

    k1max = maximum(k1)
    tail = zeros(T,M,N)
    for m in 1:M, n in 1:N
        tail[m,n] = a[m,n,k1max+1] + b[m,n,k1max+1]
    end

    for m in 1:M
        for n in 1:N
            for k in k0[m,n]:k1[m,n]
                data[k,m,n] = a[m,n,k] + b[m,n,k]
            end
        end
    end

    lgt = max(a.length, b.length)
    return ConvOp(data, tail, k0, k1, lgt)
end


text:
module MarkdownFormatTests

using Compat.Test
using Compat: @info

using Documenter

# Documenter package docs
@info("Building Documenter's docs with Markdown.")
const Documenter_root = normpath(joinpath(dirname(@__FILE__), "..", "..", "docs"))
doc = makedocs(
    debug   = true,
    root    = Documenter_root,
    modules = Documenter,
    clean   = false,
)

@testset "Markdown" begin
    @test isa(doc, Documenter.Documents.Document)

    let build_dir  = joinpath(Documenter_root, "build"),
        source_dir = joinpath(Documenter_root, "src")
        @test isdir(build_dir)
        @test isdir(joinpath(build_dir, "assets"))
        @test isdir(joinpath(build_dir, "lib"))
        @test isdir(joinpath(build_dir, "man"))

        @test isfile(joinpath(build_dir, "index.md"))
        @test isfile(joinpath(build_dir, "assets", "mathjaxhelper.js"))
        @test isfile(joinpath(build_dir, "assets", "documenter.css"))
    end

    @test doc.user.root   == Documenter_root
    @test doc.user.source == "src"
    @test doc.user.build  == "build"
    @test doc.user.clean  == false

    rm(joinpath(Documenter_root, "build"); recursive = true)

    mktempdir() do root
        let path = joinpath(root, "docs")
            Documenter.generate("DocumenterTestPackage", dir = path)
            @test isdir(path)
            @test isfile(joinpath(path, "mkdocs.yml"))
            @test isfile(joinpath(path, ".gitignore"))
            @test isfile(joinpath(path, "make.jl"))
            @test isdir(joinpath(path, "src"))
            @test isfile(joinpath(path, "src", "index.md"))
        end
    end

    @test_throws ErrorException Documenter.generate("Documenter")
    @test_throws ErrorException Documenter.generate(randstring())
end

end


text:
function update_scores(cur_player::AbstractPlayer)
  cur_board = Board()

  cur_player.game.board = cur_board

  cur_js = """
    var curBoardName, curBoardScore;
  """

  for (cur_index, cur_entry) in enumerate(cur_board.entries)
    cur_name = "$(rpad(cur_entry.name, 8, "\u00A0"))"[1:8]
    cur_score = "$(lpad(cur_entry.score, 8, "0"))"[1:8]

    cur_js *= """
      curBoardName = \$(".cs-jumbotron table tr:nth-child($(cur_index)) td:nth-child(2)");
      curBoardScore = \$(".cs-jumbotron table tr:nth-child($(cur_index)) td:nth-child(4)");

      curBoardName.text("$(cur_name)");
      curBoardScore.text("$(cur_score)");
    """
  end

  tetris_js(
    cur_player.game.scope,
    JSString(cur_js)
  )
end


text:
immutable Hyperbolic{S<:AbstractDataEnvelopment,T<:RS} <: AbstractDEA{S,T}
  Data::DEAData

  function Hyperbolic{S,T}(X,Y) where {S<:AbstractDataEnvelopment,T<:RS}
    new(DEAData(X,Y))
  end
end

function getdata(DMU::Hyperbolic)
  return DMU.Data
end

#Solve hyperbolic distance function under CRS with (Xk,Yk) as evaluation point
function (H::Hyperbolic{FreeDisposal,CRS})(Xk::Array,Yk::Array)
  X,Y = H[1:end]
  K = size(X,1)
  dompeer = NaN
  dom = find(maximum(reshape(Yk,1,:)./Y,2) .<= minimum(reshape(Xk,1,:)./X,2))
  gamma = Inf
  for k in dom
    Ixk = find(X[k,:] .> 0.0)
    Iyk = find(Y[k,:] .> 0.0)
    domgamma = sqrt(maximum(Yk[Iyk]./Y[k,Iyk])./minimum(Xk[Ixk]./X[k,Ixk]))
    if(domgamma < gamma)
      gamma = domgamma
      dompeer = k
    end
  end
  return DEAResult(gamma,[],[],[k == dompeer ? 1.0 : 0.0 for k=1:K])
end

#Solve hyperbolic distance function under NIRS/NDRS with (Xk,Yk) as evaluation point
function (H::Hyperbolic{FreeDisposal,T})(Xk::Array,Yk::Array) where T<:Union{NIRS,NDRS}
  X,Y = H[1:end]
  K = size(X,1)
  gamma = Inf
  dompeer = NaN
  dom = find(all(Y .>= reshape(Yk,1,:),2))
  for k in dom
    Ixk = find(X[k,:] .> 0.0)
    Iyk = find(Y[k,:] .> 0.0)
    alphak = sqrt(minimum(Xk[Ixk]./X[k,Ixk])*maximum(Yk[Iyk]./Y[k,Iyk]))
    if(T == NIRS)
      if(0.0 <= alphak <= 1)
        domgamma = sqrt(maximum(Yk[Iyk]./Y[k,Iyk])./minimum(Xk[Ixk]./X[k,Ixk]))
      else
        domgamma = max((1/minimum(Xk[Ixk]./X[k,Ixk])),maximum(Yk[Iyk]./Y[k,Iyk]))
      end
    elseif(T == NDRS)
      if(alphak >= 1)
        domgamma = sqrt(maximum(Yk[Iyk]./Y[k,Iyk])./minimum(Xk[Ixk]./X[k,Ixk]))
      else
        domgamma = max((1/minimum(Xk[Ixk]./X[k,Ixk])),maximum(Yk[Iyk]./Y[k,Iyk]))
      end
    end
    if(domgamma < gamma)
      gamma = domgamma
      dompeer = k
    end
  end
  return DEAResult(gamma,[],[],[k == dompeer ? 1.0 : 0.0 for k=1:K])
end

#Solve hyperbolic distance function under VRS with (Xk,Yk) as evaluation point
function (H::Hyperbolic{FreeDisposal,VRS})(Xk::Array,Yk::Array)
  X,Y = H[1:end]

  dom = find(all(Y .>= reshape(Yk,1,:),2) .& all(X .<= reshape(Xk,1,:),2))
  gamma = Inf
  dompeer = NaN
  for k in dom
    Ixk = find(X[k,:] .> 0.0)
    Iyk = find(Y[k,:] .> 0.0)
    curmin = max(maximum(X[k,Ixk]./Xk[Ixk]),1/minimum(Y[k,Iyk]./Yk[Iyk]))
    if(curmin < gamma)
      gamma = curmin
      dompeer = k
    end
  end
  return DEAResult(gamma,[],[],[k == dompeer ? 1.0 : 0.0 for k=1:size(X,1)])
end

function (H::Hyperbolic{S,T})() where {S<:AbstractDataEnvelopment,T<:RS}
  	Data = getdata(H)
  	res = Array{DEAResult}(getnrdmu(Data))

    #@sync @parallel
  	for k in eachindex(Data)
  		res[k] = H(Data[k]...)
  	end

  	return res
end

max_stars_repo_name: kepiej/NEPA


text:
function testend()
  testelement(End(),".end\n")
end
testend()

function testnode()
  testelement(Node(1,1,1),
  "N_1 x=1.000000000e+00 y=1.000000000e+00 z=1.000000000e+00\n" )
  testelement(Node(:abc,2,2,2),
  "Nabc x=2.000000000e+00 y=2.000000000e+00 z=2.000000000e+00\n")
  testelement(Node([3.0,3.0,3.0]),
  "N_1 x=3.000000000e+00 y=3.000000000e+00 z=3.000000000e+00\n")
  testelement(Node(x=4, y=4, z=4),
  "N_1 x=4.000000000e+00 y=4.000000000e+00 z=4.000000000e+00\n")
  testelement(Node(name=:abc, x=5, y=5, z=5),
  "Nabc x=5.000000000e+00 y=5.000000000e+00 z=5.000000000e+00\n")
  coordinate = [1.0,2.0,3.0]
  n1 = Node(coordinate)
  @test xyz(n1) == coordinate
end
testnode()

function testmatrices()
  result = [1.0   0.0       0.0       0.0
         0.0   0.707107  0.707107  0.0
         0.0  -0.707107  0.707107  0.0
         0.0   0.0       0.0       1.0]
  @test sum(rx(π / 4) - result) ≈ 0.0 atol=1.0e-5

  result = [0.707107  0.0  -0.707107  0.0
   0.0       1.0   0.0       0.0
   0.707107  0.0   0.707107  0.0
   0.0       0.0   0.0       1.0]
  @test sum(ry(π / 4) - result) ≈ 0.0 atol=1.0e-5

  result = [0.707107  0.707107  0.0  0.0
   -0.707107  0.707107  0.0  0.0
    0.0       0.0       1.0  0.0
    0.0       0.0       0.0  1.0]
  @test sum(rz(π / 4) - result) ≈ 0.0 atol=1.0e-5

  result = [1.0  0.0  0.0  1.0
   0.0  1.0  0.0  2.0
   0.0  0.0  1.0  3.0
   0.0  0.0  0.0  1.0]
  @test sum(txyz(1.0, 2.0, 3.0) - result) ≈ 0.0 atol=1.0e-5

  result = [ 1.0  0.0  0.0  0.0
   0.0  2.0  0.0  0.0
   0.0  0.0  3.0  0.0
   0.0  0.0  0.0  1.0]
  @test sum(scalexyz(1.0, 2.0, 3.0) - result) ≈ 0.0 atol=1.0e-5
end
testmatrices()

function testsigmarho()
  @test_throws(ArgumentError,FastHenryHelper.SigmaRho(sigma=1.0, rho=2.0))
  sr1 = FastHenryHelper.SigmaRho(sigma=1.0)
  @test sr1.sigma == 1.0
  @test isnan(sr1.rho)
  sr2 = FastHenryHelper.SigmaRho(rho=1.0)
  @test sr2.rho == 1.0
  @test isnan(sr2.sigma)
end
testsigmarho()

function testwxwywz()
  w1 = FastHenryHelper.WxWyWz(wx=1.0,wy=2.0,wz=3.0)
  @test w1.xyz == [1.0, 2.0, 3.0]
  @test w1.isdefault == false
  w2 = FastHenryHelper.WxWyWz()
  @test isnan(w2.xyz[1]) && isnan(w2.xyz[2]) && isnan(w2.xyz[3])
  @test w2.isdefault == true
end
testwxwywz()

function testwh()
  @test_throws(ArgumentError,FastHenryHelper.WH(nhinc=-1))
  @test_throws(ArgumentError,FastHenryHelper.WH(nwinc=-1))
  @test_throws(ArgumentError,FastHenryHelper.WH(rh=-1))
  @test_throws(ArgumentError,FastHenryHelper.WH(rw=-1))
  wh1 = FastHenryHelper.WH()
  @test isnan(wh1.w)
  @test isnan(wh1.h)
  @test wh1.nhinc == 0
  @test wh1.nwinc == 0
  @test isnan(wh1.rh)
  @test isnan(wh1.rw)
  wh2 = FastHenryHelper.WH(w=1.0,h=2.0,nhinc=3.0,nwinc=4.0,rh=5.0,rw=6.0)
  @test wh2.w == 1.0
  @test wh2.h == 2.0
  @test wh2.nhinc == 3.0
  @test wh2.nwinc == 4.0
  @test wh2.rh == 5.0
  @test wh2.rw == 6.0
end
testwh()

function testsegmentparameters()
  sp1 = SegmentParameters(w=1.0,h=2.0,sigma=3.0,wx=4.0,wy=5.0,wz=6.0,
                          nhinc=7.0,nwinc=8.0,rh=9,rw=10)
  sp2 = SegmentParameters(sp1,)
  sp3 = SegmentParameters(sp1,w=21.0,h=22.0,sigma=23.0,wx=24.0,wy=25.0,wz=26.0,
                          nhinc=27.0,nwinc=28.0,rh=29,rw=210)
  @test FastHenryHelper.w(sp2) == 1.0
  @test FastHenryHelper.h(sp2) == 2.0
  @test FastHenryHelper.sigma(sp2) == 3.0
  @test FastHenryHelper.wx(sp2) == 4.0
  @test FastHenryHelper.wy(sp2) == 5.0
  @test FastHenryHelper.wz(sp2) == 6.0
  @test FastHenryHelper.nhinc(sp2) == 7.0
  @test FastHenryHelper.nwinc(sp2) == 8.0
  @test FastHenryHelper.rh(sp2) == 9
  @test FastHenryHelper.rw(sp2) == 10
  @test isnan(FastHenryHelper.rho(sp2))
  @test FastHenryHelper.w(sp3) == 21.0
  @test FastHenryHelper.h(sp3) == 22.0
  @test FastHenryHelper.sigma(sp3) == 23.0
  @test FastHenryHelper.wx(sp3) == 24.0
  @test FastHenryHelper.wy(sp3) == 25.0
  @test FastHenryHelper.wz(sp3) == 26.0
  @test FastHenryHelper.nhinc(sp3) == 27.0
  @test FastHenryHelper.nwinc(sp3) == 28.0
  @test FastHenryHelper.rh(sp3) == 29
  @test FastHenryHelper.rw(sp3) == 210
  @test isnan(FastHenryHelper.rho(sp3))
end
testsegmentparameters()

function testsegment()
  sp1 = SegmentParameters(w=1.0,h=2.0,sigma=3.0,wx=4.0,wy=5.0,wz=6.0,
                          nhinc=7.0,nwinc=8.0,rh=9,rw=10)
  n1 = Node(1,2,3)
  n2 = Node(4,5,6)
  s1 = Segment(n1,n2,sp1)
  verified =
  """
  E_1 N_2 N_3
  +  w=1.000000000e+00 h=2.000000000e+00 nhinc=7 nwinc=8 rh=9.000000000e+00 rw=1.000000000e+01
  +  sigma=3.000000000e+00
  +  wx=-7.071067812e-01 wy=0.000000000e+00 wz=7.071067812e-01
  """
  @test_logs (:warn,"wx,wy,wz not perpendicular to segment length") testelement(s1,verified)
  n3 = Node(0,0,0)
  n4 = Node(1,0,0)
  s2 = Segment(n3,n4,w=.5,h=.2)
  verified =
  """
  E_1 N_2 N_3
  +  w=5.000000000e-01 h=2.000000000e-01
  """
  testelement(s2,verified)
  s3 = Segment(n3,n4,w=.5,h=.2,wx=0,wy=1,wz=0)
  verified =
  """
  E_1 N_2 N_3
  +  w=5.000000000e-01 h=2.000000000e-01
  +  wx=0.000000000e+00 wy=1.000000000e+00 wz=0.000000000e+00
  """
  testelement(s3,verified)
  s4 = Segment(n3,n4,w=.5,h=.2,wx=1,wy=1,wz=0)
  verified =
  """
  E_1 N_2 N_3
  +  w=5.000000000e-01 h=2.000000000e-01
  +  wx=0.000000000e+00 wy=1.000000000e+00 wz=0.000000000e+00
  """
  @test_logs (:warn,"wx,wy,wz not perpendicular to segment length") testelement(s4,verified)
  s5 = Segment(n3,n4,w=.5,h=.2,wx=0,wy=0,wz=1)
  verified =
  """
  E_1 N_2 N_3
  +  w=5.000000000e-01 h=2.000000000e-01
  +  wx=0.000000000e+00 wy=0.000000000e+00 wz=1.000000000e+00
  """
  testelement(s5,verified)
  s6 = Segment(n3,n4,w=.5,h=.2,wx=1,wy=1,wz=1)
  verified =
  """
  E_1 N_2 N_3
  +  w=5.000000000e-01 h=2.000000000e-01
  +  wx=-0.000000000e+00 wy=7.071067812e-01 wz=7.071067812e-01
  """
  @test_logs (:warn,"wx,wy,wz not perpendicular to segment length") testelement(s6,verified)
  @test_throws(ArgumentError,Segment(n1,n2,w=1,h=2,wx=0,wy=0,wz=0))
end
testsegment()

function testcomment()
  c = Comment("abcd")
  testelement(c,"* abcd\n")
end
testcomment()

function testdefault()
  sp1 = SegmentParameters(w=1.0,h=2.0,sigma=3.0,wx=4.0,wy=5.0,wz=6.0,
                          nhinc=7.0,nwinc=8.0,rh=9,rw=10)
  def1 = Default(sp1)
  verified =
  """
  .default
  +  w=1.000000000e+00 h=2.000000000e+00 nhinc=7 nwinc=8 rh=9.000000000e+00 rw=1.000000000e+01
  +  sigma=3.000000000e+00
  """
  testelement(def1,verified)
  def2 = Default(rho=1.0)
  verified =
  """
  .default
  +  rho=1.000000000e+00
  """
  testelement(def2,verified)
end
testdefault()

function testequiv()
  n3 = Node(:abc,3,3,3)
  n4 = Node(:def,4,4,4)
  eq1 = Equiv([n3,n4])
  testelement(eq1,".equiv Nabc Ndef\n")
  n5 = Node(5,5,5)
  n6 = Node(6,6,6)
  eq2 = Equiv([n5,n6])
  testelement(eq2,".equiv N_1 N_2\n")
end
testequiv()

function testexternal()
  n3 = Node(:abc,3,3,3)
  n4 = Node(:def,4,4,4)
  ex1 = External(n3,n4,"portname")
  testelement(ex1,".external Nabc Ndef portname\n")
  n5 = Node(5,5,5)
  n6 = Node(6,6,6)
  ex2 = External(n5,n6,"portname")
  testelement(ex2,".external N_1 N_2 portname\n")
end
testexternal()

function testfreq()
  @test_throws(ArgumentError,Freq())
  @test_throws(ArgumentError,Freq(min=100,max=10))
  @test_throws(ArgumentError,Freq(max=1000, ndec=-10))
  f1 = Freq(min=1,max=2,ndec=3)
  testelement(f1,".freq fmin=1.000000000e+00 fmax=2.000000000e+00 ndec=3.000000000e+00\n")
end
testfreq()

function testunits()
  @test_throws(ArgumentError,Units("kk"))
  u = Units("in")
  testelement(u,".units in\n")
end
testunits()

function testuniformplane()
  p = Point(x=1,y=2,z=3)
  testelement(p,"+ hole point (1.000000000e+00, 2.000000000e+00, 3.000000000e+00)\n")
  r = Rect(x1=1,y1=2,z1=3,x2=4,y2=5,z2=6)
  testelement(r,"+ hole rect (1.000000000e+00, 2.000000000e+00, 3.000000000e+00, 4.000000000e+00, 5.000000000e+00, 6.000000000e+00)\n")
  c = Circle(x=1,y=2,z=3,r=4)
  testelement(c,"+ hole circle (1.000000000e+00, 2.000000000e+00, 3.000000000e+00, 4.000000000e+00)\n")
  @test_throws(ArgumentError,UniformPlane(thick=1, seg1=2, seg2=3, nhinc=4,rh=5,sigma=6,rho=7))
  @test_throws(ArgumentError,UniformPlane(thick=1, seg1=2, seg2=3, nhinc=4,rh=-5,sigma=6))
  @test_throws(ArgumentError,UniformPlane(thick=1, seg1=2, seg2=3, nhinc=-4,rh=5,rho=7))
  @test_throws(ArgumentError,UniformPlane(thick=1, seg1=2, seg2=0, nhinc=4,rh=5,sigma=6))
  @test_throws(ArgumentError,UniformPlane(thick=1, seg1=0, seg2=3, nhinc=4,rh=5,sigma=6))
  @test_throws(ArgumentError,UniformPlane(seg1=2, seg2=3, nhinc=4,rh=5,sigma=6))
  @test_throws(ArgumentError,
    UniformPlane(
    x1=1, y1=2, z1=10,
    x2=1, y2=1, z2=10,
    x3=2, y3=1, z3=10,
    thick=2, seg1=11, seg2=12,
    segwid1 = 13, segwid2 = 14,
    sigma = 15,
    nhinc=16, rh=17,
    relx=18, rely=19, relz=20,
    nodes=[Node(1,1,11.5)]))
  n10 = Node(10,11,1)
  n11 = Node(13,14,1)
  up = UniformPlane(
  x1=1, y1=2, z1=1,
  x2=1, y2=1, z2=1,
  x3=2, y3=1, z3=1,
  thick=10, seg1=11, seg2=12,
  segwid1 = 13, segwid2 = 14,
  sigma = 15,
  nhinc=16, rh=17,
  relx=18, rely=19, relz=20,
  nodes=[n10,n11],
  holes=[p,r,c])
  verified =
  """
  G_1
  + x1=1.000000000e+00 y1=2.000000000e+00 z1=1.000000000e+00
  + x2=1.000000000e+00 y2=1.000000000e+00 z2=1.000000000e+00
  + x3=2.000000000e+00 y3=1.000000000e+00 z3=1.000000000e+00
  + thick=1.000000000e+01 seg1=11 seg2=12
  + segwid1=1.300000000e+01
  + segwid2=1.400000000e+01
  + sigma=1.500000000e+01
  + nhinc=16
  + rh=17
  + relx=1.800000000e+01
  + rely=1.900000000e+01
  + relz=2.000000000e+01
  + N_2 (1.000000000e+01,1.100000000e+01,1.000000000e+00)
  + N_3 (1.300000000e+01,1.400000000e+01,1.000000000e+00)
  + hole point (1.000000000e+00, 2.000000000e+00, 3.000000000e+00)
  + hole rect (1.000000000e+00, 2.000000000e+00, 3.000000000e+00, 4.000000000e+00, 5.000000000e+00, 6.000000000e+00)
  + hole circle (1.000000000e+00, 2.000000000e+00, 3.000000000e+00, 4.000000000e+00)
  """
  testelement(up,verified)
  up2 = deepcopy(up)
  transform!(up2,txyz(1,2,3))
  verified =
  """
  G_1
  + x1=2.000000000e+00 y1=4.000000000e+00 z1=4.000000000e+00
  + x2=2.000000000e+00 y2=3.000000000e+00 z2=4.000000000e+00
  + x3=3.000000000e+00 y3=3.000000000e+00 z3=4.000000000e+00
  + thick=1.000000000e+01 seg1=11 seg2=12
  + segwid1=1.300000000e+01
  + segwid2=1.400000000e+01
  + sigma=1.500000000e+01
  + nhinc=16
  + rh=17
  + relx=1.800000000e+01
  + rely=1.900000000e+01
  + relz=2.000000000e+01
  + N_2 (1.100000000e+01,1.300000000e+01,4.000000000e+00)
  + N_3 (1.400000000e+01,1.600000000e+01,4.000000000e+00)
  + hole point (2.000000000e+00, 4.000000000e+00, 6.000000000e+00)
  + hole rect (2.000000000e+00, 4.000000000e+00, 6.000000000e+00, 5.000000000e+00, 7.000000000e+00, 9.000000000e+00)
  + hole circle (2.000000000e+00, 4.000000000e+00, 6.000000000e+00, 4.000000000e+00)
  """
  testelement(up2,verified)
end
testuniformplane()

function testgroup()
  title = Comment("test title")
  n20 = Node(1,0,0)
  n21 = Node(10,0,0)
  sp20 = SegmentParameters(h=3,w=5,sigma=0.1)
  seg20 = Segment(n20,n21,sp20)
  g1 = Group([title,n20,n21,seg20],Dict(:a=>n20,:b=>n21))
  g2 = transform(g1,rz(π/4))
  g3 = Group([g1,g2],Dict(:c=>g1[:a], :d=>g1[:b], :e=>g2[:a], :f=>g2[:b]))
  testelement(g3[:c],"N_1 x=1.000000000e+00 y=0.000000000e+00 z=0.000000000e+00\n")
  testelement(g3[:d],"N_1 x=1.000000000e+01 y=0.000000000e+00 z=0.000000000e+00\n")
  testelement(g3[:e],"N_1 x=7.071067812e-01 y=-7.071067812e-01 z=0.000000000e+00\n")
  testelement(g3[:f],"N_1 x=7.071067812e+00 y=-7.071067812e+00 z=0.000000000e+00\n")
  verified =
  """
  * test title
  N_1 x=1.000000000e+00 y=0.000000000e+00 z=0.000000000e+00
  N_2 x=1.000000000e+01 y=0.000000000e+00 z=0.000000000e+00
  E_3 N_1 N_2
  +  w=5.000000000e+00 h=3.000000000e+00
  +  sigma=1.000000000e-01
  * test title
  N_4 x=7.071067812e-01 y=-7.071067812e-01 z=0.000000000e+00
  N_5 x=7.071067812e+00 y=-7.071067812e+00 z=0.000000000e+00
  E_6 N_4 N_5
  +  w=5.000000000e+00 h=3.000000000e+00
  +  sigma=1.000000000e-01
  +  wx=7.071067812e-01 wy=7.071067812e-01 wz=0.000000000e+00
  """
  testelement(g3,verified)
  u = Units("in")
  def1 = Default(sigma = 1.234)
  n3 = Node(:abc,3,3,3)
  n4 = Node(:def,4,4,4)
  eq1 = Equiv([n3,n4])
  n5 = Node(:geh,3,3,3)
  n6 = Node(:ijk,4,4,4)
  ex1 = External(n5,n6,"portname")
  f1 = Freq(min=1,max=2,ndec=3)
  p = FastHenryHelper.Point(x=1,y=2,z=3)
  r = FastHenryHelper.Rect(x1=1,y1=2,z1=3,x2=4,y2=5,z2=6)
  c = FastHenryHelper.Circle(x=1,y=2,z=3,r=4)
  n10 = Node(10,11,1)
  n11 = Node(13,14,1)
  up = UniformPlane(
  x1=1, y1=2, z1=1,
  x2=1, y2=1, z2=1,
  x3=2, y3=1, z3=1,
  thick=10, seg1=11, seg2=12,
  segwid1 = 13, segwid2 = 14,
  sigma = 15,
  nhinc=16, rh=17,
  relx=18, rely=19, relz=20,
  nodes=[n10,n11],
  holes=[p,r,c])
  g4 = Group([title,u,n20,n21,n3,n4,n5,n6,seg20,def1,eq1,up,ex1,f1])
  g5 = transform(g4,rx(π/2))
  verified =
  """
  * test title
  .units in
  N_1 x=1.000000000e+00 y=0.000000000e+00 z=0.000000000e+00
  N_2 x=1.000000000e+01 y=0.000000000e+00 z=0.000000000e+00
  Nabc x=3.000000000e+00 y=3.000000000e+00 z=-3.000000000e+00
  Ndef x=4.000000000e+00 y=4.000000000e+00 z=-4.000000000e+00
  Ngeh x=3.000000000e+00 y=3.000000000e+00 z=-3.000000000e+00
  Nijk x=4.000000000e+00 y=4.000000000e+00 z=-4.000000000e+00
  E_3 N_1 N_2
  +  w=5.000000000e+00 h=3.000000000e+00
  +  sigma=1.000000000e-01
  +  wx=0.000000000e+00 wy=6.123233996e-17 wz=-1.000000000e+00
  .default
  +  sigma=1.234000000e+00
  .equiv Nabc Ndef
  G_4
  + x1=1.000000000e+00 y1=1.000000000e+00 z1=-2.000000000e+00
  + x2=1.000000000e+00 y2=1.000000000e+00 z2=-1.000000000e+00
  + x3=2.000000000e+00 y3=1.000000000e+00 z3=-1.000000000e+00
  + thick=1.000000000e+01 seg1=11 seg2=12
  + segwid1=1.300000000e+01
  + segwid2=1.400000000e+01
  + sigma=1.500000000e+01
  + nhinc=16
  + rh=17
  + relx=1.800000000e+01
  + rely=1.900000000e+01
  + relz=2.000000000e+01
  + N_5 (1.000000000e+01,1.000000000e+00,-1.100000000e+01)
  + N_6 (1.300000000e+01,1.000000000e+00,-1.400000000e+01)
  + hole point (1.000000000e+00, 3.000000000e+00, -2.000000000e+00)
  + hole rect (1.000000000e+00, 3.000000000e+00, -2.000000000e+00, 4.000000000e+00, 6.000000000e+00, -5.000000000e+00)
  + hole circle (1.000000000e+00, 3.000000000e+00, -2.000000000e+00, 4.000000000e+00)
  .external Ngeh Nijk portname
  .freq fmin=1.000000000e+00 fmax=2.000000000e+00 ndec=3.000000000e+00
  """
  testelement(g5,verified)
  c1 = Comment("this is a comment")
  # Default and Freq are immutable
  g6 = Group([title,c1,u,n20,n21,n3,n4,n5,n6,seg20,eq1,up,ex1])
  deepcopyg6 = deepcopy(g6)
  for i in (4,5,6,7,8,9,10) # only verify mutable elements
    @test ~(elements(deepcopyg6)[i] === elements(g6)[i])
  end

  g7 = Group()
  elements!(g7,[title,n20,n21,seg20])
  terms!(g7,Dict(:a=>n20,:b=>n21))
  g7[:c] = n20
  @test g7[:a] === g7[:c]
  deepcopyg7 = deepcopy(g7)
  @test elements(deepcopyg7)[4].node1 === elements(deepcopyg7)[2]
  @test elements(deepcopyg7)[4].node2 === elements(deepcopyg7)[3]
  push!(g1,n20)
  @test pop!(g1) === n20
  pushfirst!(g1,n20)
  @test popfirst!(g1) === n20
  append!(g1,g2)
  prepend!(g1,g2)
  merge!(g1,g2,g3)
  verified =
  """
  * test title
  N_1 x=7.071067812e-01 y=-7.071067812e-01 z=0.000000000e+00
  N_2 x=7.071067812e+00 y=-7.071067812e+00 z=0.000000000e+00
  E_3 N_1 N_2
  +  w=5.000000000e+00 h=3.000000000e+00
  +  sigma=1.000000000e-01
  +  wx=7.071067812e-01 wy=7.071067812e-01 wz=0.000000000e+00
  * test title
  N_4 x=1.000000000e+00 y=0.000000000e+00 z=0.000000000e+00
  N_5 x=1.000000000e+01 y=0.000000000e+00 z=0.000000000e+00
  E_6 N_4 N_5
  +  w=5.000000000e+00 h=3.000000000e+00
  +  sigma=1.000000000e-01
  * test title
  N_1 x=7.071067812e-01 y=-7.071067812e-01 z=0.000000000e+00
  N_2 x=7.071067812e+00 y=-7.071067812e+00 z=0.000000000e+00
  E_3 N_1 N_2
  +  w=5.000000000e+00 h=3.000000000e+00
  +  sigma=1.000000000e-01
  +  wx=7.071067812e-01 wy=7.071067812e-01 wz=0.000000000e+00
  """
  testelement(g1,verified)
  g7 = Group([Node(1,2,3)])
  g8 = Group([g7])
  g9 = Group([g8])
  g10 = Group([g9])
  g11 = Group([g10])
  g12 = Group([g11])
  testelement(g12,"N_1 x=1.000000000e+00 y=2.000000000e+00 z=3.000000000e+00\n")
  for element in FastHenryHelper.traverseTree(g12)
    @test typeof(element) != Group
  end
end
testgroup()

function testgroup2()
  g1 = Group()
  g2 = Group([Node(1,1,1)])
  a=Node(1,1,1)
  b=Node(2,2,2)
  g3 = Group([a],Dict(:terma=>a))
  g4 = Group([a,b],Dict(:terma=>[a,b]))
  g5 = Group([a],Dict(:terma=>a),Units("mm"))
end
testgroup2()

function testrecursivedeepcopy()
  g1 = Group([Node(1,1,1)])
  e = elements(g1)
  push!(e,Node(2,2,2))
  push!(e,g1);
  g2 = deepcopy(g1);
  @test ~(g2===g1)
end
testrecursivedeepcopy()


text:
using Gen: ifelse

@gen (static) function datum(x::Float64, (grad)(inlier_std::Float64), (grad)(outlier_std::Float64),
                             (grad)(slope::Float64), (grad)(intercept::Float64))
    is_outlier = @trace(bernoulli(0.5), :z)
    std = ifelse(is_outlier, inlier_std, outlier_std)
    y = @trace(normal(x * slope + intercept, std), :y)
    return y
end

data = Map(datum)

@gen (static) function model(xs::Vector{Float64})
    n = length(xs)
    inlier_log_std = @trace(normal(0, 2), :log_inlier_std)
    outlier_log_std = @trace(normal(0, 2), :log_outlier_std)
    inlier_std = exp(inlier_log_std)
    outlier_std = exp(outlier_log_std)
    slope = @trace(normal(0, 2), :slope)
    intercept = @trace(normal(0, 2), :intercept)
    @trace(data(xs, fill(inlier_std, n), fill(outlier_std, n),
        fill(slope, n), fill(intercept, n)), :data)
end


text: println("Hello Remote World!")


text:
# PANOC problem : tiny nonconvex QP
#
# minimize      Quadratic[Q,q](x) + IndBox[l,u](x)
#
# where
#    Q = [-0.5, 0.0; 0.0, 1.0]
#    q = [0.3; 0.5]
#    l = [-1.0; -1.0]
#    u = [1.0; 1.0]

using OptiMo

# Type
mutable struct PPNCVXQP <: AbstractOptiModel
  meta::OptiModelMeta
end

# Constructor
function PPNCVXQP()
  meta = OptiModelMeta(2, 0, name="Nonconvex QP (tiny)")
  return PPNCVXQP( meta )
end

##########################################################
# Necessary methods:
# obj, grad!, prox!, objprox!
##########################################################
function OptiMo.obj(prob::PPNCVXQP, x::AbstractVector)
    # f(x)
    return 0.5 * x[2]^2 - 0.25 * x[1]^2 + 0.3 * x[1] + 0.5 * x[2]
end

function OptiMo.grad!(prob::PPNCVXQP, x::AbstractVector, dfx::AbstractVector)
    # df/dx(x)
    dfx[1] = - 0.5 * x[1] + 0.3
    dfx[2] = x[2] + 0.5
    return nothing
end

function OptiMo.prox!(prob::PPNCVXQP, x::AbstractVector, a::Real, z::AbstractVector)
    # prox_{a g}( x )
    z[1] = max(-1.0, min(x[1], 1.0))
    z[2] = max(-1.0, min(x[2], 1.0))
    return nothing
end

function OptiMo.objprox!(prob::PPNCVXQP, x::AbstractVector, a::Real, z::AbstractVector)
    # prox_{a g}( x )
    # g(z)
    z[1] = max(-1.0, min(x[1], 1.0))
    z[2] = max(-1.0, min(x[2], 1.0))
    return 0.0
end


text:
# This file is a part of EmpiricalDistributions.jl, licensed under the MIT License (MIT).

using EmpiricalDistributions
using Test

using Random, LinearAlgebra
using Distributions, StatsBase
using Adapt, ForwardDiff

@testset "uv_binned_dist" begin
    μ, σ = 1.23, 0.74
    true_dist = Normal(μ, σ)
    h = Histogram(μ-10σ:σ/10:μ+10σ)
    append!(h, rand(true_dist, 10^7))

    @test @inferred(UvBinnedDist(h)) isa UvBinnedDist
    d = UvBinnedDist(h)
    @test @inferred(UvBinnedDist{Float32}(h)) isa UvBinnedDist{Float32}
    d_32 = UvBinnedDist{Float32}(h)
    @test @inferred(UvBinnedDist{Float64}(h)) isa UvBinnedDist{Float64}
    d_64 = UvBinnedDist{Float64}(h)

    @inferred(adapt(EmpiricalDistributions._AdaptToArray(), d))._edge isa Array

    @test @inferred(convert(Histogram, d)) == normalize(h)

    @test @inferred(length(d)) == 1
    @test @inferred(size(d)) == ()
    @test @inferred(eltype(d)) == Float64

    @test all(isapprox.(mean(true_dist), @inferred(mean(d)), atol = 0.01))
    @test all(isapprox.(mode(true_dist), @inferred(mode(d)), atol = 0.05))
    @test all(isapprox.(var(true_dist), @inferred(var(d)), atol = 0.01))

    @test @inferred(minimum(d)) == μ-10σ
    @test @inferred(maximum(d)) == μ+10σ

    @test @inferred(rand(d)) isa Real
    @test @inferred(rand(d_32)) isa Float32
    @test @inferred(rand(d_64)) isa Float64
    @test @inferred(rand(d, 10)) isa Vector{<:Real}
    @test @inferred(rand(d_32, 10)) isa Vector{Float32}
    @test @inferred(rand(d_64, 10)) isa Vector{Float64}
    @test size(rand(d, 10)) == (10,)

    r = rand(d, 10^6)

    fit_dist = fit(Normal, r)
    @test isapprox(μ, fit_dist.μ, atol = 0.01)
    @test isapprox(σ, fit_dist.σ, atol = 0.01)

    @test @inferred(pdf(d, r[1])) isa Real
    @test isapprox(pdf.(Ref(d), r), pdf.(Ref(true_dist), r), rtol = 0.1)
    @test @inferred(pdf(d, 1000)) == 0
    @test @inferred(pdf(d, -1000)) == 0
    @test @inferred(pdf(d, first(h.edges[1]))) == 0
    @test @inferred(pdf(d, last(h.edges[1]))) == 0

    @test @inferred(logpdf(d, r[1])) == log(pdf(d, r[1]))
    @test @inferred(logpdf(d, 1000)) == -Inf
    @test @inferred(logpdf(d, -1000)) == -Inf

    @test @inferred(cdf(d, r[1])) isa Real
    @test isapprox(cdf.(Ref(d), r), cdf.(Ref(true_dist), r), rtol = 0.1)
    @test @inferred(cdf(d, -1000)) == 0
    @test @inferred(cdf(d, 1000)) == 1

    diff_cdf(d::Distribution, x::Real) = ForwardDiff.derivative(a -> cdf(d, a), x)
    @test isapprox(diff_cdf.(Ref(d), r), pdf.(Ref(d), r), rtol = 0.1)

    @test @inferred(logcdf(d, r[1])) isa Real
    @test isapprox(logcdf.(Ref(d), r), logcdf.(Ref(true_dist), r), rtol = 0.1)
    @test @inferred(logcdf(d, -1000)) == -Inf
    @test @inferred(logcdf(d, 1000)) == 0

    @test @inferred(quantile(d, 0.5)) isa Real
    @test @inferred(quantile(d_32, 0.5f0)) isa Float32
    @test @inferred(quantile(d_64, 0.5)) isa Float64
    @test isapprox(quantile.(Ref(d), 0.05:0.01:0.95), quantile.(Ref(true_dist), 0.05:0.01:0.95), rtol = 0.1)
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Compat.Test

struct dv_test{T} <: BAT.BATDataVector{T}
    data::Vector{T}
end

Base.append!(x::dv_test, y::dv_test) = append!(x.data, y.data)

@testset "data_vector" begin

    @testset "data_vector" begin
        x = @inferred dv_test{Float64}(Float64(1):Float64(3))
        @test typeof(x) <: DenseVector{Float64}
    end

    @testset "merge" begin
        x = @inferred dv_test{Float64}(Float64(1):Float64(3))
        y = @inferred dv_test{Float64}(Float64(4):Float64(6))
        z = @inferred dv_test{Float64}(Float64(7):Float64(9))
        @test merge(x, y, z).data == collect(Float64(1):Float64(9))
        merge!(x, y, z)
        @test x.data == collect(Float64(1):Float64(9))
    end
end


text:
# This file is part of the HORIZONS.jl package
# The HORIZONS.jl package is licensed under the MIT "Expat" License
# Copyright (c) 2017: Jorge Pérez.

#The following methods implement some functionality from JPL's Horizons telnet interface
#For more detailed info about JPL HORIZONS system, visit http://ssd.jpl.nasa.gov/?horizons

# The following code is based on the vec_tbl script, as was retrieved from:
# ftp://ssd.jpl.nasa.gov/pub/ssd/SCRIPTS/vec_tbl
# Date retrieved: Jul 19, 2017
# Credit: Jon D. Giorgini, NASA-JPL
# Jon.D.Giorgini@jpl.nasa.gov

"""
    vec_tbl(OBJECT_NAME, START_TIME, STOP_TIME, STEP_SIZE; kwargs...)

    vec_tbl(OBJECT_NAME, local_file, START_TIME, STOP_TIME, STEP_SIZE; kwargs...)

Automate the Horizons session required to produce a VECTOR table for an
object already listed in the Horizons database: a planet, natural satellite,
asteroid, comet, spacecraft, or dynamical point. If the file name `local_file`
is specified, then output is saved to that file in current folder. Otherwise,
then output is returned as a Julia string.

Generally, this script is suited for situations where the same output
format (as defined by the input file) is desired for a list of objects
specified one at a time on the script's command-line.

The current keyword arguments are:

    + `timeout=15`
    + `EMAIL_ADDR = "joe@your.domain.name"`
    + `CENTER = "@ssb"`
    + `REF_PLANE = "FRAME"`
    + `COORD_TYPE = "G"`
    + `SITE_COORD = "0,0,0"`
    + `REF_SYSTEM = "J2000"`
    + `VEC_CORR = 1`
    + `VEC_DELTA_T = false`
    + `OUT_UNITS = 1`
    + `CSV_FORMAT = false`
    + `VEC_LABELS = false`
    + `VEC_TABLE= 3`

More detailed information may be found at the HORIZONS system documentation:

https://ssd.jpl.nasa.gov/?horizons_doc

The original script vec_tbl, written by Jon D. Giorgini, may be found at src/SCRIPTS

"""
function vec_tbl(OBJECT_NAME::ObjectName, START_TIME::StartStopTime,
        STOP_TIME::StartStopTime, STEP_SIZE::StepSize; kwargs...)

    output_str, ftp_name = get_vec_tbl(OBJECT_NAME, DateTime(START_TIME), DateTime(STOP_TIME), STEP_SIZE; kwargs...)

    return output_str

    #TODO: turn output_str into a data table; possibly even a struct which saves object info + ephemeris
end

function vec_tbl(OBJECT_NAME::ObjectName, local_file::String,
        START_TIME::StartStopTime, STOP_TIME::StartStopTime,
        STEP_SIZE::StepSize; EMAIL_ADDR::String="joe@your.domain.name",
        ftp_verbose::Bool=false, kwargs...)

    output_str, ftp_name = get_vec_tbl(OBJECT_NAME, DateTime(START_TIME), DateTime(STOP_TIME), STEP_SIZE; kwargs...)

    # Retrieve file by anonymous FTP and save to file `local_file`
    ftp_init()
    # workaround `@` in email address
    ftp_email = replace(EMAIL_ADDR, "@" => "_at_")
    ftp = FTP(hostname=HORIZONS_MACHINE, username="anonymous", password=ftp_email, verbose=ftp_verbose)
    cd(ftp, HORIZONS_FTP_DIR)
    if local_file == ""
        io = download(ftp, ftp_name, ftp_name)
        close(ftp)
        ftp_cleanup()
        return ftp_name
    else
        io = download(ftp, ftp_name, local_file)
        close(ftp)
        return local_file
    end
end

function get_vec_tbl(OBJECT_NAME::ObjectName, START_TIME::DateTime,
        STOP_TIME::DateTime, STEP_SIZE::StepSize; timeout::Int=15,
        EMAIL_ADDR::String="joe@your.domain.name", CENTER::String="@ssb",
        REF_PLANE::String="FRAME", COORD_TYPE::String="G",
        SITE_COORD::String="0,0,0", REF_SYSTEM::String="J2000",
        VEC_CORR::Int=1, VEC_DELTA_T::Bool=false, OUT_UNITS::Int=1,
        CSV_FORMAT::Bool=false, VEC_LABELS::Bool=false, VEC_TABLE::VecTable=3)

    # Convert start and stop time from `DateTime`s to `String`s
    OBJECT_NAME_str = string(OBJECT_NAME)
    START_TIME_str = Dates.format(START_TIME, HORIZONS_DATE_FORMAT)
    STOP_TIME_str = Dates.format(STOP_TIME, HORIZONS_DATE_FORMAT)
    STEP_SIZE_str = string(STEP_SIZE)

    start_flag = 0

    # Connect to Horizons
    proc = ExpectProc(`telnet $HORIZONS_MACHINE 6775`, timeout)

    # Get main prompt and proceed, turning off paging, specifying I/O model,
    # and sending object look-up from command-line
    idx = expect!(proc, ["unknown host", "Horizons> "])
    if idx == 1
        warn("This system cannot find $HORIZONS_MACHINE")
        close(proc)
    elseif idx == 2
        println(proc, "PAGE")
    end

    idx = expect!(proc, ["Horizons> "])
    if idx == 1
        println(proc, "##2")
    end

    idx = expect!(proc, ["Horizons> "])
    if idx == 1
        println(proc, OBJECT_NAME_str)
    end

    # Handle object look-up confirmation
    idx = expect!(proc, [r".*Continue.*: $", r".*such object record.*", r".*Select.*<cr>: $"])
    if idx == 1
        println(proc, "yes")
        idx = expect!(proc, [r".*PK.*: $", r".*lay.*: $"])
        if idx ==1
            println(proc, "E")
        elseif idx == 2
            println(proc, "x")
            throw(ArgumentError("Cancelled -- unique object not found: $OBJECT_NAME\nObject not matched to database OR multiple matches found."))
        end
    elseif idx == 2
        # currently unable to reproduce this case on HORIZONS v4.10
        # println(proc, "x")
        # throw(ArgumentError("No such object record found: $OBJECT_NAME"))
    elseif idx == 3
        println(proc, "E")
    end

    # Request VECTOR table
    idx = expect!(proc, [r".*Observe, Elements.*: $"])
    if idx == 1
        println(proc, "V")
    end

    # Provide coordinate center
    idx = expect!(proc, [r".*Coordinate .*: $"])
    if idx == 1
        println(proc, CENTER)
    end

    # Handle coordinate center error or confirmation
    idx = expect!(proc, [r".*No site matches.*: $"s, r".*Select.*<cr>: $", r".*Coordinate center.*: $"s, r".*Confirm selected.* $", r".*Cylindrical.*: $", r".*Reference plane.*: $"])

    if idx == 1
        println(proc, "X")
        throw(ArgumentError("Cannot find CENTER = $CENTER"))
    elseif idx == 2
        # currently unable to reproduce this case on HORIZONS v4.10
        # println(proc, "X")
        # throw(ArgumentError("Non-unique CENTER = $CENTER (multiple matches); idx == $idx"))
    elseif idx == 3
        println(proc, "X")
        throw(ArgumentError("Non-unique CENTER = $CENTER (multiple matches); idx == $idx"))
    elseif idx == 4
        println(proc, "Y")
        idx = expect!(proc, [r".*Reference plane.*: $"])
        if idx == 1
            println(proc, REF_PLANE)
        end
    elseif idx == 5
        println(proc, COORD_TYPE)
        idx = expect!(proc, [r".*Unknown.*: $"s, r".*Enter c or g.*: $"s, r".*Specify.*: $"])
        if idx == 1
            println(proc, "X")
            throw(ArgumentError("Unrecognized user-input coordinate: COORD_TYPE = $COORD_TYPE"))
        elseif idx == 2
            println(proc, "X")
            throw(ArgumentError("Undefined or bad coordinate type: COORD_TYPE = $COORD_TYPE"))
        elseif idx == 3
            println(proc, SITE_COORD)
            idx = expect!(proc, [r".*Cannot read.*: $"s, r".*Specify.*: $"s, r".*Reference plane.*: $"])
            if idx == 1
                println(proc, "X")
                throw(ArgumentError("Cannot read coordinate triplet: SITE_COORD=$SITE_COORD"))
            elseif idx == 2
                println(proc, "X")
                throw(ArgumentError("Undefined site coordinate triplet: SITE_COORD = $SITE_COORD"))
            elseif idx == 3
                println(proc, REF_PLANE)
            end
        end
    elseif idx == 6
        println(proc, REF_PLANE)
    end

    # Handle reference plane error or START date
    idx = expect!(proc, [r".*Enter.*abbreviation.*: $"s, r".*Starting .*: $"])
    if idx == 1
        println(proc, "X")
        throw(ArgumentError("Enter \"ecliptic\", \"frame\" or \"body equator\" or abbreviation: REF_PLANE = $REF_PLANE\nSee Horizons documentation for available options."))
    elseif idx == 2
        start_flag = 1
        println(proc, START_TIME_str)
    end

    # Handle start date error or STOP date
    idx = expect!(proc, [r".*Cannot interpret.*: $", r".*No ephemeris.*: $"s, r".*Ending.*: $"])
    if idx == 1
        # this case is handled more "julianly" by argument typing and Base.Dates
        # println(proc, "X")
        # throw(ArgumentError("Error in date format: START_TIME = $START_TIME_str\nSee Horizons documentation for accepted formats."))
    elseif idx == 2
        println(proc, "X")
        throw(ArgumentError("START_TIME = $START_TIME_str prior to available ephemeris"))
    elseif idx == 3
        println(proc, STOP_TIME_str)
    end

    # Handle stop date error or get step size
    idx = expect!(proc, [r".*Cannot interpret.*", r".*No ephemeris.*"s, r".*Output interval.*: $"])
    if idx == 1
        # this case is handled more "julianly" by argument typing and Base.Dates
        # println(proc, "X")
        # throw(ArgumentError("Error in date format: STOP_TIME = $STOP_TIME_str\nSee Horizons documentation for accepted formats."))
    elseif idx == 2
        println(proc, "X")
        throw(ArgumentError("STOP_TIME = $STOP_TIME_str date beyond available ephemeris."))
    elseif idx == 3
        println(proc, STEP_SIZE_str)
    end

    # Handle step-size error or proceed to defaults
    idx = expect!(proc, [r".*Unknown.*: $"s, r".*Cannot use.*: $"s, r".*Accept default.*: $"s])
    if idx == 1
        println(proc, "X")
        throw(ArgumentError("STEP_SIZE = $STEP_SIZE_str error."))
    elseif idx == 2
        println(proc, "X")
        throw(ArgumentError("STEP_SIZE = $STEP_SIZE_str error."))
    elseif idx == 3
        println(proc, "N") # never accept table defaults
    end

    # Change output table defaults
    while true
        idx = expect!(proc, [r"(Cannot interpret.*\r)", r".*frame.*].*: $", r".*Corrections.*].*: $", r".*units.*].*: $", r".*CSV.*].*: $", r".*Label.*].*: $", r".*delta-T.*].*: $", r".*table type.*].*: $", r".*Select.*: $", r".*].*: $"])
        if idx == 1
            # this case is handled automatically by never accepting table defaults
            # println(proc, "X")
            # throw(ArgumentError("Error in $proc.match, $proc.before. \nSee Horizons documentation for acceptable values."))
        elseif idx == 2
            println(proc, REF_SYSTEM)
        elseif idx == 3
            println(proc, "$VEC_CORR")
        elseif idx == 4
            println(proc, "$OUT_UNITS")
        elseif idx == 5
            println(proc, yesornostring(CSV_FORMAT))
        elseif idx == 6
            println(proc, yesornostring(VEC_LABELS))
        elseif idx == 7
            println(proc, yesornostring(VEC_DELTA_T))
        elseif idx == 8
            println(proc, "$VEC_TABLE")
        elseif idx == 9
            break # Done w/default override
        # elseif idx == 10
        #     # currently unable to reproduce this case on HORIZONS v4.10
        #     # println(proc, "") # Skip unknown (new?) prompt
        end
    end
    # expect!(proc, r".*Select.*: $")

    # println(proc.before)
    # @show typeof(proc.before)
    output_str = proc.before

    # Osculating element table output has been generated. Now sitting at
    # post-output prompt. Initiate FTP file transfer.
    println(proc, "F")

    # Pick out ftp file name
    result = expect!(proc, r"File name   : (.*)\r\r\n   File type")
    proc_match = match(r"File name   : (.*)\r\r\n   File type", proc.match)
    # name of file at FTP server
    ftp_name = strip(proc_match[1]) #quit possible trailing whitespaces

    # Close telnet connection
    # println(proc, "exit")
    close(proc)

    return output_str, ftp_name
end

function vec_tbl_csv(OBJECT_NAME::ObjectName, START_TIME::StartStopTime,
        STOP_TIME::StartStopTime, STEP_SIZE::StepSize; kwargs...)

    output_str, ftp_name = get_vec_tbl(OBJECT_NAME, DateTime(START_TIME),
        DateTime(STOP_TIME), STEP_SIZE; CSV_FORMAT=true, VEC_LABELS=false,
        kwargs...)

    # get $$SOE, $$EOE offsets
    mSOE = match(r"\$\$SOE", output_str)
    mEOE = match(r"\$\$EOE", output_str)

    # get everything within SOE and EOE
    ste = output_str[mSOE.offset+7:mEOE.offset-1]
    ste = replace(ste, ",\r\n" => "\r\n") #get rid of comma at end of line

    # get and process table labels
    hdr_raw = match(r"JDTDB.*,\r\n", output_str).match
    hdr_raw = replace(hdr_raw, "Calendar Date (TDB)" => "Calendar_Date_TDB") #format calendar date string
    hdr_raw = replace(hdr_raw, ",\r\n" => "\r\n") #get rid of comma at end of line

    #string in CSV format with BOTH headers and data
    csv_str = hdr_raw*ste

    #convert data string into array
    arr = readdlm(IOBuffer(ste), ',')
    arr[:,2] = strip.( arr[:,2] ) #get rid of whitespaces

    # convert labels string into array
    hdr = convert(Array{String,2}, strip.(  readdlm( IOBuffer( hdr_raw ), ',' )  ));
    # vcat into common 2-dim array
    tbl = vcat(hdr, arr)

    #return labels and data as 2-dim array (table) and CSV-formatted string
    return tbl, csv_str
end


text:
# * Slater determinants

"""
    SlaterDeterminant(orbitals::Vector{O})

Constructs a Slater determinant from a set of spin-orbitals.

# Examples

```jldoctest
julia> SlaterDeterminant([:a, :b])
a(1)b(2) - a(2)b(1)

julia> SlaterDeterminant([:a, :b, :c])
a(1)b(2)c(3) - a(1)b(3)c(2) - a(2)b(1)c(3) + a(2)b(3)c(1) + a(3)b(1)c(2) - a(3)b(2)c(1)
```
"""
struct SlaterDeterminant{O}
    orbitals::Vector{O}
end

"""
    SlaterDeterminant(configuration::Configuration{<:SpinOrbital})

Constructs a Slater determinant from the spin-orbitals of the
spin-configuration `configuration`.

# Examples

```jldoctest
julia> SlaterDeterminant(spin_configurations(c"1s2")[1])
1s₀α(1)1s₀β(2) - 1s₀α(2)1s₀β(1)
```
"""
SlaterDeterminant(c::Configuration{<:SpinOrbital}) =
    SlaterDeterminant(c.orbitals)

"""
    length(slater_determinant)

Return the number of spin-orbitals in the Slater determinant.
"""
Base.length(sd::SlaterDeterminant) = length(sd.orbitals)

Base.show(io::IO, sd::SlaterDeterminant) =
    write(io, "|",join(sd.orbitals, " "), "|")

function Base.show(io::IO, ::MIME"text/plain", sd::SlaterDeterminant)
    noprintterms = 5
    n = length(sd)
    if n > 3
        show(io, sd)
        return
    end
    numperm = factorial(n)
    for (j,p) in enumerate(permutations(1:n))
        if j < noprintterms || j > numperm-noprintterms
            s = Combinatorics.parity(p)
            (j > 1 || s == 1) && write(io, " ", s == 0 ? "+" : "-", " ")
            for (i,c) in enumerate(p)
                write(io, "$(sd.orbitals[i])($(c))")
            end
        elseif j == noprintterms
            write(io, " + … ")
        end
    end
end

"""
    AdjointSlaterDeterminant(slater_determinant)

Representation of the Hermitian conjugate (dual vector) of a Slater
determinant. Constructed via the usual [`adjoint`](@ref) operator.
"""
struct AdjointSlaterDeterminant{O}
    sd::SlaterDeterminant{O}
end

"""
    adjoint(slater_determinant)

Construct the adjoint of `slater_determinant`

# Examples

```jldoctest
julia> SlaterDeterminant([:a, :b])'
[a(1)b(2) - a(2)b(1)]†
```
"""
Base.adjoint(sd::SlaterDeterminant) =
    AdjointSlaterDeterminant(sd)

"""
    length(adjoint_slater_determinant)

Return the number of spin-orbitals in the adjoint Slater determinant.
"""
Base.length(asd::AdjointSlaterDeterminant) = length(asd.sd)

function Base.show(io::IO, asd::AdjointSlaterDeterminant)
    write(io, "[")
    show(io, asd.sd)
    write(io, "]†")
end

export SlaterDeterminant


text:
# Driver script for all unit tests
using Jive

runtests(@__DIR__, skip=["runtests.jl"])


text:
#Hypotrochoid

using Luxor

@png begin
    origin()
    background("grey15")
    sethue("antiquewhite")
    setline(1)
    p = hypotrochoid(100, 25, 55, :stroke, stepby=0.01, vertices=true)
    for i in 0:3:15
    Luxor.poly(offsetpoly(p, i), :stroke, close=true)
    end
end

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
@testset "890.find-and-replace-pattern.jl" begin
    @test find_and_replace_pattern(["abc", "deq", "mee", "aqq", "dkd", "ccc"], "abb") ==
          ["mee", "aqq"]
end


text:
include("setting.jl");
using DelimitedFiles;
@time using OrdinaryDiffEq;

x_a_oi = readdlm("data/x_a_oi.txt");
cycles = size(x_a_oi, 1);

ti_start = 4+1;
t1 = 4.0*dT
t2 = 8.0*dT;
α = 2.5e-1 # rescaling factor

B = ones(Float64, N, N)

for d = 0:20
    coeffs = Vector();

    for i in ti_start:cycles
        # integrate time period t1, end in t
        prob = ODEProblem(lorenz96, x_a_oi[i, :], (0.0, t1));
        sol = solve(prob, Tsit5());
        xf_t1 = sol[end];

        # integrate time period t2, end in t
        prob = ODEProblem(lorenz96, x_a_oi[i-4, :], (0.0, t2));
        sol = solve(prob, Tsit5());
        xf_t2 = sol[end];

        δf = xf_t2 - xf_t1;
        pb = δf * δf';
        for j in 1:N
            pos = (j+d+40-1)%40+1;
            neg = (j-d+40-1)%40+1;
            push!(coeffs, pb[j, pos]);
            push!(coeffs, pb[j, neg]);
        end
    end

    filter!(x->!isnan(x), coeffs);
    E = sum(coeffs) / length(coeffs);
    for i in 1:N
        pos = (i+d+40-1)%40+1;
        neg = (i-d+40-1)%40+1;
        global B[i, pos] = α * E
        global B[i, neg] = α * E
    end
end

open("data/b_nmc.txt", "w") do nmc
    writedlm(nmc, B)
end

max_stars_repo_name: Typas/Data-Assimilation-Project

max_issues_repo_name: Typas/Data-Assimilation-Project

max_forks_repo_name: Typas/Data-Assimilation-Project


text:
#
#      TB - the manifold of the tangent bundle
#  Point is a Point on the n-dimensional sphere.
#
import LinearAlgebra: norm, dot
import Base: exp, log, show
export TangentBundle, TBPoint, TBTVector, show, getValue, getTangent, getBase
export distance, dot, exp, log, manifoldDimension, norm
export randomMPoint, parallelTransport, zeroTVector, typeofMPoint, typeofTVector
export validateMPoint, validateTVector
#
# Type definitions
#

@doc doc"""
    TangentBundle <: Manifold

The manifold $\mathcal M = T\mathcal N$ obtained by looking at the tangent
bundle of a [`Manifold`](@ref)s tangent spaces. The manifold obtained is of
dimension $2d$, where $d$ is the dimension of the manifold $\mathcal N$
considered.

To keep notations clear, small letters will always refer to points (`x,y`) or
tangent vectors (`ξ,η`) on the manifold $\mathcal N$, while capital letters
(`X, Y, Z` and `Ξ,Η`) will refer to points and tangent vectors in the tangent
bundle respectively.

# Abbreviation
`TB`

# Constructor
    TangentBundle(M)

generates the tangent bundle to the [`Manifold`](@ref) `M`.
"""
struct TangentBundle{M <: Manifold} <: Manifold
  name::String
  dimension::Int
  abbreviation::String
  manifold::M where {M <: Manifold}
  TangentBundle{mT}(M::mT) where {mT <: Manifold} = new{mT}("Tangent bundle of $(M.name)",
    2*manifoldDimension(M),"TB($(M.abbreviation))",M)
end
TangentBundle(M::mT) where {mT <: Manifold} = TangentBundle{mT}(M)
@doc doc"""
    getBase(M)

return the base manifold of the [`TangentBundle`](@ref) [`Manifold`](@ref) `M`.
"""
getBase(M::TangentBundle) = M.manifold

@doc doc"""
    TBPoint <: MPoint

A point $N\in \mathcal M$ on the manifold $\mathcal M = T\mathcal N$
represented by a tuple `(x,ξ)`, where $x\in\mathcal N$ is a point on the manifold
and $\xi=\xi_x\in T_x\mathcal N$ is a point in the tangent space at $x$.

Two constructors are available:
* `TBPoint(x,ξ)` to construct a tangent bundle point by specifying both
  an [`MPoint`](@ref) `x` and a [`TVector`](@ref) `ξ`.
* `TBPoint( (X) )` to construct a tangent bundle point from a tuple `X=(x,ξ)`,
i.e. the value  of another tangent bundle point.
"""
struct TBPoint{P<:MPoint, T <: TVector} <: MPoint
  value::Tuple{P,T}
  TBPoint{P,T}(x::P, ξ::T) where {P <: MPoint, T <: TVector} = new( (x,ξ) )
  TBPoint{P,T}(X::Tuple{P, T}) where {P <: MPoint, T <: TVector} = new( X )
end
TBPoint(x::P, ξ::T) where {P <: MPoint, T <: TVector} = TBPoint{P,T}(x,ξ)
TBPoint(X::Tuple{P,T}) where {P <: MPoint, T <: TVector}= TBPoint{P,T}( X )
@doc doc"""
    getValue(X)
return the value of the [`TBPoint`](@ref) `X`, i.e. the Tuple of a
[`MPoint`](@ref) and its [`TVector`](@ref).
"""
getValue(X::TBPoint{P,T}) where { P <: MPoint, T <: TVector } = X.value
@doc doc"""
    getBase(X)
return the base of the [`TBPoint`](@ref)` X`, i.e. its [`MPoint`](@ref).
"""
getBase(X::TBPoint{P,T}) where { P <: MPoint, T <: TVector } = X.value[1]
# passthrough for Extendeds
getBase(X::MPointE{TBPoint{P,T}}) where { P <: MPoint, T <: TVector } = MPointE( getBase(strip(X)) )

@doc doc"""
    getTangent(X)
return the tangent of the [`TBPoint`](@ref)` X`, i.e. the its [`TVector`](@ref).
"""
getTangent(X::TBPoint{P,T}) where { P <: MPoint, T <: TVector } = X.value[2]
#passthrough for Extendeds
getTangent(X::MPointE{TBPoint{P,T}}) where { P <: MPoint, T <: TVector } = TVectorE( getTangent(strip(X)), getBase(strip(X) ) )

@doc doc"""
    TBTVector <: TVector

A tangent vector $\Xi \in T_X\mathcal M$ on the manifold
$\mathcal M = T\mathcal N$ for the (base) manifold $\mathcal N$.
Both tangent components can be represented by elements from the base point $x$
from within $X=(x,\xi)$. Both components are from the same space since
$TT_x\mathcal N= T_x\mathcal N$, hence the tangent vector is a tuple
$(\xi\,\nu)\in T_x\mathcal N\times T_x\mathcal N$.
As for the [`TBPoint`](@ref) two constructors are available, one for stwo seperate
tangent vectors, one for a tuple of two tangent vectors.
"""
struct TBTVector{T <: TVector} <: TVector
  value::Tuple{T,T}
  TBTVector{T}(ξ::T,ν::T) where {T <: TVector} = new( (ξ,ν) )
  TBTVector{T}(Ξ::Tuple{T,T}) where {T <: TVector} = new( Ξ )
end
TBTVector(ξ::T, ν::T) where {T <: TVector} = TBTVector{T}( ξ,ν )
TBTVector(Ξ::Tuple{T,T}) where {T <:TVector}= TBTVector{T}( Ξ )
@doc doc"""
    getValue(Ξ)
return the Tuple contained in the [`TBTVector`](@ref)` Ξ`, i.e. its tuple of two
[`TVector`](@ref)s.
"""
getValue(Ξ::TBTVector) = Ξ.value;

@doc doc"""
    getBase(Ξ)
return the base of the [`TBTVector`](@ref)` Ξ`, i.e. its first [`TVector`](@ref).
"""
getBase(Ξ::TBTVector) = Ξ.value[1]
# passthrough for Extendeds
getBase(Ξ::TVectorE{TBTVector{T}}) where {T <: TVector} = TVectorE( getBase(strip(Ξ)), getBase(getBasePoint(Ξ)) )

@doc doc"""
    getTangent(Ξ)
return the tangent of the [`TBTVector`](@ref)` Ξ`, i.e. its second
[`TBTVector`](@ref).
"""
getTangent(Ξ::TBTVector) = Ξ.value[2]
getTangent(Ξ::TVectorE{TBTVector{T}}) where {T <: TVector}= TVectorE( getTangent(strip(Ξ)), getBase(getBasePoint(Ξ)) )

*(t::Number, Ξ::TBTVector) = TBTVector(t*getBase(Ξ),t*getTangent(Ξ))

+(Ξ::TBTVector, Η::TBTVector) = TBTVector(getBase(Ξ) + getBase(Η), getTangent(Ξ) + getTangent(Η) )
-(Ξ::TBTVector, Η::TBTVector) = TBTVector(getBase(Ξ) - getBase(Η), getTangent(Ξ) - getTangent(Η) )
-(Ξ::TBTVector) = TBTVector(- getBase(Ξ), -getTangent(Ξ) )
+(Ξ::TBTVector) = TBTVector( getBase(Ξ), getTangent(Ξ) )

# Functions
# ---
@doc doc"""
    distance(M,X,Y)

Compute the Riemannian distance on $\mathcal M=T\mathcal N$ by employing the
distance on the manifold for the base component and the vector norm on the
tangent space, and then take the Eucklidean Norm of the vector from
$\mathbb R^2$.
"""
distance(M::TangentBundle,X::TBPoint, Y::TBPoint) = sqrt(
  distance(M.manifold,getBase(X),getBase(Y))^2 + norm(M.manifold,getBase(X),getTangent(X)-getTangent(Y))^2
)

@doc doc"""
    dot(M,X,Ξ,Η)

Compute the Riemannian inner product for two [`TBTVector`](@ref)s `Ξ` and `Η`
from $T_X\mathcal M$ of the [`TangentBundle`](@ref)` M = TN` given by
the sum of the two inner products of the tangent vector components
"""
dot(M::TangentBundle, X::TBPoint, Ξ::TBTVector, Η::TBTVector) =
dot(M.manifold,getBase(X),getBase(Ξ),getBase(Η)) + dot(M.manifold,getBase(X),getTangent(Ξ),getTangent(Η))

@doc doc"""
    exp(M,X,Ξ[, t=1.0])

Compute the exponential map on the [`TangentBundle`](@ref) `M`$=T\mathcal N$ with
respect to the [`TBPoint`](@ref)` X=(x,ξ)` and the [`TBTVector`](@ref)` Ξ=(Ξx,Ξξ)`,
which consists of the exponential map in the first component (`exp(x,Ξx,t)` and
a (scaled) addition in the second (`ξ + tΞξ`) in the second component followed
by a parallel transport.
"""
function exp(M::TangentBundle,X::TBPoint,Ξ::TBTVector,t::Float64=1.0)
  x = exp(M.manifold,getBase(X),getBase(Ξ),t)
  ξ = parallelTransport(M.manifold, getBase(X), x,
    getTangent(X) + t*getTangent(Ξ)
  )
  return TBPoint(x,ξ)
end
@doc doc"""
    log(M,X,Y)

Compute the logarithmic map on the [`TangentBundle`](@ref)
$\mathcal M=T\mathcal N$, i.e. the `log` for the base manifold component and
a parallel transport and a minus for the tangent components.

"""
function log(M::TangentBundle,X::TBPoint,Y::TBPoint)
  Ξx = log(M.manifold, getBase(X), getBase(Y) )
  Ξξ = getTangent(X) - parallelTransport(M.manifold, getBase(X), getBase(Y), getTangent(Y))
  return TBTVector(Ξx,Ξξ)
end
@doc doc"""
    manifoldDimension(X)

returns the dimension of the [`TangentBundle`](@ref) `M`$=T\mathcal N$ to which `X`
bvelongs, which is twice the dimension of the base manifold.
"""
manifoldDimension(x::TBPoint)::Integer = 2*manifoldDimension(getBase(x))
@doc doc"""
    manifoldDimension(M)

returns the dimension of the [`TangentBundle`](@ref) `M`$=T\mathcal N$, i.e.,
twice the dimension of the base manifold `N`.
"""
manifoldDimension(M::TangentBundle)::Integer = M.dimension
@doc doc"""
    norm(M,X,Ξ)

Computes the norm of the [`TBTVector`](@ref)` Ξ` in the tangent space
$T_x\mathcal M$ at [`TBPoint`](@ref)` X` of the [`TangentBundle`](@ref) `M`.
"""
norm(M::TangentBundle, X::TBPoint, Ξ::TBTVector) = sqrt( dot(M,X,Ξ,Ξ) )
@doc doc"""
    parallelTransport(M,X,Y,Ξ)

Compute the paralllel transport of the [`TBTVector`](@ref)` Ξ` from
the tangent space $T_X\mathcal M$ at [`TBPoint`](@ref)` X` to
$T_Y\mathcal M$ at [`TBPoint`](@ref)` Y` on the [`TangentBundle`](@ref) `M`
provided that the corresponding [`geodesic`](@ref) $g(\cdot;x,y)$ is unique.
Then both components of $\Xi=(\Xi_x,\Xi_\xi)$ are parallely transported using the
parallel transport of the underlying base manifold.
"""
parallelTransport(M::TangentBundle, X::TBPoint, Y::TBPoint, Ξ::TBTVector) = TBTVector(
  parallelTransport(M.manifold, getBase(X), getBase(Y), getBase(Ξ) ),
  parallelTransport(M.manifold, getBase(X), getBase(Y), getTangent(Ξ) )
)
@doc doc"""
    randomMPoint(M)

returns a random point on the [`TangentBundle`](@ref) `M` by producing a
[`randomMPoint`](@ref)  random point on the base manifold and
[`randomTVector`](@ref) in the correspoinding tangent plane.
"""
function randomMPoint(M::TangentBundle)::TBPoint
  x = randomMPoint(M.manifold)
  ξ = randomTVector(M.manifold,x)
	return TBPoint(x,ξ)
end
@doc doc"""
    randomTVector(M,x)

returns a random tangent vector the [`TangentBundle`](@ref) `M` by producing
two [`randomTVector`](@ref)s in the correspoinding tangent plane of the
[`getBase`](@ref) of the [`TBPoint`](@ref) `x`.
"""
function randomTVector(M::TangentBundle, x::TBPoint)
  xξ = randomTVector(M.manifold,getBase(x))
  ξξ = randomTVector(M.manifold,getBase(x))
	return TBTVector(xξ,ξξ)
end
@doc doc"""
    tangentONB(M,X,Y)

constructs a tangent ONB in the tangent space of the [`TBPoint`](@ref)` X` on
the [`TangentBundle`](@ref) `M`, where $\log_XY$ is the first component.
"""
tangentONB(M::TangentBundle, X::TBPoint, Y::TBPoint) = tangentONB(M,X,log(M,X,Y))
@doc doc"""
    Η,κ = tangentONB(M,X,Ξ)

constructs a tangent ONB in $T_X\mathcal M$, i.e. in the tangent space of the
[`TBPoint`](@ref) `x` on the [`TangentBundle`](@ref) `M`
whose first vector is given by the [`TBTVector`](@ref)` Ξ`.
It is constructed by using twice the tangent ONB of the base manifold.
"""
function tangentONB(M::TangentBundle,X::TBPoint,Ξ::TBTVector)
  BaseONB,κx = tangentONB(M.manifold, getBase(X), getBase(Ξ))
  TangentONB,κξ = tangentONB(M.manifold, getBase(X), getTangent(Ξ) )
  baseTVs = [BaseONB..., [zeroTVector(M.manifold,getBase(X)) for i=1:length(TangentONB)]...]
  tangentTVs = [ [zeroTVector(M.manifold,getBase(X)) for i=1:length(BaseONB)]..., TangentONB...]
  return TBTVector.(baseTVs, tangentTVs), [κx...,κξ...]
end

typeofTVector(::Type{TBPoint{P,T}})  where {P <: MPoint, T <: TVector} = TBTVector{typeofTVector(P)}
typeofMPoint(::Type{TBTVector{T}}) where {T <: TVector} = TBPoint{typeofMPoint(T),T}

@doc doc"""
    typicalDistance(M)

returns the typical distance on the [`TangentBundle`](@ref) `M`, i.e. for
$\mathcal M = T\mathcal N$ we
obtain $t_{\mathcal M} = \sqrt{t_{\mathcal N}^2 + d_{\mathcal N}^2}$, where $d$
denotes the manifold dimension.
"""
typicalDistance(M::TangentBundle) = sqrt(
  typicalDistance(M.manifold)^2 + manifoldDimension(M.manifold)
)
@doc doc"""
    validateMPoint(M,X)

validate that the [`TBPoint`](@ref)` X` is a valid point on the
[`TangentBundle`](@ref) `M`, i.e. the first component is a point on the base
manifold and the second a tangent vector is the tangent space of the first
"""
function validateMPoint(M::TangentBundle,X::TBPoint)
  return validateMPoint(M.manifold,getBase(X)) && validateTVector(M.manifold,getBase(X),getTangent(X))
  return true
end
@doc doc"""
    validateTVector(M,X,Ξ)

validate, that the [`TBTVector`](@ref)` Ξ` is a valid tangent vector in the
tangent space of the [`TBPoint`](@ref)` X` on the [`TangentBundle`](@ref) `M`,
i.e. both components of `Ξ` are tangent vectors in the tangent space of the base
component of `X`, since the tangent space of the tangent space is represented as
the tangent space itself.
"""
function validateTVector(M::TangentBundle, X::TBPoint, Ξ::TBTVector)
  # both components have to be tangent vectors, since TM = TTM
  return validateTVector(M.manifold, getBase(X), getTangent(Ξ)) && validateTVector(M.manifold, getBase(X), getTangent(Ξ))
end
@doc doc"""
    zeroTVector(M,X)

returns a zero vector in the tangent space $T_X\mathcal M$ of the
[`TangentBundle`](@ref) $X=(x,ξ)\in T\mathcal N$ by creating two zero vectors in $T_x\mathcal M$.
"""
zeroTVector(M::TangentBundle{Mt}, X::TBPoint{P}) where {Mt <: Manifold, P <: MPoint} =
TBTVector{typeofTVector(P)}(
  zeroTVector(M.manifold,getBase(X)),
  zeroTVector(M.manifold,getBase(X))
)
# Display
# ---
show(io::IO, M::TangentBundle) = print(io,"The Tangent bundle of <$(M.manifold)>")
show(io::IO, X::TBPoint) = print(io, "TB($(getBase(X)), $(getTangent(X)))")
show(io::IO, Ξ::TBTVector) = print(io, "TBT($( getBase(Ξ) ), $( getTangent( Ξ )))")


text:
using BHAtp, Statistics

ProjDir = @__DIR__
!isdir(joinpath(ProjDir, "plots")) && mkdir(joinpath(ProjDir, "plots"))

ProjName = split(ProjDir, "/")[end]
bhaj = BHAtp.BHAJ(ProjName, ProjDir)
bhaj.ratio = 0.5

segs = [
  # Element type,  Material,    Length,     ID,   OD,    fc
  (:bit,           :steel,        0.00,   2.000, 17.500,  0.0),
  (:collar,        :steel,        1.44,   2.000,  9.500,  0.0),
  (:collar,        :steel,        3.00,   3.000,  9.500,  0.0),
  (:collar,        :steel,        1.24,   2.938,  9.500,  0.0),
  (:stabilizer,    :steel,        0.00,   2.938, 17.500,  0.0),
  (:collar,        :steel,        4.04,   2.938,  9.500,  0.0),
  (:collar,        :monel,       37.45,   2.875,  9.500,  0.0),
  (:collar,        :monel,       27.09,   3.000,  9.625,  0.0),
  (:collar,        :monel,       30.01,   2.875,  9.500,  0.0),
  (:collar,        :steel,        1.87,   3.000,  9.500,  0.0),
  (:stabilizer,    :steel,        0.00,   3.000, 17.500,  0.0),
  (:collar,        :steel,        3.73,   3.000,  9.500,  0.0),
  (:collar,        :steel,       28.73,   2.875,  9.500,  0.0),
  (:collar,        :steel,        4.22,   3.000,  9.125,  0.0),
  (:collar,        :steel,       27.98,   2.750,  7.875,  0.0),
  (:collar,        :steel,       29.45,   2.875,  7.563,  0.0),
  (:collar,        :steel,       30.87,   2.875,  7.625,  0.0),
  (:collar,        :steel,       31.10,   2.750,  8.000,  0.0),
  (:collar,        :steel,       32.18,   3.000,  8.000,  0.0),
  (:collar,        :steel,       30.88,   2.750,  7.938,  0.0),
  (:collar,        :steel,        3.50,   2.750,  8.000,  0.0)
]

traj = [
  #Heading,  Diameter
  ( 60.0,      17.5)

]

wobs = 30:5:40
incls = 30:5:60  # Or e.g. incls = [5 10 20 30 40 45 50]

@time bhaj(segs, traj, wobs, incls)

display("Fetch tp=false, wob=35, incl@bit=30 solution:")
df,df_tp = show_solution(ProjDir, 35, 30, show=false, tp=false)
df[:,[2; 5:6; 9:12]] |> display
println()

println("Fetch tp=true, wob=35, incl@bit=30 solution:")
df,df_tp = show_solution(ProjDir, 35, 30, show=false);
df[:,[2; 5:6; 9:12]] |> display
println()

tpf = create_final_tp_df(ProjDir, wobs, incls; ofu=false)
tpf |> display
println()

name = "figure_03a"

figure_03a = @pgf GroupPlot({
    group_style = {
      group_size = "1 by 2",
      vertical_sep = "0pt",
      xticklabels_at = "edge bottom"
    }
  },
  Axis({
      height = "6cm", width = "8cm",
      ylabel = "TP [°/100ft]", title = "TP and Formatioin Correction (ΔTP)"
    },
    Plot({ color="red", mark = "*" }, Table(incls, tpf.tp[1:7])),
    LegendEntry("WOB 30"),
    Plot({ color="blue", mark = "+" }, Table(incls, tpf.tp[8:14])),
    LegendEntry("WOB 35"),
    Plot({ color="green", mark = "o" }, Table(incls, tpf.tp[15:21])),
    LegendEntry("WOB 40")
  ),
  Axis({
      height = "6cm", width = "8cm",
      ylabel = "Sideforce NBS"
    },
    Plot({ color="red", mark = "*", title ="TP" }, Table(incls, tpf.ft[1:7])),
    LegendEntry("WOB 30"),
    Plot({ color="blue", mark = "+" }, Table(incls, tpf.ft[8:14])),
    LegendEntry("WOB 35"),
    Plot({ color="green", mark = "o" }, Table(incls, tpf.ft[15:21])),
    LegendEntry("WOB 40")
  )
)

savefig(joinpath(ProjDir, "plots", name), figure_03a, ProjDir)


text:
module IBAN

export
    ## methods
    is_valid_iban,
    split_iban

include("iban-constants.jl")
include("iban-validation.jl")

end # module


text:
using Oceananigans.TurbulenceClosures: AnisotropicMinimumDissipation
using Oceananigans.TimeSteppers: update_state!

function run_ocean_large_eddy_simulation_regression_test(arch, grid_type, closure)
    name = "ocean_large_eddy_simulation_" * string(typeof(closure).name.wrapper)

    spinup_steps = 10000
      test_steps = 10
              Δt = 2.0

    # Parameters
      Qᵀ = 5e-5     # Temperature flux at surface
      Qᵘ = -2e-5    # Velocity flux at surface
    ∂T∂z = 0.005    # Initial vertical temperature gradient

    # Grid
    N = L = 16
    if grid_type == :regular
        grid = RegularRectilinearGrid(size=(N, N, N), extent=(L, L, L))
    elseif grid_type == :vertically_unstretched
        zF = range(-L, 0, length=N+1)
        grid = VerticallyStretchedRectilinearGrid(architecture=arch, size=(N, N, N), x=(0, L), y=(0, L), z_faces=zF)
    end

    # Boundary conditions
    u_bcs = UVelocityBoundaryConditions(grid, top = BoundaryCondition(Flux, Qᵘ))
    T_bcs = TracerBoundaryConditions(grid, top = BoundaryCondition(Flux, Qᵀ),
                                        bottom = BoundaryCondition(Gradient, ∂T∂z))
    S_bcs = TracerBoundaryConditions(grid, top = BoundaryCondition(Flux, 5e-8))

    # Model instantiation
    model = IncompressibleModel(
             architecture = arch,
                     grid = grid,
                 coriolis = FPlane(f=1e-4),
                 buoyancy = Buoyancy(model=SeawaterBuoyancy(equation_of_state=LinearEquationOfState(α=2e-4, β=8e-4))),
                  closure = closure,
      boundary_conditions = (u=u_bcs, T=T_bcs, S=S_bcs)
    )

    # We will manually change the stop_iteration as needed.
    simulation = Simulation(model, Δt=Δt, stop_iteration=0)

    # The type of the underlying data, not the offset array.
    ArrayType = typeof(model.velocities.u.data.parent)

    ####
    #### Uncomment the block below to generate regression data.
    ####

    #=
    @warn "Generating new data for the ocean LES regression test."

    # Initialize model: random noise damped at top and bottom
    Ξ(z) = randn() * z / model.grid.Lz * (1 + z / model.grid.Lz) # noise
    T₀(x, y, z) = 20 + ∂T∂z * z + ∂T∂z * model.grid.Lz * 1e-2 * Ξ(z)
    u₀(x, y, z) = sqrt(abs(Qᵘ)) * 1e-3 * Ξ(z)
    set!(model, u=u₀, w=u₀, T=T₀, S=35)

    simulation.stop_iteration = spinup_steps-test_steps
    run!(simulation)

    checkpointer = Checkpointer(model, schedule = IterationInterval(test_steps), prefix = name,
                                dir = joinpath(dirname(@__FILE__), "data"))

    simulation.output_writers[:checkpointer] = checkpointer

    simulation.stop_iteration += 2test_steps
    run!(simulation)
    pop!(simulation.output_writers, :checkpointer)
    =#

    ####
    #### Regression test
    ####

    initial_filename = joinpath(dirname(@__FILE__), "data", name * "_iteration$spinup_steps.jld2")

    solution₀, Gⁿ₀, G⁻₀ = get_fields_from_checkpoint(initial_filename)

    Nz = grid.Nz

    model.velocities.u.data.parent .= ArrayType(solution₀.u)
    model.velocities.v.data.parent .= ArrayType(solution₀.v)
    model.velocities.w.data.parent .= ArrayType(solution₀.w)
    model.tracers.T.data.parent    .= ArrayType(solution₀.T)
    model.tracers.S.data.parent    .= ArrayType(solution₀.S)

    model.timestepper.Gⁿ.u.data.parent .= ArrayType(Gⁿ₀.u)
    model.timestepper.Gⁿ.v.data.parent .= ArrayType(Gⁿ₀.v)
    model.timestepper.Gⁿ.w.data.parent .= ArrayType(Gⁿ₀.w)
    model.timestepper.Gⁿ.T.data.parent .= ArrayType(Gⁿ₀.T)
    model.timestepper.Gⁿ.S.data.parent .= ArrayType(Gⁿ₀.S)

    model.timestepper.G⁻.u.data.parent .= ArrayType(G⁻₀.u)
    model.timestepper.G⁻.v.data.parent .= ArrayType(G⁻₀.v)
    model.timestepper.G⁻.w.data.parent .= ArrayType(G⁻₀.w)
    model.timestepper.G⁻.T.data.parent .= ArrayType(G⁻₀.T)
    model.timestepper.G⁻.S.data.parent .= ArrayType(G⁻₀.S)

    model.clock.time = spinup_steps * Δt
    model.clock.iteration = spinup_steps

    update_state!(model)

    for n in 1:test_steps
        time_step!(model, Δt, euler=false)
    end

    final_filename = joinpath(@__DIR__, "data", name * "_iteration$(spinup_steps+test_steps).jld2")

    solution₁, Gⁿ₁, G⁻₁ = get_fields_from_checkpoint(final_filename)

    test_fields = CUDA.@allowscalar (u = Array(interior(model.velocities.u)),
                                     v = Array(interior(model.velocities.v)),
                                     w = Array(interior(model.velocities.w)[:, :, 1:Nz]),
                                     T = Array(interior(model.tracers.T)),
                                     S = Array(interior(model.tracers.S)))

    correct_fields = (u = Array(interior(solution₁.u, model.grid)),
                      v = Array(interior(solution₁.v, model.grid)),
                      w = Array(interior(solution₁.w, model.grid)),
                      T = Array(interior(solution₁.T, model.grid)),
                      S = Array(interior(solution₁.S, model.grid)))

    summarize_regression_test(test_fields, correct_fields)

    @test all(test_fields.u .≈ correct_fields.u)
    @test all(test_fields.v .≈ correct_fields.v)
    @test all(test_fields.w .≈ correct_fields.w)
    @test all(test_fields.T .≈ correct_fields.T)
    @test all(test_fields.S .≈ correct_fields.S)

    return nothing
end


text:
module Temp
  function rema_!(shadow::Float64,variable::Float64,decay::Float64= 0.3)
      shadow -= (1 - decay) * (shadow - variable)
  end

  function smooth(x,decay=0.9)
    res = copy(x)
    st = x[1]
    res[1] = st
    i = 2
    for el in x[2:end]
      st = rema_!(st,el,decay)
      res[i] = st
      i+=1
    end
    res
  end

end


text:
module MultipleShootingPEM

using LossFunctions
using PyCall

const scipy_opt = PyNULL()
const scipy_sps = PyNULL()
const scipy_spslin = PyNULL()

function __init__()
    copy!(scipy_opt, pyimport_conda("scipy.optimize", "scipy"))
    copy!(scipy_sps, pyimport_conda("scipy.sparse", "scipy"))
    copy!(scipy_spslin, pyimport_conda("scipy.sparse.linalg", "scipy"))
end

include("util.jl")
include("sensitivity_equations.jl")
include("space_state_model.jl")
include("one_shoot_simulation.jl")
include("multiple_shooting.jl")
include("optimization_problem.jl")

end # module


text:
function ribbonplot(;x::AbstractVector = Int[],
                     ylow::AbstractVector = Int[],
                     yhigh::AbstractVector = Int[],
                     group::AbstractVector = Int[])

  v = barplot(x = x, y = yhigh, y2 = ylow, group = group)

  tablename = v.data[1].name

  innermark = VegaMark()
  innermark._type = "area"
  innermark.properties = VegaMarkProperties()
  innermark.properties.enter = VegaMarkPropertySet()
  innermark.properties.enter.x           = VegaValueRef(scale = "x", field = "x")
  innermark.properties.enter.y           = VegaValueRef(scale = "y", field = "y")
  innermark.properties.enter.y2          = VegaValueRef(scale = "y", field = "y2")
  innermark.properties.enter.fill        = VegaValueRef(scale = "group", field = "group")
  innermark.properties.enter.interpolate = VegaValueRef(value="monotone")

  mark = VegaMark()
  mark._type = "group"
  mark.from = VegaMarkFrom()
  mark.from.data      = tablename
  mark.from.transform = [VegaTransform(Dict{Any,Any}("type"=>"facet", "groupby" => ["group"]))]
  mark.marks = [innermark;]

  v.marks = [mark;]

  v.scales[1]._type = "linear"
  v.scales[1].zero = false
  v.scales[2].domain = VegaDataRef(data = tablename, field = ["y", "y2"])

  return v
end

max_stars_repo_name: jpthor/vega


text:
using LinearAlgebra: Symmetric, isposdef
using Distributions: MvNormal, logpdf

import ..kalman_filter_mv
import ..kalman_smoother_mv

abstract type AbstractLinearGaussian end;

struct LinearGaussian{T<:Any} <: AbstractLinearGaussian
    θ::T
end

function LinearGaussian(; kwargs...)
    LinearGaussian{typeof(kwargs.data)}(kwargs.data)
end

function Base.show(io::IO, ::LinearGaussian)
    print(io, "LinearGaussian()")
end

# Generative Model
function dpr_base(model::AbstractLinearGaussian, Zs)
    logpdf(MvNormal(model.θ.μ0, model.θ.Tau0), Zs)
end

function rpr_base(model::AbstractLinearGaussian, rng)
    rand(
        rng,
        MvNormal(model.θ.states[:, 1], Symmetric(model.θ.variances[:, :, 1])),
        model.θ.K,
    )
end

function dt_base(m::AbstractLinearGaussian, Zs, Zs_new, t)
    μ = m.θ.G(t) * Zs
    Σ = m.θ.Tau(t)
    d = logpdf(MvNormal(Σ), Zs_new - μ)
    @assert !any(isnan.(d))
    return d
end

function dm_base(m::AbstractLinearGaussian, Zs, t)
    yt = view(m.θ.Y, t, :)
    if any(ismissing.(yt))
        return zeros(eltype(Zs), size(Zs, 2))
    else
        d = logpdf(MvNormal(m.θ.Σ(t)), m.θ.F(t) * Zs .- yt)
        @assert !any(isnan.(d))
        return d
    end
end

# Inference Model
function rinit_base(m::AbstractLinearGaussian, rng)
    rand(rng, MvNormal(m.θ.states_sm[:, 1], Symmetric(m.θ.variances_sm[:, :, 1])), m.θ.K)
end

function dinit_base(m::AbstractLinearGaussian, Zs)
    logpdf(MvNormal(m.θ.states_sm[:, 1], Symmetric(m.θ.variances_sm[:, :, 1])), Zs)
end

function calc_prop(m::AbstractLinearGaussian, Zs, t)
    μ =
        m.θ.states_sm[:, t] .+
        (m.θ.Covs_sm[:, :, t-1] / m.θ.variances_sm[:, :, t-1]) *
        (Zs .- m.θ.states_sm[:, t-1])
    Σ = Symmetric(
        m.θ.variances_sm[:, :, t] -
        (m.θ.Covs_sm[:, :, t-1] / m.θ.variances_sm[:, :, t-1]) * m.θ.Covs_sm[:, :, t-1]',
    )

    @assert !any(isnan.(μ))
    @assert !any(isnan.(Σ))

    return μ, Σ
end

function rp_base(m::AbstractLinearGaussian, rng, Zs, t)
    μ, Σ = calc_prop(m, Zs, t)
    @assert isposdef(Σ)
    Zs_new = μ + rand(rng, MvNormal(Σ), size(Zs, 2))
    return Zs_new
end

function dp_base(m::AbstractLinearGaussian, Zs, Zs_new, t)
    μ, Σ = calc_prop(m, Zs, t)
    @assert isposdef(Σ)
    d = logpdf(MvNormal(Σ), Zs_new - μ)
    return d
end

function dpre_base(m::AbstractLinearGaussian, Zs, t)
    d_smooth =
        logpdf(MvNormal(m.θ.states_sm[:, t-1], Symmetric(m.θ.variances_sm[:, :, t-1])), Zs)
    d_filt = logpdf(MvNormal(m.θ.states[:, t-1], Symmetric(m.θ.variances[:, :, t-1])), Zs)

    d = (d_smooth - d_filt) .+ sum(@view(m.θ.lls[t:end]))
    @assert !any(isnan.(d))
    return d
end

function make_kalman_proposal(F0, F, G0, G, Tau, Σ, μ0, Tau0, Y)
    ll, states, variances, lls = kalman_filter_mv(F0, F, G0, G, Σ, Tau, μ0, Tau0, Y)
    states_sm, variances_sm, Covs_sm = kalman_smoother_mv(F0, F, G0, G, Σ, Tau, μ0, Tau0, Y)
    return @NT(ll, states, variances, lls, states_sm, variances_sm, Covs_sm)
end

function lg_smc_model(K, F0, F, G0, G, Tau, Σ, μ0, Tau0, Y)
    ll, states, variances, lls, states_sm, variances_sm, Covs_sm =
        make_kalman_proposal(F0, F, G0, G, Tau, Σ, μ0, Tau0, Y)

    m = LinearGaussian(
        F = F,
        G = G,
        Tau = Tau,
        Σ = Σ,
        μ0 = μ0,
        Tau0 = Tau0,
        K = K,
        states = states,
        variances = variances,
        lls = lls,
        states_sm = states_sm,
        variances_sm = variances_sm,
        Covs_sm = Covs_sm,
        Y = Y,
    )

    dpr(Zs) = dpr_base(m, Zs)
    dt(Zs, Zs_new, t) = dt_base(m, Zs, Zs_new, t)
    dm(Zs, t) = dm_base(m, Zs, t)

    rinit(rng) = rinit_base(m, rng)
    dinit(Zs) = dinit_base(m, Zs)
    rp(rng, Zs, t) = rp_base(m, rng, Zs, t)
    dp(Zs, Zs_new, t) = dp_base(m, Zs, Zs_new, t)

    dpre(Zs, t) = dpre_base(m, Zs, t)

    return @NT(rinit, rp, dinit, dpr, dp, dt, dm, dpre)
end


text:
# Переопределение сторон
Left  = West
Up    = Nord
Right = Ost
Down  = Sud

"""
    inverse(side::HorizonSide)

-- Возвращает направление, противоположное заданному
"""
inverse(side::HorizonSide) = HorizonSide(mod(Int(side)+2, 4))

"""
    left(side::HorizonSide)

-- Возвращает направление, следующее против часовой стрелке, по отношению к заданному
"""
left(side::HorizonSide) =  HorizonSide(mod(Int(side)+1, 4))

"""
    right(side::HorizonSide)

-- Возвращает направление, следующее по часовой стрелке, по отношению к заданному
"""
right(side::HorizonSide) = HorizonSide(mod(Int(side)-1, 4))

"""
    putmarkers!(r::Robot, side::HorizonSide)

-- Ставит маркеры в указанном направлении, пока не наткнется на ограду
"""
function putmarkers!(r::Robot, side::HorizonSide)
    putmarker!(r)
    while isborder(r,side)==false
        move!(r,side)
        putmarker!(r)
    end
end

"""
    putmarkers!(r::Robot, side::HorizonSide, count::Int)

-- Ставит маркеры в указанном направлении, пока не обнулится счетчик
"""
function putmarkers!(r::Robot, side::HorizonSide, count::Int)
    putmarker!(r)
    while count!=0
        move!(r,side)
        putmarker!(r)
        count -= 1
    end
end

"""
    move_by_markers(r::Robot, side::HorizonSide)

-- Двигается по маркерам в указанном направлении
"""
move_by_markers(r::Robot, side::HorizonSide) = while ismarker(r)==true move!(r,side) end

"""
    movements!(r::Robot, side::HorizonSide)

-- Перемещает Робота в заданном направлении до стенки
"""
movements!(r::Robot, side::HorizonSide) = while isborder(r,side)==false move!(r,side) end

"""
    movements!(r::Robot, side::HorizonSide, num_steps::Int)

-- Перемещает Робота в заданном направлении на заданное число шагов
"""
movements!(r::Robot, side::HorizonSide, num_steps::Int) = for _ in 1:num_steps move!(r,side) end

"""
    get_num_movements!(r::Robot, side::HorizonSide)

-- Перемещает Робота в заданном направлении до стенки и возвращает сделанное число шагов
"""
function get_num_movements!(r::Robot, side::HorizonSide)
    num_steps = 0
    while isborder(r,side)==false
        move!(r,side)
        num_steps += 1
    end
    return num_steps
end

"""
    movements!(r, sides::NTuple{2,HorizonSide}, num_steps::Vector{Any})

-- перемещает Робота по пути, представленного двумя последовательностями, sides и num_steps
-- sides - содержит последовательность направлений перемещений
-- num_steps - содержит последовательность чисел шагов в каждом из этих направлений, соответственно; при этом, если длина последовательности sides меньше длины последовательности num_steps, то предполагается, что последовательность sides должна быть продолжена периодически
"""
function movements!(r, sides::NTuple{2,HorizonSide}, num_steps::Vector{Any})
    for (i,n) in enumerate(num_steps)
        movements!(r, sides[mod(i-1, length(sides))+1], n)
    end
end

"""
    through_rectangles_into_angle(r,angle::NTuple{2,HorizonSide})

-- Перемещает Робота в заданный угол, прокладывая путь межу внутренними прямоугольными перегородками и возвращает массив, содержащий числа шагов в каждом из заданных направлений на этом пути
"""
function through_rectangles_into_angle(r,angle::NTuple{2,HorizonSide})
    num_steps = []
    while (isborder(r,angle[1]) == false || isborder(r,angle[2]) == false)
        push!(num_steps, get_num_movements!(r, angle[2]))
        push!(num_steps, get_num_movements!(r, angle[1]))
    end
    return num_steps
end

"""
    move_if_possible!(r::Robot, side::HorizonSide)

-- Двигается в заданном направлении, пока возможно, минуя внутренние перегородки
"""
function move_if_possible!(r::Robot, direct_side::HorizonSide)::Bool
    left_side = right(direct_side)
    right_side = inverse(left_side)
    num_of_steps = 0

    if isborder(r, direct_side) == false
        move!(r, direct_side)
        result = true
    else
        while isborder(r, direct_side) == true
            if isborder(r, left_side) == false
                move!(r, left_side)
                num_of_steps += 1
            else
                break
            end
        end
        if isborder(r, direct_side) == false
            move!(r, direct_side)
            while isborder(r, right_side) == true
                move!(r, direct_side)
            end
            result = true
        else
            result = false
        end
        while num_of_steps > 0
            num_of_steps -= 1
            move!(r, right_side)
        end
    end

    return result
end

"""
    movements_if_possible!(r::Robot, side::HorizonSide)

-- Перемещает Робота в заданном направлении до стенки, минуя внутренние перегородки
"""
movements_if_possible!(r::Robot, side::HorizonSide) = while isborder(r,side)==false move_if_possible!(r,side) end

"""
    movements_if_possible!(r::Robot, side::HorizonSide)

-- Перемещает Робота в заданном направлении на определенное количетво шагов, минуя внутренние перегородки
"""
movements_if_possible!(r::Robot, side::HorizonSide, num_steps::Int) = for _ in 1:num_steps move_if_possible!(r,side) end

"""
    movements_if_possible!(r, sides::NTuple{2,HorizonSide}, num_steps::Vector{Any})

-- перемещает Робота по пути, представленного двумя последовательностями, sides и num_steps
-- sides - содержит последовательность направлений перемещений
-- num_steps - содержит последовательность чисел шагов в каждом из этих направлений, соответственно; при этом, если длина последовательности sides меньше длины последовательности num_steps, то предполагается, что последовательность sides должна быть продолжена периодически
"""
function movements_if_possible!(r, sides::NTuple{2,HorizonSide}, num_steps::Vector{Any})
    for (i,n) in enumerate(num_steps)
        movements_if_possible!(r, sides[mod(i-1, length(sides))+1], n)
    end
end

"""
    putmarkers_if_possible!(r::Robot, side::HorizonSide)

-- Ставит маркеры в указанном направлении, пока не наткнется на ограду, минуя внутренние перегородки
"""
function putmarkers_if_possible!(r::Robot, side::HorizonSide)
    putmarker!(r)
    while move_if_possible!(r, side) == true
        putmarker!(r)
    end
end

"""
    putmarkers_if_possible!(r::Robot, side::HorizonSide, count::Int)

-- Ставит маркеры в указанном направлении, минуя внутренние перегородки, пока не обнулится счетчик
"""
function putmarkers_if_possible!(r::Robot, side::HorizonSide, count::Int)
    putmarker!(r)
    while count!=0
        move_if_possible!(r, side)
        putmarker!(r)
        count -= 1
    end
end

"""
    get_num_movements_if_possible!(r::Robot, side::HorizonSide)

-- Перемещает Робота в заданном направлении до стенки, минуя внутренние перегородки, и возвращает сделанное число шагов
"""
function get_num_movements_if_possible!(r::Robot, side::HorizonSide)
    num_steps = 0
    while move_if_possible!(r, side) == true
        num_steps += 1
    end
    return num_steps
end

"""
    mark_and_enumerate_if_possible!(r::Robot, side::HorizonSide)

--  Ставит маркеры в указанном направлении, пока не наткнется на ограду, не считая внутренние перегородки, а также считает количетво шагов сделанных в эту сторону и возвращает их
"""
function mark_and_enumerate_if_possible!(r::Robot, side::HorizonSide)
    num_steps = 0

    putmarker!(r)
    while move_if_possible!(r, side) == true
        putmarker!(r)
        num_steps += 1
    end
    return num_steps
end