text:
################
# AssemblyType #
################

# this type is used to steer where certain things live and assemble on
# mainly if it lives on CELLs, FACEs or BFACEs

abstract type AssemblyType end

"""
$(TYPEDEF)

causes interpolation at vertices of the grid (only for H1-conforming interpolations)
"""
abstract type AT_NODES <: AssemblyType end  # at nodes (only available for H1 conforming interpolation)

"""
$(TYPEDEF)

causes assembly/interpolation on cells of the grid
"""
abstract type ON_CELLS <: AssemblyType end  # on all cells

"""
$(TYPEDEF)

causes assembly/interpolation on faces of the grid
"""
abstract type ON_FACES <: AssemblyType end  # on all faces

"""
$(TYPEDEF)

causes assembly/interpolation on interior faces of the grid
"""
abstract type ON_IFACES <: ON_FACES end  # on interior faces

"""
$(TYPEDEF)

causes assembly/interpolation on boundary faces of the grid
"""
abstract type ON_BFACES <: AssemblyType end # on boundary faces

"""
$(TYPEDEF)

causes assembly/interpolation on edges of the grid (only in 3D)
"""
abstract type ON_EDGES <: AssemblyType end  # on all edges

"""
$(TYPEDEF)

causes assembly/interpolation on boundary edges of the grid (only in 3D)
"""
abstract type ON_BEDGES <: AssemblyType end # on boundary edges

function Base.show(io::Core.IO, ::Type{ON_CELLS})
    print(io,"ON_CELLS")
end
function Base.show(io::Core.IO, ::Type{ON_FACES})
    print(io,"ON_FACES")
end
function Base.show(io::Core.IO, ::Type{ON_BFACES})
    print(io,"ON_BFACES")
end
function Base.show(io::Core.IO, ::Type{ON_IFACES})
    print(io,"ON_IFACES")
end
function Base.show(io::Core.IO, ::Type{ON_EDGES})
    print(io,"ON_EDGES")
end
function Base.show(io::Core.IO, ::Type{ON_BEDGES})
    print(io,"ON_BEDGES")
end

ItemType4AssemblyType(::Type{ON_CELLS}) = ITEMTYPE_CELL
ItemType4AssemblyType(::Type{<:ON_FACES}) = ITEMTYPE_FACE
ItemType4AssemblyType(::Type{ON_BFACES}) = ITEMTYPE_BFACE
ItemType4AssemblyType(::Type{<:ON_EDGES}) = ITEMTYPE_EDGE
ItemType4AssemblyType(::Type{ON_BEDGES}) = ITEMTYPE_BEDGE

GridComponentNodes4AssemblyType(AT::Type{<:AssemblyType}) = GridComponent4TypeProperty(ItemType4AssemblyType(AT),PROPERTY_NODES)
GridComponentVolumes4AssemblyType(AT::Type{<:AssemblyType}) = GridComponent4TypeProperty(ItemType4AssemblyType(AT),PROPERTY_VOLUME)
GridComponentGeometries4AssemblyType(AT::Type{<:AssemblyType}) = GridComponent4TypeProperty(ItemType4AssemblyType(AT),PROPERTY_GEOMETRY)
GridComponentUniqueGeometries4AssemblyType(AT::Type{<:AssemblyType}) = GridComponent4TypeProperty(ItemType4AssemblyType(AT),PROPERTY_UNIQUEGEOMETRY)
GridComponentRegions4AssemblyType(AT::Type{<:AssemblyType}) = GridComponent4TypeProperty(ItemType4AssemblyType(AT),PROPERTY_REGION)


text:
//test 2 in julia
function a()
	x = 1
	while <x 4 do
		x += 1
	end
	print(x)
end

ext: jl


text:
export GaussianKernel, IMQKernel, RQKernel

abstract type AbstractKernel end

function kernelsum(k::AbstractKernel, x::AbstractMatrix, y::AbstractMatrix, dist::MetricOrFun)
    sum(k(dist(x,y))) / (size(x,2) * size(y,2))
end

function kernelsum(k::AbstractKernel, x::AbstractMatrix{T}, dist::MetricOrFun) where T
    l = size(x,2)
    (sum(k(dist(x,x))) - l*k(T(0)))/(l^2 - l)
end

kernelsum(k::AbstractKernel, x::AbstractVector, dist::MetricOrFun) = zero(eltype(x))

"""
	GaussianKernel(γ)

	implements the standard Gaussian kernel ``exp(-γ * x)
"""
struct GaussianKernel{T} <: AbstractKernel
	γ::T
end

(m::GaussianKernel)(x::Number) = exp(-m.γ * x)
(m::GaussianKernel)(x::AbstractArray) = exp.(-m.γ .* x)

"""
	ImqKernel(c)

Inverse polynomial kernel ``\\frac{C}{C + x}`` used in Tolstikhin, Ilya, et al.
"Wasserstein Auto-Encoders." arXiv preprint arXiv:1711.01558 (2017)

`c` being equal to 2d σ2z, which is expected squared distance between two
samples drawn from p(x).
"""
struct IMQKernel{T} <: AbstractKernel
	c::T
end

(m::IMQKernel)(x::Number) = m.c/ (m.c + x)
(m::IMQKernel)(x::AbstractArray) = m.c ./ (m.c .+ x)

"""
	RQ kernel ``\\frac{C}{C + \\|X - Y\\|}`` used
	in Tolstikhin, Ilya, et al. "Wasserstein Auto-Encoders." arXiv preprint arXiv:1711.01558 (2017)

	`c` being equal to 2d σ2z, which is expected squared distance between two samples drawn from p(x).

"""
struct RQKernel{T} <: AbstractKernel
	α::T
end

(m::RQKernel)(x::Number) = (1 + (x + eps(x))/2m.α)^-m.α
(m::RQKernel)(x::AbstractArray) = (1 .+(x .+ eps(eltype(x))) ./ (2m.α) ).^-m.α


text:
using Test

@testset "./src/strategy_conversions.jl" begin

using BellScenario

@testset "strategy conversions" begin
    @testset "DeterministicStrategy -> Vector{Int64}" begin
        s1 = DeterministicStrategy([1 1 1;0 0 0;0 0 0])
        s2 = DeterministicStrategy([1 0 0;0 1 0;0 0 1])

        norm_v1 = convert(Vector{Int64}, s1, rep="normalized")
        norm_v2 = convert(Vector{Int64}, s2, rep="normalized")
        gen_v = convert(Vector{Int64}, s1, rep="generalized")

        @test norm_v1 isa Vector{Int64}
        @test norm_v1 == [1,0,1,0,1,0]

        @test norm_v2 isa Vector{Int64}
        @test norm_v2 == [1,0,0,1,0,0]

        @test gen_v isa Vector{Int64}
        @test gen_v == [1,0,0,1,0,0,1,0,0]

        @test_throws DomainError convert(Vector{Int64}, s1, rep="non-signaling")
    end

    @testset "LocalSignaling :: Vector{Int64} -> DeterministicStrategy" begin
        norm_v1 = [1,0,0,1,0,0,1,0,0]
        gen_v1 = [1,0,0,0,1,0,0,0,1,0,0,0]

        norm_v2 = [0,0,1,0,0,0,0,0,0]

        scenario = LocalSignaling(3,4,2)

        norm_s1 = convert(DeterministicStrategy, norm_v1, scenario, rep="normalized")
        norm_s2 = convert(DeterministicStrategy, norm_v2, scenario, rep="normalized")

        gen_s1 = convert(DeterministicStrategy, gen_v1, scenario, rep="generalized")

        @test norm_s1 isa DeterministicStrategy
        @test norm_s2 isa DeterministicStrategy
        @test gen_s1 isa DeterministicStrategy

        @test norm_s1 == [1 1 1;0 0 0;0 0 0;0 0 0]
        @test norm_s2 == [0 0 0;0 0 0;1 0 0;0 1 1]
        @test gen_s1 == [1 1 1;0 0 0;0 0 0;0 0 0]

        @test_throws DomainError convert(DeterministicStrategy, norm_v1, scenario, rep="non-signaling")
    end
end

@testset "convert(DeterministicStrategy, ...)" begin
    @testset "chsh vertices non-signaling" begin
        scenario = BipartiteNonSignaling(2,2,2,2)
        vertex = [1,0,1,0,1,0,0,0]

        det_strat = convert(DeterministicStrategy, vertex, scenario)

        @test det_strat isa DeterministicStrategy
        @test det_strat == [1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1]

        strat = convert(Strategy, vertex, scenario)

        @test strat isa Strategy
        @test strat == [1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1]
    end

    @testset "chsh vertices normalized" begin
        scenario = BipartiteNonSignaling(2,2,2,2)
        vertex = [1,0,0,0,1,0,1,0,0,0,1,0]

        det_strat = convert(DeterministicStrategy, vertex, scenario, rep="normalized")

        @test det_strat isa DeterministicStrategy
        @test det_strat == [1 0 1 0;0 1 0 1;0 0 0 0;0 0 0 0]

        strat = convert(Strategy, vertex, scenario, rep="normalized")

        @test strat isa Strategy
        @test strat == [1 0 1 0;0 1 0 1;0 0 0 0;0 0 0 0]
    end

    @testset "chsh vertices generalized" begin
        scenario = BipartiteNonSignaling(2,2,2,2)
        vertex = [1,0,0,0,0,1,0,0,1,0,0,0,0,1,0,0]

        det_strat = convert(DeterministicStrategy, vertex, scenario, rep="generalized")

        @test det_strat isa DeterministicStrategy
        @test det_strat == [1 0 1 0;0 1 0 1;0 0 0 0;0 0 0 0]

        strat = convert(Strategy, vertex, scenario, rep="generalized")

        @test strat isa Strategy
        @test strat == [1 0 1 0;0 1 0 1;0 0 0 0;0 0 0 0]
    end

    @testset "non-signaling vertices" begin
        scenario = BipartiteNonSignaling(2,2,2,2)

        ns_vertex1 = [1/2,1/2,1/2,1/2,0,0,0,1/2]
        ns_strat1 = convert(Strategy, ns_vertex1, scenario)

        @test ns_strat1 isa Strategy
        @test ns_strat1 == 0.5*[0 0 0 1;1 1 1 0;1 1 1 0;0 0 0 1]

        ns_vertex2 = [1/2,1/2,1/2,1/2,1/2,0,1/2,1/2]
        ns_strat2 = convert(Strategy, ns_vertex2, scenario)

        @test ns_strat2 isa Strategy
        @test ns_strat2 == 0.5*[1 0 1 1;0 1 0 0;0 1 0 0;1 0 1 1]
    end

    @testset "DomainError strategy not deterministic" begin
        scenario = BipartiteNonSignaling(2,2,2,2)
        behavior = [0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1]

        @test_throws DomainError convert(DeterministicStrategy, behavior, scenario)
    end

    @testset "DomainError not valid rep" begin
        scenario = BipartiteNonSignaling(2,2,2,2)
        behavior = [0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1]

        @test_throws DomainError convert(Strategy, behavior, scenario, rep="fixed-direction")
    end
end

end


text:
module MotionCaptureJointCalibration

export
    PoseData,
    CalibrationProblem,
    CalibrationResult,
    solve,
    calibration_joints,
    free_joints,
    num_poses,
    num_calibration_params,
    num_markers

using Requires
using StaticArrays
using RigidBodyDynamics
using JuMP
using LinearAlgebra

using MathProgBase: SolverInterface.AbstractMathProgSolver
using RigidBodyDynamics.Graphs: TreePath

include("util.jl")
include("problem.jl")
include("result.jl")
include("deconstruct.jl")
include("residual.jl")
include("solve.jl")
include("synthetic.jl")

function __init__()
    @require DrakeVisualizer="49c7015b-b8db-5bc5-841b-fcb31c578176" begin
        @require RigidBodyTreeInspector="82daab19-8fc9-5c1e-9f69-37d6aaa0269b" begin
            @require Interact="c601a237-2ae4-5e1e-952c-7a85b0c7eef1" begin
                include("visualization.jl")
            end
        end
    end
end

end # module


text:
using RvSpectML
 using Statistics
 using Dates

make_plots = false
include("expres_1_read.jl")
order_list_timeseries = RvSpectML.make_order_list_timeseries(expres_data)
order_list_timeseries = RvSpectML.filter_bad_chunks(order_list_timeseries,verbose=true)
lambda_range_with_data = (min = maximum(d->minimum(d.λ[expres_data[1].metadata[:excalibur_mask]]),expres_data), max = minimum(d->maximum(d.λ[expres_data[1].metadata[:excalibur_mask]]),expres_data) )
RvSpectML.normalize_spectra!(order_list_timeseries,expres_data);

using StaticArrays
using Stheno, TemporalGPs
using Optim
function optimize_gp_fit_to_data(xobs::AA1, yobs::AA2, xpred::AA3;  # ; kernel::Function = matern52_sparse_kernel,
			sigmasq_obs::AA4 = 1e-16*ones(length(xobs))  #= sigmasq_cor::Real=1.0, rho::Real=1.0)  =# ) where {
				T1<:Real, AA1<:AbstractArray{T1,1}, T2<:Real, AA2<:AbstractArray{T2,1}, T3<:Real, AA3<:AbstractArray{T3,1}, T4<:Real, AA4<:AbstractArray{T4,1} }
	# global ncalls += 1
  	println("# size(xobs) = ",size(xobs), "  size(xpred) = ", size(xpred))
	#=
	kobs = make_kernel_data(xobs, kernel, sigmasq_obs=sigmasq_obs, sigmasq_cor=sigmasq_cor, rho=rho)
  	kobs_pred = make_kernel_obs_pred(xobs, xpred, kernel, sigmasq_cor=sigmasq_cor, rho=rho)
  	alpha = kobs \ yobs
	pred_mean = kobs_pred' * alpha
	=#
	function pack(;σ²=σ²,l=l)
		return [log(σ²-1e-6), log(l-1e-6)]
	end
	function unpack(θ)
	    σ² = exp(θ[1]) + 1e-6
	    l = exp(θ[2]) + 1e-6
	    return σ², l
	end

	function neg_log_marg_likel(θ)
		σ², l = unpack(θ)
	    k = σ² * stretch(Matern52(), 1 / l)
	    f_naive = GP(k, GPC())
		#f = to_sde(f_naive)
		f = to_sde(f_naive, SArrayStorage(Float64))
	    return -logpdf(f(xobs, sigmasq_obs), yobs)
	end

	#logyobs = log.(yobs)
	σ_guess = 0.5*std(yobs)
	length_guess = 0.16
	θ_guess = pack(σ²=σ_guess^2,l=length_guess)

