text:
# dqn.jl
# original: https://github.com/denizyuret/Knet.jl/blob/master/examples/reinforcement-learning/dqn/dqn_with_target.jl

module DQN

using Gym, Knet, JLD

include("replay_buffer.jl")
include("mlp.jl")
include("piecewise_schedule.jl")

function loss(w::WeightParamsR{T}, states::VecOrMat{T}, actions, targets; nh=1) where {T}
    qvals = predict_q(w, states; nh=nh)
    nrows = size(qvals, 1)
    index = actions .+ nrows .* (0:length(actions)-1)
    qpred = reshape(qvals[index], size(targets))
    mse = sum(abs2, targets .- qpred) / size(states, 2)
    return mse
end

lossgradient = gradloss(loss)

function train!(w::WeightParams{T}, prms, states, actions, targets; nh=1) where {T}
    g, mse = lossgradient(w, states, actions, targets; nh=nh)
    update!(w, g, prms)
    return mse
end

function dqn_learn(w::WeightParams{T}, opts, env, buffer, exploration, o) where {T}
    total = 0.0
    readytosave = save_interval = get(o, "save_interval", 10000)
    episode_rewards = Float32[]
    frames = Float32[]
    ob_t = reset!(env)

    n_hiddens = get(o, "n_hiddens", length(o["hiddens"]))::Int
    target_w = o["play"] || o["tupdate"] <= 1 ? w : deepcopy(w)

    for fnum = 1:o["frames"]
        o["render"] && render(env)
        #process the raw ob
        ob_t_reshaped = reshape(ob_t, size(ob_t)..., 1)
        if !o["play"] && rand() < value(exploration, fnum)
            a = sample(env.action_space)
        else
            obses_t = encode_recent(buffer, ob_t_reshaped; stack=o["stack"])
            inp = convert(o["atype"], obses_t)
            qvals = predict_q(w, inp; nh=n_hiddens)
            a = indmax(Array(qvals)) - 1
        end

        ob_t, reward, done, _ = step!(env, a)
        total += reward

        if !o["play"]
            #process the raw ob
            ob_tp1_reshaped = reshape(ob_t, size(ob_t)..., 1)
            push!(buffer, ob_t_reshaped, a + 1, reward, ob_tp1_reshaped, done)

            if can_sample(buffer, o["bs"])
                obses_t, actions, rewards, obses_tp1, dones = sample_batch(buffer, o["bs"]; stack=o["stack"])
                obses_tp1 = convert(o["atype"], obses_tp1)
                #predict next q values with the target network
                nextq = Array(predict_q(target_w, obses_tp1; nh=n_hiddens))
                maxs = maximum(nextq, 1)
                nextmax = sum(nextq .* (nextq .== maxs), 1)
                targets = reshape(rewards, 1, :) .+ (o["gamma"] .* nextmax .* dones)
                mse = train!(w, opts, convert(o["atype"], obses_t), actions, convert(o["atype"], targets); nh=n_hiddens)
            end

            if !o["no_save"] && fnum > readytosave
                model_filepath = "$(get(o, "model_prefix", ""))_$(fnum).jld"
                save_model(w, model_filepath)
                readytosave += save_interval
            end

            if o["tupdate"] > 1 && fnum % o["tupdate"] == 0
                target_w = deepcopy(w)
            end
        end

        if done
            ob_t = reset!(env)
            # o["printinfo"] && println("Frame: $fnum , Total reward: $total, Exploration Rate: $(value(exploration, fnum))")
            if o["printinfo"]
                if o["play"]
                    println("Frame: $fnum , Total reward: $total")
                else
                    println("Frame: $fnum , Total reward: $total, Exploration Rate: $(value(exploration, fnum))")
                end
            end
            push!(episode_rewards, total)
            push!(frames, fnum)
            total = 0.0
        end
    end
    return episode_rewards, frames
end

end # module

max_stars_repo_name: JuliaTagBot/JuliaGymDemo

max_issues_repo_name: JuliaTagBot/JuliaGymDemo

max_forks_repo_name: JuliaTagBot/JuliaGymDemo


text:
function test_cublas(backend::Backend, T)
  println("-- Testing CuBLAS Utilities on $T")
  dims = (2,3,5,7)
  eps = 1e-5

  println("    > memory copy between device and host")
  x = rand(T, dims)
  x2 = similar(x)
  x_blob = make_blob(backend, T, dims)
  x2_blob = make_blob(backend, T, dims)
  copy!(x_blob, x)
  copy!(x2, x_blob)
  @test all(-eps .< x2-x .< eps)

  println("    > memory copy between devices")
  fill!(x2, 0)
  copy!(x2_blob, x_blob)
  copy!(x2, x2_blob)
  @test all(-eps .< x2-x .< eps)

  println("    > axpy")
  CuBLAS.axpy(backend.cublas_ctx, length(x), convert(T,2.0), x_blob.ptr, 1, x2_blob.ptr, 1)
  copy!(x2, x2_blob)
  @test all(-eps .< 3*x - x2 .< eps)

  println("    > scal")
  copy!(x_blob, x)
  CuBLAS.scal(backend.cublas_ctx, length(x), convert(T,0.9), x_blob.ptr, 1)
  copy!(x2, x_blob)
  @test all(-eps .< 0.9*x - x2 .< eps)

  println("    > gemm")
  A = rand(T, 5,6)
  B = rand(T, 6,7)
  C = rand(T, 5,7)
  blob_A = make_blob(backend, T, (5,6, 1, 1))
  blob_B = make_blob(backend, T, (6,7, 1, 1))
  blob_C = make_blob(backend, T, (5,7, 1, 1))
  copy!(blob_A, A); copy!(blob_B, B); copy!(blob_C, C)
  CuBLAS.gemm(backend.cublas_ctx, CuBLAS.OP_N, CuBLAS.OP_N, 5, 7, 6,
      convert(T,1.0), blob_A.ptr, 5, blob_B.ptr, 6, convert(T,1.0), blob_C.ptr, 5)
  result = zeros(T, 5,7)
  copy!(result, blob_C)
  @test all(-eps .< A*B+C - result .< eps)

  println("    > blas copy")
  x = rand(T, size(x))
  copy!(x_blob, x)
  CuBLAS.copy(backend.cublas_ctx, T, length(x), x_blob.ptr.p, 1, x2_blob.ptr.p, 1)
  copy!(x2, x2_blob)
  @test all(-eps .< x2-x .< eps)

  println("    > fill!")
  fill!(x_blob, 0)
  @test all(-eps .< to_array(x_blob) .< eps)
  fill!(x_blob, 1)
  @test all(-eps .< to_array(x_blob) - 1 .< eps)
end

function test_cublas(backend::Backend)
  test_cublas(backend, Float32)
  test_cublas(backend, Float64)
end

test_cublas(backend_gpu)

ext: jl


text:
# This code calls the first-order eccentric TTV code
# compute_ttv.jl.  Please cite Agol & Deck (2015) if
# you make use of this in your research.

include("compute_ttv.jl")
#include("laplace_coefficients_initialize.jl")

