text:
module TestJoin

using FactCheck
using DataCubes
using DataFrames

facts("Join tests") do
  context("larr version tests") do
    lj_base_tbl = @larr(a=reshape(collect(1:20), 5,4), b=reshape(collect(101:120),5,4), axis1[k1=[1,2,3,4,5], k2=['a','a','b','c','d']], axis2[r1=[10,12,14,16]])
    lj_src_tbl = @larr(c=reshape(collect(51:56),3,2), axis1[k1=[5,3,1], k2=['d','a','a']], axis2[t1=["x","y"]])
    @fact leftjoin(lj_base_tbl, lj_src_tbl, 1) --> @larr(a=reshape(repmat(collect(1:20),2),5,4,2),b=reshape(repmat(collect(101:120),2),5,4,2),c=reshape(@nalift([53,NA,NA,NA,51,53,NA,NA,NA,51,53,NA,NA,NA,51,53,NA,NA,NA,51,56,NA,NA,NA,54,56,NA,NA,NA,54,56,NA,NA,NA,54,56,NA,NA,NA,54]),5,4,2),axis1[k1=collect(1:5),k2=['a','a','b','c','d']],axis2[r1=[10,12,14,16]],axis3[t1=["x","y"]])
    @fact innerjoin(lj_base_tbl, lj_src_tbl, 1) --> @larr(a=reshape(repmat([1,5,6,10,11,15,16,20],2),2,4,2),b=reshape(repmat([101,105,106,110,111,115,116,120],2),2,4,2),c=reshape([repmat([53,51],4);repmat([56,54],4)],2,4,2),axis1[k1=[1,5],k2=['a','d']],axis2[r1=[10,12,14,16]],axis3[t1=["x","y"]])
    @fact leftjoin(larr(a=[1,2,3]), larr(b=[2 3;0 1],axis1=[:x,:y])) --> larr(a=reshape(repmat([1,2,3],4),3,2,2),b=reshape([2,2,2,0,0,0,3,3,3,1,1,1],3,2,2),axis2=[:x,:y])
    @fact innerjoin(larr(a=[1,2,3]), larr(b=[2 3;0 1],axis1=[:x,:y])) --> larr(a=reshape(repmat([1,2,3],4),3,2,2),b=reshape([2,2,2,0,0,0,3,3,3,1,1,1],3,2,2),axis2=[:x,:y])
    @fact leftjoin(lj_base_tbl, lj_src_tbl, 1=>pick(lj_base_tbl,:k1,:k2)) --> leftjoin(lj_base_tbl, lj_src_tbl, 1=>pick(lj_base_tbl,[:k1,:k2]))
    @fact leftjoin(lj_base_tbl, lj_src_tbl, 1=>pick(lj_base_tbl,:k1,:k2)) --> leftjoin(lj_base_tbl, lj_src_tbl, 1)
    @fact leftjoin(lj_base_tbl, lj_src_tbl, 1=>1) --> leftjoin(lj_base_tbl, lj_src_tbl, 1)
    @fact innerjoin(lj_base_tbl, lj_src_tbl, 1=>pick(lj_base_tbl,:k1,:k2)) --> innerjoin(lj_base_tbl, lj_src_tbl, 1=>pick(lj_base_tbl,[:k1,:k2]))
    @fact innerjoin(lj_base_tbl, lj_src_tbl, 1=>pick(lj_base_tbl,:k1,:k2)) --> innerjoin(lj_base_tbl, lj_src_tbl, 1)
    @fact innerjoin(lj_base_tbl, lj_src_tbl, 1=>1) --> innerjoin(lj_base_tbl, lj_src_tbl, 1)
    @fact leftjoin(@larr(reshape(1:10,2,5),axis2[:a,:b,:c,:c,:d]), @larr(axis1[:a,:b,:c],[1,2,3]), 1=>2) --> @larr(x1=reshape(1:10,2,5), x2=[1 2 3 3 NA;1 2 3 3 NA], axis2[:a,:b,:c,:c,:d])
    @fact innerjoin(@larr("x"=>reshape(1:10,2,5),axis2[:a,:b,:c,:c,:d]), @larr(axis1[:a,:b,:c],[1,2,3]), 1=>2) --> larr("x"=>[1 3 5 7;2 4 6 8], x1=[1 2 3 3;1 2 3 3], axis2=[:a,:b,:c,:c])
    @fact leftjoin(larr(k=[:x :x :y;:z :u :v], axis1=[:x,:u], axis2=darr(r=[:x, :y, :z])), larr(axis1=darr(k=[:x,:y,:z,:m,:n,:p]), b=[1,2,3,4,5,6]), 1) --> @larr(k=[:x :x :y;:z :u :v], b=[1 1 2;3 NA NA], axis1[:x,:u], axis2[r=[:x,:y,:z]])
    @fact leftjoin(larr(k=[:x :x :y;:z :u :v], axis1=[:x,:u], axis2=darr(r=[:x, :y, :z])), larr(axis1=darr(k=[:x,:y,:z,:m,:n,:p]), b=[1,2,3,4,5,6]), 1=>1) --> @larr(k=[:x :x :y;:z :u :v], b=[1 1 1;NA NA NA], axis1[:x,:u], axis2[r=[:x,:y,:z]])
    @fact leftjoin(larr(k=[:x :x :y;:z :u :v], axis1=[:x,:u], axis2=darr(r=[:x, :y, :z])), larr(axis1=darr(k=[:x,:y,:z,:m,:n,:p]), b=[1,2,3,4,5,6]), 1=>Any[nalift([:o :x :x;:q :r :y])]) --> @larr(k=[:x :x :y;:z :u :v], b=[NA 1 1;NA NA 2], axis1[:x,:u], axis2[r=[:x,:y,:z]])
    @fact innerjoin(larr(k=[:x :x :y;:z :u :v], axis1=[:x,:u], axis2=darr(r=[:x, :y, :z])), larr(axis1=darr(k=[:x,:y,:z,:m,:n,:p]), b=[1,2,3,4,5,6]), 1) --> @larr(k=[:x :x :y;:z :u :v], b=[1 1 2;3 NA NA], axis1[:x,:u], axis2[r=[:x,:y,:z]])
    @fact innerjoin(larr(k=[:x :x :y;:z :u :v], axis1=[:x,:u], axis2=darr(r=[:x, :y, :z])), larr(axis1=darr(k=[:x,:y,:z,:m,:n,:p]), b=[1,2,3,4,5,6]), 1=>1) --> @larr(k=[:x :x :y], b=[1 1 1], axis1[[:x]], axis2[r=[:x,:y,:z]])
    @fact innerjoin(larr(k=[:x :x :y;:z :u :v], axis1=[:x,:u], axis2=darr(r=[:x, :y, :z])), larr(axis1=darr(k=[:x,:y,:z,:m,:n,:p]), b=[1,2,3,4,5,6]), 1=>Any[nalift([:o :x :x;:q :r :y])]) --> @larr(k=[:x :y;:u :v], b=[1 1;NA 2], axis1[:x,:u], axis2[r=[:y,:z]])
  end
  context("darr version tests") do
    lj_base_tbl = darr(a=reshape(collect(1:20), 5,4), b=reshape(collect(101:120),5,4), k1=repmat([1,2,3,4,5],1,4), k2=repmat(['a','a','b','c','d'],1,4), r1=repmat([10 12 14 16],5,1))
    lj_src_tbl = @larr(c=reshape(collect(51:56),3,2), axis1[k1=[5,3,1], k2=['d','a','a']], axis2[t1=["x","y"]])
    @fact leftjoin(lj_base_tbl, lj_src_tbl, 1) --> @larr(a=reshape(repmat(collect(1:20),2),5,4,2),
                                                         b=reshape(repmat(collect(101:120),2),5,4,2),
                                                         k1=reshape(repmat(1:5,8),5,4,2),
                                                         k2=reshape(repmat(['a','a','b','c','d'],8),5,4,2),
                                                         r1=reshape(repmat(reshape(repmat([10 12 14 16],5,1),20),2),5,4,2),
                                                         c=reshape(@nalift([53,NA,NA,NA,51,53,NA,NA,NA,51,53,NA,NA,NA,51,53,NA,NA,NA,51,56,NA,NA,NA,54,56,NA,NA,NA,54,56,NA,NA,NA,54,56,NA,NA,NA,54]),5,4,2),
                                                         axis3[t1=["x","y"]])
    @fact innerjoin(lj_base_tbl, lj_src_tbl, 1) --> @larr(a=reshape(repmat([1,5,6,10,11,15,16,20],2),2,4,2),
                                                          b=reshape(repmat([101,105,106,110,111,115,116,120],2),2,4,2),
                                                          k1=reshape(repmat([1,5],8),2,4,2),
                                                          k2=reshape(repmat(['a','d'],8),2,4,2),
                                                          r1=reshape(repmat([10,10,12,12,14,14,16,16],2),2,4,2),
                                                          c=reshape([repmat([53,51],4);repmat([56,54],4)],2,4,2),
                                                          axis3[t1=["x","y"]])
    @fact leftjoin(larr(a=[1,2,3]), larr(b=[2 3;0 1],axis1=[:x,:y])) --> larr(a=reshape(repmat([1,2,3],4),3,2,2),b=reshape([2,2,2,0,0,0,3,3,3,1,1,1],3,2,2),axis2=[:x,:y])
    @fact innerjoin(larr(a=[1,2,3]), larr(b=[2 3;0 1],axis1=[:x,:y])) --> larr(a=reshape(repmat([1,2,3],4),3,2,2),b=reshape([2,2,2,0,0,0,3,3,3,1,1,1],3,2,2),axis2=[:x,:y])
    @fact leftjoin(lj_base_tbl, lj_src_tbl, 1=>pick(lj_base_tbl,:k1,:k2)) --> leftjoin(lj_base_tbl, lj_src_tbl, 1=>pick(lj_base_tbl,[:k1,:k2]))
    @fact leftjoin(lj_base_tbl, lj_src_tbl, 1=>pick(lj_base_tbl,:k1,:k2)) --> leftjoin(lj_base_tbl, lj_src_tbl, 1)
    @fact innerjoin(lj_base_tbl, lj_src_tbl, 1=>pick(lj_base_tbl,:k1,:k2)) --> innerjoin(lj_base_tbl, lj_src_tbl, 1=>pick(lj_base_tbl,[:k1,:k2]))
    @fact innerjoin(lj_base_tbl, lj_src_tbl, 1=>pick(lj_base_tbl,:k1,:k2)) --> innerjoin(lj_base_tbl, lj_src_tbl, 1)
    @fact leftjoin(larr(a=[1,2,3],axis=[:a,:b,:c]),larr(b=[4,5,6],axis=[:x,:y,:z])) --> larr(a=[1 1 1;2 2 2;3 3 3],b=[4 5 6;4 5 6;4 5 6],axis=[:a,:b,:c],axis=[:x,:y,:z])
    @fact leftjoin(larr(a=[1,2,3]),larr(b=[4,5,6],axis=[:x,:y,:z])) --> larr(a=[1 1 1;2 2 2;3 3 3],b=[4 5 6;4 5 6;4 5 6],axis2=[:x,:y,:z])
    @fact leftjoin(darr(a=[1,2,3]),larr(b=[4,5,6],axis=[:x,:y,:z])) --> larr(a=[1 1 1;2 2 2;3 3 3],b=[4 5 6;4 5 6;4 5 6],axis2=[:x,:y,:z])
    @fact innerjoin(larr(a=[1,2,3],axis=[:a,:b,:c]),larr(b=[4,5,6],axis=[:x,:y,:z])) --> larr(a=[1 1 1;2 2 2;3 3 3],b=[4 5 6;4 5 6;4 5 6],axis=[:a,:b,:c],axis=[:x,:y,:z])
    @fact innerjoin(larr(a=[1,2,3]),larr(b=[4,5,6],axis=[:x,:y,:z])) --> larr(a=[1 1 1;2 2 2;3 3 3],b=[4 5 6;4 5 6;4 5 6],axis2=[:x,:y,:z])
    @fact innerjoin(darr(a=[1,2,3]),larr(b=[4,5,6],axis=[:x,:y,:z])) --> larr(a=[1 1 1;2 2 2;3 3 3],b=[4 5 6;4 5 6;4 5 6],axis2=[:x,:y,:z])
    @fact leftjoin(larr(axis=['X','Y'],s1=[1 2 1;2 1 2],s2=[:a :a :a;:b :b :b],s3=[:c :d :e;:c :e :e]), @larr(axis[s2=[:a,:b]],axis[s3=[:c,:d]],r=[1.0 2.0;3.0 4.0]), 2,1) --> larr(axis1=['X','Y'],s1=[1 2 1;2 1 2],s2=[:a :a :a;:b :b :b],s3=[:c :d :e;:c :e :e],r=[1.0 2.0 NaN;3.0 NaN NaN])
    @fact innerjoin(larr(s1=[1 2 1;2 1 2],s2=[:a :a :a;:b :b :b],s3=[:c :d :e;:c :e :e]), @larr(axis[s2=[:a,:b]],axis[s3=[:c,:d]],r=[1.0 2.0;3.0 4.0]), 2,1) --> larr(s1=[1 2;2 1],s2=[:a :a;:b :b],s3=[:c :d;:c :e],r=[1.0 2.0;3.0 NaN])
  end
end

end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export affix

"""
    affix(;kwargs...)
    affix(children::Any;kwargs...)
    affix(children_maker::Function;kwargs...)

An Affix component.
Renders Affix at given element on mouse over or any other event. For more information, see: https://mantine.dev/core/affix/
Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): Any react node that should trigger tooltip
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `className` (String; optional): Often used with CSS to style elements with common properties
- `position` (optional): Fixed position in px. position has the following type: lists containing elements 'top', 'left', 'bottom', 'right'.
Those elements have the following types:
  - `top` (Real | String; optional)
  - `left` (Real | String; optional)
  - `bottom` (Real | String; optional)
  - `right` (Real | String; optional)
- `style` (Dict; optional): Inline style override
- `zIndex` (String; optional): Root element z-index property
"""
function affix(; kwargs...)
        available_props = Symbol[:children, :id, :className, :position, :style, :zIndex]
        wild_props = Symbol[]
        return Component("affix", "Affix", "dash_mantine_components", available_props, wild_props; kwargs...)
end

affix(children::Any; kwargs...) = affix(;kwargs..., children = children)
affix(children_maker::Function; kwargs...) = affix(children_maker(); kwargs...)

ext: jl