#	return neg_log_marg_likel(θ_guess)

	results = Optim.optimize(neg_log_marg_likel, θ_guess, NelderMead()) #, show_trace=true)
	σ²_ml, l_ml = unpack(results.minimizer)
	return σ²_ml, l_ml

	#=k = σ²_ml * stretch(Matern52(), 1 / l_ml);
	f_naive = GP(k, GPC());
	#f = to_sde(f_naive)
	f = to_sde(f_naive, SArrayStorage(Float64))
	fx = f(xobs, sigmasq_obs./yobs.^2)
	f_posterior_ml = posterior(fx, log.(order_list_timeseries.chunk_list[1][ch].flux) )
	marginals(f_posterior_ml(xobs))
	#plot!(plt, f_posterior_ml(xobs); samples=10, color=:green, label="");
	#display(plt);
	=#
	#=
	@time f = GP( 0.05^2 * stretch(Matern52(), 1 / 0.01) , GPC())
	@time f_tgp = to_sde(f)
	@time fx_tgp = f_tgp(order_list_timeseries.chunk_list[1][ch].λ, order_list_timeseries.chunk_list[1][ch].var./order_list_timeseries.chunk_list[1][ch].flux.^2)
	@time f_post_tgp = posterior(fx_tgp, log.(order_list_timeseries.chunk_list[1][ch].flux) )
	@time f_pred_mean_tgp = marginals(f_post(order_list_timeseries.chunk_list[1][ch].λ))
	@time plot(order_list_timeseries.chunk_list[1][ch].λ,mean.(f_pred_mean_tgp); label="", color=:blue)
	scatter!(order_list_timeseries.chunk_list[1][ch].λ,log.(order_list_timeseries.chunk_list[1][ch].flux),markersize=1.5)
	xlims!(5329,5332)
	=#

end

function calc_marginal_gp(gp_param, xobs::AA1, yobs::AA2, xpred::AA3;  # ; kernel::Function = matern52_sparse_kernel,
				sigmasq_obs::AA4 = 1e-16*ones(length(xobs))  #= sigmasq_cor::Real=1.0, rho::Real=1.0)  =# ) where {
					T1<:Real, AA1<:AbstractArray{T1,1}, T2<:Real, AA2<:AbstractArray{T2,1}, T3<:Real, AA3<:AbstractArray{T3,1}, T4<:Real, AA4<:AbstractArray{T4,1} }
	σ², l = gp_param
	k = σ² * stretch(Matern52(), 1 / l)
	f_naive = GP(k, GPC())
	#f = to_sde(f_naive)
	f = to_sde(f_naive, SArrayStorage(Float64))
	fx = f(xobs, sigmasq_obs)
	f_posterior_ml = posterior(fx, (order_list_timeseries.chunk_list[1][ch].flux) )
	f_marg = marginals(f_posterior_ml(xpred))
	return f_marg
	#plot!(plt, f_posterior_ml(xobs); samples=10, color=:green, label="");
	#display(plt)
end

function neg_log_marg_likel(θ)
	σ², l = θ
	k = σ² * stretch(Matern52(), 1 / l)
	f_naive = GP(k, GPC())
	#f = to_sde(f_naive)
	f = to_sde(f_naive, SArrayStorage(Float64))
	xobs = order_list_timeseries.chunk_list[1][ch].λ
	yobs = (order_list_timeseries.chunk_list[1][ch].flux)
	sigmasq_obs = order_list_timeseries.chunk_list[1][ch].var
	return -logpdf(f(xobs, sigmasq_obs), yobs)
end

ch = 10
res = optimize_gp_fit_to_data(log.(order_list_timeseries.chunk_list[1][ch].λ),
  order_list_timeseries.chunk_list[1][ch].flux,log.(order_list_timeseries.chunk_list[1][ch].λ),sigmasq_obs=
  order_list_timeseries.chunk_list[1][ch].var)
  gp_param = res

@time neg_log_marg_likel(gp_param)
gp_param = [std(order_list_timeseries.chunk_list[1][ch].flux), 0.05]