function test_ttv(jmax::Integer,n1::Integer,n2::Integer,data::Vector{Float64}, WriteOutput::Bool = true, num_evals::Integer = 1, profile::Bool = false)
  @assert(jmax>=1)  # Should there be a larger minimum?
  @assert(n1>2)
  @assert(n2>2)
  @assert(length(data)==10)
  # Performs a test of the transit_times.jl routine
  # Set up planets planar-planet types for the inner and outer planets:
  p1=TTVFaster.Planet_plane_hk(data[1],data[2],data[3],data[4],data[ 5])
  p2=TTVFaster.Planet_plane_hk(data[6],data[7],data[8],data[9],data[10])
  time1 = collect(p1.trans0 .+ range(0,stop=n1-1,length=n1) .* p1.period)
  time2 = collect(p2.trans0 .+ range(0,stop=n2-1,length=n2) .* p2.period)
  alpha0=(p1.period/p2.period)^(2//3)
  # Initialize the computation of the Laplace coefficients:
  # b0=TTVFaster.LaplaceCoefficients.initialize(jmax+1,alpha0)
  # Define arrays to hold the TTVs:
  #ttv1=Array(Float64,n1)
  #ttv2=Array(Float64,n2)
  ttv_el_type = eltype(data) == Float64 ? Float64 : Number
  ttv1=zeros(ttv_el_type,n1)
  ttv2=zeros(ttv_el_type,n2)
  # Define arrays to hold the TTV coefficients and Laplace coefficients:
  f1=zeros(Float64,jmax+2,5)
  f2=zeros(Float64,jmax+2,5)
  b=zeros(Float64,jmax+2,3)
  hashsum = 0
  for i in 1:num_evals
   # Call the compute_ttv code which implements equation (33)
   TTVFaster.compute_ttv!(jmax,p1,p2,time1,time2,ttv1,ttv2)#,f1,f2,b,alpha0,b0)
    if profile
      hashsum += hash(ttv1)+hash(ttv2)
    end
  end
  if profile
    println("# Ignore this: ",hashsum) # This just makes sure optimizer doesn't optimize away important calculations.
  end

  if WriteOutput
    # Write the mean ephemeris and TTV results to two files:
    writedlm("inner_ttv.txt",[time1 ttv1])
    writedlm("outer_ttv.txt",[time2 ttv2])
  end
  return ttv1,ttv2
end

max_stars_repo_name: bmlindor/TTVFaster

max_issues_repo_name: bmlindor/TTVFaster

max_forks_repo_name: bmlindor/TTVFaster


text:
# unary.jl: Unary Array->Array operations.
# uses unary_ops from ops.jl.

using SpecialFunctions
import Base.Broadcast: broadcasted
import NNlib: relu, selu, elu

function unary_op(f, j=f, o...)
    J=Symbol(j)
    M = which(@__MODULE__, J)
    for S in (32,64)
        T = Symbol("Float$S")
        F = "$(f)_$S"
        @eval begin
            function broadcasted(::typeof($J),x::KnetArray{$T})
                y = similar(x)
                @knet8($F,(Cint,Ptr{$T},Ptr{$T}),length(y),x,y)
                return y
            end
            # Bcasted methods
            ($M).$J(x::Bcasted{<:KnetArray{$T}}) = broadcasted($J, x.value) |> Bcasted
            broadcasted(::typeof($J),x::Bcasted{<:KnetArray{$T}}) = broadcasted($J, x.value) |> Bcasted
        end
    end
    @eval begin # so we do not trigger some default Base implementation
        ($M).$J(x::Bcasted) = throw(MethodError($J,(x,)))
        broadcasted(::typeof($J),x::Bcasted) = throw(MethodError($J,(x,)))
    end
end

# Constants for selu activation from https://arxiv.org/abs/1706.02515
const λ01 = (1-erfc(1/sqrt(2))*sqrt(exp(1)))*sqrt(2pi)*(2*erfc(sqrt(2))*exp(2)+pi*erfc(1/sqrt(2))^2*exp(1)-2*(2+pi)*erfc(1/sqrt(2))*sqrt(exp(1))+pi+2)^(-0.5)
const α01 = -sqrt(2/pi)/(erfc(1/sqrt(2))*exp(1/2)-1)
const λα01 = λ01 * α01

# Define some common operations as primitives for efficiency:
# 1. Avoid creating intermediate arrays
# 2. Avoid taking derivatives of intermediate operations

for (f,g,y,dx) in
    ((:invx, :invxback, :(one(T)/xi), :(-yi*yi*dyi)),
     (:relu, :reluback, :(max(zero(T),xi)), :(ifelse(yi>0,dyi,zero(T)))),
     (:selu, :seluback, :(xi >= 0 ? T(λ01)*xi : T(λα01)*(exp(xi)-1)), :(yi >= 0 ? dyi * T(λ01) : dyi * (yi + T(λα01)))),
     (:elu,  :eluback,  :(xi >= 0 ? xi : exp(xi)-1), :(yi >= 0 ? dyi : dyi * (1+yi))),
     (:tanx, :tanhback, :(tanh(xi)), :(dyi*(one(T)-yi*yi))),
     (:sigm, :sigmback,
      # Numerically stable implementation from
      # http://timvieira.github.io/blog/post/2014/02/11/exp-normalize-trick
      :(if xi>=0; z=exp(-xi); one(T)/(one(T)+z); else; z=exp(xi); z/(one(T)+z); end),
      :(dyi*yi*(one(T)-yi))),
     )
    @eval begin
        $f(xi::T) where {T<:Number}=$y
        $g(dyi::T,yi::T) where {T<:Number}=$dx
        function broadcasted(::typeof($f),x::Array{T}) where {T<:AbstractFloat}
            y = similar(x)
            @inbounds for i=1:length(y)
                xi = x[i]
                y[i] = $y
            end
            return y
        end
        function broadcasted(::typeof($g),dy::Array{T},y::Array{T}) where {T<:AbstractFloat}
            dx = similar(dy)
            @inbounds for i=1:length(dx)
                yi = y[i]
                dyi = dy[i]
                dx[i] = $dx
            end
            return dx
        end
        @primitive $f(x),dy,y $g.(dy,y)
    end
end

"`invx(x) = 1/x`" invx
"`sigm(x) = 1/(1+exp(-x))`" sigm

"""
    relu(x)
Return `max(0,x)`.

References:
* [Nair and Hinton, 2010](https://icml.cc/Conferences/2010/abstracts.html#432). Rectified Linear Units Improve Restricted Boltzmann Machines. ICML.
* [Glorot, Bordes and Bengio, 2011](http://proceedings.mlr.press/v15/glorot11a). Deep Sparse Rectifier Neural Networks. AISTATS.
"""
relu

"""
    elu(x)

Return `(x > 0 ? x : exp(x)-1)`.

Reference: Fast and Accurate Deep Network Learning by Exponential Linear Units (ELUs) (https://arxiv.org/abs/1511.07289).
"""
elu

"""
    selu(x)

Return `λ01 * (x > 0 ? x : α01 * (exp(x)-1))` where `λ01=1.0507009873554805` and `α01=1.6732632423543778`.

Reference: Self-Normalizing Neural Networks (https://arxiv.org/abs/1706.02515).
"""
selu

# To avoid conflict with AutoGrad:
import Base: tanh
@primitive tanh(x::Array),dy,y     tanhback.(dy,y)
@primitive tanh(x::KnetArray),dy,y tanhback.(dy,y)

# For 2nd derivatives:
@primitive tanhback(dy,y),ddx  ddx.*(1 .- y.*y)  ddx.*(-2 .* dy.*y)
@primitive reluback(dy,y),ddx  ddx.*(y.>0)       nothing
@primitive invxback(dy,y),ddx  ddx.*(-y.*y)      ddx.*(-2 .* dy .* y)
@primitive sigmback(dy,y),ddx  ddx.*y.*(1 .- y)  ddx.*dy.*(1 .- 2 .* y)
@primitive seluback(dy,y),ddx  ddx.*((y.>=0).*λ01.+(y.<0).*(y.+λα01))  ddx.*dy.*(y.<0)
@primitive  eluback(dy,y),ddx  ddx.*((y.>=0).+(y.<0).*(y.+1))          ddx.*dy.*(y.<0)

# Unary plus and minus
import Base: +, -

broadcasted(::typeof(+), a::KnetArray)=a
+(a::KnetArray)=a
-(a::KnetArray)=broadcasted(-,a)

# Identity: Issue #335
broadcasted(::typeof(identity), a::KnetArray)=a

# TODO: SpecialFunctions loggamma update #486.
# SpecialFunctions 0.8: lgamma(x::Real)` is deprecated, use `(logabsgamma(x))[1]` instead.
# Other alternative is loggamma, throws a DomainError if gamma(x) is negative.
# However neither are equivalent to cuda lgamma.
# commenting this out until we implement a better solution.
# if !isdefined(SpecialFunctions, :loggamma) && isdefined(SpecialFunctions, :lgamma)
#     loggamma(x) = lgamma(x)
# end
# if isdefined(SpecialFunctions, :loggamma) && !isdefined(SpecialFunctions, :lgamma)
#     lgamma(x) = loggamma(x)
# end
# if isempty(methods(loggamma,(AutoGrad.Value,))) # Should be fixed after AutoGrad 1.1.4
#     @primitive loggamma(x),dy,y (dy.*(digamma.(x)))
# end

# Unbroadcasted zero works on arrays: this moved to karray.jl
# import Base: zero
# zero(x::KnetArray)=zero.(x)

for f in unary_ops
    if !isa(f,Tuple); f=(f,); end
    unary_op(f...)
end


text:
export leastSquares

# if Φ is a row matrix and is passed as a 1 dimensional array
# function leastSquares(Φ::AbstractVector{<:Real}, Y::AbstractVector{<:Real})
#     Φ = reshape(Φ, length(Φ), 1)
#     leastSquares(Φ, Y)
# end

# Ordinary least-squares regression
function leastSquares(Φ::AbstractMatrix{<:Real}, Y::AbstractVector{<:Real})
    @assert size(Φ, 1) == length(Y) "Dimension mismatch in matrix Φ and observation vector Y."

    # FUTURE: need to account for correlated data, if present

    coefficients = zeros(size(Y))

    condition = 1/cond(Φ) # Matrix condition number

    # Fastest, but least accurate (squares condition number)
    # if condition < 1e6
    #     ΦTΦ = Φ' * Φ
    #     coefficients = ΦTΦ \ (Φ' * Y)

    # (Much) more precise but ~ twice as slow
    if condition > eps()
        coefficients = Φ \ Y

    # Slowest, but best precision for very ill-conditioned matrices. For better conditioned matrices it is not as good as other methods.
    else
        @warn("Matrix condition number is really high. Results can become inaccurate.")
        coefficients = pinv(Φ) * Y # TODO: calibrate tolerances of pinv
    end

    return coefficients
end

max_stars_repo_name: adriangrupp/sparsePolyChaos

max_issues_repo_name: adriangrupp/sparsePolyChaos

max_forks_repo_name: adriangrupp/sparsePolyChaos


text:
###############################################################

start = 0
while start < 4
	score = (w,l,t)
	print("Welcome to a game of Rock, Paper, Scissors, Lizard, Spock!\n
	Enter game type: \n1 for Player vs. Computer \n2 for Computer vs. Computer \n3 for Player vs. Player \n4 to Quit\n")
	data = readline()
	start = parse(Int64,data)
	if start == 1
		w,l,t = 0,0,0
		while w < 99
		score = (w,l,t)
		print("\nScore: $score\n")
		v1 = "Rock"
		v2 = "Paper"
		v3 = "Scissors"
		v4 = "Lizard"
		v5 = "Spock"

		c_select = rand((v1,v2,v3,v4,v5))
		print("Please select $v1, $v2, $v3, $v4, or $v5.\n")
		h_select = string(readline()) |> lowercase |> titlecase

		print("$v1,")
		sleep(0.5)
		print(" $v2,")
		sleep(0.5)
		print(" $v3,")
		sleep(0.5)
		print(" $v4")
		sleep(0.5)
		print(" $v5,")
		sleep(0.5)
		print(" Shoot!\n")

		# check yours and computers entry for outcome: Win, Lose, Tie.
		if h_select == v1 && c_select == v1
			println("You both picked $v1. It's a tie!")
			t += 1
		elseif h_select == v1 && c_select == v2
			println("$v2 covers $v1. You lose!")
			l += 1
		elseif h_select == v1 && c_select == v3
			println("$v1 crushes $v3. You win!")
			w += 1
		elseif h_select == v1 && c_select == v4
			println("$v1 crushes $v4. You win!")
			w += 1
		elseif h_select == v1 && c_select == v5
			println("$v5 vaporizes $v1. You lose!")
			l += 1
		elseif h_select == v2 && c_select == v1
			println("$v2 covers $v1. You win!")
			w += 1
		elseif h_select == v2 && c_select == v2
			println("You both picked $v2. It's a tie!")
			t += 1
		elseif h_select == v2 && c_select == v3
			println("$v3 cuts $v2. You lose!")
			l += 1
		elseif h_select == v2 && c_select == v4
			println("$v4 eats $v2. You lose!")
			l += 1
		elseif h_select == v2 && c_select == v5
			println("$v2 disproves $v5. You win!")
			w += 1
		elseif h_select == v3 && c_select == v1
			println("$v1 smashes $v3. You lose!")
			l += 1
		elseif h_select == v3 && c_select == v2
			println("$v3 cuts $v2. You win!")
			w += 1
		elseif h_select == v3 && c_select == v3
			println("You both picked $v3. It's a tie!")
			t += 1
		elseif h_select == v3 && c_select == v4
			println("$v3 decapitates $v4. You win!")
			w += 1
		elseif h_select == v3 && c_select == v5
			println("$v5 smashes $v3. You lose!")
			l += 1
		elseif h_select == v4 && c_select == v1
			println("$v1 crushes $v4. You lose!")
			l += 1
		elseif h_select == v4 && c_select == v2
			println("$v4 eats $v2. You win!")
			w += 1
		elseif h_select == v4 && c_select == v3
			println("$v3 decapitates $v4. You lose!")
			l += 1
		elseif h_select == v4 && c_select == v4
			println("You both picked $v4. It's a tie!")
			t += 1
		elseif h_select == v4 && c_select == v5
			println("$v4 poisons $v5. You win!")
			w += 1
		elseif h_select == v5 && c_select == v1
			println("$v5 vaporizes $v1. You win!")
			w += 1
		elseif h_select == v5 && c_select == v2
			println("$v2 disproves $v5. You lose!")
			l += 1
		elseif h_select == v5 && c_select == v3
			println("$v5 smashes $v3. You win!")
			w += 1
		elseif h_select == v5 && c_select == v4
			println("$v4 poisons $v5. You lose!")
			l += 1
		elseif h_select == v5 && c_select == v5
			println("You both picked $v5. It's a tie!")
			t += 1
		else
		break
		end
		sleep(0.5)
		end
	end
	if start == 2

		cpu1_w,cpu1_l,cpu1_t = 0,0,0
		cpu2_w,cpu2_l,cpu2_t = 0,0,0

		cpu1_score = (cpu1_w,cpu1_l,cpu1_t)
		cpu2_score = (cpu2_w,cpu2_l,cpu2_t)

		while cpu1_w < 99 && cpu2_w < 99
		print("\nScores: \n$cpu1_score \n$scpu2_score \n")

		v1 = "Rock"
		v2 = "Paper"
		v3 = "Scissors"
		v4 = "Lizard"
		v5 = "Spock"
		cpu1 = "CPU1"
		cpu2 = "CPU2"

		c1_select = rand((v1,v2,v3,v4,v5))
		c2_select = rand((v1,v2,v3,v4,v5))

		print("$v1,")
		sleep(0.5)
		print(" $v2,")
		sleep(0.5)
		print(" $v3,")
		sleep(0.5)
		print(" $v4")
		sleep(0.5)
		print(" $v5,")
		sleep(0.5)
		print(" Shoot!\n")

		# check computer entries for outcome: Win, Lose, Tie.
		if c2_select == v1 && c1_select == v1
			println("You both picked $v1. It's a tie!")
			cpu1_t += 1
			cpu2_t += 1
		elseif c2_select == v1 && c1_select == v2
			println("$v2 covers $v1. $cpu1 wins!")
			cpu1_w += 1
			cpu2_l += 1
		elseif c2_select == v1 && c1_select == v3
			println("$v1 crushes $v3. $cpu2 wins!")
			cpu1_l += 1
			cpu2_w += 1
		elseif c2_select == v1 && c1_select == v4
			println("$v1 crushes $v4. $cpu2 wins!")
			cpu1_l += 1
			cpu2_w += 1
		elseif c2_select == v1 && c1_select == v5
			println("$v5 vaporizes $v1. $cpu1 wins!")
			cpu1_w += 1
			cpu2_l += 1
		elseif c2_select == v2 && c1_select == v1
			println("$v2 covers $v1. $cpu2 wins!")
			cpu1_l += 1
			cpu2_w += 1
		elseif c2_select == v2 && c1_select == v2
			println("You both picked $v2. It's a tie!")
			cpu1_t += 1
			cpu2_t += 1
		elseif c2_select == v2 && c1_select == v3
			println("$v3 cuts $v2. $cpu1 wins!")
			cpu1_w += 1
			cpu2_l += 1
		elseif c2_select == v2 && c1_select == v4
			println("$v4 eats $v2. $cpu1 wins!")
			cpu1_w += 1
			cpu2_l += 1
		elseif c2_select == v2 && c1_select == v5
			println("$v2 disproves $v5. You win!")
			w += 1
		elseif c2_select == v3 && c1_select == v1
			println("$v1 smashes $v3. $cpu1 wins!")
			cpu1_w += 1
			cpu2_l += 1
		elseif c2_select == v3 && c1_select == v2
			println("$v3 cuts $v2. $cpu2 wins!")
			cpu1_l += 1
			cpu2_w += 1
		elseif c2_select == v3 && c1_select == v3
			println("You both picked $v3. It's a tie!")
			cpu1_t += 1
			cpu2_t += 1
		elseif c2_select == v3 && c1_select == v4
			println("$v3 decapitates $v4. $cpu2 wins!")
			cpu1_l += 1
			cpu2_w += 1
		elseif c2_select == v3 && c1_select == v5
			println("$v5 smashes $v3. $cpu1 wins!")
			cpu1_w += 1
			cpu2_l += 1
		elseif c2_select == v4 && c1_select == v1
			println("$v1 crushes $v4. $cpu1 wins!")
			cpu1_w += 1
			cpu2_l += 1
		elseif c2_select == v4 && c1_select == v2
			println("$v4 eats $v2. $cpu2 wins!")
			cpu1_l += 1
			cpu2_w += 1
		elseif c2_select == v4 && c1_select == v3
			println("$v3 decapitates $v4. $cpu1 wins!")
			cpu1_w += 1
			cpu2_l += 1
		elseif c2_select == v4 && c1_select == v4
			println("You both picked $v4. It's a tie!")
			cpu1_t += 1
			cpu2_t += 1
		elseif c2_select == v4 && c1_select == v5
			println("$v4 poisons $v5. $cpu2 wins!")
			cpu1_l += 1
			cpu2_w += 1
		elseif c2_select == v5 && c1_select == v1
			println("$v5 vaporizes $v1. $cpu2 wins!")
			cpu1_l += 1
			cpu2_w += 1
		elseif c2_select == v5 && c1_select == v2
			println("$v2 disproves $v5. $cpu1 wins!")
			cpu1_w += 1
			cpu2_l += 1
		elseif c2_select == v5 && c1_select == v3
			println("$v5 smashes $v3. $cpu2 wins!")
			cpu1_l += 1
			cpu2_w += 1
		elseif c2_select == v5 && c1_select == v4
			println("$v4 poisons $v5. $cpu1 wins!")
			cpu1_w += 1
			cpu2_l += 1
		elseif c2_select == v5 && c1_select == v5
			println("You both picked $v5. It's a tie!")
			cpu1_t += 1
			cpu2_t += 1
		end
		sleep(0.5)
		end
	end
	if start == 3

		p1_w,p1_l,p1_t = 0,0,0
		p2_w,p2_l,p2_t = 0,0,0

		p1_score = (p1_w,p1_l,p1_t)
		p2_score = (p2_w,p2_l,p2_t)

		while p1_w < 99 && p2_w < 99
		print("\nScores: \n$p1_score \n$p2_score \n")

		v1 = "Rock"
		v2 = "Paper"
		v3 = "Scissors"
		v4 = "Lizard"
		v5 = "Spock"
		p1 = "Player 1"
		p2 = "Player 2"

		print("$p1, please select $v1, $v2, $v3, $v4, or $v5.\n")
		p1_select = string(readline()) |> lowercase |> titlecase
		print("$p2, please select $v1, $v2, $v3, $v4, or $v5.\n")
		p2_select = string(readline()) |> lowercase |> titlecase

		print("$v1,")
		sleep(0.5)
		print(" $v2,")
		sleep(0.5)
		print(" $v3,")
		sleep(0.5)
		print(" $v4")
		sleep(0.5)
		print(" $v5,")
		sleep(0.5)
		print(" Shoot!\n")

		# check computer entries for outcome: Win, Lose, Tie.
		if p2_select == v1 && p1_select == v1
			println("You both picked $v1. It's a tie!")
			p1_t += 1
			p2_t += 1
		elseif p2_select == v1 && p1_select == v2
			println("$v2 covers $v1. $p1 wins!")
			p1_w += 1
			p2_l += 1
		elseif p2_select == v1 && p1_select == v3
			println("$v1 crushes $v3. $p2 wins!")
			p1_l += 1
			p2_w += 1
		elseif p2_select == v1 && p1_select == v4
			println("$v1 crushes $v4. $p2 wins!")
			p1_l += 1
			p2_w += 1
		elseif p2_select == v1 && p1_select == v5
			println("$v5 vaporizes $v1. $p1 wins!")
			p1_w += 1
			p2_l += 1
		elseif p2_select == v2 && p1_select == v1
			println("$v2 covers $v1. $p2 wins!")
			p1_l += 1
			p2_w += 1
		elseif p2_select == v2 && p1_select == v2
			println("You both picked $v2. It's a tie!")
			p1_t += 1
			p2_t += 1
		elseif p2_select == v2 && p1_select == v3
			println("$v3 cuts $v2. $p1 wins!")
			p1_w += 1
			p2_l += 1
		elseif p2_select == v2 && p1_select == v4
			println("$v4 eats $v2. $p1 wins!")
			p1_w += 1
			p2_l += 1
		elseif p2_select == v2 && p1_select == v5
			println("$v2 disproves $v5. You win!")
			w += 1
		elseif p2_select == v3 && p1_select == v1
			println("$v1 smashes $v3. $p1 wins!")
			p1_w += 1
			p2_l += 1
		elseif p2_select == v3 && p1_select == v2
			println("$v3 cuts $v2. $p2 wins!")
			p1_l += 1
			p2_w += 1
		elseif p2_select == v3 && p1_select == v3
			println("You both picked $v3. It's a tie!")
			p1_t += 1
			p2_t += 1
		elseif p2_select == v3 && p1_select == v4
			println("$v3 decapitates $v4. $p2 wins!")
			p1_l += 1
			p2_w += 1
		elseif p2_select == v3 && p1_select == v5
			println("$v5 smashes $v3. $p1 wins!")
			p1_w += 1
			p2_l += 1
		elseif p2_select == v4 && p1_select == v1
			println("$v1 crushes $v4. $p1 wins!")
			p1_w += 1
			p2_l += 1
		elseif p2_select == v4 && p1_select == v2
			println("$v4 eats $v2. $p2 wins!")
			p1_l += 1
			p2_w += 1
		elseif p2_select == v4 && p1_select == v3
			println("$v3 decapitates $v4. $p1 wins!")
			p1_w += 1
			p2_l += 1
		elseif p2_select == v4 && p1_select == v4
			println("You both picked $v4. It's a tie!")
			p1_t += 1
			p2_t += 1
		elseif p2_select == v4 && p1_select == v5
			println("$v4 poisons $v5. $p2 wins!")
			p1_l += 1
			p2_w += 1
		elseif p2_select == v5 && p1_select == v1
			println("$v5 vaporizes $v1. $p2 wins!")
			p1_l += 1
			p2_w += 1
		elseif p2_select == v5 && p1_select == v2
			println("$v2 disproves $v5. $p1 wins!")
			p1_w += 1
			p2_l += 1
		elseif p2_select == v5 && p1_select == v3
			println("$v5 smashes $v3. $p2 wins!")
			p1_l += 1
			p2_w += 1
		elseif p2_select == v5 && p1_select == v4
			println("$v4 poisons $v5. $p1 wins!")
			p1_w += 1
			p2_l += 1
		elseif p2_select == v5 && p1_select == v5
			println("You both picked $v5. It's a tie!")
			p1_t += 1
			p2_t += 1
		else
		break
		end
		sleep(0.5)
		end
	end
end

max_stars_repo_name: slatteew/Julia

max_issues_repo_name: slatteew/Julia

max_forks_repo_name: slatteew/Julia


text:
"""
    transformData(flowframe, method = "asinh", cofactor = 5)

Tansforms FCS data. Currently only asinh

# Arguments:
- `flowframe`: Flowframe containing daFrame per sample
- `method`: transformation method
- `cofactor`: Cofactor for transformation
"""
function transformData(flowframe, method = "asinh", cofactor = 5)
    # loop through every file in dict
    # get the dataframe
    # convert to matrix
    # arcsinh transformation
    # convert back to dataframe
    for (k,v) in flowframe
        fcsDf = flowframe[k]
        colnames = names(fcsDf) # keep the column names
        dMatrix = Matrix(fcsDf)
        dMatrix = asinh.(dMatrix / cofactor)
        ddf = DataFrame(dMatrix)

        names!(ddf, Symbol.(colnames))
        # singleFcs["data"] = ddf
        flowframe[k] = ddf
    end
end

"""
    cleanNames!(mydata)

Replaces problematic characters in column names.
Checks if the column name contains a '-' and transforms it to and '_' and it checks if the name starts with a number.

# Arguments:
- `mydata`: dict fcsRaw or array of string
"""
function cleanNames!(mydata)
    # replace chritical characters
    # put "_" in front of colname in case it starts with a number
    # println(typeof(mydata))
    if mydata isa Dict{Any, Any}
        for (k,v) in mydata
            colnames = names(v)
            for i in eachindex(colnames)
                colnames[i] = Symbol(replace(String(colnames[i]), "-"=>"_"))
                if isnumeric(first(String(colnames[i])))
                    colnames[i] = Symbol("_" * String(colnames[i]))
                end
            end
            names!(v, colnames)
        end
    else
        for j in eachindex(mydata)
            mydata[j] = replace(mydata[j], "-"=>"_")
            if isnumeric(first(mydata[j]))
                mydata[j] = "_" * mydata[j]
            end
        end
    end

end

"""
    createDaFrame(fcsRaw, md, panel)

Creates a daFrame of type struct.
Read in the fcs raw, add sample id, subset the columns and transform

# Arguments:
- `fcsRaw`: raw FCS data
- `md`: Metadata table
- `panel`: Panel table with a column for Lineage Markers and one for Functional Markers
"""
function createDaFrame(fcsRaw, md, panel)

    # extract lineage markers
    lineageMarkers, functionalMarkers = getMarkers(panel)

    transformData(fcsRaw)

    for i in eachindex(md.file_name)
        df = fcsRaw[md.file_name[i]]
        insertcols!(df, 1, sample_id = string(md.sample_id[i]))
    end

    dfall = []
    for (k,v) in fcsRaw
        push!(dfall,v)
    end
    dfall = vcat(dfall...)
    cc = map(Symbol, vcat(lineageMarkers, functionalMarkers))
    # markers can be lineage and functional at tthe same time
    # therefore make cc unique
    unique!(cc)
    push!(cc, :sample_id)
    # reduce the dataset to lineage (and later) functional (state) markers
    dfall = dfall[:, cc]
    daf = daFrame(dfall, md, panel)
end

"""
    getMarkers(panel)

Returns the `lineageMarkers` and `functionalMarkers` on a given panel

# Arguments:
- `panel`: Panel table with a column for Lineage Markers and one for Functional Markers
"""
function getMarkers(panel)

    # extract lineage markers
    lineageMarkers = panel.Antigen[panel.Lineage .== 1, : ]
    functionalMarkers = panel.Antigen[panel.Functional .== 1, : ]

    # lineageMarkers are 2d array,
    # flatten this array by using vec:
    lineageMarkers = vec(lineageMarkers)
    functionalMarkers = vec(functionalMarkers)
    cleanNames!(lineageMarkers)
    cleanNames!(functionalMarkers)

    return lineageMarkers, functionalMarkers

end

"""
    getMetaData(f)
Collect the meta data information in a more user friendly format.

# Arguments:
- `f`: input structure with `.params` and `.data` fields
"""
function getMetaData(f)

    # declarations and initializations
    meta = f.params
    metaKeys = keys(meta)
    channel_properties = []
    defaultValue = "None"

    # determine the number of channels
    pars = parse(Int, strip(join(meta["\$PAR"])))

    # determine the range of channel numbers
    channel_numbers = 1:pars

    # determine the channel properties
    for (key,) in meta
        if key[1:3] == "\$P1"
            if !occursin(key[4], "0123456789")
                push!(channel_properties, key[4:end])
            end
        end
    end

    # define the column names
    column_names = ["\$Pn$p" for p in channel_properties]

    # create a data frame
    df = DataFrame([Vector{Any}(undef, 0) for i = 1:length(column_names)])

    # fill the data frame
    for ch in channel_numbers
        # build first each row of the datatable
        tmpV = []
        for p in channel_properties
            if "\$P$ch$p" in metaKeys
                tmp = meta["\$P$ch$p"]
            else
                tmp = defaultValue
            end
            push!(tmpV, tmp)
        end

        # push the row to the dataframe
        push!(df, tmpV)
    end

    # set the names of the df
    names!(df, Symbol.(column_names))

    return df
end


text:
function find_idvec(id::Array{UInt8, 1}, ids::Array{Array{UInt8, 1}, 1})
  for j in 1:length(ids)
    if ids[j] == id
      return j
    end
  end
  return -1
end

function scanrec!(sid::IO,
  ints::Array{Int64, 1},
  ff::Array{Array{Float64, 1}, 1},
  segs::Array{Array{Int64, 1}, 1},
  ids::Array{Array{UInt8, 1}, 1},
  fs_times::Bool,
  seg_times::Bool,
  v::Integer)

  # ===================================================================
  u16 = getfield(BUF, :uint16_buf)
  flags = getfield(BUF, :flags)

  # Fixed section of data header (48 bytes)
  pos = fastpos(sid)
  rebuffer!(sid)
  seekstart(BUF.dh_buf)
  read!(BUF.dh_buf, BUF.seq)
  read!(BUF.dh_buf, BUF.hdr)
  u16[1]          = read(BUF.dh_buf, UInt16)
  u16[2]          = read(BUF.dh_buf, UInt16)
  hh              = read(BUF.dh_buf, UInt8)
  mm              = read(BUF.dh_buf, UInt8)
  ss              = read(BUF.dh_buf, UInt8)
  skip(BUF.dh_buf, 1)
  u16[3]          = read(BUF.dh_buf, UInt16)
  BUF.n           = read(BUF.dh_buf, UInt16)
  BUF.r1          = read(BUF.dh_buf, Int16)
  BUF.r2          = read(BUF.dh_buf, Int16)
  read!(BUF.dh_buf, flags)
  BUF.tc          = read(BUF.dh_buf, Int32)
  u16[4]          = read(BUF.dh_buf, UInt16)
  u16[5]          = read(BUF.dh_buf, UInt16)

  if getfield(BUF, :swap) == true
    hdrswap!(BUF)
  end

  # ==================================================================
  # Post-read header processing

  # Check for correct byte order
  yy = u16[1]
  jj = u16[2]
  if (jj > 0x0200 || ((jj == 0x0000 || jj == 0x0100) &&
      (yy > 0x0907 || yy < 0x707)) || yy>0x0bb8)
    setfield!(BUF, :swap, !BUF.swap)
    if ((BUF.swap == true) && (BUF.wo == 0x01))
      BUF.xs = true
    end
    hdrswap!(BUF)
  end

  # Number of data points
  n = getfield(BUF, :n)

  # Output
  if v > 2
    println(stdout, String(copy(BUF.seq)), " ", String(copy(BUF.hdr)), ", fs = ", 1.0/BUF.dt, ", n = ", n)
  end
  # ==================================================================
  # Channel handling

  # Index to channel ID in ids
  j = find_idvec(BUF.hdr, ids)
  if j == -1
    append!(ints, zeros(Int64, 4))
    push!(ids, copy(BUF.hdr))
    j = length(ids)

    if fs_times
      push!(ff, Float64[])
    end
    if seg_times
      push!(segs, zeros(Int64, 2))
    end
    ints[4j-1] = one(Int64)

    if v > 0
      printstyled("New channel; SEED header = ", String(copy(BUF.hdr)), "\n", color=:red, bold=true)
    end

  elseif v == 3
    println("j = ", j)
  end

  # Check for fs changes
  if (BUF.r1 != BUF.r1_old) || (BUF.r2 != BUF.r2_old)
    update_dt!(BUF)
    if fs_times
      push!(ff[j], 1.0/BUF.dt)
      push!(ff[j], zero(Float64))
      push!(ff[j], zero(Float64))
    end
    if ints[4j] != 0
      ints[4j-1] += 1
    end
  end

  # Get fs
  fs = 1.0/getfield(BUF, :dt)

  # ===================================================================
  # Parse blockettes
  # u16[5] = position within record
  # u16[6] = number of bytes to skip to next blockette

  nsk = u16[4] - 0x0030
  u16[6] = u16[5] - 0x0030
  nblk = flags[4]
  v > 1 && println(string("Number of Blockettes = ", nblk))
  v > 2 && println(stdout, "Relative position in record = ", u16[5], " B from begin, ", u16[6], " B to next blockette, ", nsk, " B to data")

  @inbounds for i in 0x01:nblk
    fastskip(sid, u16[6])
    bt            = fastread(sid, UInt16)
    u16[6]        = fastread(sid, UInt16)
    if getfield(BUF, :swap) == true
      bt          = bswap(bt)
      u16[6]      = bswap(u16[6])
    end

    # Special handling for certain time, fs corrections
    if bt == 0x0064             # [100] Sample Rate Blockette
      BUF.dt = 1.0 / Float64(BUF.swap ? ntoh(fastread(sid, Float32)) : fastread(sid, Float32))
      setfield!(BUF, :Δ, round(Int64, sμ*BUF.dt))

    # must still be parsed (nearly) in full ... could skip fmt
    elseif bt == 0x03e8         # [1000] Data Only SEED Blockette
      BUF.fmt  = fastread(sid)
      BUF.wo   = fastread(sid)
      lx       = fastread(sid)
      fastskip(sid, 1)

      BUF.nx   = 2^lx
      BUF.xs   = ((BUF.swap == true) && (BUF.wo == 0x01))
    elseif bt == 0x03e9         # [1001] Data Extension Blockette
      fastskip(sid, 1)
      BUF.tc += signed(fastread(sid))
    end

    # Get current position relative to record begin
    u16[5]        = UInt16(fastpos(sid) - pos)

    # Update bytes to next blockette or bytes to data
    if i < nblk
      u16[6]      = u16[6] - u16[5]
    else
      nsk         = u16[4] - u16[5]
    end

    if v > 1
      printstyled(string("Position = ", fastpos(sid), ", Blockette [", bt, "]\n"), color=:light_yellow)
      println(stdout, "Relative position in record = ", u16[5], " B from begin, ", u16[6], " B to next blockette, ", nsk, " B to data")
    end
  end

  # Skip to data section
  if nsk > 0x0000
    fastskip(sid, nsk)
    if (v > 2)
      println("Skipped ", nsk, " B")
      printstyled(string("Position = ", fastpos(sid), "\n"), color=:light_green)
    end
  end

  # ===================================================================
  # Skip data
  nb = getfield(BUF, :nx) - u16[4]
  v > 1 && println("Skipping ", nb, " B")
  fastskip(sid, nb)

  # Check for time correction
  is_tc = flags[2] >> 1 & 0x01
  tc = getfield(BUF, :tc)
  if is_tc == false && tc != zero(Int32)
    δt = Int64(tc)*100
  else
    δt = zero(Int64)
  end

  # Sample rate in μs
  Δ = getfield(BUF, :Δ)

  # Start time of this record
  τ = seed_time(u16, hh, mm, ss, δt)

  # if seg_times, check to update start_time
  if seg_times
    if segs[j][end-1] == 0
      segs[j][end-1] = τ
    end
  end

  # check for gap
  if ints[4j] != 0
    gap = τ - ints[4j] - Δ

    # gap found
    if abs(gap) > div(Δ, 2)
      ints[4j-2] += 1
      if seg_times
        append!(segs[j], zeros(Int64, 2))
        segs[j][end-1] = τ
      end
    end
  elseif seg_times
    segs[j][end-1] = τ
  end

  # Set end time of current segment
  ints[4j] = τ + (n-1)*Δ

  # if seg_times tracked, always update segment end time
  if seg_times
    segs[j][end] = ints[4j]
  end

  v > 2 && printstyled(string("Position = ", fastpos(sid), "\n"), color=:light_green)

  # Append number of samples to ints[4j-3]
  ints[4j-3] = ints[4j-3] + getfield(BUF, :n)

  # Logging if fs_times == true
  if fs_times
    if ff[j][end-1] == 0.0
      ff[j][end-1] = τ*1.0e-6
    end
    ff[j][end] = ints[4j]*1.0e-6
  end

  if v > 2
    println("ints[", 4j-3, "]:ints[", 4j, "] = ", ints[4j-3:4j])
  end

  # Done
  return nothing
end

# These are nearly identical to parse_mseed, parse_mseed_file
function scan_seed!(io::IO,
  ints::Array{Int64, 1},
  ff::Array{Array{Float64, 1}, 1},
  segs::Array{Array{Int64, 1}, 1},
  ids::Array{Array{UInt8, 1}, 1},
  fs_times::Bool,
  seg_times::Bool,
  v::Integer)

  while !eof(io)
    scanrec!(io, ints, ff, segs, ids, fs_times, seg_times, v)
  end

  fill!(getfield(BUF, :hdr_old), zero(UInt8))
  setfield!(BUF, :r1_old, zero(Int16))
  setfield!(BUF, :r2_old, zero(Int16))
  return nothing
end

"""
    soh = scan_seed(fname::String[, KWs])

Scan seed file `fname` and report properties in human-readable string array `soh`.

### General Keywords
* quiet (Bool): `true` to only return compact summary strings (no stdout)
* memmap (Bool): `true` to use memory mapping
* v (Integer): `v > 0` increases scan verbosity

### Output Keywords
These are all Booleans; `false` excludes from scan.
* `npts`: Number of samples per channel (default: `true`)
* `ngaps`: Number of time gaps per channel (default: `true`)
* `nfs`: Number of unique fs values per channel (default: `true`)
* `seg_times`: Exact gap times (default: `false`)
* `fs_times`: Exact times of fs changes (default: `false`)

!!! caution

    Rarely, the number of gaps reported is off-by-one from `read_data`.
"""
function scan_seed(fname::String;
    fs_times::Bool              = false,
   seg_times::Bool              = false,
      memmap::Bool              = false,
         nfs::Bool              = true,
       ngaps::Bool              = true,
        npts::Bool              = true,
       quiet::Bool              = false,
           v::Integer           = 0
           )

#=
TO DO:
* segment time tracking
=#

  # [npts ngaps nfs end_time]
  ints = Int64[]

  # [[fs_1, te_1, fs_2, te_2, ...]_1, ...]
  ff = Array{Array{Float64, 1}, 1}(undef, 0)
  ids = Array{Array{UInt8, 1}, 1}(undef, 0)
  segs = Array{Array{Int64, 1}, 1}(undef, 0)
  io  = memmap ? IOBuffer(Mmap.mmap(fname)) : open(fname, "r")
  fastskip(io, 6)
  c = fastread(io)
  if c in (0x44, 0x52, 0x4d, 0x51)
    seekstart(io)
    scan_seed!(io, ints, ff, segs, ids, fs_times, seg_times, v)
    close(io)
  else
    close(io)
    error("Invalid file type!")
  end

  # ================================================================
  # Screen output
  nc = length(ids)

  # -----------------------------------------------------------------
  # Compact output
  soh = Array{String, 1}(undef, nc)
  for j in 1:nc
    copyto!(BUF.hdr, ids[j])
    update_hdr!(BUF)
    soh[j] = string(BUF.id_str, ", ",
            "nx = ", ints[4j-3], ", ",
            "ngaps = ", ints[4j-2], ", ",
            "nfs = ", ints[4j-1])
  end
  quiet && (return soh)

  # -----------------------------------------------------------------
  # Basic counters
  println("\n", lpad("CHANNEL", 15), " ",
          npts          ? string("| ", lpad("N PTS", 11), " ") : "",
          ngaps         ? "| N GAPS " : "",
          nfs           ? "| N FS " : "",
          "\n", "-"^16,
          npts          ? string("+", "-"^13) : "",
          ngaps         ? string("+", "-"^8)  : "",
          nfs           ? string("+", "-"^5)  : "",
          )

  for j in 1:nc
    copyto!(BUF.hdr, ids[j])
    update_hdr!(BUF)
    println(lpad(BUF.id_str, 15), " ",
            npts  ? string("| ", lpad(ints[4j-3], 11), " ") : "",
            ngaps ? string("| ", lpad(ints[4j-2],  6), " ") : "",
            nfs   ? string("| ", lpad(ints[4j-1],  4), " ") : "",
            )
  end
  println("")

  # -----------------------------------------------------------------
  # Detailed tracking
  if fs_times
    # Detailed tracking strings
    for j in 1:nc
      copyto!(BUF.hdr, ids[j])
      update_hdr!(BUF)
      println(BUF.id_str * " fs tracking:\n\n" *
        lpad("TIME WINDOW", 50) * " | " * "FS\n" * "-"^51 * "+" * "---")
      nf = div(length(ff[j]), 3)
      for i in 1:nf
        println(stdout, rpad(u2d(ff[j][3i-1]), 23),
                      " -- ",
                      rpad(u2d(ff[j][3i]), 23),
                      " | ",
                      repr(ff[j][3i-2], context=:compact=>true)
                )
      end
      println("")
    end
  end

  if seg_times
    # Detailed tracking strings
    for j in 1:nc
      copyto!(BUF.hdr, ids[j])
      update_hdr!(BUF)
      println(BUF.id_str * " segment times:\n\n")
      ns = div(length(segs[j]), 2)
      for i in 1:ns
        println(stdout, rpad(u2d(segs[j][2i-1]*1.0e-6), 23),
                      " -- ",
                      rpad(u2d(segs[j][2i]*1.0e-6), 23)
                )
      end
      println("")
    end
  end

  return soh
end


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# ╔═╡ e61b7190-bb53-11eb-0fc2-c5fb8589c6ab
using Plots, SpectralVis

# ╔═╡ 28b64763-ae5c-434a-b4fb-2aa91da6a2d7
html"""<style> main {max-width: 950px;}"""

# ╔═╡ d8844676-d76d-418d-8311-f9c6409eaa14
md"""

Available ellipse datasets:

* PGN_ellipses, field size **2°** at 48 cd/m² on white surround (24 cd/m², 42°)
* WSS_ellipses, field size **1°** on dark surround, colors shown as **short 0.063 s long flashes**
* FMC_ellipses, field size **2°** at 48 cd/m² on white surround (24 cd/m², 42°)
* WRB_ellipses, field size **2°** at varying luminances on dark surround
* DLM_ellipses, field size **2°** at varying luminances on dark surround
* 12B_ellipses, field size **10°** at varying luminances on white (8.9 cd/m²), **average of 12 observers**
* GF_ellipses, two hexagonal fields of **3°** each (12 cd/m²) on white surround: (6 cd/m², 40°), 32 years old
* AR_ellipses, two hexagonal fields of **3°** each (12 cd/m²) on white surround: (6 cd/m², 40°), 27 years old
* GW_ellipses, two hexagonal fields of **3°** each (12 cd/m²) on white surround: (6 cd/m², 40°), 42 years old
"""

# ╔═╡ 721827c8-f47e-485d-9d36-c06e861ba8ea
begin
	plot(title = "Color Matching Ellipses from various datasets, 3 times enlarged", legend = :topright, xlimits = (0, 0.75), ylimits = (0, 0.75), ratio = 1.0, size = (1200,1200))

	scatter!([0.33333],[0.33333], color = :black, shape = :cross, mswidth = 0.0, markersize = 5.0, label = "white point x, y = (1/3, 1/3)")

	φ = 3
	pgn = draw_ellipses(PGN_ellipses,φ)
	wss = draw_ellipses(WSS_ellipses,φ#=/4.5=#) # factor given by paper
	dlm = draw_ellipses(DLM_ellipses,φ)
	fmc = draw_ellipses(FMC_ellipses,φ)
	b   = draw_ellipses(B_ellipses,φ)
	ar  = draw_ellipses(AR_ellipses,φ)
	gf  = draw_ellipses(GF_ellipses,φ)
	gw  = draw_ellipses(GW_ellipses,φ)
	wrb = draw_ellipses(WRB_ellipses,φ#=*1.6=#) # factor given by paper

	for e in 1:length(b)
		plot!(b[e][1],b[e][2],color = :lightblue3, label = false)
	end
	for e in 1:length(ar)
		plot!(ar[e][1],ar[e][2],color = :blue, label = false)
		plot!(gf[e][1],gf[e][2],color = :violet, label = false)
		plot!(gw[e][1],gw[e][2],color = :green3, label = false)
	end
	for e in 1:length(pgn)
		plot!(pgn[e][1],pgn[e][2], color = :black, label= false)
		plot!(wss[e][1],wss[e][2], color = :grey, style = :dot, label= false)
		plot!(dlm[e][1],dlm[e][2], color = :red, label= false)
		plot!(fmc[e][1],fmc[e][2], color = :gold3, label= false)
		plot!(wrb[e][1],wrb[e][2], color = :orange, markersize = 1.0,label = false)
	end
	plot!([0.0,1.0],[1.0,0.0],color = :black, label = false) # limit to the right
	# labels
	plot!([0],[0],color = :black, label = "PGN, 2° (48 cd/m²) on 42° white (24 cd/m²)")
	plot!([0],[0],color = :lightblue3, label = "12B, 10° at varying luminances on white (8.9 cd/m²)")
	plot!([0],[0],color = :blue, label = "AR, 3°+3° hex (12 cd/m²) on 40° white (6 cd/m²)")
	plot!([0],[0],color = :violet, label = "GF, 3°+3° hex (12 cd/m²) on 40° white (6 cd/m²)")
	plot!([0],[0],color = :green3, label = "GW, 3°+3° hex (12 cd/m²) on 40° white (6 cd/m²)")
	plot!([0],[0],color = :grey, style = :dot, label = "WSS, 1° on dark, 0.063 s flashes")
	plot!([0],[0],color = :red, label = "DLM, 2° on dark surround, varying lumin.")
	plot!([0],[0],color = :gold3, label = "FMC, optimized for fit with PNG, 2° (48 cd/m²) on 42° (24 cd/m²)")
	plot!([0],[0],color = :orange, label = "WRB, 2° on dark surround, varying lumin.")

	plot!()
	#draw_ellipses(WRB_ellipses, :red)
	#draw_ellipses(DLM_ellipses, :cyan)
	#draw_ellipses(B_ellipses, :orange)
	#draw_ellipses(GF_ellipses, :skyblue)

end

# ╔═╡ 5a223002-a072-45da-9a6a-ad4fdcc8aefb
gw

# ╔═╡ 861a38c3-92f6-4b3a-9a19-eca15e5d7018
savefig("Color Matching Ellipses from various datasets, 3 times enlarged")

# ╔═╡ Cell order:
# ╠═e61b7190-bb53-11eb-0fc2-c5fb8589c6ab
# ╠═28b64763-ae5c-434a-b4fb-2aa91da6a2d7
# ╟─d8844676-d76d-418d-8311-f9c6409eaa14
# ╠═721827c8-f47e-485d-9d36-c06e861ba8ea
# ╠═5a223002-a072-45da-9a6a-ad4fdcc8aefb
# ╠═861a38c3-92f6-4b3a-9a19-eca15e5d7018


text:
"Check that vect is correctly specified"
function chkvect(vect::AbstractChar)
    if !(vect == 'N' || vect == 'U' || vect == 'V')
        throw(ArgumentError("vect argument must be 'N' (X is not returned), or 'U' or 'V' (X is returned)"))
    end
    vect
end

for (fname, elty) in ((:dlargv_,:Float64),
                      (:slargv_,:Float32),
                      (:zlargv_,:ComplexF64),
                      (:clargv_,:ComplexF32))
    @eval begin
        #       SUBROUTINE DLARGV( N, X, INCX, Y, INCY, C, INCC )
        #       .. Scalar Arguments ..
        #       INTEGER INCC, INCX, INCY, N
        #       .. Array Arguments ..
        #       DOUBLE PRECISION C( * ), X( * ), Y( * )
        function largv!(n::Int, x::Ptr{$elty}, incx::Int, y::Ptr{$elty}, incy::Int, c::Ptr{$elty}, incc::Int)
            ccall((@blasfunc($fname), liblapack), Nothing, (Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}), n, x, incx, y, incy, c, incc)
        end
    end
end

for (fname, elty) in ((:dlartv_,:Float64),
                      (:slartv_,:Float32),
                      (:zlartv_,:ComplexF64),
                      (:clartv_,:ComplexF32))
    @eval begin
        #       SUBROUTINE DLARTV( N, X, INCX, Y, INCY, C, S, INCC )
        #       .. Scalar Arguments ..
        #       INTEGER INCC, INCX, INCY, N
        #       .. Array Arguments ..
        #       DOUBLE PRECISION C( * ), S( * ), X( * ), Y( * )
        function lartv!(n::Int, x::Ptr{$elty}, incx::Int, y::Ptr{$elty}, incy::Int, c::Ptr{$elty}, s::Ptr{$elty}, incc::Int)
            ccall((@blasfunc($fname), liblapack), Nothing, (Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}), n, x, incx, y, incy, c, s, incc)
        end
    end
end

for (fname, elty) in ((:drot_,:Float64),
                      (:srot_,:Float32),
                      (:zrot_,:ComplexF64),
                      (:crot_,:ComplexF32))
    @eval begin
        #       SUBROUTINE DROT(N,DX,INCX,DY,INCY,C,S)
        #       .. Scalar Arguments ..
        #       DOUBLE PRECISION C,S
        #       INTEGER INCX,INCY,N
        #       .. Array Arguments ..
        #       DOUBLE PRECISION DX(*),DY(*)
        function rot!(n::Int, x::Ptr{$elty}, incx::Int, y::Ptr{$elty}, incy::Int, c::$elty, s::$elty)
            ccall((@blasfunc($fname), libblas), Nothing, (Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ref{$elty}, Ref{$elty}), n, x, incx, y, incy, c, s)
        end
    end
end

for (fname, elty) in ((:dlartg_,:Float64),
                      (:slartg_,:Float32),
                      (:zlartg_,:ComplexF64),
                      (:clartg_,:ComplexF32))
    @eval begin
        #       SUBROUTINE DLARTG( F, G, CS, SN, R )
        #       .. Scalar Arguments ..
        #       DOUBLE PRECISION CS, F, G, R, SN
        function lartg!(f::$elty, g::$elty, cs::Ref{$elty}, sn::Ref{$elty}, r::Ref{$elty})
            ccall((@blasfunc($fname), liblapack), Nothing, (Ref{$elty}, Ref{$elty}, Ref{$elty}, Ref{$elty}, Ref{$elty}), f, g, cs, sn, r)
        end
    end
end

for (fname, elty) in ((:dlar2v_,:Float64),
                      (:slar2v_,:Float32),
                      (:zlar2v_,:ComplexF64),
                      (:clar2v_,:ComplexF32))
    @eval begin
        #       SUBROUTINE DLAR2V( N, X, Y, Z, INCX, C, S, INCC )
        #       .. Scalar Arguments ..
        #       INTEGER INCC, INCX, N
        #       .. Array Arguments ..
        #       DOUBLE PRECISION C( * ), S( * ), X( * ), Y( * ), Z( * )
        function lar2v!(n::Int, x::Ptr{$elty}, y::Ptr{$elty}, z::Ptr{$elty}, incx::Int, c::Ptr{$elty}, s::Ptr{$elty}, incc::Int)
            ccall((@blasfunc($fname), liblapack), Nothing, (Ref{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}), n, x, y, z, incx, c, s, incc)
        end
    end
end

## Symmetric tridiagonalization
for (fname, elty) in ((:dsbtrd_,:Float64),
                      (:ssbtrd_,:Float32))
    @eval begin
        function sbtrd!(vect::Char, uplo::Char,
                        m::Int, k::Int, A::AbstractMatrix{$elty},
                        d::AbstractVector{$elty}, e::AbstractVector{$elty}, Q::AbstractMatrix{$elty},
                        work::AbstractVector{$elty})
            require_one_based_indexing(A)
            chkstride1(A)
            chkuplo(uplo)
            chkvect(vect)
            info = Ref{BlasInt}()
            n    = size(A,2)
            n ≠ m && throw(ArgumentError("Matrix must be square"))
            size(A,1) < k+1 && throw(ArgumentError("Not enough bands"))
            info  = Ref{BlasInt}()
            ccall((@blasfunc($fname), liblapack), Nothing,
                (Ref{UInt8}, Ref{UInt8},
                 Ref{BlasInt}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                 Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ptr{BlasInt}),
                 vect, uplo,
                 n, k, A, max(1,stride(A,2)),
                 d, e, Q, max(1,stride(Q,2)), work, info)
            LAPACK.chklapackerror(info[])
            d, e, Q
        end
    end
end

## Bidiagonalization
for (fname, elty) in ((:dgbbrd_,:Float64),
                      (:sgbbrd_,:Float32))
    @eval begin
        function gbbrd!(vect::Char, m::Int, n::Int, ncc::Int,
                        kl::Int, ku::Int, ab::AbstractMatrix{$elty},
                        d::AbstractVector{$elty}, e::AbstractVector{$elty}, Q::AbstractMatrix{$elty},
                        Pt::AbstractMatrix{$elty}, C::AbstractMatrix{$elty}, work::AbstractVector{$elty})
            info  = Ref{BlasInt}()
            ccall((@blasfunc($fname), liblapack), Nothing,
                (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
                 Ref{BlasInt}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                 Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ref{BlasInt},
                 Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ptr{BlasInt}),
                 vect, m, n, ncc,
                 kl, ku, ab, max(1,stride(ab,2)),
                 d, e, Q, max(1,stride(Q,2)),
                 Pt, max(1,stride(Pt,2)), C, max(1,stride(C,2)), work, info)
            LAPACK.chklapackerror(info[])
            d, e, Q, Pt, C
        end
    end
end

# All the eigenvalues and, optionally, eigenvectors of a real symmetric band matrix A.
for (fname, elty) in ((:dsbev_,:Float64),
                      (:ssbev_,:Float32))
    @eval begin
                # SUBROUTINE       SUBROUTINE DSBEV( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
                #     $                  INFO )
                # CHARACTER          JOBZ, UPLO
                # INTEGER            INFO, KD, LDAB, LDZ, N
                # DOUBLE PRECISION   AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )

        function sbev!(jobz::Char, uplo::Char, n::Int, kd::Int, AB::AbstractMatrix{$elty},
                       w::AbstractVector{$elty}, Z::AbstractMatrix{$elty}, work::AbstractVector{$elty})
            info  = Ref{BlasInt}()
            ccall((@blasfunc($fname), liblapack), Nothing,
                (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{$elty},
                 Ref{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ptr{BlasInt}),
                 jobz, uplo, n, kd, AB, max(stride(AB,2),1), w, Z, max(stride(Z,2),1), work, info)
            LAPACK.chklapackerror(info[])
            w, Z
        end
    end
end

# All the generalized eigenvalues and, optionally, eigenvectors of a real symmetric-definite pencil (A,B).
for (fname, elty) in ((:dsbgv_,:Float64),
                      (:ssbgv_,:Float32))
    @eval begin
                # SUBROUTINE       DSBGV( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, W, Z,
                #     $                  LDZ, WORK, INFO )
                # CHARACTER          JOBZ, UPLO
                # INTEGER            INFO, KA, KB, LDAB, LDBB, LDZ, N
                # DOUBLE PRECISION   AB( LDAB, * ), BB( LDBB, * ), W( * ), WORK( * ), Z( LDZ, * )

        function sbgv!(jobz::Char, uplo::Char, n::Int, ka::Int, kb::Int, AB::AbstractMatrix{$elty},
                       BB::AbstractMatrix{$elty}, w::AbstractVector{$elty}, Z::AbstractMatrix{$elty}, work::AbstractVector{$elty})
            info  = Ref{BlasInt}()
            ccall((@blasfunc($fname), liblapack), Nothing,
                (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
                 Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty},
                 Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ptr{BlasInt}),
                 jobz, uplo, n, ka, kb, AB, max(stride(AB,2),1), BB, max(stride(BB,2),1), w, Z, max(stride(Z,2),1), work, info)
            LAPACK.chklapackerror(info[])
            w, Z
        end
    end
end

# Symmetric/Hermitian Positive Definite banded Cholesky factorization
for (fname, elty) in ((:dpbtrf_,:Float64),
                      (:spbtrf_,:Float32),
                      (:zpbtrf_,:ComplexF64),
                      (:cpbtrf_,:ComplexF32))
    @eval begin
                # SUBROUTINE DPBTRF( UPLO, N, KD, AB, LDAB, INFO )
                # CHARACTER          UPLO
                # INTEGER            INFO, KD, LDAB, N
                # DOUBLE PRECISION   AB( LDAB, * )

        function pbtrf!(uplo::Char, m::Int, kd::Int, A::AbstractMatrix{$elty})
            require_one_based_indexing(A)
            chkstride1(A)
            chkuplo(uplo)
            info = Ref{BlasInt}()
            n    = size(A,2)
            n ≠ m && throw(ArgumentError("Matrix must be square"))
            size(A,1) < kd+1 && throw(ArgumentError("Not enough bands"))
            ccall((@blasfunc($fname), liblapack), Nothing,
                (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                 Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}),
                 uplo, n, kd,
                 A, max(1,stride(A,2)), info)
            LAPACK.chklapackerror(info[])
            A, info[]
        end
    end
end

# Symmetric/Hermitian Positive Definite banded Cholesky solver
for (fname, elty) in ((:dpbtrs_,:Float64),
                      (:spbtrs_,:Float32),
                      (:zpbtrs_,:ComplexF64),
                      (:cpbtrs_,:ComplexF32))
    @eval begin
                # SUBROUTINE DPBTRS( UPLO, N, KD, NRHS, AB, LDAB, B, LDB, INFO )
                # CHARACTER          UPLO
                # INTEGER            INFO, KD, LDAB, LDB, N
                # DOUBLE PRECISION   AB( LDAB, * ), B( LDB, * )

        function pbtrs!(uplo::Char, m::Int, kd::Int, A::AbstractMatrix{$elty},
                        B::AbstractVecOrMat{$elty})
            require_one_based_indexing(A)
            chkstride1(A)
            chkuplo(uplo)
            info = Ref{BlasInt}()
            n    = size(A,2)
            if m != n || m != size(B,1)
                throw(DimensionMismatch("matrix A has dimensions $(size(A)), but right hand side matrix B has dimensions $(size(B))"))
            end
            size(A,1) < kd+1 && throw(ArgumentError("Not enough bands"))
            ccall((@blasfunc($fname), liblapack), Nothing,
                  (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
                   Ptr{$elty}, Ref{BlasInt},
                   Ptr{$elty}, Ref{BlasInt},
                   Ptr{BlasInt}),
                  uplo, n, kd, size(B,2),
                  A, max(1,stride(A,2)),
                  B, max(1,stride(B,2)),
                  info)
            LAPACK.chklapackerror(info[])
            B
        end
    end
end

# Symmetric/Hermitian Positive Definite split-Cholesky factorization
for (fname, elty) in ((:dpbstf_,:Float64),
                      (:spbstf_,:Float32),
                      (:zpbstf_,:ComplexF64),
                      (:cpbstf_,:ComplexF32))
    @eval begin
                # SUBROUTINE DPBSTF( UPLO, N, KD, AB, LDAB, INFO )
                # CHARACTER          UPLO
                # INTEGER            INFO, KD, LDAB, N
                # DOUBLE PRECISION   AB( LDAB, * )

        function pbstf!(uplo::Char, m::Int, kd::Int, A::AbstractMatrix{$elty})
            require_one_based_indexing(A)
            chkstride1(A)
            chkuplo(uplo)
            info = Ref{BlasInt}()
            n    = size(A,2)
            n ≠ m && throw(ArgumentError("Matrix must be square"))
            size(A,1) < kd+1 && throw(ArgumentError("Not enough bands"))
            ccall((@blasfunc($fname), liblapack), Nothing,
                (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                 Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}),
                 uplo, n, kd,
                 A, max(1,stride(A,2)), info)
            LAPACK.chklapackerror(info[])
            A
        end
    end
