text:
module WignerSeitz

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

using ..Brillouin:
    AVec,
    SHOWDIGITS,
    BasisEnum, CARTESIAN, LATTICE, setting, set_setting!,
    cartesianize, latticize
import ..Brillouin:
    latticize!, cartesianize!,
    basis
using StaticArrays
using LinearAlgebra: dot, cross, norm
using PyCall
using DocStringExtensions

import Base: getindex, size, IndexStyle, show, summary

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

export Cell, wignerseitz, faces, vertices

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

const PySpatial = PyNULL()
function __init__()
    # bringing in SciPy's Spatial module (for `Voronoi` and `ConvexHull`)
    copy!(PySpatial, pyimport_conda("scipy.spatial", "scipy"))
end

# ---------------------------------------------------------------------------------------- #
# STRUCTURES

"""
$(TYPEDEF)
$(TYPEDFIELDS)
"""
struct Cell{D} <: AVec{Vector{SVector{D, Float64}}} # over of polygons
    verts   :: Vector{SVector{D, Float64}}
    faces   :: Vector{Vector{Int}}
    basis   :: SVector{D, SVector{D, Float64}}
    setting :: Ref{BasisEnum} # internal field
end

faces(c::Cell)    = c.faces
vertices(c::Cell) = c.verts

# abstract array interface
Base.@propagate_inbounds getindex(c::Cell, i::Int) = vertices(c)[faces(c)[i]]
size(c::Cell) = size(c.faces)
IndexStyle(::Type{<:Cell}) = IndexLinear()

function summary(io::IO, c::Cell{D}) where D
    print(io, "Cell{", D, "} (",
              length(c), " faces, ",
              length(vertices(c)), " vertices)")
end
function show(io::IO, ::MIME"text/plain", c::Cell{D}) where D
    summary(io, c)
    println(io, ":")

    println(io, "  verts: ", round.(first(vertices(c)), digits=SHOWDIGITS))
    foreach(v -> println(io, "         ", round.(v, digits=SHOWDIGITS)), @view vertices(c)[2:end])

    println(io, "  faces: ", first(faces(c)))
    foreach(f -> println(io, "         ", f), @view faces(c)[2:end])

    print(io, "  basis: ")
    foreach(enumerate(basis(c))) do (i,v) # dance to ensure final line isn't `\n`-terminated
        i ≠ 1 && print(io, "         ")
        if i ≠ length(basis(c))
            println(io, round.(v, digits=SHOWDIGITS))
        else
            print(io, round.(v, digits=SHOWDIGITS))
        end
    end
end
# ---------------------------------------------------------------------------------------- #
# MAIN FUNCTION

"""
    wignerseitz(basis::AbstractVector{<:SVector{D}}; merge::Bool = true, Nmax = 3)
    wignerseitz(basis::AbstractVector{AbstractVector}; merge::Bool = true, Nmax = 3)
                                                                --> Cell{D}

Given a provided basis, return a `Cell{D}` containing the vertices and associated (outward
oriented) faces of the Wigner-Seitz cell defined by `basis` in `D` dimensions. The returned
vertices are given in the the basis of `basis` (see [`cartesianize!`](@ref) for conversion).

## Keyword arguments
- `merge` (default, `true`): if `true`, co-planar faces are merged to form polygonal
  planar faces (e.g., triangles, quadrilaterals, and ngons generally). If `false`, raw
  "unprocessed" triangles (`D=3`) and segments (`D=2`) are returned instead. `merge` has no
  impact for `D=1`.
- `Nmax` (default, `3`): includes `-Nmax:Nmax` points in the initial lattice used to
  generate the underlying Voronoi tesselation. It is unwise to set this to anything lower
  than 3 without explicitly testing convergence; and probably unnecessary to increase it
  beyond 3 as well.
"""
function wignerseitz(basis::AVec{<:SVector{D,<:Real}};
            merge::Bool = true,
            Nmax::Integer = 3) where D
    # "supercell" lattice of G-vectors
    Ns = -Nmax:Nmax
    lattice = Vector{SVector{D,Float64}}(undef, length(Ns)^D)
    idx_cntr = 0
    for (idx, I) in enumerate(CartesianIndices(ntuple(_->Ns, Val(D))))
                   V =  basis[1]*I[1]
        D >= 2 && (V += basis[2]*I[2])
        D == 3 && (V += basis[3]*I[3])
        lattice[idx] = V
        iszero(I) && (idx_cntr = idx)
    end

    vor = PySpatial.Voronoi(lattice) # voronoi tesselation of lattice

    # grab all the vertices of the central voronoi region enclosing origo
    verts_cntr =  # NB: offsets by 1 due to Julia 1-based vs. Python 0-based indexing
        [vor.vertices[idx+1,:] for idx in vor.regions[vor.point_region[idx_cntr]+1]]

    # get convex hull of central vertices
    hull = PySpatial.ConvexHull(verts_cntr)
    c    = convert_to_cell(hull, basis)
    c    = reorient_normals!(c)

    # return either raw simplices or "merged" polygons
    if merge
        return latticize!(merge_coplanar!(c))
    else
        return latticize!(c)
    end
end
# overload for input as ordinary vectors; type-unstable obviously, but convenient sometimes
function wignerseitz(basis::AVec{<:AVec{<:Real}}; kwargs...)
    D = length(first(basis))
    all(V->length(V) == D, @view basis[2:end]) || error(DomainError(basis, "provided `basis` must have identical dimensions"))

    return wignerseitz(SVector{D,Float64}.(basis); kwargs...)
end

function wignerseitz(basis::AVec{<:SVector{1,<:Real}}; kwargs...)
    a = only(only(basis))
    Cell{1}([(@SVector [-a/2]), (@SVector [a/2])], [1,2], basis)
end
# ---------------------------------------------------------------------------------------- #
# UTILITIES & SUBFUNCTIONS

function convert_to_cell(hull, basis::AVec{<:SVector{D,<:Real}}) where D
    vs′ = hull.points         # vertices
    fs′ = hull.simplices .+ 1 # faces

    vs = SVector{D, Float64}.(eachrow(vs′))
    fs = Vector{Int}.(eachrow(fs′))
    return Cell(vs, fs, SVector{D, SVector{D, Float64}}(basis), Ref(CARTESIAN))
end

"""
    is_outward_oriented(c, n)

Return whether a face, specified by its center `c` and normal `n`, is outward pointing,
assuming the face is part of a convex hull centered around origo.
"""
is_outward_oriented(c::SVector{D}, n::SVector{D}) where D = dot(c,n) > zero(eltype(c))

function reorient_normals!(c::Union{D}) where D
    D == 1 && return c # no sense of "outward" in 1D
    fs = faces(c)
    for (i,f) in enumerate(fs)
        cntr = face_center(c, i)
        n    = face_normal(c, i)
        if !is_outward_oriented(cntr, n)
            reverse!(fs[i])
        end
    end
    return c
end

function merge_coplanar!(c::Cell{3})
    fs = faces(c)
    ns = face_normals(c)

    # merge the all co-planar face; step forward in a "conquering" fashion
    i = 1
    while i ≤ length(fs)
        fᵢ = fs[i]
        maybe_more_to_merge = true
        while maybe_more_to_merge
            # this outer loop is here to ensure we do the j-loop until convergence (each
            # j-loop can create new merge-options if two faces were indeed merged)
            j  = i+1
            maybe_more_to_merge = false
            while j ≤ length(fs)
                fⱼ = fs[j]
                # can be co-planar faces if they share two or more vertices; we only
                # support merging with two vertices currently, so just check for that
                if sum(∈(fᵢ), fⱼ) == 2 && is_coplanar(ns[i], ns[j])
                    fs[i] = fᵢ = merge_coplanar(fᵢ, fⱼ)
                    deleteat!(fs, j), deleteat!(ns, j)
                    was_merged = true
                    maybe_more_to_merge = true
                else
                    was_merged = false
                end
                was_merged || (j += 1)
            end
        end
        i += 1
    end

    return c
end

# compute an (oriented) face-normal from a face `f` and a list of associated vertices `vs`
function face_normal(vs::AVec{<:SVector{3,<:Real}}, f::Vector{Int})
    v₁₂ = vs[f[2]] - vs[f[1]]
    v₂₃ = vs[f[3]] - vs[f[2]]
    n = cross(v₁₂, v₂₃)

    return n/norm(n)
end
function face_normal(vs::AVec{<:SVector{2,<:Real}}, f::Vector{Int})
    # define "canonical direction" as 1->2
    v₁₂ = vs[f[2]] - vs[f[1]]
    n = SVector{2, eltype(v₁₂)}(v₁₂[2], -v₁₂[1]) # = n = [y, -x]

    return n/norm(n)
end
face_normal(c::Cell, i::Int) = face_normal(vertices(c), faces(c)[i])
face_normals(c::Cell)        = face_normal.(Ref(vertices(c)), faces(c))

# compute a *rough* center position of a polygon
function face_center(vs::Vector{SVector{D,Float64}}, f::Vector{Int}) where D
    return sum(i -> vs[i], f)/D
end
face_center(c::Cell, i::Int) = face_center(vertices(c), faces(c)[i])

# assumes normal vectors to be normalized and identically oriented
is_coplanar(n1, n2, atol::Real=1e-8) = isapprox(dot(n1, n2), 1.0, atol=atol)

function merge_coplanar(fᵢ::Vector{Int}, fⱼ::Vector{Int})
    idxs_sharedᵢ = sort!(findall(∈(fⱼ), fᵢ))
    idxs_sharedⱼ = [findfirst(==(fᵢ[i]), fⱼ)::Int for i in idxs_sharedᵢ]

    fᵢ[idxs_sharedᵢ] == fⱼ[idxs_sharedⱼ] || error("failed face identification")
    length(idxs_sharedᵢ) != 2 && throw(DomainError(length(idxs_sharedᵢ), "unexpected: faces share fewer or more than two vertices"))

    # now we do a semi-terrifying dance to "follow-along" the edges of face i, then
    # transition into face j, and finally follow any remaining edges in face j; the trick
    # is to jump at from face i to j and vice versa when we hit a shared vertex; it's made
    # tedious by the fact that the face is not "cyclical" and we have to emulate that
    # ourselves
    Dᵢ, Dⱼ = length(fᵢ), length(fⱼ)
    f = Vector{Int}(undef, Dᵢ+Dⱼ-2)
    n = 1
    # start in ngon i
    idxsᵢ, neighbourⱼ = if idxs_sharedᵢ[1] == 1 && idxs_sharedᵢ[2] == Dᵢ
        2:Dᵢ, 2
    else
        1:idxs_sharedᵢ[1], 1
    end
    for idxᵢ in idxsᵢ
        f[n] = fᵢ[idxᵢ]
        n += 1
    end

    # go to ngon j
    if idxs_sharedⱼ[neighbourⱼ] == Dⱼ
        idxsⱼ = Iterators.flatten((1:Dⱼ-1,))
    else
        if idxs_sharedⱼ[neighbourⱼ] == 1
            idxsⱼ = Iterators.flatten((idxs_sharedⱼ[neighbourⱼ]+1:Dⱼ,))
        else
            other_neighborⱼ = neighbourⱼ == 2 ? 1 : 2
            idxsⱼ = Iterators.flatten((idxs_sharedⱼ[neighbourⱼ]+1:Dⱼ, 1:idxs_sharedⱼ[other_neighborⱼ]))
        end
    end
    for idxⱼ in idxsⱼ
        f[n] = fⱼ[idxⱼ]
        n += 1
    end

    # done with ngon j; go back to ngon i
    remaining_idxsᵢ = idxs_sharedᵢ[2]+1:Dᵢ
    for idxᵢ in remaining_idxsᵢ
        f[n] = fᵢ[idxᵢ]
        n += 1
    end

    return f
end

function merge_coplanar!(c::Cell{2})
    fs′ = segments2polygon(faces(c))
    empty!(faces(c))
    push!(faces(c), fs′)
    return c
end
function segments2polygon(fs::Vector{Vector{Int}})
    # assumes there is a *single* well-defined polygon defined by `fs`
    N   = size(fs,1) # number of segments
    fs′ = Vector{Int}(undef, N)
    fs′[1:2] .= fs[1] # starting point
    for i in 3:N
        j = findfirst(f -> f[1]==(fs′[i-1]), fs)
        fs′[i] = fs[j][2]
    end
    return fs′
end

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

function latticize!(c::Cell)
    setting(c) === LATTICE && return c
    basismatrix = hcat(c.basis...)
    c.verts .= latticize.(c.verts, Ref(basismatrix))
    set_setting!(c, LATTICE)
    return c
end
function cartesianize!(c::Cell)
    setting(c) === CARTESIAN && return c
    c.verts .= cartesianize.(c.verts, Ref(c.basis))
    set_setting!(c, CARTESIAN)
    return c
end

end # module


text:
using GMT, GMT.Drawing, GMT.Gdal
using Test
using Dates, Printf

@testset "GMT" begin

try
	run(`gmt --version`)	# Will fail if GMT is not installed.
	global got_it = true
catch
	@test 1 == 1			# Let tests pass for sake of not triggering a PkgEval failure
	global got_it = false
end