text:
for JT in [judiVector, judiWeights, judiWavefield]
    @eval function eval_op(a::$JT, b::$JT, op)
        c = deepcopy(a)
        for j=1:c.nsrc
            broadcast!(op, c.data[j], a.data[j], b.data[j])
        end
        return c
    end

    @eval function eval_op(a::T, b::$JT, op) where {T<:Number}
        c = deepcopy(b)
        for j=1:c.nsrc
            broadcast!(op, c.data[j], a, b.data[j])
        end
        return c
    end

    @eval function eval_op(a::$JT, b::T, op) where {T<:Number}
        c = deepcopy(a)
        for j=1:c.nsrc
            broadcast!(op, c.data[j], a.data[j], b)
        end
        return c
    end

    @eval function eval_op_ip(a::$JT, b::T, op) where {T<:Number}
        # Needed because julia 1.1 has a different definition of these two
        op == ldiv! && (return lmul!(1f0/a, b))
        op == rdiv! && (return rmul!(a, 1f0/b))
        #apply in place op
        for j=1:a.nsrc
            op(a.data[j], b)
        end
        a
    end

    @eval function eval_op_ip(a::T, b::$JT, op) where {T<:Number}
        # Needed because julia 1.1 has a different definition of these two
        op == ldiv! && (return lmul!(1f0/a, b))
        op == rdiv! && (return rmul!(a, 1f0/b))
        #apply in place op
        for j=1:b.nsrc
            op(a, b.data[j])
        end
        a
    end

    @eval function eval_op_ip(a::$JT, b::$JT, op) where {T<:Number}
        # Needed because julia 1.1 has a different definition of these two
        op == ldiv! && (return lmul!(1f0/a, b))
        op == rdiv! && (return rmul!(a, 1f0/b))
        #apply in place op
        for j=1:b.nsrc
            op(a.data[j], b.data[j])
        end
        a
    end
end

function matmulT(a::AbstractArray{T, 2}, b) where T
    return a*vec(vcat(b.data...))
end

function getattr(o, attr::Symbol, default=o)
    try
        return getfield(o, attr)
    catch e
        return default
    end
end

function getattri(o, attr::Symbol, ind::Integer, default=o)
    try
        return getproperty(o, attr)[ind]
    catch e
        return default
    end
end

tof32(x::Number) = [Float32(x)]
tof32(x::Array{T, N}) where {N, T<:Real} = T==Float32 ? x : Float32.(x)
tof32(x::Array{Array{T, N}, 1}) where {N, T<:Real} = T==Float32 ? x : tof32.(x)
tof32(x::Array{Any, 1}) = try Float32.(x) catch e tof32.(x) end
tof32(x::StepRangeLen) = tof32.(x)
tof32(x::Array{StepRangeLen}) = tof32.(x)

# Bypass mismatch in naming and fields
Base.getproperty(obj::judiWeights, sym::Symbol) = sym == :data ? getfield(obj, :weights) : getfield(obj, sym)
Base.getproperty(W::judiWavefield, sym::Symbol) = sym == :nsrc ? length(W.data) : getfield(W, sym)

# This whole part is basically taking adavantage of julia metaprogramming to define
# arithmetic operations on our types all at once. This makes sure that all is defined
# for all types properly.

for JT in [judiVector, judiWeights, judiWavefield]
    for opo=[:+, :-, :*, :/]
        @eval begin
            $opo(a::$JT, b::T) where {T<:Number} = eval_op(a, b, $opo)
            $opo(a::T, b::$JT) where {T<:Number} = eval_op(a, b, $opo)
            $opo(a::$JT, b::$JT) = eval_op(a, b, $opo)
        end
    end

    @eval -(a::$JT) = -1*a

    for ipop=[:lmul!, :rmul!, :rdiv!, :ldiv!]
        @eval begin
            $ipop(a::$JT, b::T) where {avDT, AT, T<:Number} = eval_op_ip(a, b, $ipop)
            $ipop(a::T, b::$JT) where {avDT, AT, T<:Number} = eval_op_ip(a, b, $ipop)
            $ipop(a::$JT, b::$JT) where {avDT, AT} = eval_op_ip(a, b, $ipop)
        end
    end

    @eval begin
        @doc  """
            vcat(a::Array{$($JT)})
        Concatenate an array of $($JT) into one
        """ vcat(a::Array{$JT, 1})

        @doc  """
            vcat(a::$($JT), b::$($JT))
            Concatenate a and b into a single $($JT)
        """ vcat(a::$JT, b::$JT)
    end

    @eval *(a::AbstractArray{T, 2}, b::$JT) where T = matmulT(a, b)
end

function vcat(a::Array{T, 1}) where T<:Union{judiVector, judiWeights, judiWavefield}
    return vcat(a...)
end

##### Rebuild bad vector

function rebuild_maybe_jld(x::Vector{Any})
    try
        return tof32(x)
    catch e
        if hasproperty(x[1], :offset)
            return [Float32.(StepRangeLen(xi.ref, xi.step, xi.len, xi.offset)) for xi in x]
        end
        return x
    end
end


text:
# Autogenerated wrapper script for HDF5_jll for x86_64-linux-gnu
export libhdf5, libhdf5_hl