end

# Convert Symmetric Positive Definite generalized eigenvalue problem
# to a symmetric eigenvalue problem assuming B has been processed by
# a split-Cholesky factorization.
for (fname, elty) in ((:dsbgst_,:Float64),
                      (:ssbgst_,:Float32))
    @eval begin
                # SUBROUTINE DSBGST( VECT, UPLO, N, KA, KB, AB, LDAB,
                # BB, LDBB, X, LDX, WORK, INFO )
                # CHARACTER          UPLO, VECT
                # INTEGER            INFO, KA, KB, LDAB, LDBB, LDX, N
                # DOUBLE PRECISION   AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
                # X( LDX, * )

        function sbgst!(vect::Char, uplo::Char, n::Int, ka::Int, kb::Int,
                         AB::AbstractMatrix{$elty}, BB::AbstractMatrix{$elty},
                         X::AbstractVecOrMat{$elty}, work::AbstractVector{$elty})
            require_one_based_indexing(AB, BB, X, work)
            chkstride1(AB, BB)
            chkuplo(uplo)
            chkvect(vect)
            info  = Ref{BlasInt}()
            ccall((@blasfunc($fname), liblapack), Nothing,
                (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                 Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty},
                 Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty},
                 Ptr{BlasInt}),
                 vect, uplo, n, ka,
                 kb, AB, max(stride(AB,2),1), BB,
                 max(stride(BB,2),1), X, max(1,stride(X,2)), work,
                 info)
            LAPACK.chklapackerror(info[])
            AB
        end
    end