if (got_it)					# Otherwise go straight to end

	global dbg2 = 3			# Either 2 or 3. 3 to test the used kwargs
	global dbg0 = 0			# With 0 prints only the non-consumed options. Set to -1 to ignore this Vd

	GMT.GMT_Get_Version();
	GMT._precompile_()
	ma=[0];mi=[0];pa=[0];
	GMT.GMT_Get_Version(ma,mi,pa);
	API = GMT.GMT_Create_Session("GMT", 2, GMT.GMT_SESSION_NOEXIT + GMT.GMT_SESSION_EXTERNAL);
	GMT.GMT_Get_Ctrl(API);

	if (GMTver > v"6.1.1")
		println("		Entering: test_proj4.jl")
		include("test_proj4.jl")
		println("		Entering: test_gd_ext.jl")
		include("test_gd_ext.jl")
		println("		Entering: test_gdal.jl")
		include("test_gdal.jl")			# Fcks the automatic registering because building docs fails
		rm("point.csv")
		rm("lixo1.gmt")
		rm("lixo2.gmt")
	end
	println("		Entering: test_common_opts.jl")
	include("test_common_opts.jl")
	println("		Entering: test_B-GMTs.jl")
	include("test_B-GMTs.jl")
	println("		Entering: test_avatars.jl")
	include("test_avatars.jl")
	println("		Entering: test_GRDs.jl")
	include("test_GRDs.jl")
	println("		Entering: test_views.jl")
	include("test_views.jl")
	println("		Entering: test_PSs.jl")
	include("test_PSs.jl")
	println("		Entering: test_modern.jl")
	include("test_modern.jl")
	println("		Entering: test_P_a_T.jl")
	include("test_P_a_T.jl")
	println("		Entering: test_misc.jl")
	include("test_misc.jl")

	println("	GREENSPLINE")
	d = [0 6.44; 1820 8.61; 2542 5.24; 2889 5.73; 3460 3.81; 4586 4.05; 6020 2.95; 6841 2.57; 7232 3.37; 10903 3.84; 11098 2.86; 11922 1.22; 12530 1.09; 14065 2.36; 14937 2.24; 16244 2.05; 17632 2.23; 19002 0.42; 20860 0.87; 22471 1.26];
	greenspline(d, R="-2000/25000", I=100, S=:l, D=0, Vd=dbg2)

	println("	MAKECPT")
	makecpt(rand(10,1), E="", C=:rainbow, cptname="lixo.cpt");
	@test_throws ErrorException("E option requires that a data table is provided as well") makecpt(E="", C=:rainbow)
	cpt = makecpt(range="-1/1/0.1");
	if (GMTver > v"6.1.1")
		C = cpt4dcw("eu");
		C = cpt4dcw("PT,ES,FR", [3., 5, 8], range=[3,9,1]);
		C = cpt4dcw("PT,ES,FR", [.3, .5, .8], cmap=cpt);
	end
	@test_throws ErrorException("Unknown continent ue") cpt4dcw("ue")
	GMT.iso3to2_eu();
	GMT.iso3to2_af();
	GMT.iso3to2_na();
	GMT.iso3to2_world();
	GMT.mk_codes_values(["PRT", "ESP", "FRA"], [1.0, 2, 3], region="eu");
	@test_throws ErrorException("The region ue is invalid or has not been implemented yet.") GMT.mk_codes_values(["PRT"], [1.0], region="ue")

	println("	MAPPROJECT")
	mapproject([-10 40], J=:u29, C=true, F=true, V=:q);
	mapproject(region=(-15,35,30,48), proj=:merc, figsize=5, map_size=true);
	@test mapproject([1.0 1; 2 2], L=(line=[1.0 0; 4 3], unit=:c), Vd=dbg2) ==  "mapproject  -L+uc "
	@test mapproject([1.0 1; 2 2], L=[1.0 0; 4 3], Vd=dbg2) == "mapproject  -L "
	@test mapproject([1.0 1; 2 2], L="lixo.dat", Vd=dbg2) == "mapproject  -Llixo.dat "
	@test_throws ErrorException("Bad argument type (Tuple{}) to option L") mapproject([1.0 1; 2 2], L=(), Vd=dbg2)
	@test_throws ErrorException("line member cannot be missing") mapproject(mapproject([1.0 1; 2 2], L=(lina=[1.0 0; 4 3], unit=:c), Vd=dbg2))

	println("	PSMECA")
	meca([0.0 3.0 0.0 0 45 90 5 0 0], fill=:black, region=(-1,4,0,6), proj=:Merc, Vd=dbg2)
	meca!([0.0 3.0 0.0 0 45 90 5 0 0], fill=:black, region=(-1,4,0,6), proj=:Merc, Vd=dbg2)
	@test_throws ErrorException("Must select one convention") meca!("", [0.0 3 0 0 45 90 5 0 0 0], fill=:black, region=(-1,4,0,6), proj=:Merc)
	@test_throws ErrorException("Specifying cross-section type is mandatory") coupe([0.0 3 0 0 45 90 5 0 0], region=(-1,4,0,6))
	velo(mat2ds([0. -8 0 0 4 6 0.5; -8 5 3 3 0 0 0.5], ["4x6", "3x3"]), pen=(0.6,:red), fill_wedges=:green, outlines=true, Se="0.2/0.39/18", arrow="0.3c+p1p+e+gred", region=(-15,10,-10,10), Vd=dbg2)

	# Remove garbage
	println("	REMOVE GARBAGE")
	rm("gmt.history")
	rm("gmt.conf")
	rm("lixo.ps")
	rm("lixo.png")
	#rm("lixo.grd")
	rm("lixo.tif")
	rm("lixo.cpt")
	rm("lixo.dat")
	rm("logo.png")
	rm("lixo.eps")
	rm("lixo.jpg")
	#@static if (Sys.iswindows())  run(`rmdir /S /Q NUL`)  end

end					# End valid testing zone

end

ext: jl


text:
# auxiliary function to compute equal time correlator for su2 symmetry
function compute_eq_corr_su2(
    file    :: HDF5.File,
    Λ       :: Float64
    ;
    verbose :: Bool = true
    )       :: Float64

    # load diagonal correlations
    m, χ = read_χ(file, Λ, "diag", verbose = verbose)

    # compute equal-time on-site correlation
    eq_corr = 0.0

    for i in 1 : length(m) - 1
        eq_corr += (m[i + 1] - m[i]) * (χ[1, i + 1] + χ[1, i])
    end

    # calculate full correlator
    eq_corr *= 3.0

    # normalization from Fourier transformation
    eq_corr /= 2.0 * pi

    return eq_corr
end


text:
import RowEchelon

#=
Homework 6
Section 2.3
Problems 44
=#

# Create matrix to be inverted
A = [1 1/2 1/3 1/4 1/5 ;
     1/2 1/3 1/4 1/5 1/6 ;
     1/3 1/4 1/5 1/6 1/7 ;
     1/4 1/5 1/6 1/7 1/8 ;
     1/5 1/6 1/7 1/8 1/9]

invA = inv(A)

# print the last column
println("The last column of the inverse is:")
for i=1:5
    println(invA[i,5])
end
println("Additional decimal places are a result of floating point error.")


text:
### MCChain stores the output of a Monte Carlo iteration

type MCChain
  samples::Matrix{Float64}
  logtargets::Vector{Float64}
  gradlogtargets::Matrix{Float64}
  diagnostics::Dict
  runtime::Float64

  function MCChain(s::Matrix{Float64}, l::Vector{Float64}, g::Matrix{Float64}, d::Dict, r::Float64)
    @assert size(s, 1) == length(l) "Number of samples and logtargets values are not equal."
    if size(g) != (0, 0)
      @assert size(s) == size(g) "samples and gradlogtargets do not have the same number of rows and columns."
    end
    new(s, l, g, d, r)
  end
end

MCChain(s::Matrix{Float64}, l::Vector{Float64}, g::Matrix{Float64}, d::Dict) =
  MCChain(s, l, g, d, NaN)
MCChain(s::Matrix{Float64}, l::Vector{Float64}, g::Matrix{Float64}, r::Float64) =
  MCChain(s, l, g, Dict(), r)
MCChain(s::Matrix{Float64}, l::Vector{Float64}, d::Dict, r::Float64) = MCChain(s, l, Array(Float64, 0, 0), d, r)
MCChain(s::Matrix{Float64}, l::Vector{Float64}, g::Matrix{Float64}) = MCChain(s, l, g, Dict(), NaN)
MCChain(s::Matrix{Float64}, l::Vector{Float64}, d::Dict) = MCChain(s, l, Array(Float64, 0, 0), d, NaN)
MCChain(s::Matrix{Float64}, l::Vector{Float64}, r::Float64) = MCChain(s, l, Array(Float64, 0, 0), Dict(), r)
MCChain(s::Matrix{Float64}, l::Vector{Float64}) = MCChain(s, l, Array(Float64, 0, 0), Dict(), NaN)
MCChain() = MCChain(Array(Float64, 0, 0), Float64[], Array(Float64, 0, 0), Dict(), NaN)

function MCChain(npars::Int, nsamples::Int;
  storegradlogtarget::Bool=false, diagnostics::Dict=Dict(), runtime::Float64=NaN)
  if storegradlogtarget
    MCChain(fill(NaN, nsamples, npars), fill(NaN, nsamples), fill(NaN, nsamples, npars), diagnostics, runtime)
  else
    MCChain(fill(NaN, nsamples, npars), fill(NaN, nsamples), Array(Float64, 0, 0), diagnostics, runtime)
  end
end

function show(io::IO, c::MCChain)
  nsamples, npars = size(c.samples)
  println(io, "$npars parameters, $nsamples samples (per parameter), $(round(c.runtime, 1)) sec.")
end


text:
using PatModules
@mainmodule! example_use ("subpackage", "utils") begin
    export greetings

    function greetings()
        do_nothing()
        say_hi()
        say_bye()
    end
end


text:
using Test, TransitionalMCMC, StatsBase, Distributions, Logging, Distributed, Random, HypothesisTests

include("mcmc.jl")
include("tmcmc.jl")


text:
using Documenter, SimpleRoots

makedocs(
    sitename = "SimpleRoots.jl",
)

deploydocs(
    repo = "github.com/rafaqz/SimpleRoots.jl.git",
)


text:
# This file defines the functions pertaining to the Prerequisite type.

# It is important to note that prerequisites that define an inequality relation,
#   the relation is ALWAYS stated as
#     "prereqValue prereqRelation value_in_personnel_record".
# For example, if something requires a minimum age of 21.0, the prerequisite is
#   written as
#     "21.0 <= :age"