using Zlib_jll
using OpenSSL_jll
using LibCURL_jll
JLLWrappers.@generate_wrapper_header("HDF5")
JLLWrappers.@declare_library_product(libhdf5, "libhdf5.so.200")
JLLWrappers.@declare_library_product(libhdf5_hl, "libhdf5_hl.so.200")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, OpenSSL_jll, LibCURL_jll)
    JLLWrappers.@init_library_product(
        libhdf5,
        "lib/libhdf5.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libhdf5_hl,
        "lib/libhdf5_hl.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


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

name = "CImGui"
version = v"1.74.0"

# Collection of sources required to build CImGui
sources = [
    "https://github.com/ocornut/imgui.git" =>
    "bdce8336364595d1a446957a6164c97363349a53",

    "https://github.com/cimgui/cimgui.git" =>
    "d9e1d9a80d621cd96d9900ac340092853100416f",

    "https://github.com/Gnimuc/CImGuiBuilder.git" =>
    "3a2c95f417d2769ff5fae9beb73baf27416859cc",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
mv imgui cimgui/
mv CImGuiBuilder/wrapper/helper.c cimgui/
mv CImGuiBuilder/wrapper/helper.h cimgui/
cd cimgui/
rm CMakeLists.txt
mv ../CImGuiBuilder/wrapper/CMakeLists.txt ./
mkdir build && cd build
cmake .. -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN}
make
make install
"""

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

# The products that we will ensure are always built
products = [
    LibraryProduct("libcimgui", :libcimgui)
    LibraryProduct("libcimgui_helper", :libcimgui_helper)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: Gnimuc/CImGuiBuilder

max_issues_repo_name: Gnimuc/CImGuiBuilder

max_forks_repo_name: Gnimuc/CImGuiBuilder


text:
module CO2SystemExample1

# This is an example of the use of CO2SYS. Have a look at the code

# Changes in v1.01:
# - fixed a silly bug that prevented the resulting figures to be drawn.
#
# Steven van Heuven. svheuven@gmail.com
#
# Converted to Julia by Matthew P. Humphreys [2019-02-28]

println(" ")
println("This is an example of the use of CO2SYS.m")
println("that uses its ability to process vectors of data.")
println(" ")
println("We will generate a figure that shows the sensitivity of pH and pCO2")
println(" to changes in DIC, while keeping everything else constant")
# println(" ")
# println("(Addional info: alk=2400, si=50, po4=2, dissociation constats: Mehrbach Refit)")
println(" ")

par1type =    1 # The first parameter supplied is of type "1", which is "alkalinity"
par1     = 2400 # value of the first parameter
par2type =    2 # The first parameter supplied is of type "1", which is "DIC"
par2     = 2100:5:2300 # value of the second parameter, which is a long vector of different DIC"s!
sal      =   35 # Salinity of the sample
tempin   =   10 # Temperature at input conditions
presin   =    0 # Pressure    at input conditions
tempout  =    0 # Temperature at output conditions - doesn't matter in this example
presout  =    0 # Pressure    at output conditions - doesn't matter in this example
sil      =   50 # Concentration of silicate  in the sample (in umol/kg)
po4      =    2 # Concentration of phosphate in the sample (in umol/kg)
pHscale  =    1 # pH scale at which the input pH is reported ("1" means "Total Scale")  - doesn't matter in this example
k1k2c    =    4 # Choice of H2CO3 and HCO3- dissociation constants K1 and K2 ("4" means "Mehrbach refit")
kso4c    =    1 # Choice of HSO4- dissociation constants KSO4 ("1" means "Dickson")

# Import CO2SYS
using CO2System

# Do the calculation. See CO2SYS's help for syntax and output format
A = CO2SYS(par1,par2,par1type,par2type,sal,tempin,tempout,presin,presout,
    sil,po4,pHscale,k1k2c,kso4c)[1]

# Draw the figure
using Plots
default(show=true)

# The calculated pCO2's are in the 4th column of the output A of CO2SYS
sp1 = plot(par2, A[:,4], color=:red, label="", xlabel="DIC", ylabel="pCO2 [uatm]",
    marker=(:circle, stroke(:red)))

# The calculated pH's are in the 3th column of the output A of CO2SYS
sp2 = plot(par2, A[:,3], color=:red, label="", xlabel="DIC", ylabel="pH",
    marker=(:circle, stroke(:red)))

plot(sp1,sp2)

println("DONE!")
println(" ")
println("See CO2SystemExample1.jl to see what the syntax for this calculation was.")
println(" ")

end # module CO2SystemExample1


text:
using MzXML, Unitful
using MzXML.IntervalSets
using MzCore
using Test

@testset "MzXML" begin
    scans, info = MzXML.load("test32.mzXML")
    @test info[:msModel] == "API 3000"
    @test eltype(scans) == MzXML.MSscan{Float32,Float32}
    @test length(scans) == 2
    scan = scans[1]
    @test scan.polarity == '-'
    @test scan.msLevel == 2
    @test scan.basePeakMz == 321.4
    @test scan.totIonCurrent == 1.637827475594e06
    @test scan.retentionTime == 0.004*u"s"
    @test scans[2].retentionTime == 0.809*u"s"

    # MzCore traits
    @test intensitytype(scans[1]) === Float32
    @test mztype(scans[1]) === Float32
    axmz, axt = limits(scans)
    @test minimum(axmz.val) === 110.9f0
    @test maximum(axmz.val) === 587.6f0
    @test minimum(axt.val) === 0.004*u"s"
    @test maximum(axt.val) === 0.809*u"s"

    scans, info = MzXML.load("test32.mzXML"; timeinterval=0.2u"s" .. Inf*u"s")
    @test length(scans) == 1 && scans[1].retentionTime == 0.809*u"s"

    io = IOBuffer()
    show(io, scan)
    @test occursin("└─basePeak 321.4", String(take!(io)))

    scans, info = MzXML.load("test32.mzXML"; timeinterval=10.2u"s" .. Inf*u"s", timeshift=10u"s")
    @test length(scans) == 1 && scans[1].retentionTime == 10.809*u"s"

    idx = MzXML.index("test32.mzXML")
    @test idx == [1231, 4242]

    scans, info = MzXML.load("test64.mzXML")
    @test eltype(scans) == MzXML.MSscan{Float64,Float64}
    @test intensitytype(scans[1]) === Float64
    @test mztype(scans[1]) === Float64
end


text:
println("This is a test on the subreferencing of Associative Array")

A = ReadCSV(joinpath(Base.source_dir(),"A.csv"))

# Get rows a and b.
println("Get rows a and b:")
A1r = A["a,b,",:];
printFull(A1r)

# Get rows containing a and columns 1 thru 3.
#TODO A2r = A("a *,",1:3)
# print(A2r)

# Get rows a to b.
println("Get rows a through b:")
A3r = A["a,:,b,",:]
printFull(A3r)

# Get rows starting with a or c.
println("Get rows starting with a or c:")
A4r = A[StartsWith("a,c,"),:]
printFull(A4r)

# Get cols a and b.
println("Get cols a and b:")
A1c = A[:,"a,b,"]
printFull(A1c)

# Get rows 1 thru 3 and cols containing a.
#TODO A2c = A(1:3,"a *,")
#print(A2c)

# Get cols a to b.
println("Get cols a through b:")
A3c = A[:,"a,:,b,"]
printFull(A3c)

# Get cols starting with a or c.
println("Get cols starting with a or c:")
A4c = A[:,StartsWith("a,c,")]
printFull(A4c)

# Get all values less than b.
#TODO A1v = (A < "b,")
#print(A1v)


text:
using Test

using Random
rng = MersenneTwister(1234)

using GR
GR.__init__()

mutable struct Example
  title::AbstractString
  code::Vector{Expr}
end

const _examples = Example[

Example("Simple line plot", [:(begin
        x = 0:π/100:2π
        y = sin.(x)
        plot(x, y)
    end)]),
Example("Scatter plot", [:(begin
        x = LinRange(0, 1, 51)
        y = x .- x.^2
        scatter(x, y)
    end)]),

Example("Colored scatter plot", [:(begin
        sz = LinRange(50, 300, length(x))
        c = LinRange(0, 255, length(x))
        scatter(x, y, sz, c)
    end)]),

Example("Stem plot", [:(begin
        stem(x, y)
    end)]),

Example("Histogram plot", [:(begin
        histogram(randn(rng, 10000))
    end)]),

Example("Multi-line plot", [:(begin
        plot(randn(rng, 50))
    end)]),

Example("Overlay plot", [:(begin
        oplot(randn(rng, 50, 3))
    end)]),

Example("3-d line plot", [:(begin
        x = LinRange(0, 30, 1000)
        y = cos.(x) .* x
        z = sin.(x) .* x
        plot3(x, y, z)
    end)]),

Example("Polar plot", [:(begin
        angles = LinRange(0, 2pi, 40)
        radii = LinRange(0, 2, 40)
        polar(angles, radii)
    end)]),

Example("3-d point plot", [:(begin
        x = 2 .* rand(rng, 100) .- 1
        y = 2 .* rand(rng, 100) .- 1
        z = 2 .* rand(rng, 100) .- 1
        scatter3(x, y, z)
    end)]),

Example("Colored 3-d point plot", [:(begin
        c = 999 .* rand(rng, 100) .+ 1
        scatter3(x, y, z, c)
    end)]),

Example("Hexbin plot", [:(begin
        x = randn(rng, 100000)
        y = randn(rng, 100000)
        hexbin(x, y)
    end)]),

Example("Contour plot", [:(begin
        x = 8 .* rand(rng, 100) .- 4
        y = 8 .* rand(rng, 100) .- 4
        z = sin.(x) .+ cos.(y)
        contour(x, y, z)
    end)]),

Example("Contour plot of matrix", [:(begin
        X = LinRange(-2, 2, 40)
        Y = LinRange(0, pi, 20)
        x, y = meshgrid(X, Y)
        z = sin.(x) .+ cos.(y)
        contour(x, y, z)
    end)]),

Example("Filled contour plot", [:(begin
        x = 8 .* rand(rng, 100) .- 4
        y = 8 .* rand(rng, 100) .- 4
        z = sin.(x) .+ cos.(y)
        contourf(x, y, z)
    end)]),

Example("Filled contour plot of matrix", [:(begin
        X = LinRange(-2, 2, 40)
        Y = LinRange(0, pi, 20)
        x, y = meshgrid(X, Y)
        z = sin.(x) .+ cos.(y)
        contourf(x, y, z)
    end)]),

Example("Filled contour plot on a triangular mesh", [:(begin
        x = 8 .* rand(rng, 100) .- 4
        y = 8 .* rand(rng, 100) .- 4
        z = sin.(x) + cos.(y)
        tricont(x, y, z)
    end)]),

Example("Surface plot", [:(begin
        x = 8 .* rand(rng, 100) .- 4
        y = 8 .* rand(rng, 100) .- 4
        z = sin.(x) .+ cos.(y)
        surface(x, y, z)
    end)]),

Example("Surface plot of matrix", [:(begin
        X = LinRange(-2, 2, 40)
        Y = LinRange(0, pi, 20)
        x, y = meshgrid(X, Y)
        z = sin.(x) .+ cos.(y)
        surface(x, y, z)
    end)]),

Example("Surface plot on a triangular mesh", [:(begin
        x = 8 .* rand(rng, 100) .- 4
        y = 8 .* rand(rng, 100) .- 4
        z = sin.(x) .+ cos.(y)
        trisurf(x, y, z)
    end)]),

Example("Wireframe plot", [:(begin
        x = 8 .* rand(rng, 100) .- 4
        y = 8 .* rand(rng, 100) .- 4
        z = sin.(x) .+ cos.(y)
        wireframe(x, y, z)
    end)]),

Example("Wireframe plot of matrix", [:(begin
        X = LinRange(-2, 2, 40)
        Y = LinRange(0, pi, 20)
        x, y = meshgrid(X, Y)
        z = sin.(x) .+ cos.(y)
        wireframe(x, y, z)
    end)]),

Example("Heatmap plot", [:(begin
        X = LinRange(-2, 2, 40)
        Y = LinRange(0, pi, 20)
        x, y = meshgrid(X, Y)
        z = sin.(x) .+ cos.(y)
        heatmap(z)
    end)]),

Example("Image plot", [:(begin
        imshow(z)
    end)]),

Example("Isosurface plot", [:(begin
        s = LinRange(-1, 1, 40)
        x, y, z = meshgrid(s, s, s)
        v = 1 .- (x .^ 2 .+ y .^ 2 .+ z .^ 2) .^ 0.5
        isosurface(v, isovalue=0.2)
    end)]),

Example("Shade points", [:(begin
        GR.GR3.terminate()

        volume(randn(rng, 50, 50, 50))

        N = 1_000_000
        x = randn(rng, N)
        y = randn(rng, N)
        shade(x, y)
    end)]),

Example("Discrete plot", [:(begin
        clearws()
        xd = -2 .+ 4 * rand(rng, 100)
        yd = -2 .+ 4 * rand(rng, 100)
        zd = [xd[i] * exp(-xd[i]^2 - yd[i]^2) for i = 1:100]

        setviewport(0.1, 0.95, 0.1, 0.95)
        setwindow(-2, 2, -2, 2)
        setspace(-0.5, 0.5, 0, 90)
        setcolormap(0)
        setlinecolorind(1)
        setmarkersize(1)
        setmarkertype(-1)
        setmarkercolorind(1)
        setcharheight(0.024)
        settextalign(2, 0)
        settextfontprec(3, 0)

        x, y, z = gridit(xd, yd, zd, 200, 200)
        h = -0.5:0.05:0.5
        surface(x, y, z, 5)
        contour(x, y, h, z, 0)
        polymarker(xd, yd)
        GR.axes(0.25, 0.25, -2, -2, 2, 2, 0.01)

        updatews()
        emergencyclosegks()
    end)])
]

function basic_tests()
    @test GR.tick(1.2,3.14) == 0.5

    for ex in _examples
        @info("Testing plot: $(ex.title)")

        GR.inline("pdf")
        file_path = ENV["GKS_FILEPATH"]

        map(eval, ex.code)

        @test isfile(file_path)
        rm(file_path)
    end
end

@time basic_tests()


text:
abstract type ETMethod end
struct PriestleyTaylor <: ETMethod end
struct PenmanMonteith <: ETMethod end

"""
    potential_ET(::PriestleyTaylor, Tair, pressure, Rn; G=0.0, S=0.0, ...)
    potential_ET(::PriestleyTaylor, Tair, pressure, Rn, G, S; ...)

    potential_ET(::PenmanMonteith, Tair, pressure, Rn, VPD, Ga_h;
      G=zero(Tair),S=zero(Tair), ...)
    fpotential_ET(::PenmanMonteith, Tair, pressure, Rn, VPD, Ga_h, G, S; ...)

    potential_ET!(df, approach; ...)

Potential evapotranspiration according to Priestley & Taylor 1972 or
the Penman-Monteith equation with a prescribed surface conductance.

# Arguments
- `Tair`:      Air temperature (degC)
- `pressure`:  Atmospheric pressure (kPa)
- `Rn`:        Net radiation (W m-2)
- `VPD`:       Vapor pressure deficit (kPa)
- `Ga`:        Aerodynamic conductance to heat/water vapor (m s-1)
- `df`:        DataFrame with the above variables
- `approach`:  Approach used: Either `PriestleyTaylor()`  or `PenmanMonteith()`.
optional:
- `G=0.0`:         Ground heat flux (W m-2). Defaults to zero.
- `S=0.0`:         Sum of all storage fluxes (W m-2) . Defaults to zero.
- `constants=`[`BigleafConstants`](@ref)`()`: physical constants (cp, eps, Rd, Rgas)
for `PriestleyTaylor`:
- `alpha = 1.26`:   Priestley-Taylor coefficient
for PenmanMonteith:
- `Gs_pot = 0.6`:    Potential/maximum surface conductance (mol m-2 s-1);
- `Esat_formula`: formula used in [`Esat_from_Tair`](@ref)

# Details
Potential evapotranspiration is calculated according to Priestley & Taylor, 1972
(`approach = PriestleyTaylor()`:

``\\mathit{LE}_{pot} = (\\alpha \\, \\Delta \\, (Rn - G - S)) / (\\Delta + \\gamma)``

``\\alpha`` is the Priestley-Taylor coefficient, ``\\Delta`` is the slope
of the saturation vapor pressure curve (kPa K-1), and ``\\gamma`` is the
psychrometric constant (kPa K-1).

If `approach = PenmanMonteith()`, potential evapotranspiration is calculated according
to the Penman-Monteith equation:

``\\mathit{LE}_{pot} = (\\Delta \\, (R_n - G - S) + \\rho \\, c_p \\, \\mathit{VPD} \\; G_a) /
(\\Delta + \\gamma \\, (1 + G_a/G_{s pot})``

where ``\\Delta`` is the slope of the saturation vapor pressure curve (kPa K-1),
``\\rho`` is the air density (kg m-3),
and ``\\gamma`` is the psychrometric constant (kPa K-1).
The value of ``G_{s pot}`` is typically a maximum value of ``G_s`` observed at the site,
e.g. the ``90^{th}`` percentile of ``G_s`` within the growing season.

Ground heat flux and storage heat flux `G` or `S` are provided as optional
arguments. In the input-explicit variants, they default to zero.
In the data-frame arguments, they default to missing, which results
in assuming them to be zero which is displayed in a log-message.
Note that in difference ot the bigleaf R package, you explitly need to
care for missing values (see examples).

# Value
NamedTuple with the following entries:
- `ET_pot`: Potential evapotranspiration (kg m-2 s-1)
- `LE_pot`: Potential latent heat flux (W m-2)

# References
- Priestley, C_H_B., Taylor, R_J., 1972: On the assessment of surface heat flux
  and evaporation using large-scale parameters. Monthly Weather Review 100, 81-92.
- Allen, R_G., Pereira L_S., Raes D., Smith M., 1998: Crop evapotranspiration -
  Guidelines for computing crop water requirements - FAO Irrigation and drainage
  paper 56.
- Novick, K_A., et al. 2016: The increasing importance of atmospheric demand
  for ecosystem water and carbon fluxes. Nature Climate Change 6, 1023 - 1027.

# See also
[`surface_conductance`](@ref)

# Examples
Calculate potential ET of a surface that receives a net radiation of 500 Wm-2
using Priestley-Taylor:
```jldoctest; output=false
Tair, pressure, Rn = 30.0,100.0,500.0
ET_pot, LE_pot = potential_ET(PriestleyTaylor(), Tair, pressure, Rn)
≈(ET_pot, 0.000204; rtol = 1e-2)
# output
true
```

Calculate potential ET for a surface with known Gs (0.5 mol m-2 s-1) and Ga (0.1 m s-1)
using Penman-Monteith:
```jldoctest; output=false
Tair, pressure, Rn = 30.0,100.0,500.0
VPD, Ga_h, Gs_pot = 2.0, 0.1, 0.5
ET_pot, LE_pot = potential_ET(
  PenmanMonteith(), Tair,pressure,Rn,VPD, Ga_h; Gs_pot)
# now cross-check with the inverted equation
Gs_ms, Gs_mol = surface_conductance(
  InversePenmanMonteith(), Tair,pressure,VPD,LE_pot,Rn,Ga_h)
Gs_mol ≈ Gs_pot
# output
true
```

DataFrame variant with explicitly replacing missings using `coalesce.`:
```jldoctest; output=false
using DataFrames
df = DataFrame(
  Tair = 20.0:1.0:30.0,pressure = 100.0, Rn = 500.0, G = 105.0, VPD = 2.0,
  Ga_h = 0.1)
allowmissing!(df, Cols(:G)); df.G[1] = missing
#
# need to provide G explicitly
df_ET = potential_ET!(copy(df), PriestleyTaylor(); G = df.G)
ismissing(df_ET.ET_pot[1])
#
# use coalesce to replace missing values by zero
df_ET = potential_ET!(
  copy(df), PriestleyTaylor(); G = coalesce.(df.G, zero(df.G)))
!ismissing(df_ET.ET_pot[1])
# output
true
```
"""
function potential_ET(approach::PriestleyTaylor, Tair, pressure, Rn;
  G=zero(Tair),S=zero(Tair), kwargs...)
  #
  potential_ET(approach, Tair, pressure, Rn, G, S; kwargs...)
end
function potential_ET(::PriestleyTaylor, Tair, pressure, Rn, G, S;
  alpha=1.26,
  Esat_formula=Sonntag1990(),
  constants=BigleafConstants())
  #
  gamma  = psychrometric_constant(Tair,pressure;constants)
  Delta  = Esat_from_Tair_deriv(Tair; Esat_formula, constants)
  LE_pot = (alpha * Delta * (Rn - G - S)) / (Delta + gamma)
  ET_pot = LE_to_ET(LE_pot,Tair)
  (ET_pot = ET_pot, LE_pot = LE_pot)
end
function potential_ET!(df::AbstractDataFrame, approach::PriestleyTaylor;
  G=0.0,S=0.0, kwargs...)
  # G or S might be vectors, hence cannot use ByRow
  ft = (args...) -> potential_ET.(Ref(approach), args...,G,S; kwargs...)
  transform!(df, Cols(:Tair, :pressure, :Rn) => ft => AsTable)
end

function potential_ET(approach::PenmanMonteith, Tair, pressure, Rn, VPD, Ga_h;
  G=zero(Tair),S=zero(Tair), kwargs...)
  #
  potential_ET(approach, Tair, pressure, Rn, VPD, Ga_h, G, S; kwargs...)
end
function potential_ET(::PenmanMonteith, Tair, pressure, Rn, VPD, Ga_h, G, S;
  Gs_pot=0.6,
  Esat_formula=Sonntag1990(),
  constants=BigleafConstants()
  )
  gamma  = psychrometric_constant(Tair,pressure;constants)
  Delta  = Esat_from_Tair_deriv(Tair; Esat_formula = Esat_formula,constants)
  Gs_pot = mol_to_ms(Gs_pot,Tair,pressure;constants)
  rho    = air_density(Tair,pressure;constants)
  LE_pot = (Delta * (Rn - G - S) + rho * oftype(VPD,constants.cp) * VPD * Ga_h) /
    (Delta + gamma * (1 + Ga_h / Gs_pot))
  ET_pot = LE_to_ET(LE_pot,Tair)
  (ET_pot = ET_pot, LE_pot = LE_pot)
end
function potential_ET!(df::AbstractDataFrame, approach::PenmanMonteith;
  G=0.0,S=0.0, kwargs...
  )
  # G or S might be vectors, hence cannot use ByRow
  ft = (args...) -> potential_ET.(Ref(approach), args..., G, S; kwargs...)
  transform!(df, Cols(:Tair, :pressure, :Rn, :VPD, :Ga_h) => ft => AsTable)
end

# function potential_ET(df, approach::PriestleyTaylor;
#   G=missing,S=missing, infoGS=true, kwargs...)
#   #
#   dfGS = get_df_GS(df, G,S; infoGS)
#   f(args...) = potential_ET(args..., approach; kwargs...)
#   select(hcat(select(df,:Tair, :pressure, :Rn), dfGS; copycols = false),
#     All() => ByRow(f) => AsTable
#   )
# end
# function potential_ET(df, approach::PenmanMonteith;
#   G=missing,S=missing, infoGS=true, kwargs...)
#   #
#   dfGS = get_df_GS(df, G,S; infoGS)
#   function f(args...)
#     potential_ET(args..., approach; kwargs...)
#   end
#   select(hcat(select(df,:Tair, :pressure, :Rn, :VPD, :Ga), dfGS; copycols = false),
#     All() => ByRow(f) => AsTable
#   )
# end

# function get_df_GS(df, G, S; infoGS=true)
#   nout = nrow(df)
#   G_ = if ismissing(G)
#     infoGS && @info("Ground heat flux G is not provided and set to 0.")
#     Zeros(nout)
#   else
#     G
#   end
#   S_ = if ismissing(S)
#     infoGS && @info("Storage heat flux S is not provided and set to 0.")
#     Zeros(nout)
#   else
#     S
#   end
#   (G  = G_, S = S_)
# end

# function fill_GS_missings(df::DataFrame, G,S, missing_G_as_NA, missing_S_as_NA; infoGS=true)
#   nout = nrow(df)
#   G_ = ifelse(
#     ismissing(G),
#     (infoGS && @info("Ground heat flux G is not provided and set to 0."); Zeros(nout)),
#     @. ifelse(missing_G_as_NA, G, coalesce(G, 0.0))
#     )
#   S_ = ifelse(
#   ismissing(S),
#   (infoGS && @info("Storage heat flux S is not provided and set to 0."); Zeros(nout)),
#   @. ifelse(missing_S_as_NA, G, coalesce(S, 0.0))
#   )
#   DataFrame(G  = G_, S = S_)
# end

# function fill_vec(G; is_replace_missing = true, fillvalue = zero(G))
#   @. ifelse(is_replace_missing, coalesce(G, fillvalue), G)
# end

"""
    equilibrium_imposed_ET(Tair,pressure,VPD,Gs, Rn; ...)
    equilibrium_imposed_ET!(df; ...)

Evapotranspiration (ET) split up into imposed ET and equilibrium ET.

# Argumens
- `Tair`      : Air temperature (deg C)
- `pressure`  : Atmospheric pressure (kPa)
- `VPD`       : Air vapor pressure deficit (kPa)
- `Gs`        : surface conductance to water vapor (m s-1)
- `Rn`        : Net radiation (W m-2)
optional :
- `G=0`       : Ground heat flux (W m-2)
- `S=0`       : Sum of all storage fluxes (W m-2)
- `Esat_formula=Sonntag1990()`: formula used in [`Esat_from_Tair`](@ref)
- `constants=`[`BigleafConstants`](@ref)`()`: pysical constants (cp, eps)

# Details
Total evapotranspiration can be written in the form (Jarvis & McNaughton 6):

``ET = \\Omega \\mathit{ET}_{eq} + (1 - \\Omega) \\mathit{ET}_{imp}``

where ``\\Omega`` is the decoupling coefficient as calculated from
[`decoupling`](@ref). `ET_eq` is the equilibrium evapotranspiration i.e.,
the ET rate that would occur under uncoupled conditions, where the budget
is dominated by radiation (when Ga -> 0):

``ET_{eq} = (\\Delta \\, (R_n - G - S) \\, \\lambda) / ( \\Delta \\gamma)``

where ``\\Delta`` is the slope of the saturation vapor pressur(kPa K-1),
``\\lambda`` is the latent heat of vaporization (J kg-1), and ``\\gamma``
is the psychrometric constant (kPa K-1).
`ET_imp` is the imposed evapotranspiration rate, the ET rate
that would occur under fully coupled conditions (when Ga -> inf):

``ET_{imp} = (\\rho \\, c_p \\, \\mathit{VPD} ~ G_s \\, \\lambda) / \\gamma``

where ``\\rho`` is the air density (kg m-3).

# Value
A `NamedTuple` or `DataFrame` with the following columns:
- `ET_eq`: Equilibrium ET (kg m-2 s-1)
- `ET_imp`: Imposed ET (kg m-2 s-1)
- `LE_eq`: Equilibrium LE (W m-2)
- `LE_imp`: Imposed LE (W m-2)

# References
- Jarvis, P_G., McNaughton, K_G., 1986: Stomatal control of transpiration:
  scaling up from leaf to region. Advances in Ecological Rese1-49.
- Monteith, J_L., Unsworth, M_H., 2008: Principles of ironmPhysics.
  3rd edition. Academic Press, London.

# Examples
```jldoctest; output = false
Tair,pressure,Rn, VPD, Gs = 20.0,100.0,50.0, 0.5, 0.01
ET_eq, ET_imp, LE_eq, LE_imp = equilibrium_imposed_ET(Tair,pressure,VPD,Gs, Rn)
≈(ET_eq, 1.399424e-05; rtol = 1e-5)
# output
true
```
"""
function equilibrium_imposed_ET(Tair,pressure,VPD,Gs, Rn;
  G=zero(Tair),S=zero(Tair), kwargs...)
  equilibrium_imposed_ET(Tair,pressure,VPD,Gs, Rn, G, S; kwargs...)
end

function equilibrium_imposed_ET(Tair,pressure,VPD,Gs, Rn, G, S;
  Esat_formula=Sonntag1990(),
  constants=BigleafConstants())
  #
  rho    = air_density(Tair, pressure; constants)
  gamma  = psychrometric_constant(Tair, pressure; constants)
  Delta  = Esat_from_Tair_deriv(Tair; Esat_formula = Esat_formula, constants)
  LE_eq  = (Delta * (Rn - G - S)) / (gamma + Delta)
  LE_imp = (rho * constants.cp * Gs * VPD) / gamma
  #
  ET_imp = LE_to_ET(LE_imp,Tair)
  ET_eq  = LE_to_ET(LE_eq,Tair)
  (;ET_eq, ET_imp, LE_eq, LE_imp)
end
function equilibrium_imposed_ET!(df; G=0.0, S=0.0, kwargs...)
  f(args...) = equilibrium_imposed_ET.(args..., G, S; kwargs...)
  transform!(df, [:Tair, :pressure, :VPD, :Gs, :Rn] => f => AsTable)
end

"""
    TODO; implement decoupling.

This stub is there to satisfy links im Help-pages.
"""
function decoupling()
  error("not yet implemented.")
end


text:
"""
Start FlowAtlas using:
```julia
julia> FlowAtlas.run()
```
"""
module FlowAtlas
include_dependency("../Project.toml")

using JSServe: Session, @js_str, file_server, response_404, Routes, App, DOM, Observable, on
using JSServe, HTTP, JSON, Images, FileIO, ImageIO
using GigaScatter, ColorSchemes

using Serialization: serialize,deserialize
using FlowWorkspace: inpolygon
using Base: NamedTuple

using FlowWorkspace, StaticArrays, DataFrames, MetaGraphs, OrderedCollections
using GigaSOM, StatsBase, Impute

include("lib/colors.jl")
include("lib/selection.jl")
include("lib/embed.jl")

include("lib/tile.jl")
include("lib/gate.jl")
include("lib/count.jl")

include("lib/sys.jl")
include("lib/map.jl")
include("lib/sidebar.jl")

############################################## javascript dependencies
const ol = JSServe.Dependency( :ol, [# OpenLayers
    joinpath(@__DIR__,"assets/ol/ol.js"), joinpath(@__DIR__,"assets/ol/ol.css")
])

const d3 = JSServe.Dependency( :d3, [ # data-driven documents
    joinpath(@__DIR__, "assets/d3/d3.v6.min.js"),
])

const FileSaver = JSServe.Dependency( :FileSaver, [ # saving files
    joinpath(@__DIR__, "assets/FileSaver.min.js")
])

const style = JSServe.Dependency( :style, [ # custom styling todo(@gszep) remove remote dep
    "//cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.2/css/all.min.css", joinpath(@__DIR__, "assets/style.css"),
])

const extensions = JSServe.Dependency( :extensions, map(  extension -> joinpath(@__DIR__,extension),
    ["assets/ol/sidebar.css","assets/ol/colorbar.js","assets/Sortable.js","assets/violins.js","assets/boxplots.js","assets/utils.js"]
))

function run( path::String; files::String=joinpath(dirname(path),"*.fcs"), transform::Function=x->asinh(x/250),
        port::Int = 3141, url::String = "http://localhost:$port", cols::Symbol=:union, drop::Union{Vector{String},Nothing}=String[],
        nlevels::Int=10, nbins::Int=50, p::Real=0.1, q::Real=5e-4, channelScheme=reverse(ColorSchemes.matter), labelScheme=ColorSchemes.seaborn_colorblind,
        perplexity::Number=300, maxIter::Integer=10000, hold::Bool=true )

    indexTransform(x::AbstractVector{<:Union{Number,Missing}}) = toIndex(x; nlevels=nlevels)
    channelRange(x::AbstractVector{<:Union{Number,Missing}}) = range(quantile(skipmissing(x),(p,1-p))...,length=nlevels)

    @info "Loading FCS files..."
    data, labels, groups, gating = FlowWorkspace.load( path; files=files, transform=transform, cols=cols)

    keep = map( graph -> map( idx -> MetaGraphs.get_prop(graph,idx,:name), filter(idx -> MetaGraphs.outdegree(graph,idx) == 0, 1:MetaGraphs.nv(graph))), values(gating) )
    select!( labels, isnothing(drop) ? union(keep...,["Unlabelled"]) : Not(drop) )

    names = (
        channels = Base.names(data), populations = Base.names(labels),
        conditions = filter(name -> ~occursin(r"\d", name), Base.names(groups)),
        groups = filter(name -> occursin(r"\d", name), Base.names(groups))
    )

    ###############################################################################
    ########################################################### calculate embedding
    som_path, _ = splitext(path)
    _, embedding = embed(data, path = som_path*".som", perplexity = perplexity, maxIter = maxIter)
    embedding = (
        coordinates = map(SVector, embedding[2,:], -embedding[1,:]),
        array = embedding, extrema = extrema(embedding,dims=2)
    )

    ###############################################################################
    ######################################################### colormaps and filters
    colors = (
        labels = (
            names = OrderedDict([ name=>'#'*hex(get(labelScheme,x)) for (name,x) ∈ zip(Base.names(labels),range(0,1,length=size(labels,2))) ]),
            rows = channelview(fill(parse(RGBA{Float64},"#663F46"), size(data,1))),

        ),
        channels = (
            range =  Dict([ col => channelRange(data[!,col]) for col ∈ Base.names(data) ]),
            rows = combine(data, [ col => indexTransform => col for col ∈ Base.names(data) ]),
            colors = channelview([map(x->RGBA(get(channelScheme, x)),range(0,1,length=nlevels)); RGBA(0,0,0,0)]),
            hex = map(x->'#'*hex(x),colorview(RGB,view(channelview(map(x->RGBA(get(channelScheme, x)),range(0,1,length=nlevels))),1:3,:)))

        )
    )

    flattened_data = vcat(map(col->data[!,col],Base.names(data))...)
    channelMin,channelMax = quantile(skipmissing(flattened_data),(q,1-q))
    for (name,color) ∈ colors.labels.names colors!( Dict(["name"=>name,"color"=>color]),labels,groups,colors) end

    ##################################################### initalise filter settings
    selections = (
        codes = select(hcat(labels,groups), AsTable(:) => ByRow(encode ∘ values) => "encoding")[:,:encoding],
        names = OrderedDict([ name => true for name ∈ [Base.names(labels);Base.names(groups)] ]),
        rows = fill(true, size(data, 1))
    )

    ########################################################################## body
    app = App() do session::Session
        Map = DOM.div(id = "map", class = "sidebar-map")

        JSServe.onload(session, Map, olMap( embedding.extrema, colors; port=port))
        return DOM.div(id = "application", extensions, style, DOM.title("FlowAtlas.jl"), Map,

            HTML("""<script src="$(JSServe.Asset(joinpath(@__DIR__,"assets/ol/sidebar.js")))" ></script>"""),
            sidebar(session, names, colors; port=port)
        )
    end

    ###############################################################################
    ############################################# open app as locally hosted server
    try
        global server = JSServe.Server(app, "127.0.0.1", port;
            verbose = true, routes = Routes( "/" => app,
                r"/assetserver/" * r"[\da-f]"^40 * r"-.*" => file_server,

                r"/\d+/\d+/\d+.png" => x->tile(x,embedding,selections,colors,names),
                r"/colors" => x->colors!(x,labels,groups,colors), r"/selection" => x->selection!(x,selections,names),

                r"/gate" => x->gate(x,data,embedding,selections;channelRange=range(channelMin,channelMax,length=nbins)),
                r"/count" => x->count(x,selections),

                r"/favicon.ico" => context -> HTTP.Response(500),
                r".*" => context -> response_404() )
        )

        browser(url)
        printstyled("[FlowAtlas] ", color = :blue)

        println("Running at $url ~ Happy exploring! ")
        hold && wait()

    catch exception
        if exception isa InterruptException
            close(server)

            printstyled("[FlowAtlas] ", color = :blue)
            println("Closing... have a lovely day! 🎈")

        else
            rethrow()
        end
    end
end

@info """\n
    Welcome to FlowAtlas 🎈
    Start FlowAtlas server using:
        julia> FlowAtlas.run("workspace.wsp"; files="*.fcs")
\n"""
end

max_stars_repo_name: gszep/flow-explorer

max_issues_repo_name: gszep/flow-explorer

max_forks_repo_name: gszep/flow-explorer


text:
# using Pkg; Pkg.add("https://github.com/QuantEcon/Games.jl")

#module beautifulmind

using Games

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

# defines the simultaneous game
function init()

if uppercase("n") == uppercase(query("\n*****Simultaneous games ****\n\nDefine own game?(y/n)\n\n\tchoice: "))
    mat = Array{Float64}(undef, 2, 2, 2)
    mat[1, 1, :] = [3, 2]
    mat[1, 2, :] = [2, 2]
    mat[2, 1, :] = [2, 3]
    mat[2, 2, :] = [5, 0]

    sim = NormalFormGame(mat)
    determineStrategy(sim)
    output(sim)
else

println("Defining own game\n")
    if 1 == parse(Int,query("choose input method\n\t1.by matrix\n\t2.by vector\n\n\tchoice: "))

        matrix = Array{Int}(undef, 2, 2, 2)
        for o in 1:2
            for i in 1:2

             str = query("Input values for row $o column $i e.g(2 5): ")
            matrix[o, i, :] =  [parse(Int,split(str," ")[1]), parse(Int,split(str," ")[2])]

            end
        end
        game = NormalFormGame(matrix)
        determineStrategy(game)
        output(game)
    else

         p1 = [1 2 ; 3 4]
          str = replace(query("Input player1 strategies e.g(2 5;3 7): "),";" => "")
          str = replace(str," " => "")

          for (i,x) in enumerate(str)
             p1[i]=parse(Int,x)
          end

        p1 = Player(p1)

          p2 = [1 2 ; 3 4]
          str = replace(query("Input player2 strategies e.g(2 5;3 7): "),";" => "")
          str = replace(str," " => "")

          for (i,x) in enumerate(str)
             p2[i]=parse(Int,x)
         end

        p2 = Player(p2)

        game = NormalFormGame(p1,p2)
        determineStrategy(game)
        output(game)
    end

  end
end

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

function determineStrategy(game)

    a = query("\n\nwhat is player1's first strategy e.g(a): ")
    b = query("\nwhat is player1's second strategy e.g(b): ")

    c = query("\nwhat is player2's first strategy e.g(c): ")
    d = query("\nwhat is player2's second strategy e.g(d): ")

    act11 = game.players[1].payoff_array[1] + game.players[1].payoff_array[3]

    act12 = game.players[1].payoff_array[2] +game.players[1].payoff_array[4]

    act21 = game.players[2].payoff_array[1] + game.players[2].payoff_array[2]
    act22 = game.players[2].payoff_array[3] + game.players[2].payoff_array[4]

println("\n")

    if act11 > act12
        println("player1  dominant strategy: $a\n         dominated stratergy: $b")
    end

    if act12 > act11
        println("player1  dominant strategy: $b\n         dominated stratergy: $a")
    elseif act11 == act12
        println("player1's strategies are equal")
    end

     if(act21>act22)
        println("player2  dominant strategy: $c\n         dominated stratergy: $d")
    elseif(act22>act21)
        println("player2  dominant strategy: $d\n         dominated stratergy: $c")
    else
        println("player2's strategies are equal")
    end

   # best_responses(game.players[1], game.players[1])

    return
end

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

function output(game::NormalFormGame)
    equilibria = pure_nash(game)
    equilibriaNo = length(equilibria)

    if equilibriaNo == 0
        str = "no pure Nash equilibrium found"
    elseif equilibriaNo == 1
        str = "There is 1 equilibrium at:\n$(equilibria[1])"
    else
        str = "There is $equilibriaNo  equilibria:\n"
        for (i, E) in enumerate(equilibria)
            i < equilibriaNo ? str *= "$E, " : str *= "$E"
        end
    end
    println(join(["\n",str]))

end

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

function query(str::String)

    print(str);
   return input = readline();

end

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

init() #driver statement
#end #module


text:
###############################################################################
#Package of functions related to the Theta functions.
#
#Citation:
#Haddou, M., & Maheux, P. (2014).
#Smoothing methods for nonlinear complementarity problems.
#Journal of Optimization Theory and Applications, 160(3), 711-729.
###############################################################################

"""
Definition of Theta^1

function theta(x :: Union{Float64,Vector},
               r :: Float64)
output : Union{Float64,Vector}
"""
function theta(x :: Union{Float64,Vector},
               r :: Float64)

 if r < 0.0
  throw("DomainError: Theta undefined for r < 0.")
 elseif r==0.0 && in(0.0,x)
  throw("DomainError: Theta undefined for x=0, r=0.")
 end

 y = (x .>= 0) .* x ./ (abs.(x) .+ r) + (x .< 0) .* (x ./ r-x.^2 ./ r^2)

 return y
end

"""
Inverse of function theta
function invtheta(x :: Union{Float64,Vector},
                  r :: Float64)
"""
function invtheta(x :: Union{Float64,Vector},
                  r :: Float64)

 if in(true, x .> 1.)
  throw("DomainError: Theta^{-1} undefined for x >= 1.")
 end

 y = (x .>= 0) .* -r .* x ./ (x .- 1) + (x .< 0) .* r/2 .* (1 .- sqrt.(1 .+ 4 .* abs.(x)))

 return y
end

"""
First derivative of theta
function dtheta(x :: Union{Float64,Vector},
                r :: Float64)
"""
function dtheta(x :: Union{Float64,Vector},
                r :: Float64)
 if r < 0.0
  throw("DomainError: Theta undefined for r < 0.")
 elseif r==0.0 && in(0.0,x)
  throw("DomainError: Theta undefined for x = 0, r = 0.")
 end

 y = (x .>= 0) .* (r ./ (abs.(x) .+ r).^2) + (x .< 0) .* (1/r .-2 .* x ./ r^2)

 return y
end

"""
Second derivative of theta
function ddtheta(x :: Union{Float64,Vector},
                 r :: Float64)
"""
function ddtheta(x :: Union{Float64,Vector},
                 r :: Float64)

 if r < 0.0
  throw("DomainError: Theta undefined for r < 0.")
 elseif r==0.0 && in(0.0,x)
  throw("DomainError: Theta undefined for x = 0, r = 0.")
 end

 y = (x .>= 0) .* (-2*r ./ (abs.(x) .+ r).^3) + (x .< 0) .* (-2/(r^2))

 return y
end
#Fun fact:
# ThetaFct.ddtheta([-1.0], 0.1) != ThetaFct.ddtheta(-1.0, 0.1)

"""
alpha_theta_max : Solve the equation in alpha where psi(x,r,s,t)=s+t*theta(x,r)
In : x  (position of x)
     dx (direction of x)
     y  (position of y)
     dy (direcion of y)
     r,s,t (parameters)

output: Array with two entries
"""
function alpha_theta_max(x  :: Float64,
                         dx :: Float64,
                         y  :: Float64,
                         dy :: Float64,
                         r  :: Float64,
                         s  :: Float64,
                         t  :: Float64)

 alphaf1 = _alpha_theta_max_one_leaf(x, dx, y, dy, r, s, t)
 alphaf2 = _alpha_theta_max_one_leaf(y, dy, x, dx, r, s, t)

 #le premier fixe yg et le deuxième yh
 return [alphaf1, alphaf2]
end
function alpha_theta_max(x  :: Vector,
                         dx :: Vector,
                         y  :: Vector,
                         dy :: Vector,
                         r  :: Float64,
                         s  :: Float64,
                         t  :: Float64)

 alphaf1 = _alpha_theta_max_one_leaf(x, dx, y, dy, r, s, t)
 alphaf2 = _alpha_theta_max_one_leaf(y, dy, x, dx, r, s, t)

 #le premier fixe yg et le deuxième yh
 alpha_yg = minimum(alphaf1[alphaf1 .>= 0.0])
 alpha_yh = minimum(alphaf2[alphaf2 .>= 0.0])
 return [alpha_yg, alpha_yh]
end

"""
function _alpha_theta_max_one_leaf(x  :: Union{Float64,Vector},
                                   dx :: Union{Float64,Vector},
                                   y  :: Union{Float64,Vector},
                                   dy :: Union{Float64,Vector},
                                   r  :: Float64,
                                   s  :: Float64,
                                   t  :: Float64)

Solve the equation (a variable) :

(x+a*dx-s)*(y+a*dy-s+r)-t*(y+a*dy-s)=0 -> alphaf1p

comments: (solve the equation)
r^2[x+a*dx-s-t*((y+a*dy-s)/r-(y+a*dy-s)^2/r^2)]=0 -> alphaf1n
"""
function _alpha_theta_max_one_leaf(x  :: Union{Float64,Vector},
                                   dx :: Union{Float64,Vector},
                                   y  :: Union{Float64,Vector},
                                   dy :: Union{Float64,Vector},
                                   r  :: Float64,
                                   s  :: Float64,
                                   t  :: Float64)

 if t!=0

  a = dx .* dy
  b = dy .* (x .- s .- t) .+ dx .* (y .- s .+ r)
  c = x .* y - x .* s + x .* r - y .* s .+ s^2 .- s*r-t .* y .+ t*s

  alphaf1p = _poly_degree_2_positif(a, b, c)

  tmp1 = (alphaf1p .>= 0) .& ((x+alphaf1p .* dx .>= s) .& (y+alphaf1p .* dy .>= s))

  #si on considère la contrainte intérieure:
  #a = dy.^2
  #b = r^2.*dx-r*t.*dy+2*t*y.*dy-2*t*s.*dy
  #c = r^2.*x .- r^2*s - r*t*y .+ r*t*s+t.*y.^2.+t*s^2-t*2*s.*y
  #alphaf1n = _poly_degree_2_positif(a, b, c)
  #tmp2 = (alphaf1p .>= 0) .& (x+alphaf1p.*dx .<= s .& y+alphaf1p.*dy .<= s)

  alphaf1 = (tmp1) .* alphaf1p + (!tmp1) .* Inf

 else #cas où t=0

  tmp = (dx .!= 0) .& (((s .- x) ./ dx .> 0) .| ((s .== x) .& (dx .> 0)))
  alphaf1 = (tmp) .* (s .- x)./dx + (!tmp) .* Inf

 end #fin de condition t!=0

 return alphaf1
end

"""
function _poly_degree_2_positif(a :: Union{Float64,Vector},
                                b :: Union{Float64,Vector},
                                c :: Union{Float64,Vector})

Compute the smallest positive root of a second order polynomial
or a negative/complex value if no positive real root.
If no solution or infinite nb of solutions, return Inf.
"""
function _poly_degree_2_positif(a :: Union{Float64,Vector},
                                b :: Union{Float64,Vector},
                                c :: Union{Float64,Vector})

 disc = b.^2 .- 4 .* a .* c

 if in(true, disc.<0.0)
  return -Inf.*a #throw("FunctionError: complex roots.")
 end

 sol1 = (-b .+ sqrt.(disc)) ./ (2 .* a)
 sol2 = (-b .- sqrt.(disc)) ./ (2 .* a)

 solp  = ((sol1 .>= 0.0) .& (sol2 .>= 0.0)) .* min.(sol1,sol2)
 #solp += xor.(sol1 .< 0.0, sol2 .< 0.0) .* max.(sol1,sol2)
 solp += ((sol1 .< 0.0) .| (sol2 .< 0.0)) .* max.(sol1,sol2)

 sol  = (a .!= 0.) .* solp
 sol += ((a .== 0.) .& (b .!= 0.)) .* (-c ./ b)
 sol += ((a .== 0.) .& (b .== 0.)) .* Inf

 return sol
end


text:
using POMDPs, Pkg
using Distributions
using Parameters
using POMDPModelTools
using Random
using PyPlot
pygui(:qt5)

struct SleepState3
    complete::Int32
    productivity::Int32
    deadline::Int32
    hours_slept::Int32
    hours_since_slept::Int32
end

@with_kw mutable struct SleepPOMDP3 <: POMDP{SleepState3, Symbol, Float64}
    r_complete::Int32 = 10000
    r_missed::Int32 = -100
    sleep_multiplier::Float64 = 1.5
    work_multiplier::Float64 = 0.8
    relax_multiplier::Float64 = 1.2
    max_deadline::Int32 = 72
    discount::Float64 = 1.0
    obs_sigma::Float64 = 10
end

function sleep_plot(history, path)
    complete = []
    productivity = []
    deadline = []
    observation = []
    actions = []
    i = 0
    for (s, a, o, r) in eachstep(history, "s,a,o,r")
        push!(complete, s.complete)
        push!(productivity, s.productivity)
        push!(deadline, s.deadline)
        push!(observation, o)
        push!(actions, a)
    end
    figure(figsize=(8, 8))
    plot(collect(1:length(complete)), complete, color="green")
    plot(collect(1:length(productivity)), productivity, color="orange")
    plot(collect(1:length(deadline)), deadline, color="red")
    prev = 1
    curr = 1
    curr_seg_action = ""
    # actions = fill(:work, length(actions))
    push!(actions, :end)
    linewidth = 7.0
    for action in actions
        if action != curr_seg_action
            if curr_seg_action == :work
                axvspan(prev, curr, 0, 100, color="black", alpha=.2)
            elseif curr_seg_action == :relax
                axvspan(prev, curr, 0, 100, color="pink", alpha=.2)
            else
                axvspan(prev, curr, 0, 100, color="blue", alpha=.2)
            end
            prev = curr
            curr_seg_action = action
        end
        curr += 1
    end
    savefig(path)
end

function POMDPs.actions(m::SleepPOMDP3)
    return [:work, :sleep, :relax]
end

function POMDPs.actionindex(m::SleepPOMDP3, a::Symbol)
    if a == :work
        return 1
    elseif a == :sleep
        return 2
    elseif a == :relax
        return 3
    end
    error("Invalid Action $a")
end

function POMDPs.reward(m::SleepPOMDP3, s::SleepState3, a::Symbol)
    if s.deadline == 0 && s.complete < 100
        return -100
    elseif s.complete >= 100
        return 100
    elseif a == :relax && s.hours_since_slept < 12
        return 1
    else
        return 0
    end
end

function POMDPs.gen(m::SleepPOMDP3, s::SleepState3, a::Symbol, rng::AbstractRNG)
    deadline = s.deadline - 1
    complete = s.complete
    productivity = s.productivity
    hours_since_slept = s.hours_since_slept + 1
    hours_slept = 0
    if s.deadline == 0 || s.complete >= 100
        complete = 0
        deadline = round(rand(rng, Truncated(Normal(m.max_deadline, 5), 0.0, m.max_deadline)))
    end

    if a == :work
        if s.complete < 100 && s.deadline > 0
            complete = round(s.complete + rand(rng, Truncated(Normal(10 * (s.productivity / 100), 1), 0.0, Inf)))
        end
        productivity = round(s.productivity * rand(rng, Normal(m.work_multiplier, 0.01)))
    elseif a == :sleep
        if s.hours_slept >= 6
            productivity = min(round(s.productivity * rand(rng, Normal(m.sleep_multiplier, 0.01))), 100)
            hours_since_slept = 0
        end
        hours_slept = s.hours_slept + 1
    else
        if s.hours_since_slept < 12
            productivity = min(round(s.productivity * rand(rng, Normal(m.relax_multiplier, 0.01))), 100)
        end
    end
    return (sp=SleepState3(complete, productivity, deadline, hours_slept, hours_since_slept),)
end

function POMDPs.observation(m::SleepPOMDP3, a::Symbol, sp::SleepState3)
    return Truncated(Normal(sp.complete, m.obs_sigma), 0.0, Inf)
end

function POMDPs.discount(m::SleepPOMDP3)
    return 0.95
end

POMDPs.initialstate_distribution(m::SleepPOMDP3) = Deterministic(SleepState3(0, 100, m.max_deadline, 0, 0))

using ParticleFilters
using BasicPOMCP
using ARDESPOT
problem = SleepPOMDP3()
solver = DESPOTSolver(bounds=(DefaultPolicyLB(RandomSolver()), 2000))
planner = solve(solver, problem)
filter = SIRParticleFilter(problem, 1000)
hr = HistoryRecorder(max_steps=15*12)
history = simulate(hr, problem, planner, filter)
discounted_reward(history)

sleep_plot(history, "optimal_impossible.png")


text:
using SFML

buffer = SoundBuffer(joinpath(dirname(@__FILE__),"sound.wav"))

sound = Sound(buffer)
set_loop(sound, true)

play(sound)

while true
end


text:
struct ExploitabilitySolver{G<:Game, K}
    "Game upon which a solver is to be evaluated"
    game        :: G

    "EXPLOITING Player"
    p           :: Int

    "Child information states of being in information state `I` and taking action `a`"
    Ia_children :: Dict{Tuple{K, Int}, Set{K}}

    "Reach probability weighted utility of being in information state `I` and taking action `a`"
    utility     :: Dict{K,Vector{Float64}}

    "π⁻ⁱ(h₀ → I)"
    prob        :: Dict{K,Float64}

    "Information states where player `p` is first allowed to act"
    roots       :: Set{K}
end

function reset!(e_sol::ExploitabilitySolver)
    for U in values(e_sol.utility)
        fill!(U, 0.0)
    end
    for I in keys(e_sol.prob)
        e_sol.prob[I] = 0.0
    end
    nothing
end

function ExploitabilitySolver(game::Game, p::Int)
    I = infokeytype(game)
    e_sol = ExploitabilitySolver(
        game,
        other_player(p),
        Dict{Tuple{I,Int},Set{I}}(),
        Dict{I,Vector{Float64}}(),
        Dict{I,Float64}(),
        Set{I}()
    )
    return populate!(e_sol)
end

function ExploitabilitySolver(sol::AbstractCFRSolver, p::Int)
    return ExploitabilitySolver(sol.game, p)
end

function populate!(e_sol::ExploitabilitySolver)
    _populate!(e_sol, initialhist(e_sol.game), nothing, 0)
    return e_sol
end

function _populate!(e_sol::ExploitabilitySolver, h, last_I, last_a_idx::Int)
    game = e_sol.game
    if isterminal(game, h)
        return nothing
    elseif iszero(player(game, h))
        for a in chance_actions(game, h)
            _populate!(e_sol, next_hist(game, h, a), last_I, last_a_idx)
        end
        return nothing
    end

    A = actions(game, h)
    I = infokey(game, h)

    if player(game, h) == e_sol.p
        if isnothing(last_I)
            push!(e_sol.roots, I)
        else
            Ia_children = get(e_sol.Ia_children,(last_I, last_a_idx), nothing)
            if isnothing(Ia_children)
                e_sol.Ia_children[(last_I, last_a_idx)] = Set{typeof(I)}()
            else
                push!(Ia_children, I)
            end
        end
        if !haskey(e_sol.utility, I)
            e_sol.utility[I] = zeros(Float64, length(A))
            e_sol.prob[I] = 0.0
        end
        for (k,a) in enumerate(A)
            _populate!(e_sol, next_hist(game, h, a), I, k)
        end
    else
        for a in A
            _populate!(e_sol, next_hist(game, h, a), last_I, last_a_idx)
        end
    end
    return nothing
end

function utility!(
    e_sol::ExploitabilitySolver,
    solver::AbstractCFRSolver,
    h,
    last_I,
    last_a_idx::Int,
    π_ni::Float64)

    (;game, p) = e_sol
    if isterminal(game, h)
        # π_I doesn't represent true probability, as it's being built incrementally
        # π_I increases every time it is "hit"
        # π_I = e_sol.prob[last_I]
        # e_sol.utility[last_I][last_a_idx] += (π_ni/π_I)*utility(game, p, h)
        e_sol.utility[last_I][last_a_idx] += π_ni*utility(game, p, h)
        return nothing
    elseif iszero(player(game, h)) # chance player
        A = chance_actions(game, h)
        inv_l_a = inv(length(A))
        for a in A
            utility!(e_sol, solver, next_hist(game, h, a), last_I, last_a_idx, π_ni*inv_l_a)
        end
        return nothing
    end

    A = actions(game, h)
    I = infokey(game, h)

    if player(game, h) == p
        e_sol.prob[I] += π_ni
        for (k,a) in enumerate(A)
            h′ = next_hist(game, h, a)
            utility!(e_sol, solver, h′, I, k, π_ni)
        end
        return nothing
    else
        σ = strategy(solver, I)
        for (k,a) in enumerate(A)
            h′ = next_hist(game, h, a)
            utility!(e_sol, solver, h′, last_I, last_a_idx, π_ni*σ[k])
        end
        return nothing
    end
end

function exploitability(e_sol::ExploitabilitySolver, sol::AbstractCFRSolver)
    reset!(e_sol)
    h0 = initialhist(e_sol.game)
    utility!(e_sol, sol, h0, nothing, 0, 1.0)

    for (I,U) in e_sol.utility
        U ./= e_sol.prob[I]
    end

    u_exp = 0.0
    for I in e_sol.roots
        u_exp += e_sol.prob[I]*best_response(e_sol, I)
    end
    u_current = evaluate(sol, e_sol.p)
    return u_exp - u_current
end

function best_response(e_sol::ExploitabilitySolver, I)
    util = e_sol.utility[I]
    π_A = e_sol.prob[I]

    Q_max = -Inf
    for a_i in eachindex(util)
        children = get(e_sol.Ia_children, (I,a_i), nothing)
        Qa_i = 0.0
        Qa_i += π_A*util[a_i]
        if !isnothing(children)
            for I′ in children
                Qa_i += e_sol.prob[I′]*best_response(e_sol, I′)
            end
        end
        if Qa_i > Q_max
            Q_max = Qa_i
        end
    end
    return inv(π_A)*Q_max
end

max_stars_repo_name: WhiffleFish/HelloCFR

max_issues_repo_name: WhiffleFish/HelloCFR

max_forks_repo_name: WhiffleFish/HelloCFR


text:
#import Winston

@testset "DelaunayMeshTests" begin
  @testset "AddConstraint" begin
    local mesh = DelaunayMeshes.Mesh()
    DelaunayMeshes.setboundingbox(mesh, [-15.0, 15.0, -15.0, 15.0])
    # seed rng and store seed in file
    local seed = 4139605370
    #local seed = rand(UInt32)
    local fs = open("AddConstraint.seed", "w")
    write(fs, string(seed))
    close(fs)

    # push some random vertices
    srand(seed)
    local nvertices = 200
    local points = rand(Float64, nvertices, 2)*30. - 15.
    push!(mesh, points)

    local nringvertices = 40

    # push inner ring constraint vertices
    local θ = linspace(0., 2.*π, nringvertices + 1)[1:end-1]
    local rinner = 5.
    local router = 10.
    local innerRing = rinner * [cos(θ'); sin(θ')]'
    local outerRing = router * [cos(θ'); sin(θ')]'

    push!(mesh, innerRing)
    push!(mesh, outerRing)

    local faultyTriangles = testdelaunayness(mesh.tesselation)
    @test length(faultyTriangles) == 0
    @test length(mesh.tesselation.vertices) == (nvertices + 3 + 2*nringvertices)

    # add constraints
    local innerConstraintVertexList = [x for x in (nvertices+3+nringvertices):-1:(nvertices+3+1)]
    local outerConstraintVertexList = [x for x in (nvertices+3+nringvertices+1):(nvertices + 3 + 2*nringvertices)]

    DelaunayMeshes.addconstraint!(mesh, innerConstraintVertexList)
    DelaunayMeshes.addconstraint!(mesh, outerConstraintVertexList)

    # check that constraints are in
    local innerConstrMatrix = [innerConstraintVertexList'; circshift(innerConstraintVertexList, -1)']'
    local outerConstrMatrix = [outerConstraintVertexList'; circshift(outerConstraintVertexList, -1)']'
    local constraints = [innerConstrMatrix; outerConstrMatrix]

    validateconstraints(mesh.tesselation, constraints)
    testvertexcache(mesh.tesselation)

    # get voronoiVertices
    local vorVert = DelaunayMeshes.getvoronoivertices(mesh)

    # mark all that are between inner and outer circle
    local unScaledVorVert = DelaunayMeshes.unscalepoints(mesh, vorVert)
    local unScaledSquare = unScaledVorVert.*unScaledVorVert
    local rVor = sqrt(unScaledSquare[:,1] + unScaledSquare[:,2])
    local interiorVerticesBool = (rVor .>= rinner * ones(rVor)) & (rVor .<= router * ones(rVor))

    # plot it
    #=
    xc, yc = DelaunayMeshes.getdelaunaycoordinates(mesh.tesselation)
    p = Winston.FramedPlot(aspect_ratio=1
      #,xrange=[1.496 1.50], yrange=[1.256, 1.26]
      )
    Winston.add(p, Winston.Curve(xc, yc))
    #Winston.add(p, Winston.Points(vorVert[interiorVerticesBool,1], vorVert[interiorVerticesBool,2], kind="circle", color="green"))
    Winston.add(p, Winston.Points(
      vorVert[mesh.faceLocation,1], vorVert[mesh.faceLocation,2],
      kind="circle", color="red"))
    Winston.savefig(p, "AddConstraint.svg")
    =#

    # check that this agrees with the interior array of the triangulation
    @test length(interiorVerticesBool) == length(mesh.faceLocation)
    if !(interiorVerticesBool == mesh.faceLocation)
      for (i,e) in enumerate(interiorVerticesBool)
        if (e != mesh.faceLocation[i])
          local pos = unScaledVorVert[i, :]
          local r = sqrt(pos[1]^2 + pos[2]^2)
          println("Pos $i: expecting $e, got $(mesh.faceLocation[i]). $pos and $r (scaled:$(vorVert[i,:]))")
        end
      end
    end

    @test interiorVerticesBool == mesh.faceLocation
  end

  @testset "FindConstrainedTriangles" begin
    local mesh = DelaunayMeshes.Mesh()

    # seed rng and store seed in file
    #local seed = 2572853444
    local seed = rand(UInt32)
    local fs = open("FindConstrainedTriangles.seed", "w")
    write(fs, string(seed))
    close(fs)

    # push some random vertices
    srand(seed)
    local nvertices = 20
    local points = rand(Float64, nvertices, 2)*5.0 - 2.5
    push!(mesh, points)

   # push constraint vertices
    local constraintVertices = [
      -1.0 -1.0;
      -0.5 -1.0;
      -0.5 0.0;
      0.5 0.0;
      0.5 -1.0;
      1.0 -1.0;
      1.0 1.0;
      -1.0 1.0;
    ]
    push!(mesh, constraintVertices)

    local faultyTriangles = testdelaunayness(mesh.tesselation)
    #=
    if (length(faultyTriangles) > 0)
      for f in faultyTriangles
        local ftIndex = [f[1], f[2], f[3], f[4]]
        local vertexPoints = reduce(hcat, mesh.tesselation.vertices[ftIndex])'
        PyPlot.plot(vertexPoints[:,1], vertexPoints[:,2], "o")
      end
    end
    =#
    @test length(faultyTriangles) == 0
    @test length(mesh.tesselation.vertices) == (nvertices + 8 + 3)

    # constraints
    local n = nvertices + 3

    local constraints = [
      n+1 n+2;
      n+2 n+3;
      n+3 n+4;
      n+4 n+5;
      n+5 n+6;
      n+6 n+7;
      n+7 n+8;
      n+8 n+1;
    ]
    for i=1:size(constraints, 1)
      DelaunayMeshes.insertconstraint!(mesh.tesselation, constraints[i, 1], constraints[i,2])
    end

    validateconstraints(mesh.tesselation, constraints)
    testvertexcache(mesh.tesselation)

    # get voronoiVertices
    local vorVert = DelaunayMeshes.getvoronoivertices(mesh)

    # get faces in connected region
    local startEdge = findedgeconnectingvertices(mesh.tesselation, n+1, n+2)
    local startDualEdge = DelaunayMeshes.sym(DelaunayMeshes.rot(startEdge))
    local connectedVertices = DelaunayMeshes.getfacesinsideregion(mesh, startDualEdge)
    local connectedBool = fill(false, size(mesh.tesselation.faces, 1))
    connectedBool[connectedVertices] = true

    # get all faces inside constraint for comparison
    local insideConstraint = pointsinconstrainedregion(mesh, vorVert)

    @test connectedBool == insideConstraint

    # plot it
    #=
    xc, yc = Triangulation.GetDelaunayCoordinates(mesh.tesselation)
    PyPlot.clf()
    PyPlot.plot(xc, yc)
    PyPlot.plot(get(mesh.voronoiVertices)[insideConstraint,1], get(mesh.voronoiVertices)[insideConstraint,2], "o")
    PyPlot.savefig("FindConstrainedTriangles.svg")
    =#

  end

  @testset "push" begin
    local mesh = DelaunayMeshes.Mesh()

    # seed rng and store seed in file
    #local seed = 1985005665
    local seed = rand(UInt32)
    local fs = open("DelaunayMeshPushRandom.seed", "w")
    write(fs, string(seed))
    close(fs)

    # get dimensions of bounding box
    local testBB = [rand(), rand(), rand(), rand()] * 1e5 - 0.5e4

    # push some random vertices
    srand(seed)
    local nvertices = 200
    local points = getrandomnumbersinsideboundingbox(testBB, nvertices)
    push!(mesh, points)

    # check they are inside bounding box
    local vertexPoints = reduce(hcat, mesh.tesselation.vertices[4:end])'
    local minX = minimum(vertexPoints[:,1])
    local maxX = maximum(vertexPoints[:,1])
    local minY = minimum(vertexPoints[:,2])
    local maxY = maximum(vertexPoints[:,2])

    local dh = 0.05/2.
    @test minX >= (1.5 - dh)
    @test maxX <= (1.5 + dh)
    @test minY >= (1.25 - dh)
    @test maxY <= (1.25 + dh)

    # push some more vertices that are guaranteed to be inside
    local smallBB = [
      minimum(points[:,1]), maximum(points[:,1]),
      minimum(points[:,2]), maximum(points[:,2]),
    ]

    local morePoints = getrandomnumbersinsideboundingbox(smallBB, nvertices)
    push!(mesh, morePoints)

    # check they are still inside bounding box
    vertexPoints = reduce(hcat, mesh.tesselation.vertices[4:end])'
    minX = minimum(vertexPoints[:,1])
    maxX = maximum(vertexPoints[:,1])
    minY = minimum(vertexPoints[:,2])
    maxY = maximum(vertexPoints[:,2])

    @test minX >= (1.5 - dh)
    @test maxX <= (1.5 + dh)
    @test minY >= (1.25 - dh)
    @test maxY <= (1.25 + dh)

    # check that the vertex locations are correct
    @test length(mesh.faceLocation) == length(mesh.tesselation.faces)
    @test mesh.faceLocation[1]
    @test !mesh.faceLocation[2]
    @test all(mesh.faceLocation[3:end])
  end
end


text:
#=== elastic ===#
#=== write medium for non split pml ===#
function wMT(fid::IOStream,model::elmod2d)
    medium = model.medium
    nwf = model.nwf
    iflag = model.medium.iflag

    write(fid,"medium\n")
    write(fid,"  Type: ") # free surface or unlimited medium
    if iflag == 1
      write(fid,"Free surface 2d elastic medium\n")
    elseif iflag == 2
      write(fid,"Unlimited 2d elastic medium\n")
    end
    write(fid,"  Depth(m): $(medium.DZ)\n") # Depth
    write(fid,"  HorizonX(m): $(medium.HX)\n") # HorizonX
    write(fid,"  Single Side PML Layer Number: $(medium.ext)\n") # PML
    write(fid,"  Time Step(ms): $(round(medium.dt*1000,2))\n") # dt
    write(fid,"  Spatial Step(m): dx:$(round(medium.dx,2)), dz:$(round(medium.dz,2))\n") # dx, dy & dz
    write(fid,"  Total Simulation Time(s): $(medium.T)\n") # Total Simulation Time
    write(fid,"  Discrete Total Simulation Time: $(medium.nT)\n") # Total time sampling number
    write(fid,"  Discrete Model Size(PML) [Z,X]\n") # model size with PML
    write(fid,"     P: [$(nwf.BDntpp[1]),$(nwf.BDntpp[2])]\n")
    write(fid,"   Txz: [$(nwf.BDntxz[1]),$(nwf.BDntxz[2])]\n")
    write(fid,"    Vx: [$(nwf.BDnvx[1]),$(nwf.BDnvx[2])]\n")
    write(fid,"    Vz: [$(nwf.BDnvz[1]),$(nwf.BDnvz[2])]\n")
    write(fid,"  Discrete Model Size(No PML) [Z,X,Y]\n") # model size with no PML
    write(fid,"     P: [$(nwf.ntpp[1]),$(nwf.ntpp[2])]\n")
    write(fid,"   Txz: [$(nwf.ntxz[1]),$(nwf.ntxz[2])]\n")
    write(fid,"    Vx: [$(nwf.nvx[1]),$(nwf.nvx[2])]\n")
    write(fid,"    Vz: [$(nwf.nvz[1]),$(nwf.nvz[2])]\n\n")
end

#=== elastic ===#
#=== write medium for non split pml ===#
function wMT(fid::IOStream,model::acmod2d)
    medium = model.medium
    nwf = model.nwf
    iflag = model.medium.iflag

    write(fid,"medium\n")
    write(fid,"  Type: ") # free surface or unlimited medium
    if iflag == 1
      write(fid,"Free surface 2d acoustic medium\n")
    elseif iflag == 2
      write(fid,"Unlimited 2d acoustic medium\n")
    end
    write(fid,"  Depth(m): $(medium.DZ)\n") # Depth
    write(fid,"  HorizonX(m): $(medium.HX)\n") # HorizonX
    write(fid,"  Single Side PML Layer Number: $(medium.ext)\n") # PML
    write(fid,"  Time Step(ms): $(round(medium.dt*1000,2))\n") # dt
    write(fid,"  Spatial Step(m): $(round(medium.dx,2))\n") # dx, dy & dz
    write(fid,"  Total Simulation Time(s): $(medium.T)\n") # Total Simulation Time
    write(fid,"  Discrete Total Simulation Time: $(medium.nT)\n") # Total time sampling number
    write(fid,"  Discrete Model Size(PML) [Z,X]\n") # model size with PML
    write(fid,"     P: [$(nwf.BDntpp[1]),$(nwf.BDntpp[2])]\n")
    write(fid,"    Vx: [$(nwf.BDnvx[1]),$(nwf.BDnvx[2])]\n")
    write(fid,"    Vz: [$(nwf.BDnvz[1]),$(nwf.BDnvz[2])]\n")
    write(fid,"  Discrete Model Size(No PML) [Z,X,Y]\n") # model size with no PML
    write(fid,"     P: [$(nwf.ntpp[1]),$(nwf.ntpp[2])]\n")
    write(fid,"    Vx: [$(nwf.nvx[1]),$(nwf.nvx[2])]\n")
    write(fid,"    Vz: [$(nwf.nvz[1]),$(nwf.nvz[2])]\n\n")
end


text:
using RecyclePermanent
using RNGPool
import Statistics: mean, var
import NonUniformRandomVariateGeneration.sampleCategorical
import LinearAlgebra: mul!, det
using DataFrames

function kuznetsovLogPermanent(A::Matrix, ks::Vector{Int64}, probs::Vector{Float64})
  N, n = size(A)
  @assert N >= n
  rng = getRNG()
  v::Float64 = 0.0
  for p in 1:n
    for i in 1:N
      @inbounds probs[i] = A[i, p]
    end
    for j in 1:p-1
      @inbounds probs[ks[j]] = 0.0
    end
    sp::Float64 = sum(probs)
    if sp == 0 return -Inf end
    k::Int64 = sampleCategorical(probs, sp, rng)
    v += log(sp)
    p < n && (ks[p] = k)
  end
  return v
end

function _signedBernoulli(rng)
  return 2*(rand(rng) < 0.5) - 1
end

function GodsilGutman(A::Matrix, X::Matrix, Z::Matrix)
  m, n = size(A)
  rng = getRNG()
  for i in 1:m
    for j in 1:n
      X[i,j] = sqrt(A[i,j]) * _signedBernoulli(rng)
    end
  end
  mul!(Z, X', X)
  return det(Z)
end

function ggPermanent(A::Matrix, N::Int64)
  m, n = size(A)
  @assert m >= n
  vs::Vector{Float64} = Vector{Float64}(undef, N)
  X = Matrix{Float64}(undef, m, n)
  Z = Matrix{Float64}(undef, n, n)
  for i in 1:N
    vs[i] = GodsilGutman(A, X, Z)
  end
  return vs
end

function kmPermanent(A::Matrix, ϵ::Float64, N::Int64)
  m, n = size(A)
  @assert m >= n
  vs::Vector{Float64} = Vector{Float64}(undef, N)
  ks = Vector{Int64}(undef, n-1)
  Cs = Vector{Float64}(undef, m)
  probs = Vector{Float64}(undef, m)
  for i in 1:N
    vs[i] = exp(RecyclePermanent.kmLogPermanent(A, ks, Cs, probs, ϵ))
  end
  return vs
end

function kuznetsovPermanent(A::Matrix, N::Int64)
  m, n = size(A)
  @assert m >= n
  vs::Vector{Float64} = Vector{Float64}(undef, N)
  ks = Vector{Int64}(undef, n-1)
  probs = Vector{Float64}(undef, m)
  for i in 1:N
    vs[i] = exp(kuznetsovLogPermanent(A, ks, probs))
  end
  return vs
end

function makeBinaryMatrix(m, n)
  rng = getRNG()
  c = 1.0
  for p = 1:n
    c *= m/(m-p+1)
  end
  return (rand(rng, m, n) .< 0.5) .* (2/m * c^(1/n))
end

function experiment(m, n, trials)
  A = makeBinaryMatrix(m, n)
  x = maximum(A)
  multipliers = [1000000, 10000, 100, 1, .1]
  t = @elapsed (vs = kuznetsovPermanent(A, trials))
  df = DataFrame(Algorithm = "K", TimeMS = t/trials*1000, Mean = mean(vs), Variance = var(vs))
  t = @elapsed (vs = ggPermanent(A, trials))
  push!(df, ["GG", t/trials*1000, mean(vs), var(vs)])
  for multiplier in multipliers
    mytrials = multiplier >= 1 ? trials : trials*10
    t = @elapsed (vs = kmPermanent(A, x*multiplier, mytrials))
    push!(df, ["KM: "*string(multiplier), t/mytrials*1000, mean(vs), var(vs)])
  end
  println("Experiment: m = ", m, ", n = ", n)
  show(df)
  println()
end

setRNGs(12345)

experiment(1000, 10, 10000)
experiment(1000, 100, 10000)


text:
using YaoTensorNetwork, Yao
using Test

@testset "GraphBuilder.jl" begin
    include("GraphBuilder.jl")
end

@testset "convert.jl" begin
    include("convert.jl")
end


text:
using RedPitayaDAQServer
using PyPlot
using DataFrames
using Query
using CSV

function visualizeDeltaRead(testName, start, clusterSize, decimations=[64, 32, 16, 8])
    figure(1)
    clf()
    plotIndex = 1
    for rps = start:clusterSize
        for dec in decimations
            fileName = string(testName, "_detailed_rpc_$rps", "_dec_$dec")
            fullName = joinpath(@__DIR__, string(fileName, ".csv"))
            try
                df = DataFrame(CSV.File(fullName))
                adcBufferSize  = (1 << 25)
                subplot(length(start:clusterSize), length(decimations), plotIndex)
                hlines([adcBufferSize / adcBufferSize], xmin=-10, xmax=110)
                for i = 1:rps
                    perf = @from r in df begin
                    @where r.rp == i
                    @select {r.wpRead, r.deltaRead, r.deltaSend}
                    @collect DataFrame
                    end
                    plot(perf.deltaRead / adcBufferSize)
                end
                plotIndex += 1
            catch e
                if !isa(e, ArgumentError)
                    throw(e)
                else
                    println("$fullName does not exist")
                end
            end
        end
    end
    subplots_adjust(left=0.08, bottom=0.05, right=0.98, top=0.95, wspace=0.3, hspace=0.35)
    savefig(joinpath(@__DIR__, string(testName, "_deltaRead.png")))
    println("Created ", testName, "_deltaRead.png")
end

function visualizeWaitAndSend(testName, clusterSize, decimations=[64, 32, 16, 8])
    figure(1)
    clf()
    plotIndex = 1
    for dec in decimations
        for redPitayaIndex = 1:clusterSize
            fileName = string(testName, "_detailed_rpc_$clusterSize", "_dec_$dec")
            fullName = joinpath(@__DIR__, string(fileName, ".csv"))
            subplot(length(decimations), clusterSize, plotIndex)
            try
                df = DataFrame(CSV.File(fullName))
                perf = @from r in df begin
                @where r.rp == redPitayaIndex
                @select {r.wpRead, r.deltaRead, r.deltaSend}
                @collect DataFrame
            end
                waitTime = zeros(UInt64, length(perf.wpRead))
                for i = 2:length(perf.wpRead)
                    wpLastRequestDone = perf.wpRead[i - 1] + perf.deltaRead[i - 1] + perf.deltaSend[i - 1]
                    wait = perf.wpRead[i] + perf.deltaRead[i] - wpLastRequestDone
                    waitTime[i] = wait
                end
                plot(perf.deltaSend, label="Send")
                plot(waitTime, label="Wait")
                plotIndex += 1
            catch   e
                println(e)
            end
        end
    end
    subplots_adjust(left=0.08, bottom=0.05, right=0.98, top=0.95, wspace=0.3, hspace=0.35)
    resultFile = string(testName, "_sendAndWait_cluster_$clusterSize.png")
    savefig(joinpath(@__DIR__, resultFile))
    println("Created ", resultFile)
end

function visualizeDataRate(testName, start, clusterSize, decimations=[64, 32, 16, 8])
    figure(1)
    clf()
    plotIndex = 1
    for rps = start:clusterSize
        for dec in decimations
            fileName = string(testName, "_detailed_rpc_$rps", "_dec_$dec")
            fullName = joinpath(@__DIR__, string(fileName, ".csv"))
            try
                df = DataFrame(CSV.File(fullName))
                adcBufferSize  = (1 << 25)
                subplot(length(start:clusterSize), length(decimations), plotIndex)
                for i = 1:rps
                    perf = @from r in df begin
                    @where r.rp == i
                    @select {r.deltaSend, r.chunkSize}
                    @collect DataFrame
                    end
                    plot(dataRate.(perf.chunkSize, perf.deltaSend, dec))
                end
                plotIndex += 1
            catch e
                if !isa(e, ArgumentError)
                    throw(e)
                else
                    println("$fullName does not exist")
                end
            end
        end
    end
    subplots_adjust(left=0.08, bottom=0.05, right=0.98, top=0.95, wspace=0.3, hspace=0.35)
    resultFile = string(testName, "_dataRate.png")
    savefig(joinpath(@__DIR__, resultFile))
    println("Created ", resultFile)
end

max_stars_repo_name: tknopp/RedPitayaDAQServer

max_issues_repo_name: tknopp/RedPitayaDAQServer

max_forks_repo_name: tknopp/RedPitayaDAQServer


text:
function differentiate(obj::UnivariateMatrixFunctionDifferentiator, arg0::UnivariateMatrixFunction)
    return jcall(obj, "differentiate", UnivariateDifferentiableMatrixFunction, (UnivariateMatrixFunction,), arg0)
end


text:
using Documenter
using NSDEBase

PAGES = ["Home" => "index.md"]

makedocs(;
    sitename = "NSDEBase.jl",
    format = Documenter.HTML(),
    modules = [NSDEBase],
    pages = PAGES,
    authors = "Giancarlo A. Antonucci <giancarlo.antonucci@icloud.com>"
)

deploydocs(;
    repo = "https://github.com/giancarloantonucci/NSDEBase.jl"
)


text:
function new_struct(fields::Vector{Tuple{String,DataType}})
    name = "A" * string(hash(fields), base = 16)
    c = Expr(
        :block,
        Expr(
            :struct,
            false,
            Symbol(name),
            Expr(:block, [Expr(:(::), Symbol(f[1]), f[2]) for f in fields]...),
        ),
        Symbol(name),
    )
    eval(c)
end

myS = new_struct([("a", Int), ("b", String), ("c", Int)])
dump(myS)

code = Expr(
    :block,
    Expr(
        :struct,
        false,
        Symbol("Goma"),
        Expr(:block, Expr(:(::), Symbol("x"), Int), Expr(:(::), Symbol("y"), Int)),
    ),
    Symbol("Goma"),
)
Goma = eval(code)
println(fieldnames(Goma))
println(Goma(1, 2))

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using Test
using FlightMechanics.Models
using FlightMechanics

llh = [deg2rad(49.996908), 0.000000, 9907.31]
xyz_ecef = [4114496.258, 0.0, 4870157.031]
xyz_earth = [0., 0., -9907.31]

pllh = PositionLLH(llh...)

@test isapprox(get_llh(pllh), llh)
@test isapprox(get_xyz_earth(pllh), xyz_earth)
@test isapprox(get_xyz_ecef(pllh), xyz_ecef)
@test isapprox(get_height(pllh), llh[3])

xyz_earth = [0., 0., -9907.31]
pearth = PositionEarth(xyz_earth..., llh[1:2]...)

@test isapprox(get_llh(pearth), llh)
@test isapprox(get_xyz_earth(pearth), xyz_earth)
@test isapprox(get_xyz_ecef(pearth), xyz_ecef)
@test isapprox(get_height(pearth), -xyz_earth[3])

pecef = PositionECEF(xyz_ecef...)

@test isapprox(get_llh(pecef), llh, rtol=1e-5)
@test isapprox(get_xyz_earth(pecef), xyz_earth, rtol=1e-5)
@test isapprox(get_xyz_ecef(pecef), xyz_ecef)
@test isapprox(get_height(pecef), llh[3], atol=0.17)

pzero = Position()

@test isapprox(get_llh(pzero), [0, 0, 0])
@test isapprox(get_xyz_earth(pzero), [0, 0, 0])
@test isapprox(get_xyz_ecef(pzero), llh2ecef(0, 0, 0))
@test isapprox(get_height(pzero), 0)


text:
       # centroid
        @test centroid(tri) ≈ Point2D{T}(1//3, 1//3)

                # centroid
        @test centroid(tri) ≈ Point3D{T}(0, 1//3, 1//3)

   # centroid
        @test centroid(quad) ≈ Point2D{T}(1//2, 1//2)

                    # centroid
            p₁ = Point2D{T}(0, 0)
            p₂ = Point2D{T}(1, 0)
            p₃ = Point2D{T}(0, 1)
            p₄ = Point2D{T}(1//2, 0)
            p₅ = Point2D{T}(1//2, 1//2)
            p₆ = Point2D{T}(0, 1//2)
            tri6 = QuadraticTriangle(p₁, p₂, p₃, p₄, p₅, p₆)
            @test centroid(tri6) ≈ tri6(1//3, 1//3)

        # centroid
        p₁ = Point2D{T}(0, 0)
        p₂ = Point2D{T}(1, 0)
        p₃ = Point2D{T}(1, 1)
        p₄ = Point2D{T}(0, 1)
        p₅ = Point2D{T}(1//2,    0)
        p₆ = Point2D{T}(   1, 1//2)
        p₇ = Point2D{T}(1//2,    1)
        p₈ = Point2D{T}(   0, 1//2)
        quad8 = QuadraticQuadrilateral(p₁, p₂, p₃, p₄, p₅, p₆, p₇, p₈)
        @test centroid(quad8) ≈ quad8(1//2, 1//2)

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
module VectorAutoregressions

using DataFrames, TimeSeries
using Distributions: ccdf, TDist, Chisq, FDist
using Dates, LinearAlgebra
using StatsBase: CoefTable, autocov, autocor, crosscov, crosscor

import Base: show
import StatsBase: loglikelihood, aic, sample, coef, stderror, residuals, coeftable, fitted

export VAREstimate,
    VAR,
    LagSelection,
    StabilityCheck,
    PortmanteauTest,
    LMCorrTest,
    loglikelihood,
    aic,
    sic,
    hqc,
    coef,
    stderror,
    residuals,
    coeftable,
    fitted

include("matrix_utilities.jl")
include("VAR.jl")
include("statsbase.jl")
include("diagnostics.jl")
include("lag_selection.jl")
include("IRF.jl")

end # module


text:
module Homework

using Reactive

using Interact

using JSON
using Requests
import Requests: get

using Requires

using IJulia.CommManager

const global_config = Dict()

const state = Dict()

include("library.jl")

function reinit()

    p = joinpath(splitdir(@__FILE__)[1], "homework.js")

    # Load javascript into IJulia
    display(MIME"text/html"(),
        """<script>$(readall(p))</script>""")

    # Push some metadata to the Julia side
    comm = Comm(:HomeworkData)
    display(script("IPython.notebook.kernel.comm_manager.comms[" * JSON.json(comm.id) * """].then(
        function (comm) {
            comm.send({
                cookie: document.cookie,
                host: IPython.notebook.metadata.homework.host,
                course: IPython.notebook.metadata.homework.course,
                mode: IPython.notebook.metadata.homework.mode,
                problemset: IPython.notebook.metadata.homework.problemset
            })
        }
    )"""))

    #                      |
    #                      \
    #                      |
    #                      \
    #                      |
    #                      v

    comm.on_msg = msg -> begin
        global_config["cookie"] = msg.content["data"]["cookie"]

        if haskey(msg.content["data"], "host")
            global_config["host"] = msg.content["data"]["host"]
        end

        global_config["course"] = msg.content["data"]["course"]
        global_config["problemset"] = msg.content["data"]["problemset"]

        if haskey(msg.content["data"], "mode")
            global_config["mode"] = msg.content["data"]["mode"]
        end
    end

    # state["refresh_dashboard"] = (b = button("Refresh")) |> signal
    # display(b)
    # (state["dashboard"] = Input(Html("Loading more..."))) |> display
end

__init__() = try reinit(); catch ex; end

include("encode.jl")
include("user.jl")
include("admin.jl")

end # module


text:
"""
    scan_x(scanner, x_range, x_num_scans; kwargs..)
    scan_y(scanner, y_range, y_num_scans; kwargs..)
    scan_z(scanner, z_range, z_num_scans; kwargs..)

Move along a single axis and grab data based on a fixed increment

# Arguments
- `scanner::IntensityScan`: the structure for all controller connections.
- `x_range`: An array or tuple that contains the start and end
             point of the desired range in the x-direction.
             (Must be 1 or 2 elements long)
- `x_num_scans::Int`: The number of points/scans the user desires
                      to record along the x-axis.

# Keywords
- `verbose`::Bool: Optional, defaults to: true

# Returns
- `ScanData`

# Example
```
ps = initialize(PositionerSystem)
scope = initialize(AgilentDSOX4034A)

scanner = IntensityScan(
    xyz = ps,
    scope = scope,
    channel = 1,
    post_move_delay = 0ms,
    precapture_delay = 0µs,
)

wave_x = scan_x(scanner, [-10mm, 10mm], 5)
```
"""
function scan_x(scanner::IntensityScan, x_range, num_points; verbose=true)
    scan_data = scan_single_axis(scanner, x_range, num_points, move_x_abs; verbose=verbose)
    return scan_data
end

function scan_y(scanner, y_range, num_points; verbose=true)
    scan_data = scan_single_axis(scanner, y_range, num_points, move_y_abs; verbose=verbose)
    return scan_data
end

function scan_z(scanner, z_range, num_points; verbose=true)
    scan_data = scan_single_axis(scanner, z_range, num_points, move_z_abs; verbose=verbose)
    return scan_data
end

function scan_single_axis(scanner, axis_range, num_scans, move_func; verbose=true)
    axis = get_axis(move_func)
    check_xyz_limits(scanner, axis, axis_range)
    positions = create_positions_vector(axis_range, num_scans)

    # construct the ScanData struct
    scope_info = TcpInstruments.scope_waveform_info_get(scanner.scope, scanner.channel)
    dummy_data = TcpInstruments.scope_parse_raw_waveform(zeros(scope_info.num_points), scope_info)
    time_points = dummy_data.time
    scan_data = ScanData(axis, scope_info, time_points, scope_info.num_points, [num_scans])

    pm = Progress(num_scans; desc="Scanning $axis-direction: ", enabled=verbose)
    for scan_index in 1:num_scans
        move_func(scanner.xyz, positions[scan_index])
        sleep(raw(scanner.post_move_delay))
        ProgressMeter.update!(pm)

        scanner.trigger_function()

        sleep(raw(scanner.precapture_delay))
        data = get_data(scanner.scope, scanner.channel)
        scan_data.waveform[:, scan_index] = data.volt
        scan_data.coordinates[:, scan_index] = get_pos(scanner.xyz)

        ProgressMeter.next!(pm)
    end

    return scan_data
end

function create_positions_vector(axis_range, num_points)
    if length(axis_range) == 1
        if num_points != 1
            error(
                "Can't measure multiple points when axis_range has one position"
            )
        end
        return axis_range
    elseif length(axis_range) == 2
        return range(axis_range[begin], stop=axis_range[end], length=num_points)
    else
        error("axis_range must have 1 or 2 elements")
    end
end

function check_xyz_limits(
    scanner::IntensityScan, axis::AbstractString, axis_range::Array
)
    if axis == "x"
        check_xyz_limits(get_limits_x(scanner.xyz), axis_range)
    elseif axis == "y"
        check_xyz_limits(get_limits_y(scanner.xyz), axis_range)
    elseif axis == "z"
        check_xyz_limits(get_limits_z(scanner.xyz), axis_range)
    else
        error("Single axis must be x, y or z not: $axis")
    end
end

function check_xyz_limits(limits, axis_range::Array)
    low_limit, high_limit = limits
    for range in axis_range
        if range < low_limit || range > high_limit
            error("Scan range ($range) is outside current limits of the xyz stage ($low_limit - $high_limit)")
        end
    end
end

get_axis(move_func::Function) = split(string(move_func), '_')[2]


text:
struct Wave

    amps::Vector
    time
    framerate::Number

    signal::Signal
    #isnoise::Bool

end

function Wave(s::Signal; t = 0, framerate = 44_100)

    if t == 0
        freq = s.freq
        if !(s.freq |> typeof <: Number)
            freq = minimum(s.freq)
        end

        t = inv(freq) * 5

    end

    dt = inv(framerate)

    r = range(0, t, step = dt)

    amps = s.(r) # evaluating the signal

    if haskey(s.info, "WhiteNoise")
        amps = amps + WhiteNoise(s.info["WhiteNoise"])
    end

    if haskey(s.info, "GaussianNoise")
        amps = amps + GaussianNoise(s.info["GaussianNoise"])
    end

    if haskey(s.info, "BrownianNoise")
        amps = amps + BrownianNoise(s.info["BrownianNoise"])
    end

    if haskey(s.info, "PinkNoise")
        amps = amps + PinkNoise(s.info["PinkNoise"]...)
    end

    Wave(amps, r, framerate, s)

end

function wavplay(w::Wave)

    wavplay(w.amps, w.framerate)

end

@recipe function plot(w::Wave)

    if haskey(w.signal.info, "name")
        label --> (w.signal.info["name"] * " framerate $(w.framerate)")
    else
        label --> "Framerate $(w.framerate)"
    end
    yguide --> "Amplitude"
    xguide --> "Time in s"

    w.time, w.amps

end

function filter(w::Wave, type::Symbol, args...)

    amps = deepcopy(spec.amps)

    if type == :band
        @assert length(args) == 2
        @assert args[1] < args[2]
        amps[spec.fs .< args[1] ] .= 0
        amps[spec.fs .> args[2] ] .= 0
    end

    if type == :low
        @assert length(args) == 1
        amps[spec.fs .> args[1] ] .= 0
    end

    if type == :high
        @assert length(args) == 1
        amps[spec.fs .< args[1] ] .= 0
    end

    return Spectrum(amps, spec.fs, abs2.(amps), spec.signal, spec.framerate, spec.d)

end


text:
using Denumpyfy
using Test

@testset "Denumpyfy.jl" begin
    str = "x ** 2"
    @test replace_exponents(str)  == "x ^ 2"

    str = "range(0, self.n_grid):"
    @test replace_range(str) == "1:self.n_grid"
    str = "range(0, 1):"
    @test replace_range(str) == "1:2"
    str = "range(self.n_grid, 3):"
    @test replace_range(str) == "self.n_grid:4"
    str = "range(1, n_grid - 1):"
    @test replace_range(str) == "2:n_grid"
    str = "    for i in range(1, n_grid - 1):"
    @test replace_range(str) == "    for i in 2:n_grid"
    str = "range(1, 4):"
    @test replace_range(str) == "2:5"
    str = "range(1, 5):"
    @test replace_range(str) == "2:6"
    str = "range(1, n_grid):"
    @test replace_range(str) == "2:n_grid"
    str = "range(xs):"
    @test replace_range(str) == "1:xs"
    str = "    for i in range(1, 20):"
    @test replace_range(str) == "    for i in 2:21"
    str = "    for i in range(1, n_grid - 1):"
    @test replace_range(str) == "    for i in 2:n_grid"

    str = "len(xs)"
    @test replace_len(str) == "length(xs)"

    str = "linspace(half_delta - x_out, x_out - half_delta, n_grid)"
    @test replace_linspace(str) == "LinRange(half_delta - x_out, x_out - half_delta, n_grid)"

    str = "def construct(self, tol, it_max):"
    @test replace_def(str) == "function construct(self, tol, it_max)"

    str = "if x > 3:"
    @test replace_if(str) == "if x > 3"

    str = "while x < 4 && i == 0:"
    @test replace_while(str) == "while x < 4 && i == 0"

    str = "else:"
    @test replace_else(str) == "else"

    str = "zeros((n_grid, n_grid, n_grid))"
    @test replace_zeros(str) == "zeros(n_grid, n_grid, n_grid)"
    str = "zeros((n_grid, n_grid))"
    @test replace_zeros(str) == "zeros(n_grid, n_grid)"

    # TODO empty(()) ???

    str = "ones((n_grid, n_grid, n_grid))"
    @test replace_ones(str) == "ones(n_grid, n_grid, n_grid)"
    str = "ones((n_grid, n_grid))"
    @test replace_ones(str) == "ones(n_grid, n_grid)"

    str = "1 "
    @test handle_indexing_arg(str) == "2"
    str = "-1"
    @test handle_indexing_arg(str) == "end"
    str = " -2"
    @test handle_indexing_arg(str) == "end-1"
    strs = [" -2", "-1"]
    @test handle_indexing_arg.(strs) == ["end-1", "end"]

    strs = "alpha beta gamma sigma"
    @test replace_greek(strs) == "α β γ σ"

    str = "np.eye(5)"
    @test replace_eye(str) == "np.I(5)"

    str = "[0]"
    #@test replace_indexing(str) == "[1]"
    str = "[-1, n_grid, rho]"
    #@test replace_indexing(str) == "[end,n_grid,rho]"

    # Test multiple repeated indexes :D
    str = "ax[0] = ax[3,4] + x[3,4] + y[-1,y,5]"
    @test replace_indexing(str) == "ax[begin] = ax[4,5] + x[4,5] + y[end,y,6]"
end