end


text:
const ECOLI_MODEL_BIOMASS_IDER = "Biomass_Ecoli_core_w_GAM"
const ECOLI_MODEL_MAT_FILE = abspath(joinpath(@__DIR__, "..", "..", "data", "ecoli_core_model.mat"))

function ecoli_core_model()
    mat_model = MAT.matread(ECOLI_MODEL_MAT_FILE)["model"]
    return MetNet(mat_model; reshape = true)
end


text:
# instantiate environment
using Pkg; Pkg.instantiate()

using Gadfly
using DataFrames
using Printf
using CSV

# intersection between two circles
function intersection(c₁, r₁, c₂, r₂)
  a = r₁ * r₁
  b = r₂ * r₂
  d = sqrt(2 * (c₁ - c₂) * (c₁ - c₂))

  # trivial cases
  d ≥ r₁ + r₂      && return 0.0
  d ≤ abs(r₂ - r₁) && return π * min(a, b)

  # general case
  C1 = a * acos((d * d + a - b) / (2d*r₁))
  C2 = b * acos((d * d + b - a) / (2d*r₂))
  C3 = (1 / 2) * sqrt((-d + r₂ + r₁) * (d + r₂ - r₁) * (d - r₂ + r₁) * (d + r₂ + r₁))

  C1 + C2 - C3