# The functions of the Prerequisite type require the Personnel,
#   PersonnelDatabase, HistoryEntry, and History types.
requiredTypes = [ "personnel", "personnelDatabase", "historyEntry", "history",
    "prerequisite" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

# This function tests if the prerequisite is on the given field in the personnel
#   database.  ( necessary??? )
export isPrereqVariable
function isPrereqVariable( prereq::Prerequisite, pVar::Symbol )
    return prereq.prereqVar === pVar
end  # isPrereqVariable( prereq, pVar )

# This function tests if a person satisfies the prerequisite. If this person
#   does not have the key being checked, the result defaults to false.
export isSatisfied
function isSatisfied( prereq::Prerequisite, person::Personnel )
    # If the variable is not available, the constraint is automatically
    #   violated.
    if !hasAttribute( person, prereq.prereqVar )
        return false
    end  # if hasAttribute( person, prereq.prereqVar )

    # If the value of the attribute in the personnel record is nothing, return
    #   true if the relation is !=, and false otherwise.
    if nothing === person[ prereq.prereqVar ]
        return prereq.prereqRelation === (!=)
    end

    # The try-catch block tries to test for the constraint. If testing the
    #   constraint generates an error, for example by testing Int ∈ Int, it
    #   will automatically be violated.
    try
        return prereq.prereqRelation( prereq.prereqValue,
            person[ prereq.prereqVar ] )
    catch
        return false
    end  # try
end  # isSatisfied( prereq, person )

# This function tests if a person satisfies the prerequisite at a specific time.
#   If this person has no history of the attribute being checked, it tests on
#   the actual attribute.
function isSatisfied( prereq::Prerequisite, person::Personnel,
    timestamp::T ) where T <: Real
    # If there's no history for the attribute, check the actual attribute.
    if !hasHistory( person, prereq.prereqVar )
        return isSatisfied( prereq, person )
    end  # if !hasHistory( person, prereq.prereqVar )

    # The try-catch block tries to test for the constraint. If testing the
    #   constraint generates an error, for example by testing Int ∈ Int, it
    #   will automatically be violated.
    try
        return prereq.prereqRelation( prereq.prereqValue,
            person[ prereq.prereqVar, timestamp ] )
    catch
        return false
    end  # try
end  # isSatisfied( prereq, person, timestamp )

# This function tests if the person with the given id in the database satisfies
#   the prerequisite.
function isSatisfied( prereq::Prerequisite, dbase::PersonnelDatabase,
    index::DbIndexType )
    return isSatisfied( prereq, dbase[ index ] )
end  # isSatisfied( prereq, dbase, index )

# This function tests if the person with the given id in the database satisfies
#   the prerequisite at the given time.
function isSatisfied( prereq::Prerequisite, dbase::PersonnelDatabase,
    index::DbIndexType, timestamp::T ) where T <: Real
    return isSatisfied( prereq, dbase[ index ], timestamp )
end  # isSatisfied( prereq, dbase, index, timestamp )

# This function selects all the records in the database satisfying the specific
#   prerequisite.
export selectRecords
function selectRecords( dbase::PersonnelDatabase, prereq::Prerequisite )
    return dbase[ find( person -> isSatisfied( prereq, person ), dbase.dbase ) ]
end  # selectRecords( dbase, prereq )

# This function counts all the records in the database satisfying the specific
#   prerequisite.
export countRecords
function countRecords( dbase::PersonnelDatabase, prereq::Prerequisite )
    return count( person -> isSatisfied( prereq, person ), dbase.dbase )
end  # countRecords( dbase, prereq )

# This function counts all the records in the database satisfying the specific
#   prerequisite at the given time.
function countRecords( dbase::PersonnelDatabase, prereq::Prerequisite,
    timestamp::T ) where T <: Real
    return count( person -> isSatisfied( prereq, person, timestamp ),
        dbase.dbase )
end  # countRecords( dbase, prereq, timestamp )

function Base.show( io::IO, prereq::Prerequisite )
    showFields = [ prereq.prereqValue, ' ', prereq.prereqRelation, ' ',
        prereq.prereqVar ]

    map( x -> print( io, x ), showFields )
end  # show( io, prereq )


text:
using Pkg, LibGit2;
user_regs = joinpath(DEPOT_PATH[1],"registries");
mkpath(user_regs);
all_registries = Dict("General" => "https://github.com/JuliaRegistries/General.git",
	      "HolyLabRegistry" => "https://github.com/HolyLab/HolyLabRegistry.git");
Base.shred!(LibGit2.CachedCredentials()) do creds
  for (reg, url) in all_registries
    path = joinpath(user_regs, reg);
    LibGit2.with(Pkg.GitTools.clone(url, path; header = "registry $reg from $(repr(url))", credentials = creds)) do repo end
  end
end

max_stars_repo_name: HolyLab/ImaginePlots

max_issues_repo_name: HolyLab/ImaginePlots


text:
push!(LOAD_PATH, "../src/")
using Documenter, Batman

makedocs(
  sitename="Batman.jl",
  modules=[Batman],
  pages = Any[
    "Home" => "index.md",
    "Library" => map(s -> "lib/$(s)", sort(readdir(joinpath(@__DIR__, "src/lib"))))
  ],
  format = Documenter.HTML( prettyurls = get(ENV, "CI", nothing) == "true" )
)

deploydocs(
  repo = "github.com/MorganAskins/Batman.jl.git",
)


text:
"""
    world_age()

Get the julia world age. Borrowed from Signals.jl.

"""
world_age() = ccall( :jl_get_world_counter, Int, () )

mutable struct Scheduler
    mission              ::MissionGraph
    nm                   ::NodeManager
    worker_communications::Dict{ Int, RemoteChannel{ Channel{ WorkerCommunication } } }
    task_stats           ::Dict{ Int, JobStatistics }
    worker_channels      ::Dict{ Int, RemoteChannel{ Channel{ Any } } }
end

"""
    Scheduler( nm::NodeManager, mission::MissionGraph )

Constructs a Scheduler object from a NodeManager and MissionGraph. From a high level
the Scheduler contains mappings from which nodes pertain to which machines, and
which tasks can be assigned to which workers. This constructor initializes the DAG, by
finding all parent nodes in a DAG, and allocating as many resources to them as possible.

"""
function Scheduler( nm::NodeManager, mission::MissionGraph )
    # Find each flow from each parent to it's respective terminal
    sinks = terminalnodes( mission.g )
    #map out initial worker to task relationships - naively
    worker_task_map = initial_task_assignments( nm, mission )
    #now lets make channels for these workers to talk to the local thread!
    worker_comm_map = Dict{ Int, RemoteChannel{ Channel{ WorkerCommunication } } }()
    worker_channels = Dict{ Int, RemoteChannel{ Channel{ Any } } }()    #everyworker needs their own channel to save state of tasks in.
    task_stats      = Dict{ Int, Any }( [ node => JobStatistics() for (node, other) in mission.meta])
    @sync for (worker, metadata) in nm.computemeta
        worker_comm_map[ worker ] = RemoteChannel( () -> Channel{WorkerCommunication}(1), worker )
        worker_channels[ worker ] = RemoteChannel( () -> Channel{Any}(1), worker )
        @spawnat worker put!( worker_comm_map[ worker ], WorkerCommunication(   JobStatisticsSample(),
                                                                                worker_task_map[ worker ],
                                                                                available ) )
    end
    @info("Global states assigned to workers.")
    return Scheduler( mission, nm, worker_comm_map, task_stats, worker_channels )
end

"""
    execute_mission( sc::Scheduler )

Begins the interpretation of a mission/DAG. This is the part where work
actually get's done.

"""
function execute_mission( sc::Scheduler )
    #ToDo: Add defensive programming to ensure that all of these channels exist.
    @sync for ( worker, comm ) in sc.worker_communications
        task = take!( sc.worker_communications[ worker ] )
        if task.last_task > 0
            try
                @async begin
                    #make a single buffer to get job statistics from a called and finished fn
                    if isa( sc.mission.meta[ task.last_task ], Stash) ||
                            isa( sc.mission.meta[ task.last_task ], StashIterator )
                        #if the current task is a stash, we need to handle iteration over collections and their state in the scheduler.
                        @spawnat worker begin
                            sc.mission.meta[ task.last_task ].fn
                            dispatch_task(  sc.mission.meta[ task.last_task ].fn,
                                            sc.worker_channels[ worker ],
                                            sc.worker_communications[ worker ],
                                            task.last_task,
                                            sc.mission.meta[ task.last_task ].src )
                        end
                    else
                        @spawnat worker begin
                            sc.mission.meta[ task.last_task ].fn
                            dispatch_task(  sc.mission.meta[ task.last_task ].fn,
                                            sc.worker_channels[ worker ],
                                            sc.worker_communications[ worker ],
                                            task.last_task )
                        end
                    end
                end
            catch
                @error("Failed to spawn initial tasks.")
                #failure to do @spawnat means something funamentally bad happened :/
                #@async sc.worker_state[ worker ] = WORKER_STATE.failed
            end
        end
    end
    @info("Done distributing initial tasks. Good luck")
end

"""
    spawn_listeners(sc:Scheduler)

Starts an event listener and task dispatcher loop.

"""
function spawn_listeners(sc::Scheduler)
    #Kick off the event listener loop... This is a bit ugly but whatever, Observables, and Signals
    #follow this pattern. I'd like a more actor style or true event listener style but this is okay for now!
    @info("Worker communications established...")
    eventloop_world_age = world_age()
    while true
        #look for tasks that have completed!
        completion_listener( sc )

        sleep(0.0010)
    end #end never ending while loop
end #end function...

"""
    completion_listener(sc::Scheduler)

Updates the local storage of a scheduler when tasks are completed, or tasks
reach the end of a DAG.

"""
function completion_listener(sc::Scheduler)
    keyset = [k for k in keys(sc.worker_communications)]
    for worker in keyset
        if all( isready.( [sc.worker_communications[ worker ], sc.worker_channels[ worker ] ] ) )
            bufferworker = take!( sc.worker_communications[ worker ] )
            if ( bufferworker.last_task > 0 ) && ( bufferworker.state == ready )
                nexttask = neighbors(sc.mission.g, bufferworker.last_task)
                if length(nexttask) == 0
                    @info("Worker #$worker completed round trip of a DAG.")
                    reset_worker(   sc.worker_channels[ worker ],
                                    sc.worker_communications[ worker ] )
                else
                    nexttask = first(nexttask)
                    OnlineStats.fit!(   sc.task_stats[ bufferworker.last_task ].elapsed_time,
                                        bufferworker.task_stats.elapsed_time )
                    OnlineStats.fit!(   sc.task_stats[ bufferworker.last_task ].bytes_allocated,
                                        bufferworker.task_stats.bytes_allocated )
                    #assign next task
                    @spawnat worker begin
                        sc.mission.meta[ nexttask ].fn
                        dispatch_task(  sc.mission.meta[ nexttask ].fn,
                                        sc.worker_channels[ worker ],
                                        sc.worker_communications[ worker ],
                                        nexttask )
                    end
                    bufferworker.last_task
                end
            elseif bufferworker.state == failed
                @error("Srs error") #ToDo: handle errors
            end
        end
    end #end for workers
end

"""
    initial_task_assignments(nm::NodeManager, mission::MissionGraph)

Naively assigns all available workers to all tasks immediately available.

"""
function initial_task_assignments(nm::NodeManager, mission::MissionGraph) #:< Dict{Int,Int}
    workersavailable = total_worker_counts( nm )
    worker_queue = Dict( [ worker => 0 for (worker,tmp) in nm.computemeta  ] )

    sources = parentnodes( mission.g )
    src_dem = (sum( [ mission.meta[src].min_workers for src in sources ] ))
    if src_dem < workersavailable
        @warn("More parent node workers requested ($src_dem) then workers available ($workersavailable). Some initial tasks will be unfullfilled at initialization.")
    end
    #Distribute jobs to workers by priority
    sources_by_priority = [ [ src, mission.meta[src].min_workers, mission.meta[src].priority ] for src in sources]
    source_demand_by_priority = sort( sources_by_priority, by = x -> last(x) )
    #Satisfy minimum number of workers
    #Loop over sources, demand their minimum requirements be met
    for src_idx in 1 : length( source_demand_by_priority )
        src, demand, priority = source_demand_by_priority[ src_idx ]
        #Loop over available workers
        for worker_idx in keys( worker_queue )
            if ( source_demand_by_priority[ src_idx ][ 2 ] > 0 )
                #ensure the machine the task is bound too matches the worker
                if (nm.computemeta[ worker_idx ].address in mission.meta[ src ].machines) &&
                        (worker_queue[worker_idx] == 0)
                    worker_queue[ worker_idx ] = src
                    source_demand_by_priority[ src_idx ][2] -= 1
                    #handle if this is a stash iterator.
                    if isa.( mission.meta[ src ].src, MissionNodeIterator ) &&
                            !isnothing( mission.meta[ src ].src ) &&
                            source_demand_by_priority[ src_idx ][2] == 0
                        next!( mission.meta[ src ] )
                        source_demand_by_priority[ src_idx ][2] = source_demand_by_priority[ src_idx ][1].min_worker_constant
                    end
                end
            else
                break; #go to source loop - this source has its demand met!
            end
        end #end worker loop
    end #end source loop
    #Untested
    source_demand = [ mission.meta[src].max_workers - mission.meta[src].min_workers for src in sources ]
    if workersavailable > sum( values( worker_queue ) .> 0 )
         source_demand_by_priority = sort( sources_by_priority, by = x -> last(x) )
         #Loop over sources, demand their minimum requirements be met
         for src_idx in 1 : length( source_demand_by_priority )
             src, demand, priority = source_demand_by_priority[ src_idx ]
             #Loop over available workers
             for worker_idx in keys( worker_queue )
                 if ( source_demand_by_priority[ src_idx ][2] > 0 )
                     #ensure the machine the task is bound too matches the worker
                     if (nm.computemeta[ worker_idx].address in mission.meta[ src ].machines) &&
                             (worker_queue[worker_idx] == 0)
                         worker_queue[ worker_idx ] = src
                         source_demand_by_priority[ src_idx ][2] -= 1
                     end
                 else
                     break; #go to source loop - this source has its demand met!
                 end
             end
         end
    end
    return worker_queue
end

max_stars_repo_name: caseykneale/MrPhelps

max_issues_repo_name: caseykneale/MrPhelps

max_forks_repo_name: caseykneale/MrPhelps


text:
#!/usr/bin/julia

using ArgParse

using JLD
using Multilane
using POMDPs
using MCTS
using DataFrames
using DataFramesMeta
using PGFPlots
pushPGFPlotsPreamble("\\usepackage{siunitx}")

function curve_plot(mp::AbstractDataFrame,
                   test_key::AbstractString;
                   label=replace(test_key,"_","\\_"),
                   style=nothing,
                   )
    data = sort!(@where(mp, :test_key.==test_key), cols=[:lambda])
    # Plots.Linear(data[:time_to_lane], data[:nb_brakes], legendentry=label, style=style)
    ttl = data[:time_to_lane]
    nb = data[:nb_brakes]
    ttl_sem = data[:ttl_sem]
    nb_sem = data[:nb_sem]
    mat = cat(2, ttl, nb, ttl_sem, nb_sem, ttl_sem, nb_sem)'
    Plots.ErrorBars(mat, legendentry=label, style=style)
end

function print_tex(p::TikzPictures.TikzPicture, name)
    println("""

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        % $name
        """)
    println(p.data)
    println("""

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        """)
end

function print_tex(p::Axis, name)
    print_tex(plot(p), name)
end

function show_svg(p)
    f = tempname()
    @show f
    PGFPlots.save(string(f, ".svg"), plot(p))
    # compile = `pdflatex -batch $f.tex`
    # println("running $compile")
    # run(compile)
    @spawn run(`display -alpha off $f.svg`)
    # PGFPlots.save(string(f, ".tex"), plot(p))
end

function ttl_at_brake(mp, test_key, nb_brakes)
    data = sort!(@where(mp, :test_key.==test_key), cols=[:lambda])
    i = findfirst(data[:nb_brakes] .<= nb_brakes)
    ttl = data[:time_to_lane]
    nb = data[:nb_brakes]
    return ttl[i-1] + (ttl[i]-ttl[i-1])*(nb_brakes-nb[i-1])/(nb[i]-nb[i-1])
end

function brake_at_ttl(mp, test_key, t)
    data = sort!(@where(mp, :test_key.==test_key), cols=[:lambda])
    ttl = data[:time_to_lane]
    i = findfirst(ttl .>= t)
    nb = data[:nb_brakes]
    return nb[i-1] + (nb[i]-nb[i-1])*(t-ttl[i-1])/(ttl[i]-ttl[i-1])
end

# filenames = ["combined_results_Sep_2_11_31.jld", "combined_results_Sep_2_11_52.jld", "combined_results_Sep_8_21_57.jld"]
filenames = ["combined_results_Sep_12_11_51.jld", "combined_results_Sep_16_16_46.jld", "combined_results_Sep_16_17_51.jld"]
filenames = [joinpath(Pkg.dir("Multilane"), "data", f) for f in filenames]

if !isdefined(:results)
    results = load(first(filenames))

    for f in filenames[2:end]
        new_results = load(f)
        try
            results = merge_results!(results, new_results)
        catch ex
            warn("Careful merge of data failed, there may be inconsistent data")
            results = merge_results!(results, new_results, careful=false)
        end
    end
end

stats = results["stats"]

mp = by(stats, [:solver_key, :lambda, :problem_key, :solver_problem_key, :test_key]) do df
    DataFrame(
                  steps_in_lane=mean(df[:steps_in_lane]),
                  nb_brakes=mean(df[:nb_brakes]),
                  time_to_lane=mean(df[:time_to_lane]),
                  steps=mean(df[:steps]),
                  nb_sem=std(df[:nb_brakes])/sqrt(length(df[:nb_brakes])),
                  ttl_sem=std(df[:time_to_lane])/sqrt(length(df[:time_to_lane]))
              )
end

@show unique(mp[:test_key])

N = 1000
rng = MersenneTwister(123)
axis_args = Dict{Symbol, Any}(
    # :width=>"0.5\\columnwidth"
    :width=>"4cm",
    :style=>"xshift=4cm, yshift=3cm, axis background/.style={fill=white}",
    # :ylabel=>"\$T\$ (\\si{s})",
    # :xlabel=>"\$v_0\$ (\\si{m/s})"
)
scatter_args = Dict{Symbol, Any}(
    :markSize=>1,
    # :style=>"fill opacity=1"
)

g = GroupPlot(1,3)

bg = standard_uniform(1.0, correlation=false)
bs = Array(BehaviorModel, N)
for i in 1:N
    bs[i] = rand(rng, bg)
end
v0s_0 = Float64[bs[i].p_idm.v0 for i in 1:N]
Ts_0 = Float64[bs[i].p_idm.T for i in 1:N]

a = Axis([
    curve_plot(mp, "upper_bound_unif", label=nothing, style="dashed"),
    curve_plot(mp, "assume_normal_unif", label=nothing),
    curve_plot(mp, "mlmpc_unif", label=nothing),
    curve_plot(mp, "pomcp_unif", label=nothing)
], xmin=0,
xmax=30,
ymin=0,
ymax=1.5,
# xlabel="Time to target lane",
ylabel="Braking actions per episode",
style="grid=both",
# title = "Scenario 1: \$\\rho = 0\$"
)

s = Axis(Plots.Scatter(v0s_0, Ts_0; scatter_args...); title="Scenario 1: \$\\rho=0\$", axis_args...)

push!(g,a)
# show_svg([a,s])
# print_tex(a, "UNIFORM SOLVERS")
PGFPlots.save("/home/zach/Devel/Behavior-aware_lane_changing/acc_17_paper/fig/uniform.tex", [a,s], include_preamble=false)

@show ub = ttl_at_brake(mp, "upper_bound_unif", 0.5)
@show an = ttl_at_brake(mp, "assume_normal_unif", 0.5)
@show an - ub
@show (an - ub)/an

@show bub = brake_at_ttl(mp, "upper_bound_unif", 10.0)
@show ban = brake_at_ttl(mp, "assume_normal_unif", 10.0)
@show ban - bub
@show (ban - bub)/ban

bg = standard_uniform(1.0, correlation=true)
bs = Array(BehaviorModel, N)
for i in 1:N
    bs[i] = rand(rng, bg)
end
v0s_1 = Float64[bs[i].p_idm.v0 for i in 1:N]
Ts_1 = Float64[bs[i].p_idm.T for i in 1:N]

a = Axis([
    curve_plot(mp, "upper_bound", label=nothing, style="dashed"),
    curve_plot(mp, "assume_normal", label=nothing),
    curve_plot(mp, "mlmpc", label=nothing),
    curve_plot(mp, "pomcp", label=nothing)
], xmin=0,
xmax=30,
ymin=0,
ymax=1.5,
# xlabel="Time to target lane",
ylabel="Braking actions per episode",
style="grid=both",
# title="Scenario 2: \$\\rho = 1\$"
)

s = Axis(Plots.Scatter(v0s_1, Ts_1; scatter_args...); title="Scenario 2: \$\\rho=1\$", axis_args...)

push!(g,a)
# show_svg([a,s])
# print_tex(a, "UNIFORM SOLVERS")
PGFPlots.save("/home/zach/Devel/Behavior-aware_lane_changing/acc_17_paper/fig/correlated.tex", [a,s], include_preamble=false)

bg = standard_uniform(1.0, correlation=0.75)
bs = Array(BehaviorModel, N)
for i in 1:N
    bs[i] = rand(rng, bg)
end
v0s_75 = Float64[bs[i].p_idm.v0 for i in 1:N]
Ts_75 = Float64[bs[i].p_idm.T for i in 1:N]

a = Axis([
    curve_plot(mp, "upper_bound_075", label="Omniscient", style="dashed"),
    curve_plot(mp, "assume_normal_075", label="SAB"),
    curve_plot(mp, "mlmpc_075", label="MLMPC"),
    curve_plot(mp, "pomcp_075", label="POMCP")
], xmin=0,
xmax=30,
ymin=0,
ymax=1.5,
xlabel="Time to target lane",
ylabel="Braking actions per episode",
style="grid=both, legend columns=2, legend style={at={(0.5, -0.2)},anchor=north}",
# title="Scenario 3: \$\\rho = 0.75\$"
)

s = Axis(Plots.Scatter(v0s_75, Ts_75; scatter_args...); title="Scenario 3: \$\\rho=0.75\$", axis_args...)

push!(g,a)
# show_svg([a,s])
# print_tex(a, "UNIFORM SOLVERS")
PGFPlots.save("/home/zach/Devel/Behavior-aware_lane_changing/acc_17_paper/fig/correlated_075.tex", [a,s], include_preamble=false)

# show_svg(g)


text:
mutable struct Branch
    F_BUS::Vector{Int64}    # "f, from bus number"
    T_BUS::Vector{Int64}    # "t, to bus number"
    BR_R::Vector{Float64}   # "r, resistance (p.u.)"
    BR_X::Vector{Float64}   # "x, reactance (p.u.)"
    BR_B::Vector{Float64}   # "b, total line charging susceptance (p.u.)"
    RATE_A::Vector{Float64} # "rateA, MVA rating A (long term rating)"
    RATE_B::Vector{Float64} # "rateB, MVA rating B (short term rating)"
    RATE_C::Vector{Float64} # "rateC, MVA rating C (emergency rating)"
    TAP::Vector{Float64}    # "ratio, transformer off nominal turns ratio"
    SHIFT::Vector{Float64}  # "angle, transformer phase shift angle (degrees)"
    BR_STATUS::Vector{Bool} # "initial branch status, 1 - in service, 0 - out of service"
    ANGMIN::Vector{Float64} # "minimum angle difference, angle(Vf) - angle(Vt) (degrees)"
    ANGMAX::Vector{Float64} # "maximum angle difference, angle(Vf) - angle(Vt) (degrees)"
    PF::Vector{Float64}     # "real power injected at 'from' bus end (MW)       (not in PTI format)"
    QF::Vector{Float64}     # "reactive power injected at 'from' bus end (MVAr) (not in PTI format)"
    PT::Vector{Float64}     # "real power injected at 'to' bus end (MW)         (not in PTI format)"
    QT::Vector{Float64}     # "reactive power injected at 'to' bus end (MVAr)   (not in PTI format)"
    MU_SF::Vector{Float64}  # "Kuhn-Tucker multiplier on MVA limit at 'from' bus (u/MVA)"
    MU_ST::Vector{Float64}  # "Kuhn-Tucker multiplier on MVA limit at 'to' bus (u/MVA)"
    MU_ANGMIN::Vector{Float64}  # "Kuhn-Tucker multiplier lower angle difference limit (u/degree)"
    MU_ANGMAX::Vector{Float64}  # "Kuhn-Tucker multiplier upper angle difference limit (u/degree)"
end

"""
Extract branch information from mpc Dict. Return instance of Branch.

IN:

* mpc. Instance of MPC dict (loaded, for instance, by loadcase())

OUT:

* Instance of Branch containing all branch data.
"""
function extract_branch(mpc)
    branch = mpc["branch"]
    nc = size(branch, 2)
    nf = length(fieldnames(Branch))
    blanks = fill([NaN], nf - nc)
    return Branch([branch[:, i] for i in 1:nc]..., blanks...)
end

"""
Shortcut for returning branch info for one of the cases listed by
`casenames()`.
"""
function extract_branch(cname::AbstractString)
    mpc = loadcase(cname; describe=false)
    extract_branch(mpc)
end


text:
using CorrelationSpectroscopy
using Test

using Gtk  # https://github.com/JuliaPlots/VisualRegressionTests.jl/issues/30
using LinearAlgebra
using Plots
using Statistics
using VisualRegressionTests

"""
    demodata()

Calculate a demo dataset.

The demo dataset consists of fifteen signals at specific intervals, with peaks
at 10 and 20 rising in intensity, whereas peaks at 30 and 40 are decreasing in
intensity.
The example was inpired by this [one](https://en.wikipedia.org/wiki/Two-dimensional_correlation_analysis#Properties_of_2D_correlation_analysis>).

This function might go into the main code in the future if there is demand.

# Examples
```jldoctest
julia> x, A = demodata();

julia> length(x)
500

julia> size(A)
(15, 500)
```
"""
function demodata()
    n = 15
    x = range(0, 50, length=500)

    a₀ = 0.175
    A₁(n) = a₀ + 0.0125exp(-(n - 7)^2) + ((0.2375 - a₀) - exp(-0.8(n + 2.5)))
    A₂(n) = a₀ - 0.037exp(-0.2(n - 7)^2) + 1.025(0.2375 - a₀) / (1 + exp(-(n - 7)))
    A₃(n) = a₀*exp(-0.17(n -1)) + 0.0045(n - 1)
    A₄(n) = a₀ + 0.03exp(-0.3(n - 6)^2) - 1.6(0.2375 - a₀) / (1 + exp(-(n - 7)))

    σ = 2.0
    A = [
        A₁(t) * gaussian(xᵢ, 10, σ) +
        A₂(t) * gaussian(xᵢ, 20, σ) +
        A₃(t) * gaussian(xᵢ, 30, σ) +
        A₄(t) * gaussian(xᵢ, 40, σ) +
        0.02 * gaussian(xᵢ, 25, 10σ) +
        0.001 * randn()
        for t in 1:n, xᵢ in x
    ]

    return x, A
end
gaussian(x, μ=0.0, σ=1.0) = exp(-((x - μ) / σ)^2)

@testset "CorrelationSpectroscopy.jl" begin
    include("test_nodahilbert.jl")
    include("test_demodata.jl")
    include("test_synchronous.jl")
    include("test_asynchronous.jl")
    include("test_disrelation.jl")
end


text:
srand(123)

x = randn(10)
X = randn(10,10)+im*randn(10,10)
y0 = randn(size(x))
Y0 = randn(size(X))+im*randn(size(X))
y, Y = copy(y0), copy(Y0)

lambdas = [abs.(randn(size(x))), 0.1]
prox_col = [NormL1(lambdas[1]), NormL2(lambdas[2])]

# testing constructors
f = SeparableSum(prox_col)
y, fy = prox(f, [x, X], 1.)

y1, fy1 = prox(prox_col[1], x, 1.)
y2, fy2 = prox(prox_col[2], X, 1.)

@test norm((fy1+fy2)-fy) <= 1e-11
@test norm(y[1]-y1) <= 1e-11
@test norm(y[2]-y2) <= 1e-11

println(f)
# println(SeparableSum(repmat(prox_col,3)))


text:
# Measure execution time of a complexity, degeneracy, or redundancy function
# To run:  @time exec_time(10, complexity5, 4,4,16)  where complexity5 can be replaced by degeneracy or redundancy

maxints_for_degen = 20   #  Default
function exec_time( reps::Int64, cmplx_fcn::Function, numinputs::Int64, numoutputs::Int64, numinteriors::Int64 )
  p = Parameters(numinputs=numinputs, numoutputs=numinputs, numinteriors=numinteriors, numlevelsback=numinteriors)
  funcs=default_funcs(numinputs)
  for _  = 1:reps
    c = random_chromosome( p, funcs )
    cc = cmplx_fcn( c )
  end
end

#  julia> @time exec_time(10, complexity5, 4,4,16)
#    1.979641 seconds (31.73 M allocations: 1.084 GiB, 5.20% gc time)
#
#  julia> @time exec_time(10, degeneracy, 4,4,16)
#    4.620361 seconds (86.13 M allocations: 2.730 GiB, 4.38% gc time)
# Conclusion:  complexity5 is much faster than complexity4, complexity7
# Conclusion:  complexity5 is somewhat faster than complexity6
# Conclusion:  degeneracy is much faster than degeneracy1


text:
## Dimensions and descriptions

Base.isempty(x::RLEVector) = isempty(x.runends)
Base.size(x::RLEVector) =  (length(x),)

function Base.size(x::RLEVector, dim::Integer)
  len = length(x)
  if dim == 1
    return(len)
  else
    return( (len,1) )
  end
end

nrun(x::RLEVector) = length(x.runends)

function Base.length(x::RLEVector{T1,T2}) where {T1,T2<:Integer}
  re = x.runends
  ind = lastindex(re)
  if (ind > 0)
    @inbounds len = re[ind]
  else
    len = zero(T2)
  end
  len
end

function starts(x::RLEVector)
  re = x.runends
  rval = similar(re)
  prev = zero(eltype(re))
  @inbounds for i in eachindex(re)
        rval[i] = prev + 1
        prev = re[i]
  end
  rval
end

function starts(x::RLEVector, run::Integer)
    num_one = one(eltype(x.runends))
    run == 1 ? num_one : x.runends[run-1] + num_one
end

function widths!(x::AbstractVector)
    @inbounds for i in length(x):-1:2
        x[i] = x[i] - x[i-1]
    end
    x
end
widths(x::RLEVector) = widths!(copy(x.runends))

function widths(x::RLEVector, run::Integer)
  run == 1 ? x.runends[1] : x.runends[run] - x.runends[run-1]
end

ends(x::RLEVector) =  copy(x.runends)
_ends(x::RLEVector) =  x.runends
ends(x::RLEVector, run::Integer) =  x.runends[run]

values(x::RLEVector) = copy(x.runvalues)
_values(x::RLEVector) = x.runvalues
endtype(x::RLEVector) = eltype(rlast(x))

rfirst(x::RLEVector) = starts(x)
rfirst(x::RLEVector, run::Integer) = starts(x, run)
rwidth(x::RLEVector) = widths(x)
rwidth(x::RLEVector, run::Integer) = widths(x, run)
rlast(x::RLEVector) =  x.runends
rvalue(x::RLEVector) = x.runvalues

@deprecate rwidth widths
@deprecate rstart starts
@deprecate rlast ends
@deprecate rvalue values

@doc (@doc RLEVector) starts, widths, ends, values, rfirst, rwidth, rlast, rvalue, nrun, endtype


text:
using Onda
using Documenter

makedocs(modules=[Onda],
         sitename="Onda",
         authors="Beacon Biosignals, Inc.",
         pages=["API Documentation" => "index.md"])

# this is commented out until Onda has been open-sourced
deploydocs(repo="github.com/beacon-biosignals/Onda.jl.git")


text:
RedfieldOperator(H, R) = OpenSysOp(H, R, size(H, 1))

"""
$(TYPEDEF)

Defines RedfieldLiouvillian.

# Fields

$(FIELDS)
"""
struct RedfieldLiouvillian <: AbstractLiouvillian
    "Redfield kernels"
    kernels::Any
    "close system unitary"
    unitary::Any
    "absolute error tolerance for integration"
    atol::Float64
    "relative error tolerance for integration"
    rtol::Float64
    "cache matrix for inplace unitary"
    Ut::Union{Matrix,MMatrix}
    "cache matrix for inplace unitary"
    Uτ::Union{Matrix,MMatrix}
    "cache matrix for integration"
    Λ::Union{Matrix,MMatrix}
    "tf minus coarse grain time scale"
    Ta::Real
end

function RedfieldLiouvillian(kernels, U, Ta, atol, rtol)
    m_size = size(kernels[1][2])
    Λ = m_size[1] <= 10 ? zeros(MMatrix{m_size[1],m_size[2],ComplexF64}) :
        zeros(ComplexF64, m_size[1], m_size[2])
    # if the unitary does not in place operation, assign a pesudo inplace
    # function
    unitary = isinplace(U) ? U.func : (cache, t) -> cache .= U(t)
    RedfieldLiouvillian(kernels, unitary, atol, rtol, similar(Λ),
        similar(Λ), Λ, Ta)
end

function (R::RedfieldLiouvillian)(du, u, p, t::Real)
    s = p(t)
    for (inds, coupling, cfun) in R.kernels
        for (i, j) in inds
            function integrand(cache, x)
                R.unitary(R.Ut, t)
                R.unitary(R.Uτ, x)
                R.Ut .= R.Ut * R.Uτ'
                mul!(R.Uτ, coupling[j](p(x)), R.Ut')
                # The 5 arguments mul! will to produce NaN when it
                # should not. May switch back to it when this is fixed.
                # mul!(cache, R.Ut, R.Uτ, cfun[i, j](t, x), 0)
                mul!(cache, R.Ut, R.Uτ)
                lmul!(cfun[i, j](t, x), cache)
            end
            quadgk!(
                integrand,
                R.Λ,
                max(0.0, t - R.Ta),
                t,
                rtol=R.rtol,
                atol=R.atol,
            )
            SS = coupling[i](s)
            𝐊₂ = SS * R.Λ * u - R.Λ * u * SS
            𝐊₂ = 𝐊₂ + 𝐊₂'
            axpy!(-1.0, 𝐊₂, du)
        end
    end
end

function update_vectorized_cache!(cache, R::RedfieldLiouvillian, p, t::Real)
    iden = one(R.Λ)
    s = p(t)
    for (inds, coupling, cfun) in R.kernels
        for (i, j) in inds
            function integrand(cache, x)
                R.unitary(R.Ut, t)
                R.unitary(R.Uτ, x)
                R.Ut .= R.Ut * R.Uτ'
                mul!(R.Uτ, coupling[j](p(x)), R.Ut')
                # The 5 arguments mul! will to produce NaN when it
                # should not. May switch back to it when this is fixed.
                # mul!(cache, R.Ut, R.Uτ, cfun[i, j](t, x), 0)
                mul!(cache, R.Ut, R.Uτ)
                lmul!(cfun[i, j](t, x), cache)
            end
            quadgk!(
                integrand,
                R.Λ,
                max(0.0, t - R.Ta),
                t,
                rtol=R.rtol,
                atol=R.atol,
            )
            SS = coupling[i](s)
            SΛ = SS * R.Λ
            cache .-= (
                iden ⊗ SΛ + conj(SΛ) ⊗ iden - transpose(SS) ⊗ R.Λ -
                conj(R.Λ) ⊗ SS
            )
        end
    end
end


text:
using PyCall: pyimport
using Images: warp, colorview, channelview, RGB, imresize
using CoordinateTransformations: LinearMap, RotMatrix

np = pyimport("numpy")
cv2 = pyimport("cv2")

function arr2rgb(img)
    img = Float32.(img ./ 255.0)
    #BGR2RGB
    img = img[:, :, reverse(1:end)]
    img = permutedims(img, (3,1,2))
    img = colorview(RGB, img)
    return img
end

function rgb2arr(img)
    img = channelview(img)
    img = permutedims(img, (2,3,1))
    #RGB2BGR
    img = img[:, :, reverse(1:end)]
    img = np[:asarray](255. .* img, dtype=np[:uint8])
    return img
end

function rotate(img, degree)
    rad = degree * pi / 180
    tfm = LinearMap(RotMatrix(rad))
    return warp(img, tfm)
end

function process_something(img)
    rot_img = rotate(img, 90)
    rot_img = imresize(rot_img, (224,224))
    return rot_img
end

function video()
    cap = cv2[:VideoCapture](0)
    if !cap[:isOpened]()
        print("Error opening video stream or file")
        exit(1)
    end

    while cap[:isOpened]()
        ret_val, img = cap[:read]()
        img = arr2rgb(img)
        img = imresize(img,(224,224))
        # do something on Julia
        new_img = process_something(img)
        # concat original and processed images
        concatenated = cat([img, new_img]..., dims=2)
        out = rgb2arr(concatenated)
        cv2[:imshow]("julia app press enter esc", out)
        if cv2[:waitKey](1) == 27
            #press `esc` key to exit
            break
        end
    end
end

video()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using Gadfly

x = rand(10);y=rand(10);
plot(x=x,y=y,Guide.manual_color_key("Title", ["One", "Two"], ["green","blue"]))


text:
module XML2JSON

using LightXML
using DataStructures

export xml2json, parse_file, root

include("xml_to_json.jl")
include("xml2dict.jl")
include("utils.jl")

end # module


text:
export hamm_op

function hamm_op(d_obs::judiVector{vDT}) where {vDT}
    d_out = deepcopy(d_obs)
    nrec = length(d_obs.geometry.xloc[1])
    hamm_weights = convert(Array{vDT},reshape(hamming(nrec), 1, :))
    for i = 1:d_obs.nsrc
        d_out.data[i] = d_obs.data[i] .* (hamm_weights)
    end
    return d_out
end


text:
# This file was generated, do not modify it. # hide
sqrt.([4, 9, 16])


text:
using .Networks

using Dates
using JuMP
using Printf
using Parameters

@with_kw struct TSOOutFO <: AbstractTSO
    configs::TSOConfigs = TSOConfigs()
end

TSOInFO = TSOOutFO

function run(runnable::TSOOutFO, ech::Dates.DateTime, firmness, TS::Vector{Dates.DateTime}, context::AbstractContext)
    problem_name_l = @sprintf("tso_out_fo_%s", ech)

    tso_starts = get_starts(get_tso_schedule(context), get_generators_initial_state(context))
    market_starts = get_starts(get_market_schedule(context), get_generators_initial_state(context))
    gratis_starts = union(tso_starts, market_starts)

    runnable.configs.out_path = context.out_dir
    runnable.configs.problem_name = problem_name_l

    return tso_out_fo(get_network(context),
                    TS,
                    get_generators_initial_state(context),
                    get_scenarios(context, ech),
                    get_uncertainties(context, ech),
                    firmness,
                    get_market_schedule(context),
                    get_tso_schedule(context),
                    gratis_starts,
                    runnable.configs
                    )
end

function update_tso_schedule!(context::AbstractContext, ech, result, firmness,
                            runnable::TSOOutFO)
    tso_schedule = get_tso_schedule(context)
    tso_schedule.decider_type = DeciderType(runnable)
    tso_schedule.decision_time = ech

    for ((gen_id, ts, s), p_injected_var) in result.limitable_model.p_injected
        set_prod_value!(tso_schedule, gen_id, ts, s, value(p_injected_var))
    end
    for ((gen_id, ts, s), p_injected_var) in result.pilotable_model.p_injected
        if get_power_level_firmness(firmness, gen_id, ts) == FREE
            set_prod_value!(tso_schedule, gen_id, ts, s, value(p_injected_var))
        elseif get_power_level_firmness(firmness, gen_id, ts) == TO_DECIDE
            set_prod_definitive_value!(tso_schedule, gen_id, ts, value(p_injected_var))
        elseif get_power_level_firmness(firmness, gen_id, ts) == DECIDED
            # NOTE : decided reference may not be the TSO schedule
            @assert( ismissing(get_prod_value(tso_schedule, gen_id, ts))
                    || value(p_injected_var) ≈ get_prod_value(tso_schedule, gen_id, ts) )
            set_prod_definitive_value!(tso_schedule, gen_id, ts, value(p_injected_var))
        end
    end

    for ((gen_id, ts, s), b_on_var) in result.pilotable_model.b_on
        gen_state_value = parse(GeneratorState, value(b_on_var))
        if get_commitment_firmness(firmness, gen_id, ts) == FREE
            set_commitment_value!(tso_schedule, gen_id, ts, s, gen_state_value)
        elseif get_commitment_firmness(firmness, gen_id, ts) in [TO_DECIDE, DECIDED]
            set_commitment_definitive_value!(tso_schedule, gen_id, ts, gen_state_value)
        end
    end

    # Capping
    update_schedule_capping!(tso_schedule, context, ech, result.limitable_model)

    # loss_of_load (load-shedding)
    update_schedule_loss_of_load!(tso_schedule, context, ech, result.lol_model)

    return tso_schedule
end

function update_schedule_capping!(tso_schedule, context, ech,
                                    limitable_model::TSOLimitableModel)
    reset_capping!(tso_schedule)
    for ((gen_id,ts, s), p_injected_var) in limitable_model.p_injected
        available_prod = get_uncertainties(get_uncertainties(context, ech), gen_id, ts, s)
        injected_prod = value(p_injected_var)
        tso_schedule.capping[gen_id, ts, s] = available_prod - injected_prod
    end
end

function update_schedule_loss_of_load!(tso_schedule, context, ech, lol_model::TSOLoLModel)
    reset_loss_of_load_by_bus!(tso_schedule)

    for ((bus_id, ts, s), p_loss_of_load_var) in lol_model.p_loss_of_load
        tso_schedule.loss_of_load_by_bus[bus_id, ts, s] = value(p_loss_of_load_var)
    end
end

function update_tso_actions!(context::AbstractContext, ech, result, firmness,
                            ::TSOOutFO)
    tso_actions = get_tso_actions(context)

    # Limitations :
    # FIXME ; limit only if there is a limitation needed !?
    limitations = SortedDict{Tuple{String,DateTime}, Float64}() #TODELETE
    for ((gen_id, ts, s), p_limit_var) in result.limitable_model.p_limit
        if get_power_level_firmness(firmness, gen_id, ts) in [TO_DECIDE, DECIDED]
            @assert( value(p_limit_var) ≈ get!(limitations, (gen_id, ts), value(p_limit_var)) ) #TODELETE : checks that all values are the same across scenarios
            set_limitation_value!(tso_actions, gen_id, ts, value(p_limit_var))
        end
    end

    # Impositions
    impositions = SortedDict{Tuple{String,DateTime}, Float64}() #TODELETE
    for ((gen_id, ts, s), p_injected_var) in result.pilotable_model.p_injected
        if get_power_level_firmness(firmness, gen_id, ts) in [TO_DECIDE, DECIDED]
            @assert( value(p_injected_var) ≈ get!(impositions, (gen_id, ts), value(p_injected_var)) ) #TODELETE : checks that all values are the same across scenarios
            set_imposition_value!(tso_actions, gen_id, ts, s, value(p_injected_var), value(p_injected_var))
        end
    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.1.3
#   kernelspec:
#     display_name: Julia 1.1.1
#     language: julia
#     name: julia-1.1
# ---

# https://gist.github.com/agramfort/850437

using Statistics, LinearAlgebra

"""
    yest = lowess(x, y, f=2./3., iter=3)

    Lowess smoother: Robust locally weighted regression.
    The lowess function fits a nonparametric regression curve to a scatterplot.
    The arrays x and y contain an equal number of elements; each pair
    (x[i], y[i]) defines a data point in the scatterplot. The function returns
    the estimated (smooth) values of y.
    The smoothing span is given by f. A larger value for f will result in a
    smoother curve. The number of robustifying iterations is given by iter. The
    function will run faster with a smaller number of iterations.
"""
function lowess(x, y; f = 2/3, iter=3)

    n = length(x)
    r = Int(ceil(f * n))
    h = [sort(abs.(x .- x[i]))[r] for i in eachindex(x)]
    w = clamp.(abs.((x .- x') ./ h), 0.0, 1.0)
    w .= (1 .- w.^3) .^ 3
    yest = zeros(Float64,n)
    residuals = similar(yest)
    delta = ones(Float64,n)
    weights = similar(delta)
    for iteration in 1:iter
        for i in 1:n
            weights .= delta .* view(w, :, i)
            b = [sum(weights .* y); sum(weights .* y .* x)]
            A = SymTridiagonal([sum(weights), sum(weights .* x .* x)],
                               [sum(weights .* x)])
            beta = A \ b
            yest[i] = beta[1] + beta[2] * x[i]
        end

        residuals .= y .- yest
        s = median(abs.(residuals))
        delta .= clamp.(residuals / (6.0 * s), -1, 1)
        delta .= (1 .- delta.^2).^2
    end

    yest
end

n = 1000
x = range(0, stop=2π, length=n) |> collect
y = sin.(x) .+ 0.3 .* randn(n)
f = 0.25

using Plots, BenchmarkTools
yest = lowess(x, y, f=f, iter=1);

scatter(x, y; markersize=2, label=:noisy)
plot!(x, yest; line=(3,:red),label=:pred)


text:
using Documenter
using Latexify
using CSV

function makefigs(ext="svg")
    try
        imgdir = "doc/src/img"
        dotfiles = filter(x->endswith(x, ".dot"), readdir(imgdir))
        for dotfile in dotfiles
            run(`dot -T$ext -O doc/src/img/$dotfile`)
        end

    catch ex
        @warn "Could not update figures, perhaps dot is not installed."
        @warn ex
    end

end

function printmdtable(dir, outdir=".")
    for path in readdir(dir)
        df = CSV.read("$dir/$path")
        s = mdtable(df,latex=false)
        open("$outdir/$path.md", "w") do fp
            print(fp, string(s))
        end
    end
end

@info "Making Figures"
makefigs()
makefigs("png")

# printmdtable("examples/knowledge_graph/data", "doc/src/schema/")

@info "Loading module"
using SemanticModels
@info "Making docs"
makedocs(
modules     = [SemanticModels],
root        = "doc",
format      = Documenter.HTML(),
sitename    = "SemanticModels",
doctest     = false,
pages       = Any[
    "SemanticModels.jl"               => "index.md",
    "Intended Use Cases" => "usecases.md",
    "News" => "news.md",
    "ModelTools" =>"modeltools.md",
    "Example" => "example.md",
    "Theory" =>"theory.md",
    "Dubstep" => "dubstep.md",
    "Knowledge Graphs" => "graph.md",
    "Knowledge Extraction" => "extraction.md",
    "Validation" => "validation.md",
    "Library Reference" => "library.md",
    "Approaches" => "approach.md",
    "Slides"               => "slides.md",
    "Flu Model" => "FluModel.md",
    # "Model Types"                   => "types.md",
    # # "Reading / Writing Models"    => "persistence.md",
    # # "Plotting"                    => "plotting.md",
    # # "Parallel Algorithms"         => "parallel.md",
    "Contributing"                  => "contributing.md",
    # "License Information"           => "license.md",
    # "Citing SemanticModels"         => "citing.md"
]
)

deploydocs(
target      = "build",
deps        = nothing,
make        = nothing,
repo        = "github.com/jpfairbanks/SemanticModels.jl.git",
# julia       = "stable",
# osname      = "linux"
)


text:
function eller(n)
    result = 1
    answer,hw  = factor(n)
    for i in 1:length(answer)
        if hw[i] == 1
            result *= answer[i]-1
        else
            result *= (answer[i]^hw[i] - answer[i]^(hw[i]-1))
        end
    end
    return result
end

function factor(n)
    answer = []
    hw = []
    for i in eratosphen(n)
        if n%i==0
            push!(answer,i)
            push!(hw,1)
            n=n/i
            while n%i==0
                hw[end]+=1
                n=n/i
            end
        end
    end
    return answer,hw
end

function eratosphen(n)
    ser = fill(true,n)
    ser[1] = false
    k = 2
    while k !== nothing && k<n
        ser[k^2:k:end] .= false
        k=findnext(ser, k+1)
    end
    return findall(ser)
end


text:
# using Base.Test
# using MIToS.PDB

txt(code) = joinpath(pwd(), "data", string(uppercase(code), ".pdb"))
xml(code) = joinpath(pwd(), "data", string(uppercase(code), ".xml"))

print("""

Tests for PDB
=============
""")

print("""

Parse PDB and PDBML
-------------------
""")

print("""
2VQC => Missings
""")

let code = "2VQC"
    pdb = read(txt(code), PDBFile)
    pdbml = read(xml(code), PDBML)
    @test findfirst(x -> x.id.number == "4",  pdb) == findfirst(x -> x.id.number == "4",  pdbml)
    @test findfirst(x -> x.id.number == "73", pdb) == findfirst(x -> x.id.number == "73", pdbml)

end

print("""
Test download
""")

let code = "2VQC"
    pdb = read(txt(code), PDBFile)
    filename = downloadpdb(code)
    try
        pdbml = read(filename, PDBML)
        @test findfirst(x -> x.id.number == "4",  pdb) == findfirst(x -> x.id.number == "4",  pdbml)
        @test findfirst(x -> x.id.number == "73", pdb) == findfirst(x -> x.id.number == "73", pdbml)
    finally
        rm(filename)
    end
end

print("""
1H4A => Chain A (auth) == Chain X (label)
""")

let file = string(xml("1H4A"), ".gz"),
    auth = read(file, PDBML, label=false),
    label = read(file, PDBML)

    @test unique(ASCIIString[ res.id.chain for res in auth ]) == ["X"]
    @test unique(ASCIIString[ res.id.chain for res in label]) == ["A", "B"]

end

print("""
1SSX => Residues with insert codes: 15A 15B
""")

let code = "1SSX",
    pdb = read(txt(code), PDBFile),
    pdbml = read(xml(code), PDBML)

    @test findobjects(pdbml, Is(:number, "15A"))[1] == 1
    @test findobjects(pdbml, Is(:number, "15B"))[1] == 2

    @test findobjects(pdb, Is(:number, "15A"))[1] == 1
    @test findobjects(pdb, Is(:number, "15B"))[1] == 2

    print("""
    test @residues
    """)
    @test (@residues pdb model "*" chain "*" group "*" residue "141")[1] == collectobjects(pdb, Is(:number, "141"))[1]
    # Testing the macro in let block:
    mo = "1"
    ch = "A"
    gr = "*"
    re = "141"
    @test (@residues pdb model mo chain ch group gr residue re)[1] == (@residues pdb model "*" chain "*" group gr residue "141")[1]

    print("""
    Occupancy != 1.0 and @atom
    """)
    @test sum([ get(occ,0) for occ in  collectcaptures(collectobjects(pdbml, Is(:number, "141"))[1].atoms, :occupancy, Is(:atom, "HH22")) ]) == 1.0
    @test sum( [ atom.occupancy for atom in @atoms pdbml model "1" chain "A"  group "*" residue "141" atom "HH22" ] ) == 1.0
    # Testing the macro in let block:
    at = "HH22"
    @test sum( [ atom.occupancy for atom in @atoms pdbml model mo chain ch group gr residue re atom at ] ) == 1.0

    print("""
    Best occupancy
    """)

    atoms_141 = @atoms    pdbml model "1" chain "A"  group "*" residue "141" atom "HH22"
    resid_141 = @residues pdbml model "1" chain "A"  group "*" residue "141"

    @test bestoccupancy(atoms_141)[1].occupancy == 0.75
    @test bestoccupancy(reverse(atoms_141))[1].occupancy == 0.75
    @test bestoccupancy(PDBAtom[ atoms_141[2] ])[1].occupancy == 0.25

    @test length(resid_141[1]) == 48
    @test selectbestoccupancy(resid_141[1], collect(1:48)) == 1
    @test selectbestoccupancy(resid_141[1], [1, 2]) == 1

    @test_throws ErrorException selectbestoccupancy(resid_141[1], Int[])
    @test_throws ErrorException selectbestoccupancy(resid_141[1], collect(1:100))

    print("""
    @atom with residue wildcard
    """)

    # ATOM      2  CA  ALA A  15A     22.554  11.619   6.400  1.00  6.14           C
    @test atoms(pdb, mo, ch, "ATOM", "*", r"C.+")[1].atom == "CA"

end

print("""
`read` only atoms with the best occupancy (occupancyfilter=true)
""")

let code = "1SSX",
    pdb = read(txt(code), PDBFile, occupancyfilter=true),
    pdbml = read(xml(code), PDBML, occupancyfilter=true),

    res_pdb   = @residues pdbml model "1" chain "A"  group "*" residue "141"
    res_pdbml = @residues pdbml model "1" chain "A"  group "*" residue "141"

    atm_pdbml = @atoms pdbml model "1" chain "A"  group "*" residue "141" atom "HH22"

    @test length( atm_pdbml ) == 1
    @test atm_pdbml[1].occupancy == 0.75
    @test length( @atoms pdb   model "1" chain "A"  group "*" residue "141" atom "HH22" ) == 1

    @test length(res_pdb[1])   == 24
    @test length(res_pdbml[1]) == 24

end

print("""
1CBN => Identical PDBe ResNum for Residue 22:

        <residue dbSource="PDBe" dbCoordSys="PDBe" dbResNum="22" dbResName="SER">
          <crossRefDb dbSource="PDB" dbCoordSys="PDBresnum" dbAccessionId="1cbn" dbResNum="22" dbResName="SER" dbChainId="A"/>
          <crossRefDb dbSource="UniProt" dbCoordSys="UniProt" dbAccessionId="P01542" dbResNum="22" dbResName="P"/>
          ....
        <residue dbSource="PDBe" dbCoordSys="PDBe" dbResNum="22" dbResName="PRO">
          <crossRefDb dbSource="PDB" dbCoordSys="PDBresnum" dbAccessionId="1cbn" dbResNum="22" dbResName="PRO" dbChainId="A"/>
          <crossRefDb dbSource="UniProt" dbCoordSys="UniProt" dbAccessionId="P01542" dbResNum="22" dbResName="P"/>
          ...

""")

let code = "1CBN"
    pdb = read(txt(code), PDBFile)
    pdbml = read(xml(code), PDBML)

    @test length( @residues pdb   model "1" chain "A" group "*" residue "22" ) == 2
    @test length( @residues pdbml model "1" chain "A" group "*" residue "22" ) == 2

    @test ASCIIString[ res.id.name for res in  @residues pdb   model "1" chain "A" group "*" residue "22" ] == ["SER", "PRO"]
    @test ASCIIString[ res.id.name for res in  @residues pdbml model "1" chain "A" group "*" residue "22" ] == ["SER", "PRO"]
end

print("""
1AS5 => NMR
""")

let code = "1AS5"
    pdb = read(txt(code), PDBFile)
    pdbml = read(xml(code), PDBML)

    @test length( @residues pdbml model "1"  chain "A" group "*" residue "*" ) == 25
    @test length( @residues pdbml model "14" chain "A" group "*" residue "*" ) == 25

    @test length( @residues pdbml model "*"  chain "A" group "*" residue "*" ) == 25*14
end

print("""
1DPO => Inserted residues lack insertion letters
Single unnamed chain in 1DPO contains insertions at postions 184 (Gly, Phe), 188 (Gly, Lys), and 221 (Ala, Leu) but no insertion letters.
""")

let code = "1DPO"
    pdb = read(txt(code), PDBFile)
    pdbml = read(xml(code), PDBML)

    # Single "A" chain for PDB (auth_asym_id in PDBML)
    @test unique(ASCIIString[ res.id.chain for res in pdb  ]) == ASCIIString[ "A" ]
    # But 'A':'H' chains for PDBML (label_asym_id)
    @test unique(ASCIIString[ res.id.chain for res in pdbml]) == ASCIIString[ string(chain) for chain in 'A':'H' ]

    @test ASCIIString[ res.id.name for res in @residues pdb   model "1" chain "A" group "*" residue r"^184[A-Z]?$" ] == ASCIIString["GLY", "PHE"]
    @test ASCIIString[ res.id.name for res in @residues pdbml model "1" chain "A" group "*" residue r"^184[A-Z]?$" ] == ASCIIString["GLY", "PHE"]

    @test ASCIIString[ res.id.name for res in @residues pdb   model "1" chain "A" group "*" residue r"^188[A-Z]?$" ] == ASCIIString["GLY", "LYS"]
    @test ASCIIString[ res.id.name for res in @residues pdbml model "1" chain "A" group "*" residue r"^188[A-Z]*" ] == ASCIIString["GLY", "LYS"]

    @test ASCIIString[ res.id.name for res in @residues pdb   model "1" chain "A" group "*" residue r"^221[A-Z]?$" ] == ASCIIString["ALA", "LEU"]
    @test ASCIIString[ res.id.name for res in @residues pdbml model "1" chain "A" group "*" residue r"^221[A-Z]?$" ] == ASCIIString["ALA", "LEU"]
end

print("""
1IGY => Insertions have more than one copy of the same amino acid in a single insertion block.
For example, chain B in 1IGY contains a block of four residues inserted at sequence position 82. The block contains Leu-Ser-Ser-Leu.
""")

let code = "1IGY"
    pdb = read(txt(code), PDBFile)
    pdbml = read(xml(code), PDBML)

    @test ASCIIString[ res.id.name for res in @residues pdb   model "1" chain "B" group "*" residue r"^82[A-Z]?$" ] == ["LEU", "SER", "SER", "LEU"]
    @test ASCIIString[ res.id.name for res in @residues pdbml model "1" chain "B" group "*" residue r"^82[A-Z]?$" ] == ["LEU", "SER", "SER", "LEU"]

    @test sum( [ res.id.group for res in @residues pdb model "1" chain "D" group "*" residue "*" ] .== "HETATM" ) == length( @residues pdb model "1" chain "D" group "HETATM" residue "*" )
    @test sum( [ res.id.group for res in @residues pdb model "1" chain "D" group "*" residue "*" ] .== "ATOM" ) == length( @residues pdb model "1" chain "D" group "ATOM" residue "*" )
end

print("""
1HAG => Chain E begins with 1H, 1G, 1F, ... 1A, then 1 (in reverse alphabetic order)
""")

let code = "1HAG"
    pdb = read(txt(code), PDBFile)
    pdbml = read(xml(code), PDBML)

    @test unique(ASCIIString[ res.id.chain for res in pdb   ]) == ["E", "I"]
    @test unique(ASCIIString[ res.id.chain for res in pdbml ]) == ["A", "B", "C", "D", "E"]

    # The chain E of PDB is the chain A of PDBML
    @test ASCIIString[ res.id.number for res in @residues pdb   model "1" chain "E" group "*" residue r"^1[A-Z]?$" ] == ASCIIString[ string(1, code) for code in vcat(collect('H':-1:'A'), "") ]
    @test ASCIIString[ res.id.number for res in @residues pdbml model "1" chain "A" group "*" residue r"^1[A-Z]?$" ] == ASCIIString[ string(1, code) for code in vcat(collect('H':-1:'A'), "") ]
end

print("""
1NSA => Contains a single (unnamed) protein chain with sequence 7A-95A that continues 4-308.
""")

let code = "1NSA"
    pdb = read(txt(code), PDBFile)
    pdbml = read(xml(code), PDBML)

    # Single "A" chain for PDB (auth_asym_id in PDBML)
    @test unique(ASCIIString[ res.id.chain for res in pdb  ]) == ASCIIString[ "A" ]
    # But 'A':'F' chains for PDBML (label_asym_id)
    @test unique(ASCIIString[ res.id.chain for res in pdbml]) == ASCIIString[ string(chain) for chain in 'A':'F' ]

    ind = findobjects(pdbml, Is(:number, "95A"))[1]
    @test pdbml[ind + 1].id.number == "4"

    ind = findobjects(pdb, Is(:number, "95A"))[1]
    @test pdb[ind + 1].id.number == "4"
end

print("""
1IAO => Contains in chain B (in this order) 1S, 2S, 323P-334P, 6-94, 94A, 95-188, 1T, 2T
""")

let code = "1IAO"
    pdb = read(txt(code), PDBFile)
    pdbml = read(xml(code), PDBML)

    pdb_B   = @residues pdb   model "1" chain "B" group "*" residue "*"
    pdbml_B = @residues pdbml model "1" chain "B" group "*" residue "*"

    @test pdb_B[   findobjects(pdb_B,   Is(:number, "2S"))[1]   + 1].id.number == "323P"
    @test pdb_B[   findobjects(pdb_B,   Is(:number, "334P"))[1] + 1].id.number == "6"
    @test pdb_B[   findobjects(pdb_B,   Is(:number, "94"))[1]   + 1].id.number == "94A"
    @test pdb_B[   findobjects(pdb_B,   Is(:number, "94A"))[1]  + 1].id.number == "95"
    @test pdb_B[   findobjects(pdb_B,   Is(:number, "188"))[1]  + 1].id.number == "1T"
    @test pdb_B[   findobjects(pdb_B,   Is(:number, "1T"))[1]   + 1].id.number == "2T"

    @test pdbml_B[ findobjects(pdbml_B, Is(:number, "2S"))[1]   + 1].id.number == "323P"
    @test pdbml_B[ findobjects(pdbml_B, Is(:number, "334P"))[1] + 1].id.number == "6"
    @test pdbml_B[ findobjects(pdbml_B, Is(:number, "94"))[1]   + 1].id.number == "94A"
    @test pdbml_B[ findobjects(pdbml_B, Is(:number, "94A"))[1]  + 1].id.number == "95"
    @test pdbml_B[ findobjects(pdbml_B, Is(:number, "188"))[1]  + 1].id.number == "1T"
    @test pdbml_B[ findobjects(pdbml_B, Is(:number, "1T"))[1]   + 1].id.number == "2T"
end

print("""

RESTful PDB Interface
=====================
""")

@test getpdbdescription("4HHB")["resolution"] == "1.74"
@test_throws KeyError getpdbdescription("104D")["resolution"] # NMR

print("""

Write PDB files
===============
""")

let code = "2VQC", io = IOBuffer()
    pdb = read(txt(code), PDBFile)
    print(io, pdb, PDBFile)
    printed = split(takebuf_string(io), '\n')

    @test length(printed) == 609 # Only ATOM, HETATM & END + 1 because the trailing \n
    @test printed[1]   == "ATOM      1  N   THR A   4       2.431  19.617   6.520  1.00 24.37           N  "
    @test printed[607] == "HETATM  607  O   HOH A2025      13.807  38.993   2.453  1.00 33.00           O  "
end

let code = "1AS5", # NMR
    io = IOBuffer()

    pdb = read(txt(code), PDBFile)
    print(io, pdb, PDBFile)
    printed = split(takebuf_string(io), '\n')

    @test sum(map(x -> startswith(x, "MODEL "), printed)) == 14 # 14 models
    @test sum(map(x -> x == "ENDMDL", printed)) == 14 # 14 models
end

let code = "1IAO", # 2 Chains
    io = IOBuffer()

    pdb = read(txt(code), PDBFile)
    print(io, pdb, PDBFile)
    printed = split(takebuf_string(io), '\n')

    # MIToS only prints TER for the ATOM group if the chain changes.
    # Some modified residues are annotated as HETATM in the middle of the ATOM chain:
    # TER can not be printed from ATOM to HETATM if the chain doesn’t change.

    # Only prints TER between chain A and B
    @test sum(map(x -> startswith(x, "TER "), printed)) == 1

    @test filter(r"TER ", printed)[1] == "TER    1418      TRP A 178 "
end

print("""

read/write consistency
----------------------
""")

let io = IOBuffer()
    for code in ["2VQC", "1IAO", "1NSA", "1HAG", "1IGY", "1DPO", "1AS5", "1CBN", "1SSX"]
        println(code)

        readed = read(txt(code), PDBFile)
        print(io, readed, PDBFile)
        readed_writed_readed = parse(takebuf_string(io), PDBFile)

        @test readed_writed_readed == readed
    end
end


text:
# this file is to run JWAS
using DataFrames,CSV,Random,Distributions,JWAS,DelimitedFiles, LinearAlgebra
Random.seed!(123)

method            = ARGS[1]
rep               = parse(Int, ARGS[2])  #rep1,...,rep10
pct               = parse(Float64, ARGS[3])

@show method,rep,pct

mainpath="/group/qtlchenggrp/tianjing/omics_sample_all/nonlinear/varmi_2/random_missing/"
cd(mainpath*"$pct/$method/rep$rep/")

############ READ DATA ##########
#read phenotype
phenofile="/group/qtlchenggrp/tianjing/omics_data/nonlinear/varmi_2/data1/y.csv"
phenotypes = CSV.read(phenofile,DataFrame,delim = ',',header=true,missingstrings=["NA"])

#read omics
omicsfile  = "/group/qtlchenggrp/tianjing/omics_data/nonlinear/varmi_2/data1/partial_missing_data/omics_partialmissing.pct$pct.rep$rep.csv"
omics  = CSV.read(omicsfile,DataFrame,delim = ',',header=true,missingstrings=["NA"]);

insertcols!(omics,2,:y => phenotypes[:,:y]);
omics[!,:y]= convert(Vector{Union{Missing,Float64}}, omics[!,:y]);
omics[956:1055,:y] .= missing  #testing dataset

@show [count(ismissing,col) for col in eachcol(omics)]

omics_name=["gene$i" for i in 1:1200]

@show size(omics)
@show omics[1:5,1:5];
@show omics[end-4:end,1:5];

#read genotype
if occursin("gblup",method)  #NN-GBLUP
	genofile   = "/group/qtlchenggrp/tianjing/omics_data/andrew1055/data1/GRM_geno.csv"
	genotypes  = get_genotypes(genofile;separator=',',method="GBLUP",header=false);
end

model_equation = "y = intercept + genotypes";

model = build_model(model_equation,
                    num_hidden_nodes=1200,
                    latent_traits=omics_name,
                    nonlinear_function="sigmoid");  #"linear" for linear activation function
out   = runMCMC(model,omics,chain_length=5000,printout_model_info=false);

max_stars_repo_name: zhaotianjing/NN-LMM

max_issues_repo_name: zhaotianjing/NN-LMM

max_forks_repo_name: zhaotianjing/NN-LMM


text:
function get_limit_vars(cur_scan::AbstractScan, cur_limit::Symbol, cur_field::Symbol; cur_first::Union{Void,Int}=nothing, cur_last::Union{Void,Int}=nothing)
  cur_field_symbol = Symbol("$(cur_limit)_reactors")

  cur_reactor_list = getfield(cur_scan, cur_field_symbol)

  cur_values = []

  for cur_reactor in cur_reactor_list
    cur_reactor.is_good || continue

    cur_value = nothing

    if isdefined(cur_reactor, cur_field)
      cur_value = getfield(cur_reactor, cur_field)
    elseif isdefined(Fussy, cur_field)
      cur_value = getfield(Fussy, cur_field)(cur_reactor)
    else
      getfield(cur_reactor, cur_field) # throw error
    end

    push!(cur_values, cur_value)
  end

  cur_length = length(cur_values)

  if cur_first != nothing && cur_first < cur_length
    return cur_values[1:cur_first]
  elseif cur_last != nothing && cur_last < cur_length
    return cur_values[end-(cur_last-1):end]
  end

  cur_values
end


text:
#!/bin/bash
#=
JULIA="${JULIA:-julia}"
JULIA_CMD="${JULIA_CMD:-$JULIA --color=yes --startup-file=no}"
# below gets the directory name of the script, even if there is a symlink involved
# from https://stackoverflow.com/questions/59895/get-the-source-directory-of-a-bash-script-from-within-the-script-itself
SOURCE="${BASH_SOURCE[0]}"
while [ -h "$SOURCE" ]; do # resolve $SOURCE until the file is no longer a symlink
  DIR="$( cd -P "$( dirname "$SOURCE" )" >/dev/null 2>&1 && pwd )"
  SOURCE="$(readlink "$SOURCE")"
  [[ $SOURCE != /* ]] && SOURCE="$DIR/$SOURCE" # if $SOURCE was a relative symlink, we need to resolve it relative to the path where the symlink file was located
done
DIR="$( cd -P "$( dirname "$SOURCE" )" >/dev/null 2>&1 && pwd )"
export JULIA_PROJECT=$DIR/..
export JULIA_LOAD_PATH=@:@stdlib  # exclude default environment
exec $JULIA_CMD -e 'include(popfirst!(ARGS))' "$SOURCE" "$@"
=#
using ArgParse
s = ArgParseSettings()
@add_arg_table s begin
    "pulse_file"
        help = "name of the pulse containing ljh file to use to make basis"
        required = true
        arg_type = String
    "noise_file"
        help = "name of the noise_analysis hdf5 file to use to make basis"
        required = true
        arg_type = String
    "--outputfile", "-o"
        arg_type = String
        help="specify the path of the outputfile, otherwise it will make one up based on pulse_file"
    "--replaceoutput", "-r"
        help = "if this is present, it will overwrite any existing output files with the same name as it plans to use"
        action = :store_true
    "--n_pulses_for_train"
        arg_type = Int
        default = 3000
        help = "number of pulses for training"
    "--n_basis"
        arg_type = Int
        default = 6
        help = "number of basis vectors to calculate"
    "--n_loop"
        arg_type = Int
        default = 5
        help = "number of training loops"
    "--frac_keep"
        arg_type = Float64
        default = 0.8
        help ="""The fraction of training pulses to keep. Each loop cuts a fraction of pulses
        with the highest residuals, after all loops, this fraction of pulses remain uncut."""
    "--tsvd_method"
        default = "noisemass3"
        help = """which truncated SVD method to use, supports `TSVDmass3`, `noisemass3`, `TSVD`, and `full`.
        The default is `noisemass3`, in which the usual MASS method is applied (find the pulse
        average, a time-correction, and a constant, and use SVD on the residuals to find any
        additional components). The `TVSDmass3` is like `noisemass3` but uses approximate ARMA models.
        The other 2 methods just apply (truncated) SVD to the training data.
        The results of the last two should be nearly identical, and `TSVD` is faster.
        But you can try `full`, which computes the full SVD and retains the leading `n_basis` vectors,
        as a sanity check if the basis vectors look weird.
        Note that for methods TSVDmass3 and noisemass3, the value of `n_basis` must be at least 3."""
    "--maxchannels"
        default = 1000000
        help = "process at most this many channels, in order from lowest channel number"
        arg_type = Int
end
parsed_args = parse_args(ARGS, s)
using Pope.NoiseAnalysis
using Pope.LJH
using Pope
using HDF5
if parsed_args["outputfile"]==nothing
    parsed_args["outputfile"] = Pope.outputname(parsed_args["pulse_file"],"model")
end
display(parsed_args);println()
if !parsed_args["replaceoutput"] && isfile(parsed_args["outputfile"])
    println("intended output file $(parsed_args["outputfile"]) exists, pass --replaceoutput if you would like to overwrite it")
    exit(1) # anything other than 0 indicated process unsuccesful
end

ljhdict = LJH.allchannels(parsed_args["pulse_file"],parsed_args["maxchannels"]) # ordered dict mapping channel number to filename
outputh5 = h5open(parsed_args["outputfile"],"w")
Pope.make_basis_all_channel(outputh5, ljhdict, parsed_args["noise_file"],
    parsed_args["frac_keep"],
    parsed_args["n_loop"],
    parsed_args["n_pulses_for_train"],
    parsed_args["n_basis"],
    parsed_args["tsvd_method"])
close(outputh5)


text:
@testset "Fill NaNs" begin
	# NaN will be replaced by interpolated value
	data = DataFrame(Temp=[10.,11.,12.,14.,NaN, 14.,15.,16],
		Grav=[1.,2.,NaN,NaN,NaN,6.,NaN,8.],
		datetime=collect(DateTime(2010,1,1,1):Dates.Hour(1):DateTime(2010,1,1,8)));
	out = fillnans(data[:Temp],2);
	@test out[5] ≈ 14.
	@test all(.!isnan.(out)) == true
	# To find indices corresponding to corrected values just search for difference
	# `fillnans` was called:
	corrindex = findall(isnan.(data[:Temp]) .& .!isnan.(out)); # will return [5]
	@test corrindex == [5]

	# NaNs will NOT be replaced as the missing window is too long (3>2), except
	# for the last NaN
	out2 = fillnans(data[:Grav],2);
	@test findall(isnan.(data[:Grav]) .& .!isnan.(out2)) == [7]
	@test out2[7] ≈ 7.0

	# NaNs will be replace as window is longer (or equal to missing gap)
	out3 = fillnans(data[:Grav],3);
	@test findall(isnan.(data[:Grav]) .& .!isnan.(out3)) == [3,4,5,7]
	@test sum(out3[[3,4,5,7]]) ≈ 3+4+5+7.0
end

@testset "Replace NaNs and NAs" begin
	data1 = DataFrame(Temp=[10.,NaN,30.,40.],
					  pres=[1.0,2.0,NaN,4.0],
		   	datetime=[DateTime(2010,1,1,0),DateTime(2010,1,1,1),
			 DateTime(2010,1,1,2),DateTime(2010,1,1,4)]);
 	replacenans!(data1,0.0)
	@test data1[:Temp] == [10.,0.0,30.,40.]
	@test data1[:pres] == [1.0,2.0,0.0,4.0]

	data1 = DataFrame(Temp=[10,missing,30,40],
					  pres=[1.0,2.0,NaN,4.0],
		   	datetime=[DateTime(2010,1,1,0),DateTime(2010,1,1,1),
			 DateTime(2010,1,1,2),DateTime(2010,1,1,4)]);
 	na2nan!(data1)
	@test data1[:Temp][[1,3,4]] == [10,30,40]
	@test ismissing(data1[:Temp][2])
	@test data1[:pres][[1,2,4]] == [1.,2.,4.]
	@test !ismissing(data1[:pres][3])
	@test isnan(data1[:pres][3])
end

@testset "Convert missing to NaNs" begin
	df = DataFrame(datetime=[DateTime(2000),DateTime(2001)],
				datetime2 = [DateTime(2002),missing],
				grav=[missing,3.],
				pres=[1,missing],
				temp=[missing,missing]);
	df[:datetime3] = convert(Array{Union{Missing, DateTime},1},df[:datetime])

	@test missing2nan(df[:datetime]) == [DateTime(2000),DateTime(2001)]
	@test missing2nan(df[:datetime3]) == [DateTime(2000),DateTime(2001)]
	@test missing2nan(df[:datetime2])[1] == DateTime(2002)
	@test ismissing(missing2nan(df[:datetime2])[2])
	@test missing2nan(df[:grav])[2] == 3.
	@test missing2nan(df[:pres])[1] == 1.
	@test isnan(missing2nan(df[:grav])[1])
	@test isnan(missing2nan(df[:pres])[2])
	@test isnan(missing2nan(df[:temp])[1])
	@test isnan(missing2nan(df[:temp])[1])

	# input should stay unchanged
	@test ismissing(df[:grav][1])
	@test df[:pres][1] == 1
	@test eltype(typeof(df[:temp])) == Missing

	# test for dataframe input
	dfout = missing2nan(df)

	@test missing2nan(dfout[:datetime]) == [DateTime(2000),DateTime(2001)]
	@test missing2nan(dfout[:datetime3]) == [DateTime(2000),DateTime(2001)]
	@test missing2nan(dfout[:datetime2])[1] == DateTime(2002)
	@test ismissing(missing2nan(dfout[:datetime2])[2])
	@test missing2nan(dfout[:grav])[2] == 3.
	@test missing2nan(dfout[:pres])[1] == 1.
	@test isnan(missing2nan(dfout[:grav])[1])
	@test isnan(missing2nan(dfout[:pres])[2])
	@test isnan(missing2nan(dfout[:temp])[1])
	@test isnan(missing2nan(dfout[:temp])[1])

	# input should stay unchanged
	@test ismissing(df[:grav][1])
	@test df[:pres][1] == 1
	@test eltype(typeof(df[:temp])) == Missing
end


text:
type MarginDataView <: Grid
  handle :: Ptr{GObject}
  ItemName :: Label
  MaxBuy :: Label
  MinSell :: Label
  Margin :: Label
  function MarginDataView()
    grid = @Grid()
    self = Gtk.gobject_move_ref(new(grid.handle),grid)
    setproperty!(self, "orientation", 1)
    self.ItemName = @Label("No Item")
    self.MaxBuy = @Label("Max Buy: -")
    self.MinSell = @Label("Min Sell: -")
    self.Margin = @Label("Margin: -")
    push!(self, self.ItemName)
    push!(self, self.MaxBuy)
    push!(self, self.MinSell)
    push!(self, self.Margin)
    self
  end
end
function UpdateMarginDataView(view, data :: ItemSummery)
  setproperty!(view.ItemName, "label", string(data.itemID))
  setproperty!(view.MaxBuy, "label", data.maxbuy)
  setproperty!(view.MinSell, "label", data.minsell)
  setproperty!(view.Margin, "label", data.margin)
end

max_stars_repo_name: barcharcraz/eve-julia-gui

max_issues_repo_name: barcharcraz/eve-julia-gui

max_forks_repo_name: barcharcraz/eve-julia-gui


text:
"""
    rope_pos_tension(dF,Ls, Ns, ds, diameter_rope, Cs, ρ_air, v, β₀)
                    -> Vector{PositionTuple} , Vector{Force}
Rope positions and tensions, starting at the unknown position anchor,
ending at the known position anchor. Should resemble a catenary.

TODO: Currently, rope mass is neglected. Make wind + inerta an input function df(ds, β) -> ForceTuple. Check against catenary.

Ls       scalar length of rope.
Ns       Number of rope sections (one less than rope nodes).
diameter_rope
Cs      shape coefficient, for wind in the section plane.
ρ_air   density of air at design temperature.
v       scalar velocity of relative wind referred to the rope, not reduced by projection
α_w_rel direction of relative wind referred to the rope section, measured from global x axis around z.
Fk      Force acting on start of rope, a ForceTuple like (10.0N, 100.0N) referred to the chord coordinate system
α_chord Angle from global x axis to the chord x axis, positive from leading edge to trailing edge
"""
function rope_pos_tension(Ls, Ns, diameter_rope, Cs, ρ_air, v, α_w_rel, Fk, α_chord)

    # Length of rope section
    ds = Ls / Ns
    # Most parameters are constant, make a short form function
    dFthis(β) = dF(ds, diameter_rope, Cs, ρ_air, v, α_w_rel, β)
    # Rope node positions, origo at start
    ps = Vector{typeof((0.0m, 0.0m))}()

    # Rope tension at nodes
    Ts = Vector{typeof(0.0N)}()

    # Initial values
    β₀  = α_chord + atan(Fk[2], Fk[1]) - π |> °
    x₀, y₀ = 0.0m, 0.0m
    β₋₁ = β₀
    Ts₀ = hypot(Fk) |> N
    curL = 0.0m
    while curL < Ls
        push!(ps, (x₀, y₀))
        push!(Ts, Ts₀)
        # Estimate outgoing angle
        βe =  β₀ + β₀ - β₋₁
        # External force vector on this piece
        df = 0.5dFthis(β₀) + 0.5dFthis(βe)
        # Force vector, kite end
        F₀ = -Ts₀ .* (cos(β₀), sin(β₀))
        # Force vector, buggy end, from force equilibrium
        F₁ = -(F₀ + df)
        # Check that this is actually tension
        @assert hypot(F₀ + F₁) < hypot(F₁)  "Failed simple tension check, hypot(F₀ + F₁) < hypot(F₁) at $curL"
        # Tension, buggy end
        Ts₁ = hypot(F₁)
        # Angle, buggy end
        β₁ = atan(F₁[2], F₁[1]) |> °
        # Averaged angle over span
        βm = (2β₀ + β₁) / 3
        # Position, buggy end (moment equilibrium would be more exact, tension could be considered)
        (x₁, y₁) = (x₀, y₀) + ds .* (cos(βm), sin(βm))
        # Prepare next iteration
        β₋₁ = β₀
        x₀, y₀ = x₁, y₁
        Ts₀ = Ts₁
        curL += ds
        β₀ = β₁
    end
    push!(ps, (x₀, y₀))
    push!(Ts, Ts₀)
    ps, Ts
end

"
projection(α_w_rel, β)

Wind projection factor for a straight section of rope. Covers
    1) Shorter projected length of section.
    2) Slower projected wind speed. This effect is squared, as in velocity normal to the projection is squared.

α_w_rel   direction of relative wind referred to the rope section, from x axis around z.
β         axial orientation of rope section, from x axis around z.
"
function projection(α_w_rel, β)
    angle = α_w_rel  - (β + 90°)
    (abs(cos(angle)))^3
end

"
    dq(diameter_rope, Cs, ρ_air, v, α_w_rel, β)
Scalar load per length due to wind acting normal to a straight rope section

diameter_rope
Cs      shape coefficient, for wind in the section plane
ρ_air   density of air at design temperature
v       scalar velocity of relative wind referret to the rope
α_w_rel direction of relative wind referred to the rope section, from x axis around z.
β       axial orientation of rope section, from x axis around z.
"
dq(diameter_rope, Cs, ρ_air, v, α_w_rel, β) =  0.5ρ_air * Cs * diameter_rope * v^2 * projection(α_w_rel, β) |> N*mm^-1

"
    dF(ds, diameter_rope, Cs, ρ_air, v, β)
    --> (force_x, force_y)::ForceTuple

Wind load, globally oriented force vector per length of straight rope section

ds scalar length of rope section
diameter_rope
Cs      shape coefficient, for wind in the section plane
ρ_air   density of air at design temperature
v       scalar velocity of relative wind referret to the rope
α_w_rel direction of relative wind referred to the rope section, from x axis around z.
β  axial orientation of rope section relative to x axis
"
function dF(ds, diameter_rope, Cs, ρ_air, v, α_w_rel, β)
    dQ = dq(diameter_rope, Cs, ρ_air, v, α_w_rel, β)
    ds * dQ  .* (cos(β + π / 2 ), sin(β+ π / 2)) .|> N
end

const ρ_SK75 = 0.975 * 0.001kg/(10mm)^3
const σ_TS_SK75 = 3600N/mm^2

"
Fill factors reduce the fibre area compared to the circumscribed circle, see 'test_17.jl'.
"
fill_factor(d) = (-0.003*d/mm + 0.853)

"""
Spin factors reduce the fibre breaking strength, see 'test_17.jl'.
"""
spin_factor(d) = 5*10^-5 * (d/mm)^2 + -0.0059*(d/mm) + 0.5504

"""
Circumscribed area as a function of diameter
"""
circle_area(d) = π / 4 * d^2

"""
Area filled with fibres as a function of diameter, see 'test17.jl'.
"""
area_filled(d) = fill_factor(d)* circle_area(d)

"""
Assuming no fibre breaking strenght reduction, but including the fibre fill factor
"""
ideal_spin_strength(d) = area_filled(d) * σ_TS_SK75 |> kN

"""
Rope breaking strength as a function of diameter, see 'test17.jl'
"""
rope_breaking_strength(d) = spin_factor(d) * ideal_spin_strength(d)

"""
Rope weight per length as a function of diameter, see 'test17.jl'
"""
rope_weight(d) = ρ_SK75 * area_filled(d) |> kg/m


text:
import Base: @deprecate, @deprecate_binding, depwarn

###
### v0.15 deprecations
###

### Add empty exists method for JLD,MAT to extend to smooth over deprecation process PR#790
export exists
function exists end

### Changed in PR#776
@deprecate create_dataset(parent::Union{File,Group}, path::AbstractString, dtype::Datatype, dspace::Dataspace,
    lcpl::Properties, dcpl::Properties, dapl::Properties, dxpl::Properties) HDF5.Dataset(HDF5.h5d_create(parent, path, dtype, dspace, lcpl, dcpl, dapl), HDF5.file(parent), dxpl) false

### Changed in PR#798
@deprecate get_dims(dspace::Union{Dataspace,Dataset,Attribute}) get_extent_dims(dspace) false
@deprecate set_dims!(dset::Dataspace) set_extent_dims(dset) false


text:
# This file is a part of SimilaritySearch.jl

using SimilaritySearch
using Test

@testset "Result set" begin
    k = 4
    V = rand(Float32, 50)
    Vsorted = sort!([i => v for (i, v) in enumerate(V)], by=x->x[2])[1:k]
    res = KnnResult(k)

    for i in eachindex(V)
        push!(res, i, V[i])
    end

    @show res Vsorted
    @test minimum(res) == first(Vsorted)[end]
    @test maximum(res) == last(Vsorted)[end]
    @test argmin(res) == first(Vsorted)[1]
    @test argmax(res) == last(Vsorted)[1]

    @test collect(keys(res)) == first.(Vsorted)
    @test collect(values(res)) == last.(Vsorted)
end