using Plots
ch = 17
 #gp_param = [.1639394167390819*4, 0.14584100679829712/5615*4]
 gp_param = [ 0.4890909216856761, 5.800274590507981e-5 ]  .* 2
 println("gp_param = ", gp_param, " -log L = ", neg_log_marg_likel(gp_param))
 #xplt = range(log.(minimum(order_list_timeseries.chunk_list[1][ch].λ)),stop=log.(maximum(order_list_timeseries.chunk_list[1][ch].λ)),length=2*length(order_list_timeseries.chunk_list[1][ch].λ))
 xplt = log.(order_list_timeseries.chunk_list[1][ch].λ)
 marg = calc_marginal_gp(gp_param,log.(order_list_timeseries.chunk_list[1][ch].λ),
  order_list_timeseries.chunk_list[1][ch].flux,xplt,sigmasq_obs=
  order_list_timeseries.chunk_list[1][ch].var)
 plot(exp.(xplt), mean.(marg), #= yerr=std.(marg).*mean.(marg) =# )
 scatter!(order_list_timeseries.chunk_list[1][ch].λ, order_list_timeseries.chunk_list[1][ch].flux .- 0.0*mean.(marg), markersize=1.5)
 #xlims!(5295,5300)
 #xlims!(5625,5627)
 #ylims!(0.9,1.01)

nxplt = length(xplt)
 σ², l = gp_param
 k = σ² * stretch(Matern52(), 1 / l)
 #f_naive = GP(k, GPC())
 f_naive = GP(k, GPC())
 #f = to_sde(f_naive)
 f = to_sde(f_naive, SArrayStorage(Float64))
 fx = f(log.(order_list_timeseries.chunk_list[1][ch].λ), order_list_timeseries.chunk_list[1][ch].var./order_list_timeseries.chunk_list[1][ch].flux.^2)
 f_posterior_ml = posterior(fx, log.(order_list_timeseries.chunk_list[1][ch].flux) )

m = (mean.(marginals(f_posterior_ml((xplt)))))
 dmdlnl = zeros(size(m))
 dmdlnl[2:end-1] .= (m[3:end].-m[1:end-2])./(xplt[3]-xplt[1])
 d2mdlnl2[1] = dmdlnl[2]
 d2mdlnl2[end] = dmdlnl[end-1]
 d2mdlnl2 = zeros(size(m))
 d2mdlnl2[2:end-1] .= (m[3:end].+m[1:end-2].-2.0.*m[2:end-1])./(xplt[2]-xplt[1]).^2
 d2mdlnl2[1] = d2mdlnl2[2]
 d2mdlnl2[end] = d2mdlnl2[end-1]
 plot(exp.(xplt),exp.(m),color=:black)
 plot!(exp.(xplt),dmdlnl./std(dmdlnl),color=:blue)
 plot!(exp.(xplt),d2mdlnl2./std(d2mdlnl2),color=:red)
 #idx_line_cand = findall(x->x>0.5,d2mdlnl2./std(d2mdlnl2))
 idx_line_cand = findall(x->x>2e8,d2mdlnl2)
 mean_idx = convert_boundaries_to_means(idx_line_cand)
 line_ctrs = exp.(xplt[mean_idx])
 idx_good_line_ctrs = findall(x->abs(x)<1e3,dmdlnl[convert_boundaries_to_means(idx_line_cand)])
 good_line_ctrs = exp.(xplt[mean_idx[idx_good_line_ctrs]])
 scatter!(line_ctrs,zeros(length(line_ctrs)))
 scatter!(good_line_ctrs,0.2.+zeros(length(good_line_ctrs)))
 scatter!(order_list_timeseries.chunk_list[1][ch].λ,order_list_timeseries.chunk_list[1][ch].flux,markersize=1.5)
ylims!(-0.05,1.25)
xlims!(5635.15,5635.65)

log(5635.65/5635.15)
log(5635.65/5635.15)/log(order_list_timeseries.chunk_list[1][ch].λ[2941]/order_list_timeseries.chunk_list[1][ch].λ[2940])

std(d2mdlnl2)
xlims!(5601,5603.5)
xlims!(5588,5594)
xlims!(5594,5600)
xlims!(5600,5606)
ylims!(0.8,1.25)
xlims!(5640,5644)
 #ylims!(-1,2)
xlims!(5625,5627)
xlims!(5585,5595)
xlims!(5591,5593)
xlims!(5615,5620)
ylims!(0.4,0.65)
xlims!(5635,5645)
xlims!(5635,5636)

histogram(log10.(abs.(dmdlnl[convert_boundaries_to_means(idx_line_cand)])),nbins=60)

findall(x->abs(x)<3e3,dmdlnl[convert_boundaries_to_means(idx_line_cand)])
findall(x->abs(x)<1e3,dmdlnl[convert_boundaries_to_means(idx_line_cand)])

function how_many_lines_found(threshold)
	idx_line_cand = findall(x->x>threshold,d2mdlnl2)
	line_ctrs = exp.(xplt[convert_boundaries_to_means(idx_line_cand)])
	length(line_ctrs)
end

thresholds = [1e10, 3e9, 2e9, 1e9, 5e8, 3e8, 2e8, 1.5e8, 1e8, 8e7, 5e7, 3e7, 2e7, 1e7, 3e6, 1e6]
scatter(log10.(thresholds), how_many_lines_found.(thresholds))

function convert_boundaries_to_means(x::AbstractArray{T,1}) where { T<:Real }
	idx_start = 1
	idx_stop = 1
	centers = Int64[]
	for i in 2:length(x)
		if x[i]-x[i-1] != 1
			idx_stop = i-1
			if idx_start!=1
				push!(centers,floor(Int,(x[idx_start]+x[idx_stop])/2 ))
			end
			idx_start = i
		end
	end
	return centers
end

idx_boundaries = findall(x->x!=0.75,idx_line_cand[2:end].-idx_line_cand[1:end-1])

line_ctrs = map(i->(exp(xplt[floor(Int,idx_boundaries[i-1])])+exp(xplt[floor(Int,idx_boundaries[i])]))/2, 2:length(idx_boundaries))

line_ctrs

exp(xplt[idx_boundaries][end])
exp(maximum(xplt))
size(m)

#=

@time f = GP(0.5^2 * stretch(Matern52(), 1 / 0.1) + 0.5^2 * stretch(Matern52(), 1 / 50) , GPC())
ch = 4
 @time fx = f(order_list_timeseries.chunk_list[1][ch].λ, order_list_timeseries.chunk_list[1][ch].var./order_list_timeseries.chunk_list[1][ch].flux.^2)
 @time f_posterior = f | Obs(fx, log.(order_list_timeseries.chunk_list[1][ch].flux) )
 @time f_poster_on_grid = f_posterior(order_list_timeseries.chunk_list[1][ch].λ)
 @time fpred_mean = marginals(f_poster_on_grid);
 @time plot!(f_poster_on_grid; label="", color=:red)
 scatter!(order_list_timeseries.chunk_list[1][ch].λ,log.(order_list_timeseries.chunk_list[1][ch].flux))

using TemporalGPs
#@time f = GP(0.1^2 * stretch(Matern52(), 1 / 0.1) + 0.1^2 * stretch(Matern52(), 1 / 50) , GPC())
std(log.(order_list_timeseries.chunk_list[1][ch].flux))
@time f = GP( 0.05^2 * stretch(Matern52(), 1 / 0.01) , GPC())
@time f_tgp = to_sde(f)
@time fx_tgp = f_tgp(order_list_timeseries.chunk_list[1][ch].λ, order_list_timeseries.chunk_list[1][ch].var./order_list_timeseries.chunk_list[1][ch].flux.^2)
@time f_post_tgp = posterior(fx_tgp, log.(order_list_timeseries.chunk_list[1][ch].flux) )
@time f_pred_mean_tgp = marginals(f_post(order_list_timeseries.chunk_list[1][ch].λ))
@time plot(order_list_timeseries.chunk_list[1][ch].λ,mean.(f_pred_mean_tgp); label="", color=:blue)
scatter!(order_list_timeseries.chunk_list[1][ch].λ,log.(order_list_timeseries.chunk_list[1][ch].flux),markersize=1.5)
xlims!(5329,5332)

@time logpdf(f_post_tgp(order_list_timeseries.chunk_list[1][ch].λ), log.(order_list_timeseries.chunk_list[1][ch].flux))

=#


text:
using Test
using FASTX
using FormatSpecimens
using BioSymbols
import BioGenerics
import BioGenerics.Testing: intempdir
import BioSequences:
    @dna_str,
    @rna_str,
    @aa_str,
    LongDNASeq,
    LongAminoAcidSeq,
    LongSequence,
    AminoAcidAlphabet,
    DNAAlphabet,
    DNA_N,
    DNA_A,
    DNA_G,
    Alphabet

@testset "FASTA" begin
    @testset "Record" begin
        record = FASTA.Record()
        @test !BioGenerics.isfilled(record)
        @test_throws ArgumentError FASTA.identifier(record)
        @test_throws ArgumentError FASTA.description(record)
        @test_throws ArgumentError FASTA.sequence(record)

        record = FASTA.Record(">foo\nACGT\n")
        @test BioGenerics.isfilled(record)
        @test BioGenerics.hasseqname(record)
        @test FASTA.hasidentifier(record)
        @test BioGenerics.seqname(record) == FASTA.identifier(record) == "foo"
        @test !FASTA.hasdescription(record)
        @test FASTA.description(record) === nothing
        @test BioGenerics.hassequence(record)
        @test FASTA.hassequence(record)
        @test FASTA.sequence(record) == dna"ACGT"
        @test collect(FASTA.sequence_iter(DNA, record)) == [DNA_A, DNA_C, DNA_G, DNA_T]
        @test FASTA.sequence(record, 2:3) == LongDNASeq(collect(FASTA.sequence_iter(DNA, record, 2:3))) == dna"CG"
        @test FASTA.sequence(String, record) == "ACGT"
        @test FASTA.sequence(String, record, 2:3) == "CG"

        record1 = FASTA.Record("id", "desc", "AGCT")
        record2 = FASTA.Record("id", "desc", "AGCT")
        @test record1 == record2
        @test hash(record1) == hash(record2)
        @test unique([record1, record1, record2, record2]) == [record1] == [record2]

        record1 = FASTA.Record("id", "AGCT")
        record2 = FASTA.Record("id", "AGCT")
        @test record1 == record2
        @test FASTA.Record() == FASTA.Record()
        @test FASTA.Record() != record1
        @test hash(record1) == hash(record2)
        @test unique([record1, record1, record2, record2]) == [record1] == [record2]

        @test FASTA.Record("id", "AGCT") != FASTA.Record("id2", "AGCT")
        @test FASTA.Record("id", "AGCT") != FASTA.Record("id", "TAGC")
        @test FASTA.Record("id", "desc", "AGCT") != FASTA.Record("id", "AGCT")
        @test FASTA.Record("id", "desc", "AGCT") != FASTA.Record("id", "desc", "TAGC")
        @test FASTA.Record("id", "desc", "AGCT") != FASTA.Record("id", "desc2", "AGCT")

        @test hash(FASTA.Record("id", "AGCT")) != hash(FASTA.Record("id2", "AGCT"))
        @test hash(FASTA.Record("id", "AGCT")) != hash(FASTA.Record("id", "TAGC"))
        @test hash(FASTA.Record("id", "desc", "AGCT")) != hash(FASTA.Record("id", "AGCT"))
        @test hash(FASTA.Record("id", "desc", "AGCT")) != hash(FASTA.Record("id", "desc", "TAGC"))
        @test hash(FASTA.Record("id", "desc", "AGCT")) != hash(FASTA.Record("id", "desc2", "AGCT"))

        record = FASTA.Record("""
        >CYS1_DICDI fragment
        SCWSFSTTGNVEGQHFISQNKL
        VSLSEQNLVDCDHECMEYEGE
        """)
        @test BioGenerics.isfilled(record)
        @test FASTA.identifier(record) == "CYS1_DICDI"
        @test FASTA.description(record) == "fragment"
        @test FASTA.header(record) == "CYS1_DICDI fragment"
        @test FASTA.sequence(record) == aa"SCWSFSTTGNVEGQHFISQNKLVSLSEQNLVDCDHECMEYEGE"
        @test FASTA.sequence(record, 10:15) == aa"NVEGQH"

        # PR 37
        s = ">A \nTAG\n"
        rec = first(iterate(FASTA.Reader(IOBuffer(s))))
        @test isempty(rec.description)
    end

    output = IOBuffer()
    writer = FASTA.Writer(output, 5)
    write(writer, FASTA.Record("seq1", dna"TTA"))
    write(writer, FASTA.Record("seq2", "some description", dna"ACGTNN"))
    flush(writer)
    @test String(take!(output)) == """
    >seq1
    TTA
    >seq2 some description
    ACGTN
    N
    """

    reader = FASTA.Reader(IOBuffer(
    """
    >seqA some description
    QIKDLLVSSSTDLDTTLKMK
    ILELPFASGDLSM
    >seqB
    VLMALGMTDLFIPSANLTG*
    """))
    record = FASTA.Record()
    @test read!(reader, record) === record
    @test FASTA.identifier(record) == "seqA"
    @test FASTA.description(record) == "some description"
    @test FASTA.header(record) == "seqA some description"
    @test FASTA.sequence(record) == aa"QIKDLLVSSSTDLDTTLKMKILELPFASGDLSM"
    @test copyto!(LongAminoAcidSeq(FASTA.seqlen(record)), record) == aa"QIKDLLVSSSTDLDTTLKMKILELPFASGDLSM"
    @test read!(reader, record) === record
    @test FASTA.identifier(record) == "seqB"
    @test !FASTA.hasdescription(record)
    @test FASTA.sequence(record) == aa"VLMALGMTDLFIPSANLTG*"
    @test copyto!(LongAminoAcidSeq(FASTA.seqlen(record)), record) == aa"VLMALGMTDLFIPSANLTG*"
    @test_throws ArgumentError copyto!(LongAminoAcidSeq(10), FASTA.Record())
    @test_throws ArgumentError FASTA.extract(reader, AminoAcidAlphabet(), "seqA", 2:3)

    function test_fasta_parse(filename, valid)
        filepath = joinpath(path_of_format("FASTA"), filename)
        # Reading from a stream
        stream = open(FASTA.Reader, filepath)
        @test eltype(stream) == FASTA.Record
        if valid
            for seqrec in stream end
            @test true  # no error
            @test close(stream) === nothing
        else
            @test_throws Exception begin
                for seqrec in stream end
            end
            return
        end

        # in-place parsing
        stream = open(FASTA.Reader, filepath)
        entry = eltype(stream)()
        while !eof(stream)
            read!(stream, entry)
        end

        # Check round trip
        output = IOBuffer()
        writer = FASTA.Writer(output, width = 60)
        outputB = IOBuffer()
        writerB = FASTA.Writer(outputB, width = -1)
        expected_entries = Any[]
        for seqrec in open(FASTA.Reader, filepath)
            write(writer, seqrec)
            write(writerB, seqrec)
            push!(expected_entries, seqrec)
        end
        flush(writer)
        flush(writerB)

        seekstart(output)
        seekstart(outputB)
        read_entries = FASTA.Record[]
        read_entriesB = FASTA.Record[]
        for seqrec in FASTA.Reader(output)
            push!(read_entries, seqrec)
        end
        for seqrec in FASTA.Reader(outputB)
            push!(read_entriesB, seqrec)
        end
        @test expected_entries == read_entries
        @test expected_entries == read_entriesB
    end

    valid_specimens = list_valid_specimens("FASTA") do specimen
        !hastag(specimen, "comments")
    end
    invalid_specimens = list_invalid_specimens("FASTA")
    for specimen in valid_specimens
        test_fasta_parse(filename(specimen), true)
    end
    for specimen in invalid_specimens
        test_fasta_parse(filename(specimen), false)
    end

    # Test trailing whitespace
    reader = FASTA.Reader(IOBuffer(">A\nTA \n>B C\nCA \nCC \n"))
    @test length(collect(reader)) == 2
    close(reader)

    @testset "Faidx" begin
        # Need to make sure it can handle mixed Unix and Windows newlines
        fastastr = """
        >chr1
        CCACACCACACCCACACACC
        >chr2
        ATGCATGCATGCAT\r
        GCATGCATGCATGC\r
        >chr3
        AAATAGCCCTCATGTACGTCTCCTCCAAGCCCTGTTGTCTCTTACCCGGA
        TGTTCAACCAAAAGCTACTTACTACCTTTATTTTATGTTTACTTTTTATA
        >chr4
        TACTT
        """
        # generated with `samtools faidx`
        faistr = """
        chr1	20	6	20	21
        chr2	28	33	14	16
        chr3	100	71	50	51
        chr4	5	179	5	6
        """
        mktempdir() do dir
            filepath = joinpath(dir, "test.fa")
            write(filepath, fastastr)
            write(filepath * ".fai", faistr)
            open(FASTA.Reader, filepath, index=filepath * ".fai") do reader
                chr3 = reader["chr3"]
                @test FASTA.identifier(chr3) == "chr3"
                @test FASTA.sequence(chr3) == dna"""
                AAATAGCCCTCATGTACGTCTCCTCCAAGCCCTGTTGTCTCTTACCCGGA
                TGTTCAACCAAAAGCTACTTACTACCTTTATTTTATGTTTACTTTTTATA
                """

                seq = FASTA.extract(reader, DNAAlphabet{2}(), "chr2", 10:20)
                @test seq == dna"TGCATGCATGC"
                @test Alphabet(seq) == DNAAlphabet{2}()

                chr2 = reader["chr2"]
                @test FASTA.identifier(chr2) == "chr2"
                @test FASTA.sequence(chr2) == dna"""
                ATGCATGCATGCAT
                GCATGCATGCATGC
                """

                chr4 = reader["chr4"]
                @test FASTA.identifier(chr4) == "chr4"
                @test FASTA.sequence(chr4) == dna"""
                TACTT
                """

                chr1 = reader["chr1"]
                @test FASTA.identifier(chr1) == "chr1"
                @test FASTA.sequence(chr1) == dna"""
                CCACACCACACCCACACACC
                """

                @test_throws ArgumentError reader["chr5"]
            end
        end

        # invalid index
        @test_throws ArgumentError FASTA.Reader(IOBuffer(fastastr), index=π)
    end

    @testset "append" begin
        intempdir() do
            filepath = "test.fa"
            writer = open(FASTA.Writer, filepath)
            write(writer, FASTA.Record("seq1", dna"AAA"))
            close(writer)
            writer = open(FASTA.Writer, filepath, append=true)
            write(writer, FASTA.Record("seq2", dna"CCC"))
            close(writer)
            seqs = open(collect, FASTA.Reader, filepath)
            @test length(seqs) == 2
        end
    end
end

@testset "FASTQ" begin
    @test isa(FASTQ.Record("1", dna"AA", UInt8[10, 11]), FASTQ.Record)
    @test isa(FASTQ.Record("1", "desc.", dna"AA", UInt8[10, 11]), FASTQ.Record)
    @test_throws ArgumentError FASTQ.Record("1", dna"AA", UInt8[10])

    output = IOBuffer()
    writer = FASTQ.Writer(output)
    write(writer, FASTQ.Record("1", dna"AN", UInt8[11, 25]))
    write(writer, FASTQ.Record("2", "high quality", dna"TGA", UInt8[40, 41, 45]))
    flush(writer)
    @test String(take!(output)) == """
    @1
    AN
    +
    ,:
    @2 high quality
    TGA
    +
    IJN
    """

    output = IOBuffer()
    writer = FASTQ.Writer(output, quality_header=true)
    write(writer, FASTQ.Record("1", dna"AN", UInt8[11, 25]))
    write(writer, FASTQ.Record("2", "high quality", dna"TGA", UInt8[40, 41, 45]))
    flush(writer)
    @test String(take!(output)) == """
    @1
    AN
    +1
    ,:
    @2 high quality
    TGA
    +2 high quality
    IJN
    """

    # PR 37
    s = "@A \nTAG\n+\nJJJ\n"
    rec = first(iterate(FASTQ.Reader(IOBuffer(s))))
    @test isempty(rec.description)

    # Test issue 25
    lines = ["@A", "A", "+", "F", "@B", "CA", "+", "CF"]
    io = reader = FASTQ.Reader(IOBuffer(join(lines, "\r\n") * "\r\n"))
    record = FASTQ.Record()
    read!(reader, record)
    read!(reader, record)
    @test eof(reader)

    # Test issue # 48 - allow trailing whitespace
    reader = FASTQ.Reader(IOBuffer("@A\nTA\n+\nFF\n\n\n"))
    @test length(collect(reader)) == 1
    close(reader)

    @testset "Record" begin
        record = FASTQ.Record()
        @test !BioGenerics.isfilled(record)

        record = FASTQ.Record("""
        @SRR1238088.1.1 HWI-ST499:111:D0G94ACXX:1:1101:1173:2105
        AAGCTCATGACCCGTCTTACCTACACCCTTGACGAGATCGAAGGA
        +SRR1238088.1.1 HWI-ST499:111:D0G94ACXX:1:1101:1173:2105
        @BCFFFDFHHHHHJJJIJIJJIJJJJJJJJIJJJJIIIJJJIJJJ
        """)

        seq = dna"AAGCTCATGACCCGTCTTACCTACACCCTTGACGAGATCGAAGGA"

        @test BioGenerics.isfilled(record)
        @test FASTQ.hasidentifier(record) == BioGenerics.hasseqname(record) == true
        @test FASTQ.identifier(record) == BioGenerics.seqname(record) == "SRR1238088.1.1"
        @test FASTQ.hasdescription(record)
        @test FASTQ.description(record) == "HWI-ST499:111:D0G94ACXX:1:1101:1173:2105"
        @test FASTQ.header(record) == FASTQ.identifier(record) * " " * FASTA.description(record)
        @test FASTQ.hassequence(record) == BioGenerics.hassequence(record) == true
        @test FASTQ.sequence(LongDNASeq, record) == seq
        @test LongDNASeq(collect(FASTQ.sequence_iter(DNA, record))) == seq
        @test LongDNASeq(collect(FASTQ.sequence_iter(DNA, record, 3:7))) == dna"GCTCA"
        @test copyto!(LongDNASeq(FASTQ.seqlen(record)), record) == seq
        @test_throws ArgumentError copyto!(LongDNASeq(10), FASTQ.Record())
        @test FASTQ.sequence(record) == BioGenerics.sequence(record) == seq
        @test FASTQ.sequence(String, record) == "AAGCTCATGACCCGTCTTACCTACACCCTTGACGAGATCGAAGGA"
        @test FASTQ.hasquality(record)
        @test FASTQ.quality(record) == b"@BCFFFDFHHHHHJJJIJIJJIJJJJJJJJIJJJJIIIJJJIJJJ" .- 33

        record1 = FASTQ.Record("id", "desc", "AAGCT", collect("@BCFF"))
        record2 = FASTQ.Record("id", "desc", "AAGCT", collect("@BCFF"))
        @test record1 == record2
        @test hash(record1) == hash(record2)
        @test unique([record1, record1, record2, record2]) == [record1] == [record2]

        record1 = FASTQ.Record("id", "AAGCT", collect("@BCFF"))
        record2 = FASTQ.Record("id", "AAGCT", collect("@BCFF"))
        @test record1 == record2
        @test FASTQ.Record() == FASTQ.Record()
        @test FASTQ.Record() != record1
        @test hash(record1) == hash(record2)
        @test unique([record1, record1, record2, record2]) == [record1] == [record2]

        @test FASTQ.Record("id", "AAGCT", collect("@BCFF")) != FASTQ.Record("id2", "AAGCT", collect("@BCFF"))
        @test FASTQ.Record("id", "AAGCT", collect("@BCFF")) != FASTQ.Record("id", "AGTCA", collect("@BCFF"))
        @test FASTQ.Record("id", "AAGCT", collect("@BCFF")) != FASTQ.Record("id", "AAGCT", collect("@BCFG"))
        @test FASTQ.Record("id", "desc", "AAGCT", collect("@BCFF")) != FASTQ.Record("id", "AAGCT", collect("@BCFF"))
        @test FASTQ.Record("id", "desc", "AAGCT", collect("@BCFF")) != FASTQ.Record("id", "desc", "AGTCA", collect("@BCFF"))
        @test FASTQ.Record("id", "desc", "AAGCT", collect("@BCFF")) != FASTQ.Record("id", "desc2", "AAGCT", collect("@BCFF"))

        @test hash(FASTQ.Record("id", "AAGCT", collect("@BCFF"))) != hash(FASTQ.Record("id2", "AAGCT", collect("@BCFF")))
        @test hash(FASTQ.Record("id", "AAGCT", collect("@BCFF"))) != hash(FASTQ.Record("id", "AGTCA", collect("@BCFF")))
        @test hash(FASTQ.Record("id", "AAGCT", collect("@BCFF"))) != hash(FASTQ.Record("id", "AAGCT", collect("@BCFG")))
        @test hash(FASTQ.Record("id", "desc", "AAGCT", collect("@BCFF"))) != hash(FASTQ.Record("id", "AAGCT", collect("@BCFF")))
        @test hash(FASTQ.Record("id", "desc", "AAGCT", collect("@BCFF"))) != hash(FASTQ.Record("id", "desc", "AGTCA", collect("@BCFF")))
        @test hash(FASTQ.Record("id", "desc", "AAGCT", collect("@BCFF"))) != hash(FASTQ.Record("id", "desc2", "AAGCT", collect("@BCFF")))

        record = FASTQ.Record("""
        @SRR1238088.1.1
        AAGCTCATGACCCGTCTTACCTACACCCTTGACGAGATCGAAGGA
        +
        @BCFFFDFHHHHHJJJIJIJJIJJJJJJJJIJJJJIIIJJJIJJJ
        """)
        @test BioGenerics.isfilled(record)
        @test !FASTQ.hasdescription(record)
    end

    function test_records(rs1, rs2)
        if length(rs1) != length(rs2)
            return false
        end
        for (r1, r2) in zip(rs1, rs2)
            if FASTQ.identifier(r1) != FASTQ.identifier(r2) ||
               FASTQ.sequence(r1)   != FASTQ.sequence(r2)   ||
               FASTQ.quality(r1)    != FASTQ.quality(r2)
                return false
            end
        end
        return true
    end

    function test_fastq_parse(filename, valid)
        # Reading from a reader
        filepath = joinpath(path_of_format("FASTQ"), filename)
        reader = open(FASTQ.Reader, filepath)
        @test eltype(reader) == FASTQ.Record
        if valid
            for record in reader end
            @test true  # no error
            @test close(reader) === nothing
        else
            @test_throws Exception begin
                for record in reader end
            end
            return
        end

        # in-place parsing
        reader = open(FASTQ.Reader, filepath)
        record = eltype(reader)()
        try
            while true
                read!(reader, record)
            end
        catch ex
            close(reader)
            if !isa(ex, EOFError)
                rethrow()
            end
        end

        # Check round trip
        output = IOBuffer()
        writer = FASTQ.Writer(output)
        expected_entries = FASTQ.Record[]
        for record in open(FASTQ.Reader, filepath)
            write(writer, record)
            push!(expected_entries, record)
        end
        flush(writer)
        seekstart(output)
        read_entries = FASTQ.Record[]
        for record in FASTQ.Reader(output)
            push!(read_entries, record)
        end

        return test_records(expected_entries, read_entries)
    end

    function test_fastq_fasta_translation(filename)
        ST = LongSequence{DNAAlphabet{4}}

        filepath = joinpath(path_of_format("FASTQ"), filename)
        reader = open(FASTQ.Reader, filepath)
        tmpfile = tempname()
        writer = open(FASTA.Writer, tmpfile)
        transcribe(reader, writer)
        close(reader)
        close(writer)
        fqrecords = open(rdr -> collect(rdr), FASTQ.Reader, filepath)
        farecords = open(rdr -> collect(rdr), FASTA.Reader, tmpfile)

        return all(identifier.(fqrecords) .== identifier.(farecords)) &&
               all(sequence.(ST, fqrecords) .== sequence.(ST, farecords))
    end

    valid_specimens = list_valid_specimens("FASTQ") do specimen
        spec_tags = hastags(specimen) ? sort(tags(specimen)) : String[]
        invalid_tags = sort!(["gaps", "rna", "comments", "linewrap"])
        return isempty(intersect(spec_tags, invalid_tags))
    end
    invalid_specimens = list_invalid_specimens("FASTQ")
    for specimen in valid_specimens
        test_fastq_parse(filename(specimen), true)
        if filename(specimen) != "zero_length.fastq"
            @test test_fastq_fasta_translation(filename(specimen))
        end
    end
    for specimen in invalid_specimens
        test_fastq_parse(filename(specimen), false)
    end

    @testset "invalid quality encoding" begin
        # Sanger full range (note escape characters before '$' and '\')
        record = FASTQ.Record("""
        @FAKE0001 Original version has PHRED scores from 0 to 93 inclusive (in that order)
        ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTAC
        +
        !"#\$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\\]^_`abcdefghijklmnopqrstuvwxyz{|}~
        """)

        # the range is not enough in these encodings
        for encoding in (:solexa, :illumina13, :illumina15)
            @test_throws ErrorException FASTQ.quality(record, encoding)
        end

        # the range is enough in these encodings
        for encoding in (:sanger, :illumina18)
            @test FASTQ.quality(record, encoding) == collect(0:93)
        end
    end

    @testset "fill ambiguous nucleotides" begin
        input = IOBuffer("""
        @seq1
        ACGTNRacgtnr
        +
        BBBB##AAAA##
        """)
        @test FASTQ.sequence(first(FASTQ.Reader(input, fill_ambiguous=nothing))) == dna"ACGTNRACGTNR"
        seekstart(input)
        @test FASTQ.sequence(first(FASTQ.Reader(input, fill_ambiguous=DNA_A)))   == dna"ACGTAAACGTAA"
        seekstart(input)
        @test FASTQ.sequence(first(FASTQ.Reader(input, fill_ambiguous=DNA_G)))   == dna"ACGTGGACGTGG"
        seekstart(input)
        @test FASTQ.sequence(first(FASTQ.Reader(input, fill_ambiguous=DNA_N)))   == dna"ACGTNNACGTNN"
        seekstart(input)
        @test FASTQ.sequence(LongSequence{DNAAlphabet{2}}, first(FASTQ.Reader(input, fill_ambiguous=DNA_A))) == dna"ACGTAAACGTAA"
    end

    @testset "FASTQRead" begin
        record = FASTQ.Record("""
                   @SRR1238088.1.1 HWI-ST499:111:D0G94ACXX:1:1101:1173:2105
                   AAGCTCATGACCCGTCTTACCTACACCCTTGACGAGATCGAAGGA
                   +SRR1238088.1.1 HWI-ST499:111:D0G94ACXX:1:1101:1173:2105
                   @BCFFFDFHHHHHJJJIJIJJIJJJJJJJJIJJJJIIIJJJIJJJ
                   """)
        r1 = FASTQ.Read(record)
        r2 = FASTQ.Read(record,:sanger)
        @test r1.quality == r2.quality
        @test String(r1.sequence) == "AAGCTCATGACCCGTCTTACCTACACCCTTGACGAGATCGAAGGA"
        @test Int64.(r1.quality)[1:5] == [31, 33, 34, 37, 37]
        @test length(r1) == 45
        @test String(sequence(r1)) == "AAGCTCATGACCCGTCTTACCTACACCCTTGACGAGATCGAAGGA"
        @test Int64.(quality(r1))[1:5] == [31, 33, 34, 37, 37]
        @test String(r1[3:6].sequence) == "GCTC"
        @test Int64.(r1[3:6].quality) == [34,37,37,37]
        @test String(r1[3].sequence) == "G"
        @test sprint(show, r1) == "FASTX.FASTQ.FASTQRead{DNAAlphabet{4}}:\n   identifier: SRR1238088.1.1\n  description: HWI-ST499:111:D0G94ACXX:1:1101:1173:2105\n     sequence: AAGCTCATGACCCGTCTTACCTACACCCTTGACGAGATCGAAGGA\n      quality: [31, 33, 34, 37, 37, 37, 35, 37, 39, 39, 39, 39, 39, 41, 41, 41, 40, 41, 40, 41, 41, 40, 41, 41, 41, 41, 41, 41, 41, 41, 40, 41, 41, 41, 41, 40, 40, 40, 41, 41, 41, 40, 41, 41, 41]"
    end

    @testset "Conversion to FASTA" begin
        fqrecord = FASTQ.Record("""
                   @SRR1238088.1.1 HWI-ST499:111:D0G94ACXX:1:1101:1173:2105
                   AAGCTCATGACCCGTCTTACCTACACCCTTGACGAGATCGAAGGA
                   +SRR1238088.1.1 HWI-ST499:111:D0G94ACXX:1:1101:1173:2105
                   @BCFFFDFHHHHHJJJIJIJJIJJJJJJJJIJJJJIIIJJJIJJJ
                   """)
        @test sequence(FASTA.Record(fqrecord)) == sequence(fqrecord)
        @test identifier(FASTA.Record(fqrecord)) == identifier(fqrecord)
    end

end

@testset "Quality scores" begin
    @testset "Decoding base quality scores" begin
        function test_decode(encoding, values, expected)
            result = Array{Int8}(undef, length(expected))
            FASTQ.decode_quality_string!(encoding, values, result, 1, length(result))
            @test result == expected
        end

        test_decode(FASTQ.SANGER_QUAL_ENCODING,
                    UInt8['!', '#', '$', '%', '&', 'I', '~'],
                    Int8[0, 2, 3, 4, 5, 40, 93])

        test_decode(FASTQ.SOLEXA_QUAL_ENCODING,
                    UInt8[';', 'B', 'C', 'D', 'E', 'h', '~'],
                    Int8[-5, 2, 3, 4, 5, 40, 62])

        test_decode(FASTQ.ILLUMINA13_QUAL_ENCODING,
                    UInt8['@', 'B', 'C', 'D', 'E', 'h', '~'],
                    Int8[0, 2, 3, 4, 5, 40, 62])

        test_decode(FASTQ.ILLUMINA15_QUAL_ENCODING,
                    UInt8['C', 'D', 'E', 'h', '~'],
                    Int8[3, 4, 5, 40, 62])

        test_decode(FASTQ.ILLUMINA18_QUAL_ENCODING,
                    UInt8['!', '#', '$', '%', '&', 'I', '~'],
                    Int8[0, 2, 3, 4, 5, 40, 93])
    end

    @testset "Encoding base quality scores" begin
        function test_encode(encoding, values, expected)
            result = Array{UInt8}(undef, length(expected))
            FASTQ.encode_quality_string!(encoding, values, result, 1, length(result))
            @test result == expected
        end

        test_encode(FASTQ.SANGER_QUAL_ENCODING,
                    Int8[0, 2, 3, 4, 5, 40, 93],
                    UInt8['!', '#', '$', '%', '&', 'I', '~'])

        test_encode(FASTQ.SOLEXA_QUAL_ENCODING,
                    Int8[-5, 2, 3, 4, 5, 40, 62],
                    UInt8[';', 'B', 'C', 'D', 'E', 'h', '~'])

        test_encode(FASTQ.ILLUMINA13_QUAL_ENCODING,
                    Int8[0, 2, 3, 4, 5, 40, 62],
                    UInt8['@', 'B', 'C', 'D', 'E', 'h', '~'])

        test_encode(FASTQ.ILLUMINA15_QUAL_ENCODING,
                    Int8[3, 4, 5, 40, 62],
                    UInt8['C', 'D', 'E', 'h', '~'])

        test_encode(FASTQ.ILLUMINA18_QUAL_ENCODING,
                    Int8[0, 2, 3, 4, 5, 40, 93],
                    UInt8['!', '#', '$', '%', '&', 'I', '~'])
    end

end


text:
###############################################################
######################GAMMA####################################
###############################################################

immutable PGamma <: ContinuousUnivariateDistribution
    params::Array{Float64}
end

PGamma() = PGamma([1.])

for op in [:pdf, :logpdf, :cdf, :quantile]
    @eval ($op)(pdist::PGamma, t::Real) = ($op)(Gamma(exp(pdist.params[1]),exp(-pdist.params[1])), t)
end

rand(pdist::PGamma) = rand(Gamma(exp(pdist.params[1]),exp(-pdist.params[1])))

ds(pdist::PGamma, t, int_coefs, i)    =  exp(pdist.params[1])*(log(t)-t) + int_coefs.aux[1]
d²s(pdist::PGamma, t, int_coefs,i, j)  = ds(pdist, t, int_coefs,i)+ int_coefs.aux[2]

function dl!(gl, pdist::PGamma, t, int_coefs)
    gl[2] = exp(pdist.params[1])
    gl[1] = muladd(gl[2],(log(t)-t), int_coefs.aux[1])
    gl[2] *= (t-1)
end

function d²l!(gl, hl, pdist::PGamma, t, int_coefs)
    hl[2,2] = gl[2] = exp(pdist.params[1])
    gl[1] = muladd(gl[2],(log(t)-t), int_coefs.aux[1])
    gl[2] *= (t-1)
    hl[1,1] = gl[1] + int_coefs.aux[2]
    hl[1,2] = hl[2,1] = gl[2]
    hl[2,2] *= -t
end

auxvec(pdist::PGamma) = [exp(pdist.params[1])*(-polygamma(0,exp(pdist.params[1]))+1+pdist.params[1]),
-exp(2*pdist.params[1])*polygamma(1,exp(pdist.params[1]))+exp(pdist.params[1])]

function initialize_aft(S::AbstractVector,X::AbstractArray, pdist::PGamma, N)
    St = [(s.t₀+s.t₁)/2 for s in S]
    valid = isfinite.(St)
    res = glm(X[valid,:], St[valid], Gamma(), LogLink())
    T = 1/GLM.dispersion(res, true)
    return vcat(log(T),coef(res))
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    WeightedCrossValidation(weighter, partitioner;
                            lambda=1.0, loss=Dict())

Weighted cross-validation in which samples are split into folds with
`partitioner` method and are weighted with `weighter` method. Weights
are raised to `lambda` power in `[0,1]`. Optionally, specify `loss`
function from `LossFunctions.jl` for some of the variables.

    WeightedCrossValidation(weighter, k; shuffle=true,
                            lambda=1.0, loss=Dict())

Alternatively, specify the desired number of folds `k` and `shuffle`
options for a `UniformPartitioner`.

## References

* Sugiyama et al. 2006. [Importance-weighted cross-validation for
  covariate shift](https://link.springer.com/chapter/10.1007/11861898_36)
* Sugiyama et al. 2007. [Covariate shift adaptation by importance weighted
  cross validation](http://www.jmlr.org/papers/volume8/sugiyama07a/sugiyama07a.pdf)
"""
struct WeightedCrossValidation{W<:AbstractWeighter,
                               P<:AbstractPartitioner,
                               T<:Real} <: AbstractErrorEstimator
  weighter::W
  partitioner::P
  lambda::T
  loss::Dict{Symbol,SupervisedLoss}

  function WeightedCrossValidation{W,P,T}(weighter, partitioner, lambda, loss) where {W,P,T}
    @assert 0 ≤ lambda ≤ 1 "lambda must lie in [0,1]"
    new(weighter, partitioner, lambda, loss)
  end
end

WeightedCrossValidation(weighter::W, partitioner::P;
                        lambda=1.0, loss=Dict()) where {W<:AbstractWeighter,
                                                        P<:AbstractPartitioner} =
  WeightedCrossValidation{W,P,typeof(lambda)}(weighter, partitioner, lambda, loss)

WeightedCrossValidation(weighter::W, k::Int; shuffle=true,
                        lambda=1.0, loss=Dict()) where {W<:AbstractWeighter} =
  WeightedCrossValidation(weighter, UniformPartitioner(k, shuffle), lambda=lambda, loss=loss)

function Base.error(solver::AbstractLearningSolver,
                    problem::LearningProblem,
                    eestimator::WeightedCrossValidation)
  # retrieve problem info
  sdata = sourcedata(problem)
  ovars = outputvars(task(problem))
  loss  = eestimator.loss
  for var in ovars
    if var ∉ keys(loss)
      loss[var] = defaultloss(sdata[1,var])
    end
  end

  # folds for cross-validation
  folds  = subsets(partition(sdata, eestimator.partitioner))
  nfolds = length(folds)

  solutions = pmap(1:nfolds) do k
    # source and target indices
    sinds = [ind for i in vcat(1:k-1, k+1:nfolds) for ind in folds[i]]
    tinds = folds[k]

    # source and target data
    train = view(sdata, sinds)
    hold  = view(sdata, tinds)

    # setup and solve sub-problem
    subproblem = LearningProblem(train, hold, task(problem))
    solve(subproblem, solver)
  end

  # weight all samples
  weights = weight(sdata, eestimator.weighter) .^ eestimator.lambda

  result = pmap(ovars) do var
    losses = map(1:nfolds) do k
      hold = view(sdata, folds[k])
      w = view(weights, folds[k])
      y = hold[var]
      ŷ = solutions[k][var]
      value(loss[var], y, ŷ, AggMode.WeightedSum(w)) / length(y)
    end
    var => mean(losses)
  end

  Dict(result)
end


text:
#=##############################################################################
# DESCRIPTION
    Vehicle and components performance montiro.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Apr 2020
  * License   : MIT
=###############################################################################

function generate_monitor_vahana(vehicle, rho, RPMref, nsteps, save_path, Vinf)

    # Reference parameters for calculating coefficients
    # NOTE: make b, ar, and qinf equals to 1.0 to obtain dimensional force
    b_ref, ar_ref = 1.0, 1.0
    qinf = 1.0
    Jref = 1.0

    # Force axis labels
    CL_lbl = "Lift (N)"
    CD_lbl = "Induced drag (N)"

    # Directions of force components
    L_dir = [0, 0, 1]
    D_dir = [-1, 0, 0]

    # Rotor monitors
    rotors_monitors = []

    for (si, rotors) in enumerate(vehicle.rotor_systems)
        monitor_name = "rotorsys$(si)"
        rotors_monitor = uns.generate_monitor_rotors(rotors, Jref, rho, RPMref,
                                                        nsteps;
                                                        save_path=save_path,
                                                        run_name=monitor_name,
                                                        figname=monitor_name)
        push!(rotors_monitors, rotors_monitor)
    end

    # Main wing monitor
    monitor_name = "wing_main"
    mainwing_system = vlm.get_wing(vehicle.vlm_system, "MWing")
    mainwing_monitor = uns.generate_monitor_wing(mainwing_system, Vinf, b_ref, ar_ref,
                                                    rho, qinf, nsteps;
                                                    save_path=save_path,
                                                    run_name=monitor_name,
                                                    figname=monitor_name,
                                                    CL_lbl=CL_lbl,
                                                    CD_lbl=CD_lbl,
                                                    L_dir=L_dir,
                                                    D_dir=D_dir)

    # Tandem wing monitor
    monitor_name = "wing_tandem"
    tandemwing_system = vlm.get_wing(vehicle.vlm_system, "TWing")
    tandemwing_monitor = uns.generate_monitor_wing(tandemwing_system, Vinf, b_ref, ar_ref,
                                                    rho, qinf, nsteps;
                                                    save_path=save_path,
                                                    run_name=monitor_name,
                                                    figname=monitor_name,
                                                    CL_lbl=CL_lbl,
                                                    CD_lbl=CD_lbl,
                                                    L_dir=L_dir,
                                                    D_dir=D_dir)

    # Stitch the simulation monitor together
    function monitor(args...)

        for rotors_monitor in rotors_monitors
            rotors_monitor(args...)
        end

        mainwing_monitor(args...)
        tandemwing_monitor(args...)

        return false
    end

    return monitor
end


text:
using MyExample, Test

@testset "MyExample.jl" begin
    @test MyExample.func(2,2) == 4
    @test MyExample.func(3,3) == 9
end

@testset "Derivative tests" begin
    @test MyExample.dfunc(2,2) == 2
end


text:
export Bins, count!, hellinger, data_discrepancy, hellinger_stat_simulations

struct Bins
    Z::Array{Function} # lenght(z) = n List of membership functions
    q::Array{Float64} # lenght(q) = n empirical probabilities associated with each interval
    count::Array{Int64} # simulation counts in each bin
    Bins(Z,q,count) = length(Z) == length(q) ==length(count) ? new(Z,q,count) : error("length mismatch")
end

function Bins(data; min = -Inf, max = +Inf) # censoring cutoffs of continuous measurement and number of measurements outside cutoff
    lows = sum(data .< min)
    highs = sum(data .>= max)
    datapoints = sort(filter(x -> min <= x < max, data)) # get things in order
    midpoints = [(datapoints[ii] + datapoints[ii+1])/2 for ii in 1:length(datapoints)-1]
    boundaries = cat(-Inf, min, midpoints, max, Inf; dims =1)
    Z = [x -> boundaries[ii] <= x <= boundaries[ii+1] for ii in 1:length(boundaries)-1]
    q = cat(lows, ones(length(datapoints)), highs; dims=1)
    q = q/sum(q)
    count = zeros(Int64,length(q))
    Bins(Z,q,count)
end

function count!(x,b::Bins)
    ind = findfirst(z -> b.Z[z](x) ,eachindex(b.Z))
    b.count[ind] += 1
end

function hellinger(b::Bins)
    total = sum(b.count)
    sum( (sqrt(q) - sqrt(c/total))^2 for (c,q) in zip(b.count,b.q))
end

function observed_discrepancy(reference_data, simulated_data; min = -Inf, max = +Inf)
    bins = Bins(reference_data)
    for x in simulated_data
        count!(x,bins)
    end
    hellinger(bins)
end

function data_discrepancy(bins::Bins, θ_vec, ab_profiles... ; samples = 10^4)
    # the data are expected to be minus-log-mutant-frequency objects.
    # The cutoffs are the minimum and maximum data observed.
    b = bins
    b.count .= 0
    for θ in θ_vec
        sampler_array = map(x->WrightSampler(x...),
            hcat([θ .*  profile for profile in ab_profiles]...))
        for _ in 1:samples
            x = -log(escape_freq(sample!.(sampler_array)))
            count!(x, b)
        end
    end
    hellinger(b)
end

function hellinger_stat_simulations(θ_vec, ab_profiles... ;
    data_size=16, min=0.1, max=56/3, # parameters of the data
    samples = 10^3, # parameters of the simulation
    trials =10^3) # * length of θ_vec = generated samples
    datastream = []
    for θ in θ_vec
        sampler_array = map(x->WrightSampler(x...),
                hcat(θ .* [ ab_profiles...]...))
        for ii in 1:trials*data_size
        push!(datastream, -log(escape_freq(sample!.(sampler_array))))
        end
    end
    out =[]
    while length(datastream) >= data_size
        data = [];
        for jj in 1:data_size
            ind = rand(1:length(datastream))
            push!(data, splice!(datastream, ind))
        end
        b = Bins(data; min=min,max=max)
        push!(out, data_discrepancy(b, θ_vec, ab_profiles... ; samples = samples))
    end
    return out
end

max_stars_repo_name: StatPhysBio/Tomoko

max_issues_repo_name: StatPhysBio/Tomoko

max_forks_repo_name: StatPhysBio/Tomoko


text:
struct Interaction
    response::Response
    type::Set{AnalyticProperty}
    function Interaction(response, type)
        return new(response, Set(type))
    end
    function Interaction(response, type::AnalyticProperty)
        return new(response, Set([type,]))
    end
end

Base.isequal(a::Interaction, b::Interaction) = (a.response == b.response) && issetequal(a.type, b.type)
Base.:(==)(a::Interaction, b::Interaction) = Base.isequal(a, b)

function short(name::Response)
    if name == ChargeCharge
        return "cc"
    elseif name == SpinSpin
        return "σσ"
    elseif name == UpUp
        return "↑↑"
    elseif name == UpDown
        return "↑↓"
    else
        @error("$name is not implemented!")
    end
end

function short(type::AnalyticProperty)
    if type == Instant
        return "Ins"
    elseif type == Dynamic
        return "Dyn"
    elseif type == D_Instant
        return "dIns"
    elseif type == D_Dynamic
        return "dDyn"
    else
        @error("$type is not implemented!")
    end
end

function symbol(name::Response, type::AnalyticProperty, addition = nothing)
    if isnothing(addition)
        return Symbol("$(short(name))$(short(type))")
    else
        return Symbol("$(short(name))$(short(type))$(addition)")
    end

end

@with_kw struct GenericPara
    diagType::DiagramType
    innerLoopNum::Int

    isFermi::Bool = true
    spin::Int = 2
    loopDim::Int = 3
    interaction::Vector{Interaction} = [Interaction(ChargeCharge, [Instant,]),] # :ChargeCharge, :SpinSpin, ...
    weightType::DataType = Float64

    firstLoopIdx::Int = firstLoopIdx(diagType)
    totalLoopNum::Int = firstLoopIdx + innerLoopNum - 1

    #### turn the following parameters on if there is tau variables ########
    hasTau::Bool = false
    firstTauIdx::Int = firstTauIdx(diagType)
    totalTauNum::Int = firstTauIdx + innerTauNum(diagType, innerLoopNum, interactionTauNum(hasTau, interaction)) - 1
    #if there is no imaginary-time at all, then set this number to zero!
    ########################################################################
    filter::Vector{Filter} = [NoHatree,] #usually, the Hatree subdiagram should be removed
    transferLoop = [] #Set it to be the transfer momentum/frequency if you want to check the diagrams are proper or not
    extra::Any = Nothing
end

# function Base.show(io::IO, para::GenericPara)

# end

function Base.getproperty(obj::GenericPara, sym::Symbol)
    # if sym === :hasTau
    #     return obj.totalTauNum > 0
    if sym == :interactionTauNum
        return interactionTauNum(obj.hasTau, obj.interaction)
    elseif sym == :innerTauNum
        # println(innerTauNum(obj.diagType, obj.innerLoopNum, obj.interactionTauNum))
        return innerTauNum(obj.diagType, obj.innerLoopNum, obj.interactionTauNum)
    else # fallback to getfield
        return getfield(obj, sym)
    end
end

function Base.isequal(p::GenericPara, q::GenericPara)
    for field in fieldnames(typeof(p)) #fieldnames doesn't include user-defined entries in Base.getproperty
        if field == :filter
            if Set(p.filter) != Set(q.filter)
                return false
            end
        elseif field == :transferLoop
            if (isempty(p.transferLoop) && isempty(q.transferLoop) == false) || (isempty(p.transferLoop) == false && isempty(q.transferLoop))
                return false
            elseif isempty(p.transferLoop) == false && isempty(q.transferLoop) == false
                if (p.transferLoop ≈ q.transferLoop) == false
                    return false
                end
            end
        elseif field == :interaction
            if (p.interaction ⊆ q.interaction) == false || (q.interaction ⊆ p.interaction) == false
                return false
            end
        else
            if Base.getproperty(p, field) != Base.getproperty(q, field)
                return false
            end
        end
    end
    return true
end

Base.:(==)(a::GenericPara, b::GenericPara) = Base.isequal(a, b)

"""
    function innerTauNum(diagType::DiagramType, innerLoopNum, interactionTauNum)

    internal imaginary-time degrees of freedom for a given diagram type and internal loop number.
    For the vertex functions (self-energy, polarization, vertex3, and vertex4), innerTauNum is equivalent to tauNum.
    For the Green function, tauNum = innerTauNum + external tauNum
"""
function innerTauNum(diagType::DiagramType, innerLoopNum, interactionTauNum)
    if diagType == Ver4Diag
        return (innerLoopNum + 1) * interactionTauNum
    elseif diagType == SigmaDiag
        return innerLoopNum * interactionTauNum
    elseif diagType == GreenDiag
        return innerLoopNum * interactionTauNum
    elseif diagType == PolarDiag
        return 1 + innerTauNum(Ver3Diag, innerLoopNum - 1, interactionTauNum)
    elseif diagType == Ver3Diag
        return 1 + innerTauNum(Ver4Diag, innerLoopNum - 1, interactionTauNum)
    else
        error("not implemented!")
    end
end

function interactionTauNum(hasTau::Bool, interactionSet)
    if hasTau == false
        return 0
    end
    for interaction in interactionSet
        if Dynamic in interaction.type || D_Dynamic in interaction.type
            return 2
        end
    end
    return 1
end

function firstTauIdx(diagType::DiagramType, offset::Int = 0)
    if diagType == GreenDiag
        return 3 + offset
    elseif diagType == Ver3Diag
        return 1 + offset
    elseif diagType == PolarDiag
        return 1 + offset
    else
        return 1 + offset
    end
end

function firstLoopIdx(diagType::DiagramType, offset::Int = 0)
    if diagType == Ver4Diag #three extK
        return 4 + offset
    elseif diagType == SigmaDiag #one extK
        return 2 + offset
    elseif diagType == GreenDiag #one extK
        return 2 + offset
    elseif diagType == PolarDiag #one extK
        return 2 + offset
    elseif diagType == Ver3Diag #two extK
        return 3 + offset
    else
        error("not implemented!")
    end
end

function totalTauNum(diagType::DiagramType, innerLoopNum, interactionTauNum, offset::Int = 0)
    return firstTauIdx(diagType, offset) + innerTauNum(diagType, innerLoopNum, interactionTauNum) - 1
end

function totalLoopNum(diagType::DiagramType, innerLoopNum, offset::Int = 0)
    return firstLoopIdx(diagType, offset) + innerLoopNum - 1
end

function totalTauNum(para, diagType::Symbol = :none)
    return para.totalTauNum
    # if diagType == :Ver4
    #     return (para.internalLoopNum + 1) * para.interactionTauNum
    # else
    #     error("not implemented!")
    # end
end

function totalLoopNum(para, diagType::Symbol = :none)
    return para.totalLoopNum
end


text:
##==============================================================================
## OLD Victory Park Example code, don't delete until replacements are coded
##==============================================================================

# function projNewLandmPoints(vps::Graphs.ExVertex, br::Array{Float64,1}, cov::Array{Float64,2})
#   # TODO -- convert to use Distributions and common projection function
#   Xps = getVal(vps)
#   lmPts = zeros(2,size(Xps,2))
#   for i in 1:size(Xps,2)
#       ent = [cov[1,1]*randn(); cov[2,2]*randn()]
#       init = vec(Xps[1:2,i])+randn(2)
#       lmPts[:,i] = solveLandm(br + ent, vec(Xps[:,i]), init)
#   end
#   return lmPts
# end

# function calcIntersectVols( fgl::G, predLm::BallTreeDensity;
#                             currage=0, maxdeltaage=Inf) where G <: AbstractDFG
#     # TODO upgrade to using MMD test
#     # all landmarks of interest
#     ll = ls(fgl, r"l\d")
#     # output result
#     rr = Dict{String, RemoteRef}()
#     fetchlist = String[]
#     iv = Dict{String, Float64}()
#     for l in ll
#       pvlm = getVariable(fgl,l)
#       # TODO -- can be improved via query in DB case
#       if currage - pvlm.attributes["maxage"] < maxdeltaage
#         p = getBelief(fgl, l)
#         rr[l] = remotecall(uppA(), intersIntgAppxIS, p,predLm)
#         push!(fetchlist, l)
#       else
#         println("calcIntersectVols -- ignoring $(l) because maxdeltaage exceeded")
#         iv[l] = 0
#       end
#     end
#     max = 0
#     maxl = String("")
#     for l in fetchlist #ll
#       # p = getBelief(fgl, l)
#       # tv = intersIntgAppxIS(p,predLm)
#       # iv[l] = tv
#       tv = fetch(rr[l])
#       iv[l] = tv
#       if max < tv  max = tv; maxl = l; end
#     end
#     return iv, maxl
# end

# function maxIvWithoutID(ivs::Dict{String, Float64}, l::T) where {T <: AbstractString}
#   max = 0
#   maxl = String("")
#   for i in ivs
#     if max < i[2] && i[1] != l;  max = i[2]; maxl = i[1]; end
#   end
#   return maxl
# end

# # binary tests to distinguish how to automatically add a landmark to the existing factor graph
# function doAutoEvalTests(fgl::G, ivs::Dict{T, Float64}, maxl::T, lmid::Int, lmindx::Int) where {G <: AbstractDFG, T <: AbstractString}
#   maxAnyval = maxl != String("") ? ivs[maxl] : 0.0
#   # maxid = fgl.IDs[maxl]
#   lmidSugg = lmid != -1 # a landmark ID has been suggested
#   maxAnyExists = maxAnyval > 0.03 # there is notable intersection with previous landm
#   lmIDExists = haskey(fgl.v, lmid) # suggested lmid already in fgl
#   newlmindx = lmindx
#   if lmIDExists
#     lmSuggLbl = String(getVert(fgl,lmid).label) # TODO -- wasteful
#   else
#     newlmindx = lmindx + 1
#     lmSuggLbl = String(string('l',newlmindx))
#   end
#   maxl2 = lmIDExists ? maxIvWithoutID(ivs, lmSuggLbl) : String("")
#   maxl2Exists = lmIDExists ? (maxl2 != "" ? ivs[maxl2] > 0.03 : false) : false # there is notable intersection with previous landm
#   intgLmIDExists = lmIDExists ? ivs[lmSuggLbl] > 0.03 : false

#   return lmidSugg, maxAnyExists, maxl2Exists, maxl2, lmIDExists, intgLmIDExists, lmSuggLbl, newlmindx
# end

# function evalAutoCases!(fgl::G, lmid::Int, ivs::Dict{T, Float64}, maxl::T,
#                         pose::T, lmPts::Array{Float64,2}, br::Array{Float64,1}, cov::Array{Float64,2}, lmindx::Int;
#                         N::Int=100, solvable::Int=1 ) where {G <: AbstractDFG, T <: AbstractString}
#   lmidSugg, maxAnyExists, maxl2Exists, maxl2, lmIDExists, intgLmIDExists, lmSuggLbl, newlmindx = doAutoEvalTests(fgl,ivs,maxl,lmid, lmindx)

#   println("evalAutoCases -- found=$(lmidSugg), $(maxAnyExists), $(maxl2Exists), $(lmIDExists), $(intgLmIDExists)")

#   vlm = Union{}; fbr = Union{};
#   if (!lmidSugg && !maxAnyExists)
#     #new landmark and UniBR constraint
#     v,L,lm = getLastLandm2D(fgl)
#     vlm = newLandm!(fgl, lm, lmPts, cov, N=N,solvable=solvable)
#     fbr = addBRFG!(fgl, pose, lm, br, cov, solvable=solvable)
#   elseif !lmidSugg && maxAnyExists
#     # add UniBR to best match maxl
#     vlm = getVariable(fgl,maxl)
#     fbr = addBRFG!(fgl, pose, maxl, br, cov, solvable=solvable)
#   elseif lmidSugg && !maxl2Exists && !lmIDExists
#     #add new landmark and add UniBR to suggested lmid
#     vlm = newLandm!(fgl, lmSuggLbl, lmPts, cov, N=N, solvable=solvable)
#     fbr = addBRFG!(fgl, pose, lmSuggLbl, br, cov, solvable=solvable)
#   elseif lmidSugg && !maxl2Exists && lmIDExists && intgLmIDExists
#     # doesn't self intesect with existing lmid, add UniBR to lmid
#     vlm = getVariable(fgl, lmid)
#     fbr = addBRFG!(fgl, pose, lmSuggLbl, br, cov, solvable=solvable)
#   elseif lmidSugg && maxl2Exists && !lmIDExists
#     # add new landmark and add MMBR to both maxl and lmid
#     vlm = newLandm!(fgl, lmSuggLbl, lmPts, cov, N=N, solvable=solvable)
#     addMMBRFG!(fgl, pose, [maxl2;lmSuggLbl], br, cov, solvable=solvable)
#   elseif lmidSugg && maxl2Exists && lmIDExists && intgLmIDExists
#     # obvious case, add MMBR to both maxl and lmid. Double intersect might be the same thing
#     println("evalAutoCases! -- obvious case is happening")
#     addMMBRFG!(fgl, pose, [maxl2;lmSuggLbl], br, cov, solvable=solvable)
#     vlm = getVariable(fgl,lmSuggLbl)
#   elseif lmidSugg && maxl2Exists && lmIDExists && !intgLmIDExists
#     # odd case, does not intersect with suggestion, but does with some previous landm
#     # add MMBR
#     @warn "evalAutoCases! -- no self intersect with suggested $(lmSuggLbl) detected"
#     addMMBRFG!(fgl, pose, [maxl;lmSuggLbl], br, cov, solvable=solvable)
#     vlm = getVariable(fgl,lmSuggLbl)
#   elseif lmidSugg && !maxl2Exists && lmIDExists && !intgLmIDExists
#   #   # landm exists but no intersection with existing or suggested lmid
#   #   # may suggest some error
#     @warn "evalAutoCases! -- no intersect with suggested $(lmSuggLbl) or map detected, adding  new landmark MM constraint incase"
#     v,L,lm = getLastLandm2D(fgl)
#     vlm = newLandm!(fgl, lm, lmPts, cov, N=N, solvable=solvable)
#     addMMBRFG!(fgl, pose, [lm; lmSuggLbl], br, cov, solvable=solvable)
#   else
#     error("evalAutoCases! -- unknown case encountered, can reduce to this error to a warning and ignore user request")
#   end

#   return vlm, fbr, newlmindx
# end

# function addAutoLandmBR!(fgl::G,
#                          pose::T,
#                          lmid::Int,
#                          br::Array{Float64,1},
#                          cov::Array{Float64,2},
#                          lmindx::Int;
#                          N::Int=100,
#                          solvable::Int=1  ) where {G <: AbstractDFG, T <: AbstractString}
#   #
#   vps = getVariable(fgl, pose)
#   lmPts = projNewLandmPoints(vps, br, cov)
#   lmkde = kde!(lmPts)
#   currage = parse(Int, pose[2:end])
#   ivs, maxl = calcIntersectVols(fgl, lmkde, currage=currage,maxdeltaage=10)

#   # There are 8 cases of interest
#   vlm, fbr, newlmindx = evalAutoCases!(fgl, lmid, ivs, maxl,pose,lmPts, br,cov,lmindx,N=N,solvable=solvable)

#   return vlm, fbr, newlmindx
# end

# function newLandm!(fg::AbstractDFG, lm::T, wPos::Array{Float64,2}, sig::Array{Float64,2};
#                   N::Int=100, solvable::Int=1, labels::Vector{T}=String[]) where {T <: AbstractString}

#     vert=addVariable!(fg, Symbol(lm), Point2, N=N, solvable=solvable, tags=union(["LANDMARK";], labels))
#     # TODO -- need to confirm this function is updating the correct memory location. v should be pointing into graph
#     # vert=addVariable!(fg, Symbol(lm), wPos, sig, N=N, solvable=solvable, tags=labels)

#     vert.attributes["age"] = 0
#     vert.attributes["maxage"] = 0
#     vert.attributes["numposes"] = 0
#     updateFullVert!(fg, vert)

#     println("newLandm! -- added $(lm)")
#     return vert
# end

# function addBRFG!(fg::G,
#                   pose::T,
#                   lm::T,
#                   br::Array{Float64,1},
#                   cov::Array{Float64,2};
#                   solvable::Int=1  ) where {G <: AbstractDFG, T <: AbstractString}
#   #
#   vps = getVert(fg,pose)
#   vlm = getVert(fg,lm)
#   testlbl = vps.label*vlm.label
#   for nei in getOutNeighbors(fg, vlm)
#     if nei.label == testlbl
#       # TODO -- makes function call brittle
#       @warn "We already have $(testlbl), skipping this constraint"
#       return nothing
#     end
#   end
#   @show keys(vlm.attributes)
#   np = vlm.attributes["numposes"]
#   la = vlm.attributes["age"]
#   nage = parse(Int,pose[2:end])
#   vlm.attributes["numposes"] = np+1
#   vlm.attributes["age"] = ((la*np)+nage)/(np+1)
#   vlm.attributes["maxage"] = nage
#   updateFullVert!(fg, vlm)

#   pbr = Pose2Point2BearingRange(Normal(br[1], cov[1,1]), Normal(br[2],  cov[2,2]))  #{Normal, Normal}
#   @show vps, vlm
#   f = addFactor!(fg, [vps;vlm], pbr, solvable=solvable, graphinit=true ) #[vps;vlm],

#   # only used for max likelihood unimodal tests.
#   u, P = pol2cart(br[[2;1]], diag(cov))
#   infor = inv(P^2)
#   # addLandmMeasRemote(vps.index,vlm.index,u,infor) # for iSAM1 remote solution as reference
#   return f
# end

# function addMMBRFG!(fg::G,
#                     syms::Array{Symbol,1}, br::Array{Float64,1},
#                     cov::Array{Float64,2}; w::Vector{Float64}=Float64[0.5;0.5],
#                     solvable::Int=1) where G <: AbstractDFG
#     #
#     # vps = getVert(fg,pose)
#     # vlm1 = getVert(fg,lm[1])
#     # vlm2 = getVert(fg,lm[2])

#     pbr = Pose2Point2BearingRange(Normal(br[1],cov[1,1]),  Normal(br[2],cov[2,2]))
#     syms = Symbol.([pose;lm...])
#     f = addFactor!(fg, syms, pbr, multihypo=[1.0; w...], solvable=solvable, graphinit=true )
#     return f
# end

# function projNewLandm!(fg::G,
#                        pose::T,
#                        lm::T,
#                        br::Array{Float64,1},
#                        cov::Array{Float64,2};
#                        addfactor=true,
#                        N::Int=100,
#                        solvable::Int=1,
#                        labels::Vector{T}=String[]  ) where {G <: AbstractDFG, T <: AbstractString}
#     #
#     vps = getVariable(fg, pose)

#     lmPts = projNewLandmPoints(vps, br, cov)
#     vlm = newLandm!(fg, lm, lmPts, cov, N=N, solvable=solvable, tags=labels) # cov should not be required here
#     if addfactor
#       fbr = addBRFG!(fg, pose, lm, br, cov, solvable=solvable)
#       return vlm, fbr
#     end
#     return vlm
# end

##==============================================================================
## OLD Victory Park Example code, don't delete until replacements are coded
##==============================================================================


text:
# addprocs(4)

using Caesar, RoME
using TransformUtils
using Rotations, CoordinateTransformations
using PyLCM
using PyCall
using Distributions

using CloudGraphs

using KernelDensityEstimate

using JLD, HDF5

@show lcmtpath = joinpath(dirname(@__FILE__),"lcmtypes")

# println("Compile LCM types")
run(`lcm-gen -p --ppath $(lcmtpath) $(lcmtpath)/rome_pose_node_t.lcm`)
run(`lcm-gen -p --ppath $(lcmtpath) $(lcmtpath)/rome_pose_pose_nh_t.lcm`)
# run(`lcm-gen -p --ppath $(lcmtpath) $(lcmtpath)/rome_point_cloud_t.lcm`)
run(`lcm-gen -p --ppath $(lcmtpath) $(lcmtpath)/hauv_submap_points_t.lcm`)

println("Adding lcmtypes dir to Python path: $(lcmtpath)")
unshift!(PyVector(pyimport("sys")["path"]),lcmtpath)

@pyimport rome
@pyimport hauv

# random test and viewing function
# function whatsthere(channel, msgdata)
#   odomsg = rome.pose_pose_nh_t[:decode](msgdata)
#   @show odomsg[:utime]
#   nothing
# end

# ! = non-const function (eyecandy only)
# . = dot-notation (element-wise execution)
function lcmposenodes!(vc, slaml::SLAMWrapper, msgdata,
        lastposeRPZ::Vector{Float64}, MSGDATA;
        record=false,
        model=nothing )
  #
  println("lcmposenodes! start")
  !record ? nothing : push!(MSGDATA[:pose_node_t], deepcopy(msgdata))

  posenode = rome.pose_node_t[:decode](msgdata)
  @show posenode[:id]
  tA = [posenode[:mean][1];posenode[:mean][2];posenode[:mean][3]]
  EA = Euler(posenode[:mean][6], posenode[:mean][5], posenode[:mean][4])
  q = TransformUtils.convert(Quaternion, EA)
  newpose = SE3(tA, q)
  am = Translation(tA...) ∘ LinearMap(Quat(q.s,q.v...))
  # updaterealtime!(vc,:body, am)
  # visualizerealtime(vc)

  # Dx = slaml.lastpose\newpose # get delta transform between poses
  # copy!(slaml.lastpose.t, newpose.t)
  # copy!(slaml.lastpose.R.R, newpose.R.R)
  # println("lcmposenodes! mid")

  lastposeRPZ[1], lastposeRPZ[2] = posenode[:mean][6], posenode[:mean][5]
  lastposeRPZ[3] = posenode[:mean][3]

  if posenode[:id] == 0
    # first pose prior
    dcovar = [posenode[:covar][1];posenode[:covar][2];posenode[:covar][3];posenode[:covar][4];posenode[:covar][5];posenode[:covar][6]]
    # TODO -- this wont do, must add to existing slaml.fg
    # initialize fg by setting the first pose w/ prior
    slaml.fg = identitypose6fg(fg=slaml.fg, initpose=newpose, initCov=diagm(dcovar))

# ???
    lbls, = ls(slaml.fg)
    slaml.usrid2lbl[posenode[:id]] = lbls[1]
    slaml.lbl2usrid[lbls[1]] = posenode[:id]
  end

  # update rov pose in visuzalizer
  # TODO: add rovt as argument to fcn
  model( vc, am )

  # change to using pose prior and interpose translation only
  # println("Added node $(v.label)")
  # println("lcmposenodes! end")
  nothing
end

function lcmodomsg!(vc, slaml::SLAMWrapper, msgdata,
        lastposeRPZ::Vector{Float64}, MSGDATA;
        record=false  )
  #
  println("lcmodomsg! start")
  @show lastposeRPZ
  !record ? nothing : push!(MSGDATA[:pose_pose_nh_t], deepcopy(msgdata))

  odomsg = rome.pose_pose_nh_t[:decode](msgdata)
  # @show odomsg[:mean_dim]
  @show odomsg[:node_1_id], odomsg[:node_2_id]
  # @show odomsg[:confidence]
  @show tA = [odomsg[:mean][1];odomsg[:mean][2];odomsg[:mean][3]]
  roll,pitch,yaw = odomsg[:mean][6], odomsg[:mean][5], odomsg[:mean][4]
  println("dXYZ=$(round(tA, 2)), dRPY=$(round(collect((roll,pitch,yaw)),3))")
  EA = Euler(roll,pitch,yaw)
  Dx = SE3(tA, EA)
  dcovar = [odomsg[:covar][1];odomsg[:covar][2];odomsg[:covar][3];odomsg[:covar][6];odomsg[:covar][5];odomsg[:covar][4]]
  dcovar = dcovar

  # println("lcmodomsg! mid")

  if abs(odomsg[:confidence] - 1.0) < 1e-10
    constrs = IncrementalInference.FunctorInferenceType[] # factor super type
    println("adding PartialPriorRollPitchZ")
    prpz = PartialPriorRollPitchZ(
      MvNormal(  [lastposeRPZ[1];lastposeRPZ[2]],diagm([odomsg[:covar][6];odomsg[:covar][5]])  ),
      Normal(  lastposeRPZ[3],odomsg[:covar][3]  )
    )
    # println("adding PartialPose3XYYaw")
    # xyy = PartialPose3XYYaw(
    #   MvNormal([odomsg[:mean][1];odomsg[:mean][2];odomsg[:mean][4]],
    #            [odomsg[:covar][1];odomsg[:covar][2];odomsg[:covar][4]] )
    # )

    xyy = Pose3Pose3( MvNormal([tA; EA], diagm(dcovar)^2 ) )

    push!( constrs, prpz )
    push!( constrs, xyy )

    # may or may not add a new pose node, depending on message order
    v,f = addposeFG!(slaml, constrs, saveusrid=odomsg[:node_2_id], ready=0 )
    setDBAllReady!(slaml.fg)

  else
    # loop closure case
    println("LOOP CLOSURE")
    odoc3 = Pose3Pose3NH( MvNormal(veeEuler(Dx), diagm(1.0./dcovar)), [0.5;0.5]) # define 50/50% hypothesis
    between_verts = Symbol[Symbol("x$(odomsg[:node_1_id]+1)");Symbol("x$(odomsg[:node_2_id]+1)")]
    # order in between_verts must match the order assumed in the factor functor
    addFactor!(slaml.fg,between_verts,odoc3)
  end

  # update viewer
  visualizeallposes!(vc, slaml.fg)

  # println("lcmodomsg! end")
  nothing
end

# call back for point cloud messages
function lcmpointcloudmsg!(slaml::SLAMWrapper,
                          msgdata,
                          MSGDATA;
                          record::Bool=false  )
  #
  !record ? nothing : push!(MSGDATA[:submap_points_t], deepcopy(msgdata))

  msg = hauv.submap_points_t[:decode](msgdata)
  @show msg[:id]
  @show typeof(msg[:points_local])

  node_name =  Symbol("x$(odomsg[:node_1_id]+1)")
  # function appendvertbigdata!(
  vert = getVert(slaml.fg, )
  appendvertbigdata!(slaml.fg, )

  nothing
end

function runlistener!(
      fl::Vector{Bool},
      slam::SLAMWrapper,
      vc;
      MSGDATA=nothing,
      until::Symbol=:x999,
      record=false,
      model = nothing  )
  #
  lastposeRPZ = zeros(3)

  if record
    MSGDATA[:pose_node_t] = Vector{Vector{UInt8}}()
    MSGDATA[:pose_pose_nh_t] = Vector{Vector{UInt8}}()
    MSGDATA[:submap_points_t] = Vector{Vector{UInt8}}()
  end

  # closure
  gg = (channel, msg_data) -> lcmposenodes!(vc, slam, msg_data, lastposeRPZ, MSGDATA, record=record, model=model )
  odohdl = (channel, msg_data) -> lcmodomsg!(vc, slam, msg_data, lastposeRPZ, MSGDATA, record=record )
  pchdl = (channel, msg_data) -> lcmpointcloudmsg!(slam, msg_data, MSGDATA, record=record )

  print("Preparing LCM...")
  lc = LCM()
  # subscribe(lc, "ROME_NODES", gg)
  # subscribe(lc, "ROME_FACTORS", odohdl)
  # subscribe(lc, "ROME_POINTCLOUD", pchdl)
  subscribe(lc, "SUBMAPS", pchdl)

  while fl[1]
    handle(lc)
    if slam.lastposesym == until
      break
    end
  end
  println("done")
  nothing
end

function lcmsendpose(lc, vert::Graphs.ExVertex, usrid::Int)
  msg = rome.pose_node_t()
  # p = getVertKDE(vert)
  p = getKDE(vert)
  mv = getKDEfit(p)
  # pme = getKDEMean(vert)
  msg[:utime] = 0
  msg[:id] = usrid
  msg[:mean_dim] = 6
  msg[:mean] = mv.μ
  msg[:covar_dim] = 6
  msg[:covar] = diag(mv.Σ)

  publish(lc, "MMISAM_POSE_EST_MEAN", msg)
  nothing
end
function lcmsendallposes(lc, slaml::SLAMWrapper)
  for lbluid in slaml.lbl2usrid
    @show lbluid
    lcmsendpose(lc, getVert(slaml.fg, lbluid[1]), lbluid[2])
  end
  nothing
end

function setupSLAMinDB(;cloudGraph=nothing, addrdict=nothing)
  if cloudGraph != nothing
    return SLAMWrapper(initfg(sessionname=addrdict["session"], cloudgraph=cloudGraph), nothing, 0)
  else
    return SLAMWrapper(RoME.initfg(), nothing, 0)
  end
end

vc = startdefaultvisualization()
sleep(3.0)
rovt = loadmodel(:rov)
rovt(vc)

Nparticles = 100

# for local dictionary usage
# slam = setupSLAMinDB()

0

# for remote db usage
include(joinpath(dirname(@__FILE__),"..","database","blandauthremote.jl"))
addrdict["session"] = "SESSHAUVDEV"
# addrdict are the server credentials and session id
# if none are given, will prompt user
cloudGraph, addrdict = standardcloudgraphsetup(addrdict=addrdict)
# # cloudGraph, addrdict = standardcloudgraphsetup()
slam = setupSLAMinDB(cloudGraph=cloudGraph, addrdict=addrdict)

flags = Bool[true]

MSGDATA = Dict{Symbol, Vector{Vector{UInt8}}}()

println("LCM listener running...")
# can call @async runlistener!)
@async runlistener!(flags, slam, vc, until=:x1, MSGDATA=MSGDATA, record=false, model=rovt)
flags[1] = false

####

MSGDATA[:submap_points_t]
msg = hauv.submap_points_t[:decode](MSGDATA[:submap_points_t][1])

hauv.submap_points_t

# slamdb10 = deepcopy(slam)

savejld(slam.fg, file="hauvpedro.jld")

fullLocalGraphCopy!(slam.fg)

writeGraphPdf(slam.fg)
run(`evince fg.pdf`)

#

# # visualizeDensityMesh!(vc, slam.fg, :x30)
visualizeallposes!(vc, slam.fg, drawtype=:fit)

f1 = getfnctype(slam.fg, slam.fg.fIDs[:x35x57])
f1.Zij

# loopclosures
@show lcfncs = lsf(slam.fg, Pose3Pose3NH)

session = addrdict["session"]
@show IDs = getPoseExVertexNeoIDs(cloudGraph.neo4j.connection, session=session, reqbackendset=false);

using IncrementalInference

plotKDEresiduals(slam.fg, :x15x61, marg=[1;2], N=Nparticles)

using Gadfly

plotMCMC(tree, :x60, dims=[1;2], levels=1)

tree = wipeBuildNewTree!(slam.fg, drawpdf=true)
@async run(`evince bt.pdf`)

inferOverTree!(slam.fg,tree,N=Nparticles, dbg=true)

pl = spyCliqMat(tree, :x61)
Gadfly.draw(PDF("spycliqmatx61.pdf",15cm,20cm),pl)
run(`evince spycliqmatx61.pdf`)
run(`cp hauv67.jld /home/dehann/Pictures/hauv/`)

[solve(slam.fg) for i in 1:1] #, densitymeshes=[:x1;:x33;:x60])
visualize(slam.fg, vc, drawtype=:fit)

using IncrementalInference

fg = slam.fg

fg = slamdict10.fg

getVal(fg, :x1)

fullLocalGraphCopy!(slam.fg)

sp = plotKDEresiduals(fg, :x1x2, marg=[1;2], N=Nparticles)
sp = plotKDEresiduals(fg, :x1x2, marg=[3;6], N=Nparticles)

# Gadfly.draw(PNG("/home/dehann/Pictures/sp3.png",20cm,15cm),sp)

savejld(fg) # file=tempfg.jld

using RoME

# Juno.breakpoint("/home/dehann/.julia/v0.5/RoME/src/fgos.jl",19)

fgu = loadjld(file="test.jld")

writeGraphPdf(fgu)
run(`evince fg.pdf`)

# end

1
# # loopclosures
# @show lcfncs = lsf(fg, Pose3Pose3NH)
# #

#
#
# # plotKDEofnc(fg, :x1x2, marg=[1;2], N=Nparticles)
#
# donum = 5
# plotKDEofnc(fg, lcfncs[donum], marg=[1;2], N=Nparticles)
# plotKDEofnc(fg, lcfncs[donum], marg=[6;3], N=Nparticles)
# plotKDEofnc(fg, lcfncs[donum], marg=[4;5], N=Nparticles)
#
#

fgu = loadjld(file="hauvpedro.jld")

plotKDE(fgu, :x60, marg=[1;2])
#

IncrementalInference.plotKDEresiduals(fgu, :x20x21, marg=[1;2], N=100)

IncrementalInference.lsf(fgu, PartialPriorRollPitchZ)

# #
# #
# # println("AD-HOC DEBUGGING")
# #
# # slam = setupSLAMinDB()
# #
# #
# # lcmposenodes!(vc, slam, MSGDATA[:pose_node_t][1], record=false)
# # lcmposenodes!(vc, slam, MSGDATA[:pose_node_t][2], record=false)
# #
# #
# # lcmodomsg!(vc, slam, MSGDATA[:pose_pose_nh_t][1], record=false)
# #

#


text:
"""
Find each mesh color.

# Notes
`meshes` is typically the output from [`meshes_from_topology`].

# Exemples
opf= read_opf("simple_OPF_shapes.opf")
mesh_data= meshes_from_topology(opf)
shapes= extract_opf_shapes(opf)
materials= extract_opf_materials(opf)
cols= meshes_color(opf,mesh_data)
mesh(merge_meshes(mesh_data["mesh"]), color= cols)
"""
function meshes_color(opf,meshes::Dict{String,Dict}= meshes_from_topology(opf))

    mesh_ids= collect(keys(meshes["mesh"]))
    materials= extract_opf_materials(opf)

    # n_triangles= sum([length(i.second.faces) for i in meshes["mesh"]])
    n_triangles= Dict{Int32,Int32}()
    for i in meshes["mesh"]
        push!(n_triangles, i.first => length(i.second.vertices))
    end
    n_triangles_tot= sum(values(n_triangles))
    prev_max= [0]
    # mesh_color= Dict{Int32, RGBA{Float64}}() # To get only the color per mesh ID
    triangle_color= Array{RGBA{Float64}}(undef, n_triangles_tot)

    for i in 1:length(mesh_ids)
        material_id= meshes["attributes"][mesh_ids[i]]["materialIndex"]
        col= RGBA(materials[material_id].diffuse...)
        # push!(mesh_color, mesh_ids[i] => col)
        for j in ((1:n_triangles[mesh_ids[i]]) .+ prev_max[1])
            triangle_color[j]= col
        end
        prev_max[1] += n_triangles[mesh_ids[i]]
    end

    return triangle_color
end


text:
# Define the dynamics model of the model.
struct QuadrotorGame{N,M,P,SVu,SVx,SVz,T} <: AbstractGameModel
    n::Int  # Number of states
    m::Int  # Number of controls
	p::Int  # Number of players
	ni::Vector{Int}  # Number of states for each player
	mi::Vector{Int}  # Number of controls for each player
	pu::SVu # Indices of the each player's controls
	px::SVx # Indices of the each player's x and y positions
	pz::SVz # Indices of the each player's states
	mass::T # mass of the quadrotor, in kg (default = 0.5)
    J::Diagonal{T,SVector{3,T}} # inertia of the quadrotor, in kg⋅m² (default = `Diagonal([0.0023, 0.0023, 0.004])`)
    Jinv::Diagonal{T,SVector{3,T}} #
    gravity::SVector{3,T} # gravity vector, in kg/m² (default = [0,0,-9.81])
    motor_dist::T # distance between the motors, in m (default = 0.1750)
    kf::T # motor torque constant (default = 0.0245)
    km::T # motor force constant (default = 1.0)
end

function QuadrotorGame(;p::Int=2, mass::T=0.5) where {T}
	@assert p <= 4 # if p > 4, you need to define the dynamics function associated with it.
	J = Diagonal(SVector{3,T}([0.0023, 0.0023, 0.004]))
	J_inv = Diagonal(SVector{3,T}(1 ./ [0.0023, 0.0023, 0.004]))
	gravity = SVector{3,T}([0,0,-9.81])
	motor_dist = 0.1750
	# kf = 0.0245
	# kf = 0.245
	kf = 1.245
	km = 1.0

	# p = number of players
	n = 12p
	m = 4p
	pu = [SVector{4,Int}([i + (j-1)*p for j=1:4]) for i=1:p]
	px = [SVector{2,Int}([i + (j-1)*p for j=1:2]) for i=1:p]
	pz = [SVector{12,Int}([i + (j-1)*p for j=1:12]) for i=1:p]
	TYPE = typeof.((pu,px,pz))
	ni = 12*ones(Int,p)
	mi = 4*ones(Int,p)
	return QuadrotorGame{n,m,p,TYPE...,T}(
		n,m,p,ni,mi,pu,px,pz,
		mass, J, J_inv, gravity, motor_dist, kf, km)
end

inertia(model::QuadrotorGame) = model.J
inertia_inv(model::QuadrotorGame) = model.Jinv
mass(model::QuadrotorGame) = model.mass

function forces(model::QuadrotorGame{N,M,P,SVu,SVx,SVz,T}, x, u, i::Int) where {N,M,P,SVu,SVx,SVz,T}
	R = MRP(x[3P+i], x[4P+i], x[5P+i]) # MRP indices for player i
    kf = model.kf
    g = model.gravity
    m = model.mass

    w1 = u[0P+i]
    w2 = u[1P+i]
    w3 = u[2P+i]
    w4 = u[3P+i]

    F1 = max(0,kf*w1);
    F2 = max(0,kf*w2);
    F3 = max(0,kf*w3);
    F4 = max(0,kf*w4);
    F = @SVector [0., 0., F1+F2+F3+F4] #total rotor force in body frame

    f = m*g + R*F # forces in world frame
    return f
end

function moments(model::QuadrotorGame{N,M,P,SVu,SVx,SVz,T}, x, u, i::Int) where {N,M,P,SVu,SVx,SVz,T}
    kf, km = model.kf, model.km
    L = model.motor_dist

    w1 = u[0P+i]
    w2 = u[1P+i]
    w3 = u[2P+i]
    w4 = u[3P+i]

    F1 = max(0,kf*w1);
    F2 = max(0,kf*w2);
    F3 = max(0,kf*w3);
    F4 = max(0,kf*w4);

    M1 = km*w1;
    M2 = km*w2;
    M3 = km*w3;
    M4 = km*w4;
    tau = @SVector [L*(F2-F4), L*(F3-F1), (M1-M2+M3-M4)] #total rotor torque in body frame
    return tau
end

function wrenches(model::QuadrotorGame, x::SVector, u::SVector, i::Int)
    F = forces(model, x, u, i)
    M = moments(model, x, u, i)
    return F, M
end

function RobotDynamics.dynamics(model::QuadrotorGame{N,M,P,SVu,SVx,SVz,T},
		x::SVector, u::SVector, i::Int) where {N,M,P,SVu,SVx,SVz,T}
	#
	q = @SVector [x[3P+i], x[4P+i], x[5P+i]]
	v = @SVector [x[6P+i], x[7P+i], x[8P+i]]
	ω = @SVector [x[9P+i], x[10P+i], x[11P+i]]
	R = MRP(q)

    # Original dynamics
    F,τ = wrenches(model, x, u, i)
    m = mass(model)
    J = inertia(model)
    Jinv = inertia_inv(model)

    xdot = v
    qdot = Rotations.kinematics(R,ω)
    vdot = F ./ m
    ωdot = Jinv*(τ - ω × (J*ω))

	return xdot, qdot, vdot, ωdot
end

function RobotDynamics.dynamics(model::QuadrotorGame{N,M,1,SVu,SVx,SVz,T},
		x::SVector, u::SVector) where {N,M,SVu,SVx,SVz,T}

	xdot1, qdot1, vdot1, ωdot1 = RobotDynamics.dynamics(model, x, u, 1)
	return @SVector [
		xdot1[1],
		xdot1[2],
		xdot1[3],
		qdot1[1],
		qdot1[2],
		qdot1[3],
		vdot1[1],
		vdot1[2],
		vdot1[3],
		ωdot1[1],
		ωdot1[2],
		ωdot1[3],
		]
end

function RobotDynamics.dynamics(model::QuadrotorGame{N,M,2,SVu,SVx,SVz,T},
		x::SVector, u::SVector) where {N,M,SVu,SVx,SVz,T}

	xdot1, qdot1, vdot1, ωdot1 = RobotDynamics.dynamics(model, x, u, 1)
	xdot2, qdot2, vdot2, ωdot2 = RobotDynamics.dynamics(model, x, u, 2)
	return @SVector [
		xdot1[1], xdot2[1],
		xdot1[2], xdot2[2],
		xdot1[3], xdot2[3],
		qdot1[1], qdot2[1],
		qdot1[2], qdot2[2],
		qdot1[3], qdot2[3],
		vdot1[1], vdot2[1],
		vdot1[2], vdot2[2],
		vdot1[3], vdot2[3],
		ωdot1[1], ωdot2[1],
		ωdot1[2], ωdot2[2],
		ωdot1[3], ωdot2[3],
		]
end

function RobotDynamics.dynamics(model::QuadrotorGame{N,M,3,SVu,SVx,SVz,T},
		x::SVector, u::SVector) where {N,M,SVu,SVx,SVz,T}

	xdot1, qdot1, vdot1, ωdot1 = RobotDynamics.dynamics(model, x, u, 1)
	xdot2, qdot2, vdot2, ωdot2 = RobotDynamics.dynamics(model, x, u, 2)
	xdot3, qdot3, vdot3, ωdot3 = RobotDynamics.dynamics(model, x, u, 3)
	return @SVector [
		xdot1[1], xdot2[1], xdot3[1],
		xdot1[2], xdot2[2], xdot3[2],
		xdot1[3], xdot2[3], xdot3[3],
		qdot1[1], qdot2[1], qdot3[1],
		qdot1[2], qdot2[2], qdot3[2],
		qdot1[3], qdot2[3], qdot3[3],
		vdot1[1], vdot2[1], vdot3[1],
		vdot1[2], vdot2[2], vdot3[2],
		vdot1[3], vdot2[3], vdot3[3],
		ωdot1[1], ωdot2[1], ωdot3[1],
		ωdot1[2], ωdot2[2], ωdot3[2],
		ωdot1[3], ωdot2[3], ωdot3[3],
		]
end

function RobotDynamics.dynamics(model::QuadrotorGame{N,M,4,SVu,SVx,SVz,T},
		x::SVector, u::SVector) where {N,M,SVu,SVx,SVz,T}

	xdot1, qdot1, vdot1, ωdot1 = RobotDynamics.dynamics(model, x, u, 1)
	xdot2, qdot2, vdot2, ωdot2 = RobotDynamics.dynamics(model, x, u, 2)
	xdot3, qdot3, vdot3, ωdot3 = RobotDynamics.dynamics(model, x, u, 3)
	xdot4, qdot4, vdot4, ωdot4 = RobotDynamics.dynamics(model, x, u, 4)
	return @SVector [
		xdot1[1], xdot2[1], xdot3[1], xdot4[1],
		xdot1[2], xdot2[2], xdot3[2], xdot4[2],
		xdot1[3], xdot2[3], xdot3[3], xdot4[3],
		qdot1[1], qdot2[1], qdot3[1], qdot4[1],
		qdot1[2], qdot2[2], qdot3[2], qdot4[2],
		qdot1[3], qdot2[3], qdot3[3], qdot4[3],
		vdot1[1], vdot2[1], vdot3[1], vdot4[1],
		vdot1[2], vdot2[2], vdot3[2], vdot4[2],
		vdot1[3], vdot2[3], vdot3[3], vdot4[3],
		ωdot1[1], ωdot2[1], ωdot3[1], ωdot4[1],
		ωdot1[2], ωdot2[2], ωdot3[2], ωdot4[2],
		ωdot1[3], ωdot2[3], ωdot3[3], ωdot4[3],
		]
end

dim(model::QuadrotorGame) = 3


text:
function pr_feedback(history)
    return history[end].a == 1
end

sim = simulate(mdl1; feedback=pr_feedback, init_θ = (α = 1.0, β=100.0, Values = reshape([1.0, 0.0], 2, 1)))

@test mean(sim.data.r) == 1.0

max_stars_repo_name: sqwayer/AnimalBehavior

max_issues_repo_name: sqwayer/AnimalBehavior


text:
function writeattributes(file, p::Parameters)
    attrs(file)["Integrator"] = p.integrator.name
    attrs(file)["Thermostat"] = p.thermo.name
    attrs(file)["Potential"] = p.potential.name
    attrs(file)["nsteps"] = p.nsteps
    attrs(file)["dtsampling"] = p.dtsampling
    attrs(file)["dtintegration"] = p.dt
    attrs(file)["Q"] = p.Q
    attrs(file)["T"] = p.T
    attrs(file)["nsimulations"] = p.nsimulations
end

"""
A .hdf5 file is generated where the lyapunov spectrum for each initial condition (nsimulations = number of initial conditions) is saved.
"""
function lyapunov_exponents(p::Parameters)

    data = Matrix{Float64}(p.nsimulations,6)

    for i in 1:p.nsimulations
        data[i,:] = lyapunov_simulation(p)
        #println("Simulation$i done")
    end

    return data
end

"""
    A long trajectory is generated by parts (coincident with nsimulations) and saved in an .hdf5  file
"""
function trajectory(p::Parameters)

    beta = 1./p.T
    r0 = initialcondition(beta, p.Q)
    tx = trajectory_simulation(p, r0)
    r0 = tx[end,:][2:end]
    println("Part 1 done. ")

    for i in 2:p.nsimulations
        traj = trajectory_simulation(p, r0)
        traj = traj[2:end,:]
        traj[:,1] += (i-1)*p.nsteps*p.dt
        tx = vcat(tx, traj)
        r0 = traj[end,:][2:end]
        println("Part $i done.")
    end

    return tx
end

"""
Creates the file "filename".hdf5 with attributes associated to the passed parameters and
whose dataset is the Matrix generated by the function lyapunov_exponents
"""
function save_lyapunov(filename::String, data::Matrix{Float64}, p::Parameters)
    try
        mkdir("../data")
    end

    file = h5open("../data/filename.hdf5", "w")
    writeattributes(file,p)

    file["data"] = data
    close(file)

    println("The file $filename.hdf5 was saved in the folder `../data/`")
end

"""
Creates the file "filename".hdf5 with attributes associated to the passed parameters and
whose dataset is the Matrix generated by the function lyapunov_exponents
"""
function save_trajectory(filename::String, data::Matrix{Float64}, p::Parameters)
    try
        mkdir("../data")
    end

    file = h5open("../data/filename.hdf5", "w")
    writeattributes(file,p)

    file["data"] = data
    close(file)

    println("The file $filename.hdf5 was saved in the folder `../data/`")

end

max_stars_repo_name: dapias/STODs

max_issues_repo_name: dapias/STODs

max_forks_repo_name: dapias/STODs


text:
using ProjectMoana
using ColorTypes
## TESTING RAY METHODS ###############################################

@testset "Cameras: Ray" begin
    ray1 = Ray(Point(1.0, 2.0, 3.0), Vec(5.0, 4.0, -1.0))
    ray2 = Ray(Point(1.0, 2.0, 3.0), Vec(5.0, 4.0, -1.0))
    ray3 = Ray(Point(5.0, 1.0, 4.0), Vec(3.0, 9.0, 4.0))
    ray4 = Ray(Point(1.0, 2.0, 4.0), Vec(4.0, 2.0, 1.0))

    @test isapprox(ray1, ray2)
    @test isapprox(ray1, ray3) == false
    @test isapprox(at(ray4, 0.0), ray4.origin)
    @test isapprox(at(ray4, 1.0), Point(5.0, 4.0, 5.0))
    @test isapprox(at(ray4, 2.0), Point(9.0, 6.0, 6.0))

end

@testset "Cameras: Camera" begin

    ray = Ray( Point(1.0, 2.0, 3.0), Vec(6.0, 5.0, 4.0))
    transformation = translation( Vec(10.0, 11.0, 12.0) ) * rotation_x(pi/2)
    transformed = transformation * ray

    @test isapprox( Point(11.0, 8.0, 14.0), transformed.origin)
    @test isapprox( Vec(6.0, -4.0, 5.0), transformed.dir)

    camera = OrthogonalCamera(2.0)
    ray1 = fire_ray(camera, 0.0, 0.0)
    ray2 = fire_ray(camera, 1.0, 0.0)
    ray3 = fire_ray(camera, 0.0, 1.0)
    ray4 = fire_ray(camera, 1.0, 1.0)

    @test isapprox( 0.0, squared_norm(cross(ray1.dir,ray2.dir)) )
    @test isapprox( 0.0, squared_norm(cross(ray1.dir,ray3.dir)) )
    @test isapprox( 0.0, squared_norm(cross(ray1.dir,ray4.dir)) )

    @test isapprox( at(ray1,1.0), Point(0.0, 2.0, -1.0) )
    @test isapprox( at(ray2,1.0), Point(0.0, -2.0, -1.0) )
    @test isapprox( at(ray3,1.0), Point(0.0, 2.0, 1.0) )
    @test isapprox( at(ray4,1.0), Point(0.0, -2.0, 1.0 ))

    cam = OrthogonalCamera(1.0, translation(Vec(0.0,-1.0,0.0)*2.0)*rotation_z(pi/2))
    ray = fire_ray(cam, 0.5, 0.5)

    @test isapprox( at(ray, 1.0), Point(0.0, -2.0, 0.0))
end

## TESTING ImageTracer METHODS ###############################################

image  = HdrImage(4, 2)
camera = PerspectiveCamera(2)
tracer = ImageTracer(image, camera)

ray1 = fire_ray(tracer, 0, 0, 2.5, 1.5)
ray2 = fire_ray(tracer, 2, 1, 0.5, 0.5)

fire_all_rays(tracer, ray -> RGB(1.0, 2.0, 3.0 ))

top_left_ray = fire_ray(tracer, 1, 1, 0.0, 0.0)
bottom_right_ray = fire_ray(tracer, 4, 2, 1.0, 1.0)

@testset "Cameras: ImageTracer" begin
    # uv submapping
    @test isapprox(ray1[1], ray2[1])

    # orientation
    @test at(bottom_right_ray[1], 1.0) ≈ Point(0.0, -2.0, -1.0)
    @test at(top_left_ray[1], 1.0) ≈ Point(0.0, 2.0, 1.0)

    # image coverage
    for row ∈ 1:image.height
        for col ∈ 1:image.width
            @test image.pixels[get_pixel(image, col, row)] == RGB(1.0, 2.0, 3.0)
        end
    end

end


text:
include("parameters.jl")
include("transmatrix.jl")
include("expressions.jl")
include("steadystate.jl")
include("classify.jl")

function Shocks(names::Array{String})
    dictionary = Dict(zip(names,1:length(names)))
    syms = createsymversions(names)
    return Shocks(names,dictionary,syms)
end

function parseequationsandvars(vars::Array{String},equations::Array{String})

    if length(equations)!=length(vars)
        error("You have $(length(equations)) for $(length(vars)) endogenous variables.")
    end
    var_dictionary = Dict(zip(vars,1:length(vars)))
    # Define arrays for all the symbolic objects, except states, because we don't know them yet
    leads_sym = createsymversions(vars,"_lead")
    contemp_sym = createsymversions(vars)

    # Check the equations have zero or one equals signs. If one, convert to a one-sided expressions
    zeroedequations = rearrangeequations(equations)

    # Sympify the expressions and detect state variables while doing so
    eqsyms,var_isstate,states_sym,states_indices = sympifyallexpressions(zeroedequations,leads_sym,vars)

    # Split the equations into 'exogenous' and everything else. The exogenous ones are equations that depend entirely on
    # states and shocks - e.g. stochastic equations.
    eq_isexo,var_isexo = findexoequations(eqsyms,leads_sym,contemp_sym)
    # Now go through and mark which equations and variables are static- this gets used later to sub out non-transition equations
    eq_isstatic,var_isstatic = findstaticexpressions(eqsyms,leads_sym,contemp_sym,states_sym,var_isstate,eq_isexo,var_isexo)

    equations = Equations(equations,eq_isstatic,eq_isexo,eqsyms)
    vars = Variables(vars,var_dictionary,var_isstate,var_isstatic,var_isexo,states_indices,states_sym,leads_sym,contemp_sym)
    return equations,vars
end

function RSDSGEModel(var_names::Array{String},shock_names::Array{String},parameter_names::Array{String},equations::Array{String},parameter_values::Array,regime_transmatrix::Array{Float64,2},ss_guess::Array{Float64,1}=zeros(Float64,length(var_names)))
    # Some useful constants first
    numvars = length(var_names)
    numshocks = length(shock_names)
    numparameters = length(parameter_names)
    numregimes = size(regime_transmatrix,1)

    # Check all the names are ok
    checknames(var_names)
    checknames(shock_names)
    checknames(parameter_names)

    transmatrix = TransitionMatrix(regime_transmatrix)
    shocks = Shocks(shock_names)
    equations,vars = parseequationsandvars(var_names,equations)
    meta = Meta(countnz(vars.isstate),numvars,numshocks,numregimes,numparameters,countnz(vars.isstatic),countnz(vars.isexo))
    parameters = Parameters!(meta,equations,vars,shocks,transmatrix,parameter_names,parameter_values,ss_guess)
    # Check for unknown symbols in any of the equations
    checksymbols(equations,vars,shocks,parameters)

    # Create the steady state system
    steadystate = SteadyState(meta,equations,vars,parameters,shocks,ss_guess)

    # Create the perturbation system
    perturbationsystem = PerturbationSystem(meta,parameters,vars,shocks,equations,steadystate,transmatrix)

    return RSDSGEModel(meta,parameters,vars,shocks,equations,steadystate,transmatrix,perturbationsystem)
end