end

# -----------------
# SHIFT FUNCTIONS
# -----------------
kldiv(δ, τ) = δ^2 + τ^2 - log(τ^4) - 1

function jaccard(c₁, r₁, c₂, r₂)
  I = intersection(c₁, r₁, c₂, r₂)
  U = π * r₁^2 + π * r₂^2 - I

  1 - I / U
end

jaccard(δ, τ) = jaccard(0, 3, 3*√2δ + 0, 3τ)

function novelty(c₁, r₁, c₂, r₂)
  I = intersection(c₁, r₁, c₂, r₂)
  O = π * r₂^2 - I
  B = π * r₂^2
  R = (O - I) / B

  (R + 1) / 2
end

novelty(δ, τ) = novelty(0, 3, 3*√2δ + 0, 3τ)

# covariate shift configuration
shiftconfig(δ, τ) = 2δ ≤ 1 - τ ?
                    "inside" :
                    2δ ≥ 1 + τ ?
                    "outside" :
                    "partial"

# -------------
# MAIN SCRIPT
# -------------

# load results table
df = CSV.read(joinpath(@__DIR__,"gaussian.csv"), missingstring="NaN")

# drop outliers due to numerical instability
df = filter(row -> row[:DRV] ≤ 0.5 && row[:ACTUAL] ≤ 0.5, df)
df = dropmissing(df)

# shift functions
df[!,:KLDivergence]    = kldiv.(df[!,:δ], df[!,:τ])
df[!,:JaccardDistance] = jaccard.(df[!,:δ], df[!,:τ])
df[!,:NoveltyFactor]   = novelty.(df[!,:δ], df[!,:τ])

# shift configuration
df[!,:config] = shiftconfig.(df[!,:δ], df[!,:τ])

# correlation length as factor
df[!,:rfactor] = map(df[!,:r]) do r
  @sprintf "r=%.1f" r
end

# filter "inside" configuration
ff = filter(row -> row[:config] == "inside", df)

# set plotting theme
theme = Gadfly.get_theme(Val(:default))
Gadfly.push_theme(theme)
theme = style(point_size=2.5px,
              key_position=:top,
              default_color=colorant"black")
colors = ("#1b9e77","#7570b3","#d95f02")

# generalization error vs. shift function
Gadfly.with_theme(theme) do
  xcols = (:KLDivergence,:JaccardDistance,:NoveltyFactor)
  set_default_plot_size(24cm, 14cm)
  p = plot(df, x=Col.value(xcols...), y=:ACTUAL,
       xgroup=Col.index(xcols...), ygroup=:MODEL,
       color=:config,
       Guide.xlabel("Covariate shift"),
       Guide.ylabel("Error by models"),
       Guide.title("Error vs. covariate shift"),
       Guide.colorkey(title="Configuration"),
       Scale.color_discrete_manual(colors...),
       Geom.subplot_grid(layer(Geom.point), free_x_axis=true,
                         layer(yintercept=[0.0], Geom.hline(color="gray", style=:dash)),
                         layer(yintercept=[0.5], Geom.hline(color="gray", style=:dash)),
                         Coord.cartesian(xmin=0.0,ymax=0.5),
                         Guide.ylabel(orientation=:vertical)))
  p |> SVG(joinpath(@__DIR__,"gaussian-plot1.svg"))
  p
end

# generalization error by different methods
Gadfly.with_theme(theme) do
  ycols = (:CV,:BCV,:DRV,:ACTUAL)
  set_default_plot_size(24cm, 18cm)
  p1 = plot(df, x=:NoveltyFactor, y=Col.value(ycols...),
       xgroup=Col.index(ycols...), color=:rfactor,
       Guide.xlabel("Covariate shift by methods"),
       Guide.ylabel("Error"),
       Guide.title("Error vs. covariate shift"),
       Guide.colorkey(title="Correlation length"),
       Scale.color_discrete_manual(colors...),
       Geom.subplot_grid(layer(Geom.point),
                         layer(yintercept=[0.0], Geom.hline(color="gray", style=:dash)),
                         layer(yintercept=[0.5], Geom.hline(color="gray", style=:dash)),
                         Coord.cartesian(xmin=0.0,xmax=1.0,ymax=0.5)))
  p2 = plot(ff, x=:rfactor, y=Col.value(ycols...),
            xgroup=Col.index(ycols...), color=:rfactor,
            Guide.xlabel("Correlation length by methods"),
            Guide.ylabel("Error"),
            Guide.title("Error vs. correlation length for inside configuration"),
            Guide.colorkey(title="Correlation length"),
            Scale.color_discrete_manual(colors...),
            Geom.subplot_grid(Geom.boxplot))
  p = vstack(p1, p2)
  p |> SVG(joinpath(@__DIR__,"gaussian-plot2.svg"))
  p
end

# generalization error by different correlation lengths
theme2 = style(point_size=1.8px,
               key_position=:top,
               default_color=colorant"black")
Gadfly.with_theme(theme2) do
  set_default_plot_size(24cm, 18cm)
  ycols = (:CV,:BCV,:DRV,:ACTUAL)
  p1 = plot(df, x=:NoveltyFactor, y=Col.value(ycols...),
       xgroup=:rfactor, color=Col.index(ycols...),
       Guide.xlabel("Covariate shift by correlation length"),
       Guide.ylabel("Error"),
       Guide.title("Error vs. covariate shift"),
       Guide.colorkey(title="Method"),
       Scale.color_discrete_manual(colors...),
       Geom.subplot_grid(layer(Geom.point), layer(Geom.line, Stat.smooth(smoothing=1.0)),
                         layer(yintercept=[0.0], Geom.hline(color="gray", style=:dash)),
                         layer(yintercept=[0.5], Geom.hline(color="gray", style=:dash)),
                         Coord.cartesian(xmin=0.0,xmax=1.0,ymax=0.5)))
  xcols = (:CV,:BCV,:DRV)
  p2 = plot(ff, x=Col.value(xcols...), y=:ACTUAL, xgroup=Col.index(xcols...),
            Guide.xlabel("Estimated Error"), Guide.ylabel("Actual Error"),
            Guide.title("Q-Q plot for inside configuration"),
            Geom.subplot_grid(layer(Geom.point,Stat.qq),
                              layer(Geom.abline(color="grey",style=:dot))))
  p = vstack(p1, p2)
  p |> SVG(joinpath(@__DIR__,"gaussian-plot3.svg"))
  p
end

max_issues_repo_name: IBM/geostats-gen-error

max_forks_repo_name: IBM/geostats-gen-error


text:
## BEST BASIS TYPES
"""
    BestBasisType

Abstract type for best basis. Current available types are:
- [`LSDB`](@ref)
- [`JBB`](@ref)
- [`BB`](@ref)
- [`SIBB`](@ref)
"""
abstract type BestBasisType end

"""
    LSDB([; cost, redundant])

Least Statistically Dependent Basis (LSDB).

# Keyword Arguments
- `cost::LSDBCost`: (Default: `DifferentialEntropyCost()`) Cost function for LSDB.
- `redundant::Bool`: (Default: `false`) Whether the performed wavelet transform is
  redundant. Set `redundant=true` when running LSDB with redundant wavelet transforms such
  as SWT or ACWT.

**See also:** [`BestBasisType`](@ref), [`JBB`](@ref), [`BB`](@ref), [`SIBB`](@ref)
"""
@with_kw struct LSDB <: BestBasisType
    cost::LSDBCost = DifferentialEntropyCost()
    redundant::Bool = false
end

"""
    JBB([; cost, redundant])

Joint Best Basis (JBB).

# Keyword Arguments
- `cost::JBBCost`: (Default: `LoglpCost(2)`) Cost function for JBB.
- `redundant::Bool`: (Default: `false`) Whether the performed wavelet transform is
  redundant. Set `redundant=true` when running LSDB with redundant wavelet transforms such
  as SWT or ACWT.

**See also:** [`BestBasisType`](@ref), [`LSDB`](@ref), [`BB`](@ref), [`SIBB`](@ref)
"""
@with_kw struct JBB <: BestBasisType    # Joint Best Basis
    cost::JBBCost = LoglpCost(2)
    redundant::Bool = false
end

"""
    BB([; cost, redundant])

Standard Best Basis (BB).

# Keyword Arguments
- `cost::BBCost`: (Default: `ShannonEntropyCost()`) Cost function for BB.
- `redundant::Bool`: (Default: `false`) Whether the performed wavelet transform is
  redundant. Set `redundant=true` when running LSDB with redundant wavelet transforms such
  as SWT or ACWT.

**See also:** [`BestBasisType`](@ref), [`LSDB`](@ref), [`JBB`](@ref), [`SIBB`](@ref)
"""
@with_kw struct BB <: BestBasisType     # Individual Best Basis
    cost::BBCost = ShannonEntropyCost()
    redundant::Bool = false
end

"""
    SIBB([; cost])

Shift Invariant Best Basis (SIBB).

# Keyword Arguments
- `cost::BBCost`: (Default: `ShannonEntropyCost()`) Cost function for SIBB.

**See also:** [`BestBasisType`](@ref), [`LSDB`](@ref), [`JBB`](@ref),
    [`BB`](@ref)
"""
@with_kw struct SIBB <: BestBasisType   # Shift invariant best basis
    cost::BBCost = ShannonEntropyCost()
end

## ----- TREE COST -----
# LSDB for 1D signals
"""
    tree_costs(X, method)

Returns the cost of each node in a binary tree in order to find the best basis.

# Arguments
- `X::AbstractArray{T} where T<:AbstractFloat`: A set of decomposed signals, of sizes
  `(n,L,k)` for 1D signals or `(n,m,L,k)` for 2D signals, where:
    - `n`: Length of signal (1D) or vertical length of signal (2D).
    - `m`: Horizontal length of signal (2D).
    - `L`: Number of decomposition levels plus 1 (for standard wavelet decomposition) or
        number of nodes in the tree (for redundant transforms such as ACWT and SWT).
    - `k`: Number of signals.
- `method::BestBasisType`: Type of best basis, ie. `BB()`, `JBB()` or `LSDB()`.

!!! note
    For standard best basis (`BB()`), only one signal is processed each time, and therefore
    the inputs `X` should have dimensions `(n,L)` or `(n,m,L)` instead.

# Returns
- `Vector{T}`: A vector containing the costs at each node.

# Examples
```julia
using Wavelets, WaveletsExt

X = generatesignals(:heavisine, 6) |> x -> duplicatesignals(x, 5, 2, true)
wt = wavelet(WT.db4)
Xw = wpdall(X, wt)

tree_costs(Xw, JBB())
tree_costs(Xw, LSDB())
```

**See also:** [`bestbasistree`](@ref), [`bestbasis_treeselection`](@ref)
"""
function tree_costs(X::AbstractArray{T,3}, method::LSDB) where T<:AbstractFloat
    n, L, _ = size(X)

    if method.redundant
        costs = Vector{T}(undef, L)
        for i in eachindex(costs)
            j = getdepth(i,:binary)
            costs[i] = coefcost(X[:,i,:], method.cost) / (1<<j)
        end
    else
        costs = Vector{T}(undef, gettreelength(1<<L))
        i = 1
        for d in 0:(L-1)
            n₀ = nodelength(n, d)
            for node in 0:(1<<d-1)
                rng = (node*n₀+1):((node+1)*n₀)
                costs[i] = coefcost(X[rng, d+1, :], method.cost)
                i += 1
            end
        end
    end
    return costs
end
# LSDB for 2D signals
function tree_costs(X::AbstractArray{T,4}, method::LSDB) where T<:AbstractFloat
    n, m, L, _ = size(X)

    if method.redundant
        costs = Vector{T}(undef,L)
        for i in eachindex(costs)
            d = getdepth(i,:quad)
            costs[i] = coefcost(X[:,:,i,:], method.cost) / (1<<(2*d))
        end
    else
        costs = Vector{T}(undef, gettreelength(1<<L,1<<L))
        for i in eachindex(costs)
            d = getdepth(i,:quad)
            rng₁ = getrowrange(n,i)
            rng₂ = getcolrange(m,i)
            costs[i] = coefcost(X[rng₁,rng₂,d+1,:], method.cost)
        end
    end
    return costs
end

# JBB for 1D signals
function tree_costs(X::AbstractArray{T,3}, method::JBB) where T<:AbstractFloat
    # For each level, compute mean, sum of squares, and variance
    (n, L, N) = size(X)                     # L = levels if wpd, nodes if swpd
    EX = sum(X, dims = 3) ./ N              # calculate E(X)
    EX² = sum(X .^ 2, dims = 3) ./ N        # calculate E(X²)
    VarX = EX² - (EX) .^ 2                  # calculate Var(X)
    VarX = reshape(VarX, (n, L))
    σ = VarX .^ 0.5                         # calculate σ = √Var(X)
    @assert all(σ .>= 0)

    if method.redundant
        costs = Vector{T}(undef, L)
        for i in eachindex(costs)
            j = floor(Integer, log2(i))
            costs[i] = coefcost(σ[:, i], method.cost) / (1<<j)
        end
    else
        costs = Vector{T}(undef, gettreelength(1<<L))
        i = 1   # iterates over the nodes for the costs variable
        for lvl in 0:(L-1)
            n₀ = nodelength(n, lvl)
            for node in 0:(2^lvl-1)
                rng = (node * n₀ + 1):((node + 1) * n₀)
                coef = σ[rng, lvl+1]
                costs[i] = coefcost(coef, method.cost)
                i += 1
            end
        end
    end
    return costs
end
# JBB for 2D signals
function tree_costs(X::AbstractArray{T,4}, method::JBB) where T<:AbstractFloat
    # For each level, compute mean, sum of squares, and variance
    n,m,L,N = size(X)
    EX = sum(X, dims=4) / N
    EX² = sum(X.^2, dims=4) / N
    VarX = EX² - EX.^2 |> y -> dropdims(y, dims=4)
    σ = sqrt.(VarX)
    @assert all(σ .≥ 0)

    if method.redundant
        costs = Vector{T}(undef,L)
        for i in eachindex(costs)
            d = getdepth(i,:quad)
            costs[i] = coefcost(σ[:,:,i], method.cost) / (1<<(2*d))
        end
    else
        costs = Vector{T}(undef,gettreelength(1<<L,1<<L))
        for i in eachindex(costs)
            d = getdepth(i,:quad)
            rng₁ = getrowrange(n,i)
            rng₂ = getcolrange(m,i)
            costs[i] = coefcost(σ[rng₁,rng₂,d+1], method.cost)
        end
    end
    return costs
end

# Standard Best Basis for 1D signals
function tree_costs(X::AbstractArray{T,2}, method::BB) where T<:AbstractFloat
    nrm = norm(X[:,1])
    L = size(X, 2)                      # count of levels if wpd, nodes if swpd
    n = size(X,1)

    if method.redundant                # swpd
        costs = Vector{T}(undef, L)
        for i in eachindex(costs)
            j = floor(Integer, log2(i))
            costs[i] = coefcost(X[:, i], method.cost, nrm) / (1<<j)
        end
    else
        costs = Vector{T}(undef, 2^L - 1)
        i = 1
        for lvl in 0:(L-1)
            n₀ = nodelength(n, lvl)
            for node in 0:(2^lvl-1)
                rng = (node * n₀ + 1):((node + 1) * n₀)
                costs[i] = coefcost(X[rng, lvl+1, :], method.cost, nrm)
                i += 1
            end
        end
    end
    return costs
end
# Standard Best Basis for 2D signals
function tree_costs(X::AbstractArray{T,3}, method::BB) where T<:AbstractFloat
    nrm = norm(X[:,:,1])
    n,m,L = size(X)

    if method.redundant
        costs = Vector{T}(undef,L)
        for i in eachindex(costs)
            d = getdepth(i,:quad)
            costs[i] = coefcost(X[:,:,i], method.cost, nrm) / (1<<(2*d))
        end
    else
        costs = Vector{T}(undef, gettreelength(1<<L,1<<L))
        for i in eachindex(costs)
            d = getdepth(i,:quad)
            rng₁ = getrowrange(n,i)
            rng₂ = getcolrange(m,i)
            costs[i] = coefcost(X[rng₁,rng₂,d+1], method.cost)
        end
    end
    return costs
end

# SIWPD for 1D signals
"""
    tree_costs(y, tree, method)

Computes the cost for each node from the SIWPD decomposition.

# Arguments
- `y::AbstractArray{T,2} where T<:Number`: A SIWPD decomposed signal.
- `tree::AbstractVector{BitVector}`: The full SIWPD tree.
- `method::SIBB`: The `SIBB()` method.

# Returns
- `Vector{Vector{Union{T,Nothing}}}`: SIWPD best basis tree.

!!! warning
    Current implementation works but is unstable, ie. we are still working on better
    syntax/more optimized computations/better data structure.
"""
function tree_costs(y::AbstractArray{T,2}, tree::AbstractVector{BitVector},
        method::SIBB) where T<:Number

    nn = length(tree)
    ns = size(y,1)
    @assert size(y,2) == nn
    tree_costs = Vector{Vector{Union{T,Nothing}}}(undef, nn)
    nrm = norm(y[:,1])

    for i in eachindex(tree)
        level = floor(Int, log2(i))
        len = nodelength(ns, level)
        # number of nodes corresponding to Ω(i,j)
        costs = Vector{Union{AbstractFloat,Nothing}}(nothing, length(tree[i]))
        for j in eachindex(tree[i])
            if tree[i][j]
                shift = j-1                     # current shift
                nstart = shift*len + 1
                nend = (shift+1) * len
                costs[j] = coefcost(y[nstart:nend,i], method.cost, nrm)
            end
        end
        tree_costs[i] = costs
    end
    return tree_costs
end


text:
using Test
using PyCall
import Ghost: play!, Tape, Input
import ONNX: ONNXCtx, push_call!, julia2onnx, onnx2julia

function ort_run(path, args...)
    onnxruntime = pyimport("onnxruntime")
    ort_session = onnxruntime.InferenceSession(path)
    ort_inputs = PyDict(Dict([ort_session.get_inputs()[i].name => args[i] for i=1:length(args)]))
    return ort_session.run(nothing, ort_inputs)
end

function ort_test(tape::Tape, args...)
    mktemp() do path, _
        r1 = play!(tape, args...)
        ONNX.save(path, tape)
        r2 = ort_run(path, julia2onnx.(args)...)[1] |> onnx2julia
        tape2 = ONNX.load(path, args...; exec=true)
        r3 = tape2[tape2.result].val
        @test isapprox(r1, r2)
        @test isapprox(r1, r3)
    end
end

function ort_test(fn::Function, args...)
    tape = Tape(ONNXCtx())
    inp = [push!(tape, Input(arg)) for arg in args]
    res = push_call!(tape, fn, inp...)
    tape.result = res
    ort_test(tape, args...)
end


text:
"""
    FunctionTransform(f)

Take a function or object `f` as an argument which is going to act on each vector individually.
Make sure that `f` is supposed to act on a vector.
For example replace `f(x)=sin(x)` by `f(x)=sin.(x)`
```
    f(x) = abs.(x)
    tr = FunctionTransform(f)
```
"""
struct FunctionTransform{F} <: Transform
    f::F
end

(t::FunctionTransform)(x) = t.f(x)

_map(t::FunctionTransform, x::AbstractVector{<:Real}) = map(t.f, x)
_map(t::FunctionTransform, x::ColVecs) = ColVecs(mapslices(t.f, x.X; dims=1))
_map(t::FunctionTransform, x::RowVecs) = RowVecs(mapslices(t.f, x.X; dims=2))

duplicate(t::FunctionTransform,f) = FunctionTransform(f)

Base.show(io::IO, t::FunctionTransform) = print(io, "Function Transform: ", t.f)


text:
export QNVal,
       QN,
       name,
       val,
       modulus,
       isactive,
       isfermionic

struct QNVal
  name::SmallString
  val::Int
  modulus::Int

  QNVal() = new(SmallString(),0,0)

  function QNVal(name,v::Int,m::Int=1)
    am = abs(m)
    if am > 1
      new(SmallString(name),mod(v,am),m)
    end
    new(SmallString(name),v,m)
  end

end

name(qv::QNVal) = qv.name
val(qv::QNVal) = qv.val
modulus(qv::QNVal) = qv.modulus
isactive(qv::QNVal) = (modulus(qv) != 0)
isfermionic(qv::QNVal) = (modulus(qv) < 0)
Base.:<(qv1::QNVal,qv2::QNVal) = (name(qv1) < name(qv2))
Base.:-(qv::QNVal) =  QNVal(name(qv),-val(qv),modulus(qv))

function pm(qv1::QNVal,qv2::QNVal,fac::Int)
  if name(qv1) != name(qv2)
    error("Cannot add QNVals with different names \"$(name(qv1))\", \"$(name(qv2))\"")
  end
  if modulus(qv1) != modulus(qv2)
    error("QNVals with matching name \"$(name(qv1))\" cannot have different modulus values ")
  end
  m1 = modulus(qv1)
  if m1 == 1
    return QNVal(name(qv1),val(qv1)+fac*val(qv2),1)
  end
  return QNVal(name(qv1),Base.mod(val(qv1)+fac*val(qv2),abs(m1)),m1)
end

Base.:+(qv1::QNVal,qv2::QNVal) = pm(qv1,qv2,+1)
Base.:-(qv1::QNVal,qv2::QNVal) = pm(qv1,qv2,-1)

const ZeroVal = QNVal()

const maxQNs = 4
const QNStorage = SVector{maxQNs,QNVal}
const MQNStorage = MVector{maxQNs,QNVal}

struct QN
  store::QNStorage

  function QN()
    s = QNStorage(ntuple(_ ->ZeroVal,Val(maxQNs)))
    new(s)
  end

  QN(s::QNStorage) = new(s)
end

function QN(qvs...)
  m = MQNStorage(ntuple(_->ZeroVal,Val(maxQNs)))
  for (n,qv) in enumerate(qvs)
    m[n] = QNVal(qv...)
  end
  Nvals = length(qvs)
  sort!(m[1:Nvals];by=name,alg=InsertionSort)
  for n=1:(length(qvs)-1)
    if name(m[n])==name(m[n+1])
      error("Duplicate name \"$(name(m[n]))\" in QN")
    end
  end
  return QN(QNStorage(m))
end

QN(name,val::Int,modulus::Int=1) = QN((name,val,modulus))

Base.getindex(q::QN,n::Int) = getindex(q.store,n)

function val(q::QN,name_)
  sname = SmallString(name_)
  for n=1:maxQNs
    name(q[n]) == sname && return val(q[n])
  end
  return 0
end

function modulus(q::QN,name_)
  sname = SmallString(name_)
  for n=1:maxQNs
    name(q[n]) == sname && return modulus(q[n])
  end
  return 0
end

function combineQNs(a::QN,b::QN,operation)
  !isactive(b[1]) && return a

  ma = MQNStorage(a.store)

  for nb=1:maxQNs
    !isactive(b[nb]) && break
    bname = name(b[nb])
    for na=1:maxQNs
      aname = name(a[na])
      if !isactive(ma[na])
        ma[na] = b[nb]
        break
      elseif name(ma[na]) == bname
        ma[na] = operation(ma[na],b[nb])
        break
      elseif (bname < aname) && (na==1 || bname > name(ma[na-1]))
        for j=maxQNs:-1:(na+1)
          ma[j] = ma[j-1]
        end
        ma[na] = b[nb]
        break
      end
    end
  end
  return QN(QNStorage(ma))
end

function Base.:+(a::QN,b::QN)
  return combineQNs(a,b,+)
end

function Base.:-(a::QN,b::QN)
  return combineQNs(a,b,-)
end

function Base.:(==)(a::QN,b::QN)
  function valsMatch(x::QN,y::QN)
    for xv in x.store
      val(xv) == 0 && continue
      found = false
      for yv in y.store
        name(yv)!=name(xv) && continue
        val(yv)!=val(xv) && return false
        found = true
      end
      found || return false
    end
    return true
  end

  return valsMatch(a,b) && valsMatch(b,a)
end

function Base.:(<)(qa::QN,qb::QN)
  a = 1
  b = 1
  while a<=maxQNs && b<=maxQNs && (isactive(qa[a])||isactive(qb[b]))
    aval = val(qa[a])
    bval = val(qb[b])
    if !isactive(qa[a])
      if 0 == bval
        b += 1
        continue
      end
      return 0 < bval
    elseif !isactive(qb[b])
      if 0 == aval
        a += 1
        continue
      end
      return aval < 0
    else # both are active
      aname = name(qa[a])
      bname = name(qb[b])
      if aname < bname
        if aval == 0
          a += 1
          continue
        end
        return aval < 0
      elseif bname < aname
        if 0 == bval
          b += 1
          continue
        end
        return 0 < bval
      else  # aname == bname
        if aval == bval
          a += 1
          b += 1
          continue
        end
        return aval < bval
      end
    end
  end
  return false
end

function Base.show(io::IO,q::QN)
  print(io,"QN(")
  for n=1:maxQNs
    v = q[n]
    !isactive(v) && break
    n > 1 && print(io,",")
    print(io,"(\"$(name(v))\",$(val(v))")
    if modulus(v) != 1
      print(io,",$(modulus(v)))")
    else
      print(io,")")
    end
  end
  print(io,")")
end


text:
#=
some code taken and modified from PiecewiseLinearOpt.jl
https://github.com/joehuchette/PiecewiseLinearOpt.jl
Copyright (c) 2016: Joey Huchette

We homogenize using a given perspective variable v. σ ≥ 0 and sum(σ) = v are already imposed.
In the "bounded" case, we assume v ∈ [0,1] (MIP formulation can be linear/big-M).
In the "nonneg" case, we assume v ≥ 0 (MIP formulation needs conic constraints).
=#

import PiecewiseLinearOpt
const PLO = PiecewiseLinearOpt

const FloatOrVar = Union{Float64, JuMP.VariableRef}

# SOS2 formulation types
abstract type PWLSOS2 end

struct SOS2 <: PWLSOS2 end
struct CCBounded <: PWLSOS2 end
struct CCNonneg <: PWLSOS2 end
struct LogIBBounded <: PWLSOS2 end
struct LogIBNonneg <: PWLSOS2 end

function PWL_SOS2(::SOS2, model::JuMP.Model, σ::Vector{JuMP.VariableRef}, ::FloatOrVar)
    JuMP.@constraint(model, σ in JuMP.SOS2())
    return
end

function PWL_SOS2(
    ::CCBounded,
    model::JuMP.Model,
    σ::Vector{JuMP.VariableRef},
    binper::FloatOrVar, # binary or 1
)
    n = length(σ)
    y = JuMP.@variable(model, [1:(n - 1)], Bin)
    JuMP.@constraints(model, begin
        sum(y) <= binper
        σ[1] <= y[1]
        [i in 2:(n - 1)], σ[i] <= y[i - 1] + y[i]
        σ[n] <= y[n - 1]
    end)
    return
end

function PWL_SOS2(
    ::CCNonneg,
    model::JuMP.Model,
    σ::Vector{JuMP.VariableRef},
    binper::FloatOrVar, # binary or 1
)
    n = length(σ)
    y = JuMP.@variable(model, [1:(n - 1)], Bin)
    z = JuMP.@variable(model) # TODO or add n variables?
    RSOC = JuMP.RotatedSecondOrderCone()
    JuMP.@constraints(model, begin
        sum(y) <= binper
        [z, y[1], σ[1]] in RSOC
        [i in 2:(n - 1)], [z, y[i - 1] + y[i], σ[i]] in RSOC
        [z, y[n - 1], σ[n]] in RSOC
    end)
    return
end

function _logIB_H(k::Int)
    _H = PLO.reflected_gray_codes(k)
    d = length(_H)
    H = Dict(i => _H[i] for i in 1:d)
    H[0] = H[1]
    H[d + 1] = H[d]
    return H
end

function PWL_SOS2(
    ::LogIBBounded,
    model::JuMP.Model,
    σ::Vector{JuMP.VariableRef},
    binper::FloatOrVar, # binary or 1
)
    n = length(σ)
    k = ceil(Int, log2(n - 1))
    y = JuMP.@variable(model, [1:k], Bin)
    H = _logIB_H(k)
    for j in 1:k
        H1 = JuMP.@expression(model, sum(σ[i] for i in 1:n if H[i - 1][j] == H[i][j] == 1))
        H0 = JuMP.@expression(model, sum(σ[i] for i in 1:n if H[i - 1][j] == H[i][j] == 0))
        JuMP.@constraint(model, H1 <= y[j])
        JuMP.@constraint(model, H0 <= binper - y[j])
    end
    return
end

function PWL_SOS2(
    ::LogIBNonneg,
    model::JuMP.Model,
    σ::Vector{JuMP.VariableRef},
    binper::FloatOrVar,
)
    n = length(σ)
    k = ceil(Int, log2(n - 1))
    y = JuMP.@variable(model, [1:k], Bin)
    z = JuMP.@variable(model) # TODO or add k variables?
    RSOC = JuMP.RotatedSecondOrderCone()
    H = _logIB_H(k)
    for j in 1:k
        H1 = JuMP.@expression(model, sum(σ[i] for i in 1:n if H[i - 1][j] == H[i][j] == 1))
        H0 = JuMP.@expression(model, sum(σ[i] for i in 1:n if H[i - 1][j] == H[i][j] == 0))
        JuMP.@constraint(model, [z, y[j], H1] in RSOC)
        JuMP.@constraint(model, [z, binper - y[j], H0] in RSOC)
    end
    return
end


text:
include("dictionaries.jl")

function generate_genotype(model::GeneModel, instruction::RandomGenome,
                           id::Int)
    """
    returns a geneotype, based on the specified size, connectivity, and
    interaction type.

    interaction_types:
    binary - all interacting genes have interaction of 1 or -1
    normal - normal distribution mean 0, sd 1
    uniform - uniform distribution over [-1,1]
    """
    genes = instruction.genomesize
    grn = rand(interaction_dict[instruction.interactiontype],genes,genes)

    num_unconnected = convert(Int64, floor((1-instruction.connectivity)*genes^2))
    locs = sample(1:genes^2,num_unconnected,replace=false,ordered=true)
    for loc in locs
        grn[id] = 0.0
    end
    model.individuals[id].grn = grn
    return 0
end


text:
module MutableNamedTuples

export MutableNamedTuple

struct MutableNamedTuple{N,T <: Tuple{Vararg{<:Ref}}}
    nt::NamedTuple{N, T}
end

MutableNamedTuple(; kwargs...) = MutableNamedTuple(NamedTuple{keys(kwargs)}(Ref.(values(values(kwargs)))))

function MutableNamedTuple{names}(tuple::Tuple) where names
    MutableNamedTuple(NamedTuple{names}(Ref.(tuple)))
end

Base.keys(::MutableNamedTuple{names}) where {names} = names
Base.values(mnt::MutableNamedTuple) = getindex.(values(getfield(mnt, :nt)))
refvalues(mnt::MutableNamedTuple) = values(getfield(mnt, :nt))

Base.NamedTuple(mnt::MutableNamedTuple) = NamedTuple{keys(mnt)}(values(mnt))
Base.Tuple(mnt::MutableNamedTuple) = values(mnt)

function Base.show(io::IO, mnt::MutableNamedTuple{names}) where {names}
    print(io, "MutableNamedTuple", NamedTuple(mnt))
end

Base.getproperty(mnt::MutableNamedTuple, s::Symbol) = getproperty(NamedTuple(mnt), s)

function Base.setproperty!(mnt::MutableNamedTuple, s::Symbol, x)
    i = findfirst(x -> x === s, keys(mnt))
    refvalues(mnt)[i][] = x
end

Base.length(mnt::MutableNamedTuple) = length(getfield(mnt, :nt))
Base.iterate(mnt::MutableNamedTuple, iter=1) = iterate(NamedTuple(mnt), iter)
Base.firstindex(mnt::MutableNamedTuple) = 1
Base.lastindex(mnt::MutableNamedTuple) = lastindex(NamedTuple(mnt))
Base.getindex(mnt::MutableNamedTuple, i::Int) = getfield(NamedTuple(mnt), i)
Base.getindex(mnt::MutableNamedTuple, i::Symbol) = getfield(NamedTuple(mnt), i)
function Base.indexed_iterate(mnt::MutableNamedTuple, i::Int, state=1)
    Base.indexed_iterate(NamedTuple(mnt), i, state)
end

end # module


text:
import StatsFuns: log2π
import SpecialFunctions: loggamma

@node GammaShapeRate Stochastic [ out, (α, aliases = [ shape ]), (β, aliases = [ rate ]) ]

@average_energy GammaShapeRate (q_out::Any, q_α::PointMass, q_β::Any) = begin
    mean(loggamma, q_α) - mean(q_α) * mean(log, q_β) - (mean(q_α) - 1.0) * mean(log, q_out) + mean(q_β) * mean(q_out)
end

@average_energy GammaShapeRate (q_out::Any, q_α::GammaDistributionsFamily, q_β::Any) = begin
    mean(loggamma, q_α) - mean(q_α) * mean(log, q_β) - (mean(q_α) - 1.0) * mean(log, q_out) + mean(q_β) * mean(q_out)
end


text:
# ------------------------------------------------------------------
# Copyright (c) 2017, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    AbstractPath

A path on a spatial domain of type `D`.
"""
abstract type AbstractPath{D<:AbstractDomain} end

"""
    Base.length(path)

Return the length of a `path`.
"""
Base.length(p::AbstractPath) = npoints(p.domain)

"""
    Base.eltype(::Type{P})

Return element type of path type `P`.
"""
Base.eltype(::Type{P}) where {P<:AbstractPath} = Int

#------------------
# IMPLEMENTATIONS
#------------------
include("paths/simple_path.jl")
include("paths/random_path.jl")


text:
"Generate read (unpack) method for structs."
function generate_read(T::Type)
    fc = fieldcount(T)
    types = [fieldtype(T, i) for i = 1:fc]

    # Create unpack function expression
    function_expression = :(function Base.read(io::IO, t::Type{$T}) end)

    # Create Type call expression and add parameters
    type_expression = :(($T)())
    for t in types
        read_expression = :(read(io, $t))
        append!(type_expression.args, 0)  # dummy with known length
        type_expression.args[end] = read_expression
    end

    # Replace empty function body with Type call
    function_expression.args[2] = type_expression

    eval(function_expression)
end

"Generate write (pack) method for structs."
function generate_write(T::Type)
    # Create pack function expression
    function_expression = :(function Base.write(io::IO, T::$T) end)

    body_expression = quote end
    for t in fieldnames(T)
        append!(body_expression.args, 0)  # dummy with known length
        write_expression = :(write(io, T.$t))
        body_expression.args[end] = write_expression
    end

    # Return nothing at the end
    append!(body_expression.args, 0)  # dummy with known length
    body_expression.args[end] = :(nothing)

    # Replace empty function body with write calls
    function_expression.args[2] = body_expression

    eval(function_expression)
end

function generate_io(T::Type)
    generate_read(T)
    generate_write(T)
end

"""Generate IO expressions macro."""
macro gen_io(typ::Expr)
    T = typ.args[2]
    if isexpr(T, :(<:))
        T = T.args[1]
    end
    if isexpr(T, :curly)
        T = T.args[1]
    end

    ret = Expr(:toplevel, :(Base.@__doc__ $(typ)))
    push!(ret.args, :(generate_io($T)))
    return esc(ret)
end


text:
"""
### Kinetic matrix elements
"""
function kinetic(a::PGBF,b::PGBF)
  return a.norm*b.norm*kinetic(a.expn,a.x,a.y,a.z,a.I,a.J,a.K,
                b.expn,b.x,b.y,b.z,b.I,b.J,b.K)
end

function kinetic(aexpn::Float64,ax::Float64,ay::Float64,az::Float64,aI::Int64,aJ::Int64,aK::Int64,
      bexpn::Float64,bx::Float64,by::Float64,bz::Float64,bI::Int64,bJ::Int64,bK::Int64)
  overlap0 = overlap(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK)
  overlapx1 = overlap(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI+2,bJ,bK)
  overlapy1 = overlap(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ+2,bK)
  overlapz1 = overlap(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK+2)
  overlapx2 = overlap(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI-2,bJ,bK)
  overlapy2 = overlap(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ-2,bK)
  overlapz2 = overlap(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK-2)
  term0 = bexpn*(2*(bI+bJ+bK)+3)*overlap0
  term1 = -2*(bexpn^2)*(overlapx1+overlapy1+overlapz1)
  term2 = -0.5*(bI*(bI-1)*overlapx2+bJ*(bJ-1)*overlapy2+bK*(bK-1)*overlapz2)
  return term0+term1+term2
end

kinetic(a::CGBF,b::CGBF) = contract(kinetic,a,b)

function test_kinetic()
  s = pgbf(1.0)
  c = cgbf(0.0,0.0,0.0)
  push!(c,1,1)
  @assert isapprox(amplitude(c,0,0,0),0.71270547)
  @assert isapprox(kinetic(1.,0.,0.,0.,0,0,0,1.,0.,0.,0.,0,0,0),2.9530518648229536)
  @assert isapprox(kinetic(s,s),1.5)
  @assert isapprox(kinetic(c,c),1.5)
  @printf "test_kinetic is passed\n"
end


text:
# Conditional Inference
export CondBool, CondInt

struct CondBool <: Dist{Bool}
    b1::DistBool
    b2::DistBool
end

function infer(c::CondBool)
    infer(c.b1 & c.b2)/infer(c.b2)
end

struct CondInt <: Dist{Int}
    b1::DistInt
    b2::DistBool
end

function infer(c::CondInt)
    d = c.b1
    b = c.b2
    mb = max_bits(d)
    ans = Vector(undef, 2^mb)
    non_zero_index = 0
    for i=0:2^mb - 1
        a = infer(CondBool(prob_equals(d, i), b))
        if !(a ≈ 0)
            non_zero_index = i+1
        end
        ans[i + 1] = a
    end
    ans[1:non_zero_index]
end