text:
# solve the ground state of Lieb-Liniger model
# also demonstrates the usage of preconditioners in Optim.jl

using TensorKit
using TensorOperations
using TensorKitAD
using KrylovKit
using Zygote
using KrylovKit
using LinearAlgebra
using Plots
using Optim
using ChainRules
using ChainRulesCore
using Printf

using Revise
using statmech_tm_solver

function Base.show(io::IO, t::OptimizationState)
    @printf io "%4d   %14.14e   %14.14e\n" t.iteration t.value t.g_norm
end
#########################################################

# parameters
c = 2
μ = 1
L = 50
d = 1
nsteps1 = 200 # ordinary LBFGS steps
nsteps2 = 1000 # preconditioned LBFGS steps

# cost function (energy) and its gradient
function f(ψ_arr::Array{ComplexF64, 3})
    ψ = convert_to_cmps(ψ_arr)
    return energy_lieb_liniger(ψ, c, L, μ)
end

function g!(ψg_arr::Array{ComplexF64, 3}, ψ_arr::Array{ComplexF64, 3})
    copyto!(ψg_arr, f'(ψ_arr))
end

# construct the preconditioner
function _precondprep!(P::preconditioner, ψ_arr::Array{ComplexF64, 3})
    ψ = convert_to_cmps(ψ_arr)
    ψgrad = convert_to_tensormap(f'(ψ_arr), 2)

    P1 = preconditioner(ψ, ψgrad, L; gauge=:periodic)
    P.map = P1.map
    P.invmap = P1.invmap
    P.proj = P1.proj
end

ψm = cmps(rand, 2, d)
#ψm = quickload("lieb_liniger_c$(c)_mu$(μ)_L$(L)_chi4") |> convert_to_cmps

for χ in [2; 4; 8; 12; 16; 20]
    global ψm

    # important: random init or start with existing results
    #ψ0 = expand(ψm, χ, 0.01)
    ψ0 = quickload("lieb_liniger_c$(c)_mu$(μ)_L$(L)_chi$(χ)") |> convert_to_cmps

    ψ_arr0 = convert_to_array(ψ0)
    E0 = energy_lieb_liniger(ψ0, c, L, μ)

    res = optimize(f, g!, ψ_arr0, LBFGS(), Optim.Options(show_trace=true, iterations=nsteps1))
    ψ_arr = Optim.minimizer(res);
    ψm = convert_to_cmps(ψ_arr);
    ψm_arr = convert_to_tensormap(f'(ψ_arr), 2);

    P0 = preconditioner(ψm, ψm_arr, L);
    res = optimize(f, g!, ψ_arr, LBFGS(P = P0, precondprep=_precondprep!), Optim.Options(show_trace=true, iterations=nsteps2))
    ψm = convert_to_cmps(Optim.minimizer(res));

    # calculate observables. total, kinetic and potential energies
    E = energy_lieb_liniger(ψm, c, L, μ)
    env = finite_env(K_mat(ψm, ψm), L)
    env = permute(env, (2, 3), (4, 1))

    op_ρ = particle_density(ψm)
    op_k = kinetic(ψm)
    ρ_value = tr(env * op_ρ)
    k_value = tr(env * op_k)

    @show χ, E, ρ_value, k_value

    quicksave("lieb_liniger_c$(c)_mu$(μ)_L$(L)_chi$(χ)", convert_to_tensormap(ψm))
end


text:
import IntervalArithmetic: radius, mid, mag, mig, sup, inf, diam, bisect, dist
export radius, mid, mag, mig, sup, inf, diam, bisect, dist

for f in (:radius, :mid, :mag, :mig, :sup, :inf, :diam)
    @eval $f(a::NumberInterval) = $f(Interval(a))
end

mid(a::NumberInterval, α) = mid(Interval(a), α)

function bisect(a::NumberInterval, α = IntervalArithmetic.where_bisect)
    c, d = bisect(Interval(a), α)
    return NumberInterval(c), NumberInterval(d)
end

function dist(a::NumberInterval, b::NumberInterval)
    return dist(Interval(a), Interval(b))
end

ext: jl


text:
using SpikeStats
using Test

testset = poissonSpikes(10,100,0.001)

println(names(SpikeStats))

@testset "SpikeStats.jl" begin
    @test typeof(testset) == Array{Float64,1}
    @test round(computeCV(diff(testset)), digits=2) == 0.93
    @test round(computeFano(testset,0.001), digits=2) == 0.90
end


text:
# samples Bernoully trials based on the simple randomn number generator

function BernoulliTrial(n, seed)

    rng = MersenneTwister(seed)
    A = [0, 1]

    return rand(rng, A, n)

end

# samples a coin toss, either win or lose

function coinTrial(n, seed)

    rng = MersenneTwister(seed)
    A = [-1, 1]

    return rand(rng, A, n)

end


text:
# Autogenerated wrapper script for libidn_jll for i686-linux-musl
export idn, libidn

JLLWrappers.@generate_wrapper_header("libidn")
JLLWrappers.@declare_executable_product(idn)
JLLWrappers.@declare_library_product(libidn, "libidn.so.12")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        idn,
        "bin/idn",
    )

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

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using LinearSpace
using Test

struct MT i::Int end
#scalartype(mt::MT) = Float64

@linearspace MT #Float64

x, y = MT(2), MT(3)

z = x + 3y

@test z.a == [1, 3]
@test z.v == [x, y]

z = 1x
@test z.a == [1]
@test z.v == [x]

@test lctype(MT) == MT_LC

ts = MT[]
as = Float64[]
z = x + 3y
for (a,t) in terms(z)
    @show a, t
    push!(as,a)
    push!(ts,t)
end

@test as == [1.0, 3.0]
@test ts == [x,y]

ts = MT[]
as = Float64[]
z = x
for (a,t) in terms(z)
    @show a, t
    push!(as,a)
    push!(ts,t)
end

@test as == [1.0]
@test ts == [x]


text:
module Vlasiator

using Requires
using StaticArrays
using Printf: @sprintf
using UnPack
using LinearAlgebra: ×, dot, ⋅, norm, I
using Statistics: mean
using EzXML
using FLoops
using Mmap: mmap
using WriteVTK
using LazyGrids: ndgrid
using LaTeXStrings
using Dates

include("utility/rotation.jl")
include("utility/log.jl")
include("utility/curvature.jl")
include("vlsv/vlsvreader.jl")
include("vlsv/vlsvutility.jl")
include("utility/plot.jl")

export
   # vlsvreader
   MetaVLSV, VarInfo,
   load, readvariable, readparameter, readvariablemeta, readvcells,
   hasvariable, hasparameter, hasname,
   # vlsvutility
   getcell, getslicecell, getlevel, refineslice, getcellcoordinates, getvcellcoordinates,
   getchildren, getparent, isparent, getsiblings,
   getcellinline, getnearestcellwithvdf, getcellwithvdf, write_vtk, issame,
   # plot helper
   SI, RE, Log, Linear, SymLog,
   # log
   readlog

precompile(load, (String,))
precompile(readvariable, (MetaVLSV, String, Bool))

function __init__()
   @require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" begin
      include("plot/pyplot.jl")
   end
   @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
      include("plot/plots.jl")
   end
end

end


text:
export unary_operators, binary_operators, nary_operators, compute
export AbstractExpression

# Define abstract base type
abstract type AbstractExpression end

# Define Algebraic Operators
include("abstract_operations.jl")

# Define Domains
include("abstract_domains.jl")

# Define Fields
include("abstract_fields.jl")

# Define Data
include("abstract_data.jl")

# Define equations and systems
include("abstract_equations.jl")

# Define Wrapper functions and derivative utils
include("abstract_utils.jl")

# compute function
compute(a::AbstractExpression) = throw(error("compute not defined for $(typeof(a))"))

# Include Generic Evaluation Rules and Output Format and Exports
for unary_operator in unary_operators
    b_name, b_symbol = Meta.parse.(unary_operator)
    @eval compute(a::$b_name{𝒮}) where {𝒮} = $b_symbol(compute(a.term))
    @eval function Base.show(io::IO, operation::$b_name{𝒮}) where {𝒮}
        print(io, $b_symbol, "(", operation.term, ")")
    end
    @eval export $b_name
end

for binary_operator in binary_operators
    b_name, b_symbol = Meta.parse.(binary_operator)
    @eval compute(a::$b_name{𝒮, 𝒯}) where {𝒮, 𝒯} = $b_symbol(compute(a.term1), compute(a.term2))
    @eval function Base.show(io::IO, operation::$b_name{𝒮, 𝒯}) where {𝒮, 𝒯}
        # print(io, "(", operation.term1, $b_symbol , operation.term2, ")")
        # clearly a great option
        color_numbers = [30:33, 65:69, 136:142, 202:207]
        choices = collect(Iterators.flatten(color_numbers))
        color = 226 # rand(choices)
        printstyled(io, "(", color = color)
        print(io, operation.term1)
        printstyled(io, $b_symbol, color = color )
        print(io,  operation.term2)
        printstyled(io, ")", color = color)
    end
    @eval export $b_name
end

for nary_operator in nary_operators
    b_name, b_symbol = Meta.parse.(nary_operator)
    @eval compute(a::$b_name{𝒮}) where {𝒮} = $b_symbol(compute.(a.terms)...)
    @eval function Base.show(io::IO, operation::$b_name{𝒮}) where {𝒮}
        print(io, $b_symbol, "(" )
        for term in operation.terms
            print(io, term, ",")
        end
        print(io, "0)")
    end
    @eval export $b_name
end


text:
export Arc

mobius(a,b,c,d,z) = (a*z+b)/(c*z+d)
mobius(a,b,c,d,z::Number) = isinf(z) ? a/c : (a*z+b)/(c*z+d)
mobiusinv(a,b,c,d,z) = mobius(d,-b,-c,a,z)

mobiusD(a,b,c,d,z) = (a*d-b*c)/(d+c*z)^2
mobiusinvD(a,b,c,d,z) = mobiusD(d,-b,-c,a,z)
"""
    Arc(c,r,(θ₁,θ₂))

represents the arc centred at `c` with radius `r` from angle `θ₁` to `θ₂`.
"""
struct Arc{T,V<:Real,TT} <: SegmentDomain{TT}
    center::T
    radius::V
    angles::Tuple{V,V}
    Arc{T,V,TT}(c,r,a) where {T,V,TT} = new{T,V,TT}(convert(T,c),convert(V,r),Tuple{V,V}(a))
end

Arc(c::T,r::V,t::Tuple{V1,V2}) where {T<:Number,V<:Real,V1<:Real,V2<:Real} =
    Arc{promote_type(T,V,V1,V2),
        promote_type(real(T),V,V1,V2),
        Complex{promote_type(real(T),V,V1,V2)}}(c,r,t)
Arc(c::Vec{2,T},r::V,t::Tuple{V,V}) where {T<:Number,V<:Real} =
    Arc{Vec{2,promote_type(T,V)},
        promote_type(real(T),V),
        Vec{2,promote_type(real(T),V)}}(c,r,t)
Arc(c::Vec{2,T},r::V,t::Tuple{V1,V2}) where {T<:Number,V<:Real,V1<:Real,V2<:Real} =
    Arc{Vec{2,promote_type(T,V,V1,V2)},
        promote_type(real(T),V,V1,V2),
        Vec{2,promote_type(T,V,V1,V2)}}(c,r,t)

Arc(c::Tuple,r,t) = Arc(Vec(c...),r,t)
Arc(c,r,t0,t1) = Arc(c,r,(t0,t1))

complex(a::Arc{V}) where {V<:Vec} = Arc(complex(a.center...),a.radius,a.angles)

isambiguous(d::Arc) =
    isnan(d.center) && isnan(d.radius) && isnan(d.angles[1]) && isnan(d.angles[2])
convert(::Type{Arc{T,V}},::AnyDomain) where {T<:Number,V<:Real} =
    Arc{T,V}(NaN,NaN,(NaN,NaN))
convert(::Type{IT},::AnyDomain) where {IT<:Arc} =
    Arc(NaN,NaN,(NaN,NaN))

isempty(d::Arc) = false

reverseorientation(a::Arc) = Arc(a.center,a.radius,reverse(a.angles))

arclength(d::Arc) = d.radius*(d.angles[2]-d.angles[1])

function mobiuspars(z0,r,t0,t1)
    c=exp(im*t0/2)+exp(im*t1/2)
    f=exp(im*t0/2)-exp(im*t1/2)
    d=exp(im*t0/2+im*t1/2)*r
    -c,c*(z0+d),f,f*(d-z0)
end

mobiuspars(a::Arc) = mobiuspars(a.center,a.radius,a.angles...)

for OP in (:mobius,:mobiusinv,:mobiusD,:mobiusinvD)
    @eval $OP(a::Arc,z) = $OP(mobiuspars(a)...,z)
end

tocanonical(a::Arc{T},x) where {T<:Number} = real(mobius(a,x))
tocanonicalD(a::Arc{T},x) where {T<:Number} = mobiusD(a,x)
fromcanonical(a::Arc{T,V,TT},x) where {T<:Number,V<:Real,TT<:Complex} =
    mobiusinv(a,x)
fromcanonicalD(a::Arc{T},x) where {T<:Number} = mobiusinvD(a,x)

tocanonical(a::Arc{V},x::Vec) where {V<:Vec} =
    tocanonical(complex(a),complex(x...))
fromcanonical(a::Arc{V},x::Number) where {V<:Vec} =
    Vec(reim(fromcanonical(complex(a),x))...)
fromcanonicalD(a::Arc{V},x::Number) where {V<:Vec} =
    Vec(reim(fromcanonicalD(complex(a),x))...)

## information

issubset(a::Arc,b::Arc) =
    a.center==b.center && a.radius==b.radius &&
                    (b.angles[1]≤a.angles[1]≤a.angles[2]≤b.angles[2] ||
                            b.angles[1]≥a.angles[1]≥a.angles[2]≥b.angles[2])

# Algebra

*(c::Real,d::Arc) =
    Arc(c*d.center,c*d.radius,(sign(c)*d.angles[1],sign(c)*d.angles[2]))
*(d::Arc,c::Real) = *(c,d)

for op in (:+,:-)
    @eval begin
        $op(c::Number,d::Arc) = Arc($op(c,d.center),d.radius,d.angles)
        $op(d::Arc,c::Number) = Arc($op(d.center,c),d.radius,d.angles)
    end
end

# allow exp(im*Segment(0,1)) for constructing arc
function exp(d::IntervalOrSegment{<:Complex})
    @assert isapprox(real(leftendpoint(d)),0) && isapprox(real(rightendpoint(d)),0)
    Arc(0,1,(imag(leftendpoint(d)),imag(rightendpoint(d))))
end


text:
function module_init()
    # Locate the directory in which the data files are located
    this_filename_dir = dirname(@__FILE__())
    data_dir = joinpath(this_filename_dir, "..", "data")

    # Verify that the CORRCOEF file exists
    corrcoef_filename = joinpath(data_dir, "CORRCOEF")
    @assert isfile(corrcoef_filename)

    # Verify that the EGM96 file exists
    egm96_filename = joinpath(data_dir, "EGM96")
    @assert isfile(egm96_filename)

    # Initialize the f477 C library
    ccall((:init_arrays, libf477),
        Void,
        (Cstring, Cstring),
        corrcoef_filename,
        egm96_filename)
end

# Initialize the module, reading the model data from disk.
module_init()

function undulation(lat::Number, lon::Number)
    return ccall((:undulation, libf477),
        Float64,
        (Float64, Float64, Int, Int),
        deg2rad(lat), deg2rad(lon), 360, 361)
end


text:
using Test
using Logic: unify, @var, @predicate, Variable, Success, Failure, Clause, Conj, Disj

@test Success() == Success()
@test unify(1, 1) == Success()
@test unify(1, 2) == Failure(1, 2)

@var X::Int
@var Y::String
@var Z::String
@var T::Any
@test unify(X, 1) == unify(1, X) == Success(Dict(X => 1))
@test unify(X, X) == Success(Dict(X => X))
@test unify(X, Y) == Failure(X, Y)
@test unify(Y, Z) == Success(Dict(Y => Z))

@predicate parentp(p::Variable{String}, c::Variable{String})

@assert unify(parentp("Abraham", "Isaac"), parentp(Y, Z)) == Success(Dict(Y => "Abraham", Z => "Isaac"))
@assert unify(parentp("Abraham", Z), parentp(Y, "Isaac")) == Success(Dict(Y => "Abraham", Z => "Isaac"))

@predicate can(verb)
@predicate source(node, port)
@test unify(can(source("flask1", X)), can(T)) == Success(Dict(T => source("flask1", X)))
@test unify(can(source("flask1", "0")), can(source(Y, Z))) == Success(Dict(Y => "flask1", Z => "0"))

@predicate male(name::Variable{String})
@predicate father(f::Variable{String}, c::Variable{String})
@show unify(father("Abraham", "Isaac"), Clause(father(Y, Z), Conj(parentp(Y, Z), male(Y))))


text:
"""
This example uses the ProfileData struct and GP.
"""

include("GP1.jl")

log_γs = -3.0:0.1:3.0

V_name = Dict("T" =>"Temperature [°C]", "wT"=>"Temperature flux [°C⋅m/s]")

v_str = "T"
# v_str = "wT"

N = 4
# N = 2

##

p = me_file_comparison2(2, l2_norm, 4.8:0.01:6.0, v_str; N=2)
savefig(pwd() * "/hyperparameter_landscapes/rescaled_$(v_str)/compare_sim_SE_γ_landscapes_l2norm_k2_neighborhood_cscrapr_N2.png")

##
#  *--*--*--*--*--*--*
#  | Animations      |
#  *--*--*--*--*--*--*

function animate_profile_trained_on_gs16(filename, k::Int64, γ, d, D)

    mydata = construct_profile_data(filename, v_str, D; N=N)

    # find the minimizing gamma value then animate
    # min_gamma, min_error = get_min_gamma(2, data, normalize, l2_norm, log_γs)
    kernel = get_kernel(k, γ, 0.0)

    data16 = construct_profile_data("general_strat_16_profiles.jld2", v_str, D; N=N, verbose=true)
    𝒢 = construct_gpr(data16.x_train, data16.y_train, kernel; distance_fn=d, z=data16.zavg);

    gpr_prediction = get_gpr_pred(𝒢, mydata)

    animation_set = 1:10:(mydata.Nt-2)
    anim = @animate for i in animation_set
        exact = mydata.v[:,i]
        day_string = string(floor(Int, mydata.t[i]/86400))
        p1 = scatter(gpr_prediction[i], mydata.zavg, label = "GP")
        # xlims=(minimum(data.v[:,1]),maximum(data.v[:,1]))
        xlims=(18,20)

        if i<data16.Nt
            exact16 = data16.v[:,i]
        else
            exact16 = data16.v[:,data16.Nt]
        end

        plot!(exact, mydata.z, legend = :topleft, label = "LES", xlabel = "$(V_name["T"])", ylabel = "depth", title = "i = $(i)", xlims=xlims)
        plot!(exact16, data16.z, legend = :topleft, label = "LES gs 16", xlabel = "$(V_name["T"])", ylabel = "depth", xlims=xlims)

    end

    return anim
end

filename="general_strat_32_profiles"
anim = animate_profile_trained_on_gs16("$(filename).jld2", 2, 0.3, l2norm_strat_penalty, 16)
gif(anim, pwd() * "/../les/data_sandreza/$(filename)/gp_γ4.1_M12_l2norm_strat_penalty_trainedOnGs8_16_32_scaled.gif", fps = 20)

##

filename="general_strat_32_profiles"
anim = animate_profile("$(filename).jld2", 2, 5.50, l2_norm, 16, "T"; N=2)
gif(anim, pwd() * "/../les/data_sandreza/$(filename)/gp_γ4.1_M12_l2norm_γ5.5_N2_scaled.gif", fps = 20)

##
#  *--*--*--*--*--*--*--*--*--*
#  | Error metric comparison  |
#  *--*--*--*--*--*--*--*--*--*

filename = "general_strat_16_profiles.jld2"
data = construct_profile_data(filename, v_str, 16; N=N)
log_γs = -3.0:0.1:3.0 # hyperparameter slider range

for k in 1:5
    #with normalization
    p = error_metric_comparison(k, data, l2_norm, log_γs, true)
    savefig(pwd() * "/hyperparameter_landscapes/rescaled_$(v_str)/kernel$(k)_γ_gs16_l2norm.png")
    # #no normalization
    # p = error_metric_comparison(k, data, l2_norm, log_γs, false)
    # savefig(pwd() * "/hyperparameter_landscapes/no_normalization_$(v_str)/kernel$(k)_γ_gs16_l2norm.png")
end

##
#  *--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*
#  | ME on true check for gs8 vs gs16 vs gs32: fixed range of γs  |
#  *--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*

# ------- L2 norm -------
for k in 1:5

    # if k==2; log_γs = 0.0:0.1:6.0
    # else; log_γs = -3.0:0.1:3.0
    # end

    log_γs = -3.0:0.1:3.0

    #with normalization
    p = me_file_comparison(k, l2_norm, log_γs, v_str)
    savefig(pwd() * "/hyperparameter_landscapes/rescaled_$(v_str)/compare_sims_kernel$(k)_γ_l2norm.png")
    # #no normalization
    # p = me_file_comparison(k, l2_norm, log_γs, v_str; normalize=false)
    # savefig(pwd() * "/hyperparameter_landscapes/rescaled_$(v_str)/compare_sims_kernel$(k)_γ_l2norm.png")
end
# ------- H1 norm -------
for k in 1:5

    # if k==2; log_γs = 0.0:0.1:6.0
    # else; log_γs = -3.0:0.1:3.0
    # end

    log_γs = -3.0:0.1:3.0

    #with normalization
    p = me_file_comparison(k, h1_norm, log_γs, v_str)
    savefig(pwd() * "/hyperparameter_landscapes/rescaled_$(v_str)/compare_sims_kernel$(k)_γ_h1norm.png")
    # #no normalization
    # p = me_file_comparison(k, h1_norm, log_γs, v_str; normalize=false)
    # savefig(pwd() * "/hyperparameter_landscapes/no_normalization_$(v_str)/compare_sims_kernel$(k)_γ_h1norm.png")
end

# ------- H^-1 norm -------
for k in 1:5

    # if k==2; log_γs = 0.0:0.1:6.0
    # else; log_γs = -3.0:0.1:3.0
    # end

    log_γs = -3.0:0.1:3.0

    #with normalization
    p = me_file_comparison(k, hm1_norm, log_γs, v_str)
    savefig(pwd() * "/hyperparameter_landscapes/rescaled_$(v_str)/compare_sims_kernel$(k)_γ_hm1norm.png")
    # #no normalization
    # p = me_file_comparison(k, hm1_norm, log_γs, v_str; normalize=false)
    # savefig(pwd() * "/hyperparameter_landscapes/no_normalization_$(v_str)/compare_sims_kernel$(k)_γ_hm1norm.png")
end

##
#  *--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*
#  | ME on true check for gs8 vs gs16 vs gs32: neighborhood of the min  |
#  *--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*
# custom γ ranges
# l2_norm

p = me_file_comparison(1, l2_norm, 3.8:0.001:4.3, v_str; normalize=true)
savefig(pwd() * "/plots/hyperparameter_landscapes_constlims_with_normalization_$(v_str)/compare_les_SE_γ_landscapes_l2norm.png")

p = me_file_comparison(2, l2_norm, 3.8:0.001:4.3, v_str; normalize=true)
savefig(pwd() * "/plots/hyperparameter_landscapes_constlims_with_normalization_$(v_str)/compare_les_M12_γ_landscapes_l2norm.png")

p = me_file_comparison(3, l2_norm, -0.3:0.001:0.3, v_str; normalize=true)
savefig(pwd() * "/plots/hyperparameter_landscapes_constlims_with_normalization_$(v_str)/compare_les_M32_γ_landscapes_l2norm.png")

p = me_file_comparison(4, l2_norm, -0.3:0.01:0.3, v_str; normalize=true)
savefig(pwd() * "/plots/hyperparameter_landscapes_constlims_with_normalization_$(v_str)/compare_les_M52_γ_landscapes_l2norm.png")

p = me_file_comparison(5, l2_norm, -0.2:0.001:0.2, v_str; normalize=true)
savefig(pwd() * "/plots/hyperparameter_landscapes_constlims_with_normalization_$(v_str)/compare_les_RQ_α1_γ_landscapes_l2norm.png")

## compare distance metrics

# function mes_plot_file_comparison(k::Int64, filename, normalize, d, γs)
#     # Mean error on greedy check correlated with mean error on true check?
#
#     results = Dict()
#     for dm in [l2_norm, h1_norm, hm1_norm]
#
#         𝒟 = construct_profile_data(file, "T", 16; N=N)
#         mets  = zeros(length(γs)) # mean error (true check)
#
#         for i in 1:length(γs)
#             kernel = get_kernel(k, γs[i], 0.0)
#             𝒢 = construct_gpr(𝒟.x_train, 𝒟.y_train, kernel; distance_fn=dm, z=𝒟.zavg, normalize=normalize);
#             # -----compute mean error for true check----
#             total_error = 0.0
#             gpr_prediction = get_gpr_pred(𝒢, 𝒟)
#             for i in 1:𝒟.Nt-2
#                 exact    = 𝒟.y[i+1]
#                 predi    = gpr_prediction[i+1]
#                 total_error += l2_norm(exact, predi) # euclidean distance
#             end
#             mets[i] = total_error/(𝒟.Nt-2)
#         end
#
#         results[dm]=mets
#     end
#
#     r1 = results[l2_norm]
#     γ=γs[argmin(r1)]
#     p1 = plot(γs, r1, xlabel="log(γ)", ylabel="ME, true check", title="l² norm, log(γ)=$(γ), min = $(round(minimum(r1);digits=5))", legend=false, yscale=:log10) # 1D plot: mean log marginal loss vs. γ
#     vline!([γ])
#
#     r2 = results[h1_norm]
#     γ=γs[argmin(r2)]
#     p2  = plot(γs, r2,  xlabel="log(γ)", ylabel="ME, true check", title="H¹ norm, log(γ)=$(γ), min = $(round(minimum(r2);digits=5))", legend=false, yscale=:log10)  # 1D plot: mean error vs. γ
#     vline!([γ])
#
#     r3 = results[hm1_norm]
#     γ=γs[argmin(r3)]
#     p3  = plot(γs, r3,  xlabel="log(γ)", ylabel="ME, true check", title="H⁻¹ norm, log(γ)=$(γ), min = $(round(minimum(r3);digits=5))", legend=false, yscale=:log10)  # 1D plot: mean error vs. γ
#     vline!([γ])
#
#     return plot(p1, p2, p3, layout = @layout [a ; b; c])
# end

## automate finding gamma that minimizes the mean error

# filename = "general_strat_16_profiles.jld2"
v_str = "wT"
N = 2

filename = "general_strat_16_profiles.jld2"
data = construct_profile_data(filename, v_str, D; N=N)

normalize = true

for kern in 1:5
    for d in [l2_norm, h1_norm, hm1_norm]

        # find the minimizing gamma value then animate
        min_gamma, min_error = get_min_gamma(kern, data, normalize, d, -0.3:0.1:0.3)
        kernel = get_kernel(kern, min_gamma, 0.0)
        𝒢 = construct_gpr(data.x_train, data.y_train, kernel; distance_fn=d, z=data.zavg, normalize=normalize);
        gpr_prediction = get_gpr_pred(𝒢, data)

        animation_set = 1:30:(data.Nt-2)
        anim = @animate for i in animation_set

            exact = data.v[:,i+1]
            day_string = string(floor(Int, data.t[i]/86400))
            p1 = scatter(gpr_prediction[i+1], data.zavg, label = "GP")
            plot!(exact, data.z, legend = :topleft, label = "LES", xlabel = "$(V_name[v_str])", ylabel = "depth", title="gs32, k=$(kern), $(d), log(γ)=$(min_gamma), error=$(round(min_error;digits=7)) day " * day_string, xlims = (-1e-5,4e-5))
            display(p1)

        end

        # save animation
        gif(anim, pwd() * "/../les/data_sandreza/$(filename[1:end-5])/gp_min_gamma_wnormalize_γ$(min_gamma)_k$(kern)_d$(d)_$(v_str)_N$(N).gif", fps = 15)

    end
end

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
"""
Cosine taper a N-dimensional array along given dimensions.

# Arguments
* `x::Array{Float64,N}` :
* `perc::Float64` : taper percentage
"""
function taper(x::AbstractArray,perc::Float64;dims=[1])
	xout=copy(x);
	taper!(xout,perc,dims=dims)
	return xout
end

"""
# no allocations
* `dims=[1]` : taper along these dimensions (default only first)
"""
function taper!(x, perc=0.0; bperc=perc, eperc=perc, dims=[1])
	for dim in dims
		nt=size(x,dim)
		nttb=min(round(Int,nt*bperc*0.01), nt)
		ntte=min(round(Int,nt*eperc*0.01), nt)
		kb=inv(2.0*round(Int,nt*bperc*0.01)-1)*pi
		ke=inv(2.0*round(Int,nt*eperc*0.01)-1)*pi
		for i in CartesianIndices(size(x))
			if(1≤i[dim]≤nttb)
				x[i] *= sin((i[dim]-1)*kb)
			end
			if(nt-ntte+1≤i[dim]≤nt)
				x[i] *= sin((-i[dim]+nt)*ke)
			end
		end
	end
end

"""
Tapering is necessary to be able to input random signal into finite-difference code
Filtering tapering are applied only if the length of the time series is greater than 10
"""
function get_tapered_random_tmax_signal(tgrid;
					fmin=nothing,
					fmax=nothing,
					tmaxfrac::Float64=1.0,
					dist=Uniform(-2.0, 2.0),
					sparsep=1.0,
					taperperc=20.
					)
	filt_flag=(length(tgrid) > 5) && (!(fmin === nothing)) && (!(fmax===nothing))

	fs = inv(step(tgrid));
	if(filt_flag)
		designmethod = Butterworth(6);
		filtsource = Bandpass(fmin, fmax; fs=fs);
	end

	itind = argmin(abs.(tgrid.-abs(tmaxfrac)*tgrid[end]))
	if(tmaxfrac>0.0)
		its=1:itind
	elseif(tmaxfrac<0.0)
		its=itind:length(tgrid)
	end
	# 20% taper window
	twin = taper(ones(length(its)),taperperc)
	X = zeros(length(its))
	wavsrc = zeros(length(tgrid))
	if(filt_flag)
		X[:] = rand(dist, length(its)) .* twin
	else
		X[:] = rand(dist, length(its))
	end
	if(sparsep ≠ 1.0)
		Xs=sprandn(length(X), sparsep)
		X[findall(Xs.==0.0)].=0.0
	end
	# band limit
	(filt_flag) && (filt!(X, digitalfilter(filtsource, designmethod), X))

	(length(X) ≠ 1) && normalize!(X)
	wavsrc[its] = X
	return wavsrc
end


text:
using PyPlot
using DelimitedFiles

R = readdlm("result.txt")
n = [1, 4, 9, 16, 25]

close("all")
loglog(n[end:-1:1], R[1:5, 3], "or-", label = "1 Core")
loglog(n[4:-1:1], R[6:9, 3], "xg-", label = "2 Cores")
loglog(n[3:-1:1], R[10:12, 3], "^b-", label = "3 Cores")
loglog(n[end:-1:1], R[1:5, 5], "or--")
loglog(n[4:-1:1], R[6:9, 5], "xg--")
loglog(n[3:-1:1], R[10:12, 5], "^b--")
loglog(n[end-2:-1:1], 8 ./n[end-2:-1:1], "k--", label = "\$\\mathcal{O}(1/n_{{proc}})\$")
legend()
title("Strong Scaling for Forward Computation")
xlabel("MPI Processors")
ylabel("Time (sec)")
savefig("forward_strong.png")

close("all")
loglog(n[end:-1:1], R[1:5, 4], "or-", label = "1 Core")
loglog(n[4:-1:1], R[6:9, 4], "xg-", label = "2 Cores")
loglog(n[3:-1:1], R[10:12, 4], "^b-", label = "3 Cores")
loglog(n[end:-1:1], R[1:5, 6], "or--")
loglog(n[4:-1:1], R[6:9, 6], "xg--")
loglog(n[3:-1:1], R[10:12, 6], "^b--")
loglog(n[end-2:-1:1], 12 ./n[end-2:-1:1], "k--", label = "\$\\mathcal{O}(1/n_{{proc}})\$")
legend()
title("Strong Scaling for Gradient Backpropagation")
xlabel("MPI Processors")
ylabel("Time (sec)")
savefig("backward_strong.png")


text:
# Functions that are useful for testing
#using Test, JLD2, D4M, PyPlot, SparseArrays
using Test, D4M, JLD, PyPlot, SparseArrays
#include("loadD4M.jl")

import Base.isequal

UnionArray = Array{Union{AbstractString,Number}}

# Check whether two associative arrays are "equal"
function isequal(A::Assoc,B::Assoc)
    pass = true
    why = ""
    rA,cA,vA = find(A)
    rB,cB,vB = find(B)
    adjA = A.A
    adjB = B.A

    # Check sizes of r,c,v match
    if ~(size(rA,1)==size(rB,1) && size(rA,2)==size(rB,2))
        pass = false
        why = why*"Row key sizes don't match:\n\tsize(rA): "*string(size(rA))*"\n\tsize(rB):"*string(size(rB))*"\n"
    end
    if ~(size(cA)==size(cB) && size(rA,2)==size(rB,2))
        pass = false
        why = why*"Column key sizes don't match:\n\tsize(cA): "*string(size(cA))*"\n\tsize(cB):"*string(size(cB))*"\n"
    end
    if ~(size(vA)==size(vB) && size(rA,2)==size(rB,2))
        pass = false
        why = why*"Value sizes don't match:\n\tsize(vA): "*string(size(vA))*"\n\tsize(vB):"*string(size(vB))*"\n"
    end

    # Check r,c,v match
    if ~all(rA .== rB)
        pass = false
        why = why*"Row keys don't match\n"
    end
    if ~all(cA .== cB)
        pass = false
        why = why*"Column keys don't match\n"
    end
    if ~all(vA .== vB)
        pass = false
        why = why*"Values don't match\n"
    end

    # Check sparse matrices match
    if ~(adjA == adjB)
        pass = false
        why = why*"Sparse matrices don't match\n"
    end

    if ~pass
        println(why)
    end

    return pass
end

# Type and size testing of Assoc and its components
function testassoc(A)
    pass = true
    why = ""

    # Check is A is correct data type
    if ~isa(A,Assoc)
        pass = false
        why = why*"Assoc wrong type: "*string(typeof(A))*"\n"
    end

    # Check row,col,val are correct data type
    if ~isa(A.row,UnionArray)
        pass = false
        why = why*"Row wrong type: "*string(typeof(A.row))*"\n"
    end
    if ~isa(A.col,UnionArray)
        pass = false
        why = why*"Col wrong type: "*string(typeof(A.col))*"\n"
    end
    if ~isa(A.val,UnionArray)
        pass = false
        why = why*"Val wrong type: "*string(typeof(A.val))*"\n"
    end

    # Check sparse matrix is right data type
    if ~isa(A.A,SparseMatrixCSC)
        pass = false
        why = why*"Sparse matrix wrong type: "*string(typeof(A.A))*"\n"
    end

    # Check size of row and col arrays match dimensions of sparse matrix
    if ~length(A.row) == size(A.A,1)
        pass = false
        why = why*"Row length does not match sparse matrix dimension 1:\n\tlength(A.row): "*string(length(A.row))*"\n\tsize(A.A,1): "*string(size(A.A,1))*"\n"
    end
    if ~length(A.col) == size(A.A,2)
        pass = false
        why = why*"Col length does not match sparse matrix dimension 2:\n\tlength(A.col): "*string(length(A.col))*"\n\tsize(A.A,2): "*string(size(A.A,2))*"\n"
    end

    # Check that numeric values are in sparse matrix, not in A.val
    if ~isempty(A) && ~(!all(isa.(A.val,Number)) || A.val == [1.0])
        pass = false
        why = why*"Numeric values in wrong place"
    end

    if ~pass
        println(why)
    end

    return pass
end

function testassoc(A,B)
    pass = true
    pass = pass && testassoc(A)
    pass = pass && isequal(A,B)
end

# Run Tests
@testset "All Tests" begin
    #include("loadD4M.jl")
    @testset "Intro Tests" begin
        intro_res = load("test_assocs/1intro.jld")["intro_res"]
        @testset "Assoc Intro" begin
            include("../examples/1Intro/1AssocIntro/AI1_SetupTEST.jl")
            @test testassoc(A,intro_res["AI1_setup"])

            include("../examples/1Intro/1AssocIntro/AI2_SubsrefTEST.jl")
            @test testassoc(A1r,intro_res["AI2_subsref"]["A1r"])
            @test testassoc(A3r,intro_res["AI2_subsref"]["A3r"])
            @test testassoc(A4r,intro_res["AI2_subsref"]["A4r"])
            @test testassoc(A1c,intro_res["AI2_subsref"]["A1c"])
            @test testassoc(A3c,intro_res["AI2_subsref"]["A3c"])
            @test testassoc(A4c,intro_res["AI2_subsref"]["A4c"])

            include("../examples/1Intro/1AssocIntro/AI3_MathTEST.jl")
            @test testassoc(A,intro_res["AI3_math"]["A"])
            @test testassoc(sum(A,1),intro_res["AI3_math"]["sum1"])
            @test testassoc(sum(A,2),intro_res["AI3_math"]["sum2"])
            @test testassoc(Aab,intro_res["AI3_math"]["Aab"])
            @test testassoc(F,intro_res["AI3_math"]["F"])
            @test testassoc(Fn,intro_res["AI3_math"]["Fn"])
            @test testassoc(AtA,intro_res["AI3_math"]["AtA"])

            include("../examples/1Intro/1AssocIntro/AI4_AdvConstructTEST.jl")
            @test testassoc(A00,intro_res["AI4_construct"]["A00"])
            @test testassoc(A01,intro_res["AI4_construct"]["A01"])
            @test testassoc(A02,intro_res["AI4_construct"]["A02"])
            @test testassoc(A03,intro_res["AI4_construct"]["A03"])
            @test testassoc(A04,intro_res["AI4_construct"]["A04"])
            @test testassoc(A05,intro_res["AI4_construct"]["A05"])
            @test testassoc(A06,intro_res["AI4_construct"]["A06"])
            @test testassoc(A11,intro_res["AI4_construct"]["A11"])
            @test testassoc(A12,intro_res["AI4_construct"]["A12"])
            @test testassoc(A13,intro_res["AI4_construct"]["A13"])
            @test testassoc(A14,intro_res["AI4_construct"]["A14"])
            @test testassoc(A15,intro_res["AI4_construct"]["A15"])
            @test testassoc(A16,intro_res["AI4_construct"]["A16"])

        end
        @testset "Edge Art" begin
            include("../examples/1Intro/2EdgeArt/EA1_GraphTEST.jl")
            @test testassoc(Ev,intro_res["EA1_graph"]["Ev"])
            @test testassoc(Av,intro_res["EA1_graph"]["Av"])
            @test testassoc(Ae,intro_res["EA1_graph"]["Ae"])

            include("../examples/1Intro/2EdgeArt/EA2_SubsrefTEST.jl")
            @test testassoc(Eo,intro_res["EA2_subsref"]["Eo"])
            @test testassoc(Eog,intro_res["EA2_subsref"]["Eog"])

            include("../examples/1Intro/2EdgeArt/EA3_SubGraphTEST.jl")
            @test testassoc(EvO,intro_res["EA3_subgraph"]["EvO"])
            @test testassoc(EvG,intro_res["EA3_subgraph"]["EvG"])
            @test testassoc(AvOG,intro_res["EA3_subgraph"]["AvOG"])
            #@test testassoc(AeOG2)
            #@test testassoc(AeOG2,intro_res["EA3_subgraph"]["AeOG2"])

        end
    end
    @testset "Apps" begin
        @testset "Entity Analysis" begin
            include("../examples/2Apps/1EntityAnalysis/EA1_ReadTEST.jl")
            apps_EA1 = load("test_assocs/2apps_1EA1.jld")["apps_EA1"]
            @test testassoc(E,apps_EA1["E"])
            include("../examples/2Apps/1EntityAnalysis/EA2_StatTEST.jl")
            apps_EA2 = load("test_assocs/2apps_1EA2.jld")["apps_EA2"]
            @test testassoc(sum(logical(col2type(E,"/")),1),apps_EA2["Ent"])
            @test testassoc(En,apps_EA2["En"])
            @test testassoc(An,apps_EA2["An"])
            @test sum(getadj(An[:,StartsWith("LOCATION/,")]),dims=2) == apps_EA2["plot"]
            include("../examples/2Apps/1EntityAnalysis/EA3_FacetTEST.jl")
            apps_EA3 = load("test_assocs/2apps_1EA3.jld")["apps_EA3"]
            @test testassoc(F,apps_EA3["F"])
            @test testassoc(Fn,apps_EA3["Fn"])
            @test testassoc(F' > 1,apps_EA3["F2"])
            @test testassoc(Fn' > 0.02,apps_EA3["Fn2"])
            include("../examples/2Apps/1EntityAnalysis/EA4_GraphTEST.jl")
            apps_EA4 = load("test_assocs/2apps_1EA4.jld")["apps_EA4"]
            @test testassoc(Ae,apps_EA4["Ae"])
            @test testassoc(Ep,apps_EA4["Ep"])
            @test testassoc(Ap,apps_EA4["Ap"])
            @test testassoc(Ad,apps_EA4["Ad"])
            include("../examples/2Apps/1EntityAnalysis/EA5_GraphQueryTEST.jl")
            apps_EA5 = load("test_assocs/2apps_1EA5.jld")["apps_EA5"]
            @test testassoc(A,apps_EA5["A"])
            @test testassoc(An,apps_EA5["An"])
            @test testassoc((A[p,x] > 4) & (An[p,x] > 0.3),apps_EA5["Amf"])
            @test testassoc(A[p1,p1],apps_EA5["Ap1"])
            @test testassoc(A[p2,p2] > 1,apps_EA5["Ap2"])
        end
        PyPlot.close()
        PyPlot.close()
        PyPlot.close()
        PyPlot.close()
        PyPlot.close()
        PyPlot.close()
        @testset "Track Analysis" begin
            include("../examples/2Apps/2TrackAnalysis/findtrackgraph.jl")
            include("../examples/2Apps/2TrackAnalysis/findtracks.jl")
            include("../examples/2Apps/2TrackAnalysis/TA1_BuildTEST.jl")
            apps_TA1 = load("test_assocs/2apps_2TA1.jld")["apps_TA1"]
            @test testassoc(E3,apps_TA1["E3"])
            @test testassoc(Etx,apps_TA1["Etx"])
            @test testassoc(Ext,apps_TA1["Ext"])
            @test testassoc(At,apps_TA1["At"])
            @test testassoc(Ax,apps_TA1["Ax"])
            include("../examples/2Apps/2TrackAnalysis/TA2_QueryTEST.jl")
            apps_TA2 = load("test_assocs/2apps_2TA2.jld")["apps_TA2"]
            @test testassoc(A,apps_TA2["A"])
            @test testassoc(A[:,p1*p2],apps_TA2["Ap1p2"])
            @test testassoc(A[t,:] == x,apps_TA2["At"])
            include("../examples/2Apps/2TrackAnalysis/TA3_GraphTEST.jl")
            apps_TA3 = load("test_assocs/2apps_2TA3.jld")["apps_TA3"]
            @test testassoc(A,apps_TA3["A"])
            @test testassoc(G,apps_TA3["G"])
            @test testassoc(G > 5,apps_TA3["G5"])
            @test testassoc(Go,apps_TA3["Go"])
            @test testassoc((Go > 2) & ((Go ./ G) > 0.2),apps_TA3["GoG"])
        end
    end
    PyPlot.close()
    PyPlot.close()
    PyPlot.close()
    @testset "Scaling" begin
        @testset "Parallel Database" begin
            include("../examples/3Scaling/2ParallelDatabase/KronGraph500NoPerm.jl")
            println("pDB01_DataTEST.jl")
            include("../examples/3Scaling/2ParallelDatabase/pDB01_DataTEST.jl")
            PyPlot.close()
            PyPlot.close()
            println("pDB02_FileTEST.jl")
            include("../examples/3Scaling/2ParallelDatabase/pDB02_FileTEST.jl")
            println("pDB03_AssocTEST.jl")
            include("../examples/3Scaling/2ParallelDatabase/pDB03_AssocTEST.jl")
            @test testassoc(testA)
            @test testassoc(testE)
            println("pDB04_DegreeTEST.jl")
            include("../examples/3Scaling/2ParallelDatabase/pDB04_DegreeTEST.jl")
            PyPlot.close()
            PyPlot.close()
            @test testassoc(Aall)
            #=
            @testset "Using Database" begin
                println("pDB05_DBsetupTEST.jl")
                include("../examples/3Scaling/2ParallelDatabase/pDB05_DBsetupTEST.jl")
                println("pDB06_AdjInsertTEST.jl")
                include("../examples/3Scaling/2ParallelDatabase/pDB06_AdjInsertTEST.jl")
                println("pDB07_AdjQueryTEST.jl")
                include("../examples/3Scaling/2ParallelDatabase/pDB07_AdjQueryTEST.jl")
                PyPlot.close()
                @test testassoc(Adeg)
                @test testassoc(A)
                println("pDB08_AdjQueryItTEST.jl")
                include("../examples/3Scaling/2ParallelDatabase/pDB08_AdjQueryItTEST.jl")
                println("pDB09_AdjJoinTEST.jl")
                include("../examples/3Scaling/2ParallelDatabase/pDB09_AdjJoinTEST.jl")
                println("pDB10_EdgeInsertTEST.jl")
                include("../examples/3Scaling/2ParallelDatabase/pDB10_EdgeInsertTEST.jl")
                println("pDB11_EdgeQueryTEST.jl")
                include("../examples/3Scaling/2ParallelDatabase/pDB11_EdgeQueryTEST.jl")
                println("pDB12_EdgeQueryItTEST.jl")
                include("../examples/3Scaling/2ParallelDatabase/pDB12_EdgeQueryItTEST.jl")
                println("pDB13_EdgeJoinTEST.jl")
                include("../examples/3Scaling/2ParallelDatabase/pDB13_EdgeJoinTEST.jl")
                PyPlot.close()
            end
            =#
        end

        PyPlot.close()
        PyPlot.close()
        PyPlot.close()
        @testset "Matrix Performance" begin
            println("MP1_DenseTEST.jl")
            include("../examples/3Scaling/3MatrixPerformance/MP1_DenseTEST.jl")
            println("MP2_SparseTEST.jl")
            include("../examples/3Scaling/3MatrixPerformance/MP2_SparseTEST.jl")
            println("MP3_AssocTEST.jl")
            include("../examples/3Scaling/3MatrixPerformance/MP3_AssocTEST.jl")
            println("MP4_AssocCatKeyTEST.jl")
            include("../examples/3Scaling/3MatrixPerformance/MP4_AssocCatKeyTEST.jl")
            println("MP5_AssocCatValKeyTEST.jl")
            include("../examples/3Scaling/3MatrixPerformance/MP5_AssocCatValKeyTEST.jl")
            println("MP6_AssocPlus.jl")
            include("../examples/3Scaling/3MatrixPerformance/MP6_AssocPlus.jl")
        end
    end
end


text:
@testset "Xor" begin
@testset "Xor basic all solutions" begin
    m = Model(optimizer_with_attributes(
        CS.Optimizer,
        "all_solutions" => true,
        "logging" => [],
    ))
    @variable(m, 1 <= x <= 5, Int)
    @variable(m, 1 <= y <= 5, Int)
    @constraint(m, (x <= 2) ⊻ (y <= 2))
    optimize!(m)

    num_sols = MOI.get(m, MOI.ResultCount())
    @test num_sols == count((i <= 2) ⊻ (j <= 2) for i=1:5, j=1:5)
    for i in 1:num_sols
        xval = convert(Int, JuMP.value(x; result=i))
        yval = convert(Int, JuMP.value(y; result=i))
        @test (xval <= 2) ⊻ (yval <= 2)
    end
end
end


text:
using DynamicHMCModels

ProjDir = @__DIR__

df = DataFrame(CSV.read(joinpath(ProjDir, "..", "..", "data", "chimpanzees.csv"), delim=';'))

Base.@kwdef struct Chimpanzees_02
    "Number of actors"
    N_actors::Int
    pulled_left::Vector{Int}
    prosoc_left::Vector{Int}
    condition::Vector{Int}
    actor::Vector{Int}
end

function make_transformation(model::Chimpanzees_02)
    as((a = as(Vector, model.N_actors), bp = asℝ, bpC = asℝ))
end

model = Chimpanzees_02(; N_actors = maximum(df.actor), pulled_left = df.pulled_left,
                    prosoc_left = df.prosoc_left, condition = df.condition,
                    actor = df.actor)

function (model::Chimpanzees_02)(θ)
    @unpack pulled_left, prosoc_left, condition, actor = model
    @unpack a, bp, bpC = θ
    ℓ_likelihood = mapreduce(+, actor, condition, prosoc_left,
       pulled_left) do actor, condition, prosoc_left, pulled_left
           p = logistic(a[actor] + (bp + bpC * condition) * prosoc_left)
           logpdf(Bernoulli(p), pulled_left)
       end
    P = Normal(0, 10)
    ℓ_prior = logpdf(P, bpC) + logpdf(P, bp) + sum(a -> logpdf(P, a), a)
    ℓ_prior + ℓ_likelihood
end

P = TransformedLogDensity(make_transformation(model), model)
∇P = ADgradient(:ForwardDiff, P)
results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000)
posterior = P.transformation.(results.chain)

DynamicHMC.Diagnostics.EBFMI(results.tree_statistics)

DynamicHMC.Diagnostics.summarize_tree_statistics(results.tree_statistics)

a3d = Array{Float64, 3}(undef, 1000, 9, 1);
for j in 1:1
  for i in 1:1000
    a3d[i, 1, j] = values(posterior[i].bp)
    a3d[i, 2, j] = values(posterior[i].bpC)
    a3d[i, 3:9, j] = values(posterior[i].a)
  end
end

# Create MCMCChains object

parameter_names = ["bp", "bpC"]
pooled_parameter_names = ["a[$i]" for i in 1:7]
sections =   Dict(
  :parameters => parameter_names,
  :pooled => pooled_parameter_names
)
cnames = vcat(parameter_names, pooled_parameter_names)
chns = create_mcmcchains(a3d, cnames, sections, start=1)

# Result rethinking

rethinking = "
Iterations = 1:1000
Thinning interval = 1
Chains = 1,2,3,4
Samples per chain = 1000

Empirical Posterior Estimates:
        Mean        SD       Naive SE       MCSE      ESS
a.1 -0.74503184 0.26613979 0.0042080396 0.0060183398 1000
a.2 10.77955494 5.32538998 0.0842018089 0.1269148045 1000
a.3 -1.04982353 0.28535997 0.0045119373 0.0049074219 1000
a.4 -1.04898135 0.28129307 0.0044476339 0.0056325117 1000
a.5 -0.74390933 0.26949936 0.0042611590 0.0052178124 1000
a.6  0.21599365 0.26307574 0.0041595927 0.0045153523 1000
a.7  1.81090866 0.39318577 0.0062168129 0.0071483527 1000
 bp  0.83979926 0.26284676 0.0041559722 0.0059795826 1000
bpC -0.12913322 0.29935741 0.0047332562 0.0049519863 1000
";

# Describe draws

describe(chns) |> display

# Describe pooled draws

describe(chns, sections=[:pooled])


text:
using Revise # reduce need for recompile
using Plots
using LinearAlgebra

push!(LOAD_PATH, "./src") # user defined modules
using CommonUtils, Basis1D

"Approximation parameters"
N   = 3 # The order of approximation
K1D = 16
CFL = 2
T   = 10 # endtime

"Mesh related variables"
VX = LinRange(-1,1,K1D+1)
EToV = transpose(reshape(sort([1:K1D; 2:K1D+1]),2,K1D))

"Construct matrices on reference elements"
r,w = gauss_lobatto_quad(0,0,N)
V = vandermonde_1D(N, r)
Dr = grad_vandermonde_1D(N, r)/V
M = inv(V*V')

"Nodes on faces, and face node coordinate"
wf = [1;1]
Vf = vandermonde_1D(N,[-1;1])/V
LIFT = M\(transpose(Vf)*diagm(wf)) # lift matrix

"Construct global coordinates"
V1 = vandermonde_1D(1,r)/vandermonde_1D(1,[-1;1])
x = V1*VX[transpose(EToV)]

"Connectivity maps"
xf = Vf*x
mapM = reshape(1:2*K1D,2,K1D)
mapP = copy(mapM)
mapP[1,2:end] .= mapM[2,1:end-1]
mapP[2,1:end-1] .= mapM[1,2:end]

"Make periodic"
mapP[1] = mapM[end]
mapP[end] = mapM[1]

"Geometric factors and surface normals"
J = repeat(transpose(diff(VX)/2),N+1,1)
nxJ = repeat([-1;1],1,K1D)
rxJ = 1

"initial conditions"
u = @. exp(-25*x^2)

"Time integration"
rk4a,rk4b,rk4c = rk45_coeffs()
CN = (N+1)*(N+2)/2  # estimated trace constant
dt = CFL * 2 / (CN*K1D)
Nsteps = convert(Int,ceil(T/dt))
dt = T/Nsteps

"pack arguments into tuples - will "
ops = (Dr,LIFT,Vf)
vgeo = (rxJ,J)
fgeo = (nxJ)

function rhs(u,ops,vgeo,fgeo,mapP)
    # unpack args
    Dr,LIFT,Vf = ops
    rxJ,J = vgeo
    nxJ = fgeo

    uM = Vf*u # can replace with nodal extraction
    du = uM[mapP]-uM

    ux = rxJ*(Dr*u)
    tau = 1 # upwind penalty parameter
    rhsu = ux + .5*LIFT*(@. du*nxJ - tau*abs(nxJ)*du)

    return -rhsu./J
end

"plotting nodes"
Vp = vandermonde_1D(N,LinRange(-1,1,10))/V
gr(size=(300,300),legend=false,markerstrokewidth=1,markersize=2)
plt = plot(Vp*x,Vp*u)

resu = zeros(size(x))
@gif for i = 1:Nsteps
    for INTRK = 1:5
        rhsu = rhs(u,ops,vgeo,fgeo,mapP)
        @. resu = rk4a[INTRK]*resu + dt*rhsu
        @. u   += rk4b[INTRK]*resu
    end

    if i%10==0 || i==Nsteps
        println("Number of time steps $i out of $Nsteps")
        # display(plot(Vp*x,Vp*u,ylims=(-.1,1.1)))
        # push!(plt, Vp*x,Vp*u,ylims=(-.1,1.1))
        plot(Vp*x,Vp*u,ylims=(-.1,1.1),title="Timestep $i out of $Nsteps",lw=2)
        scatter!(x,u)
        # sleep(.0)
    end
end every 5

# plot(Vp*x,Vp*u,ylims=(-.1,1.1))

max_stars_repo_name: jlchan/ESDG-CNS

max_issues_repo_name: jlchan/ESDG-CNS

max_forks_repo_name: jlchan/ESDG-CNS


text:
struct XYRobot
    robot::Robot
    coord::Coord
end

max_stars_repo_path: functions/struct XYRobot.jl

max_stars_repo_name: Igetsuga/AlgorithmicsHomework

max_issues_repo_path: functions/struct XYRobot.jl

max_issues_repo_name: Igetsuga/AlgorithmicsHomework

max_forks_repo_path: functions/struct XYRobot.jl

max_forks_repo_name: Igetsuga/AlgorithmicsHomework


text:
>> ethopian_multiplication(17, 34)
578
>> ethopian_multiplication_arr(44, 71)
3124

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Test

@testset "Testing solution to Exercise 1" begin

@testset "Testing that variables exist" begin
   @test isdefined(Main,:student_name)  # Main refers to the Main namespace, i.e. where
   @test isdefined(Main,:student_dept)  # variables go by default if they're not in a module.
   @test isdefined(Main,:student_year)  # Most packages have
end;

@testset "Test value of student_name is not obviously wrong" begin
   @test typeof(student_name) == String
   @test length(student_name) > 3
   @test occursin(' ',student_name)
   @test student_name != instructor
end;

@testset "Test value of student_dept is not obviously wrong" begin
   @test typeof(student_dept) == String
   @test length(student_dept) >= 3
end;

@testset "Test value of student_year is not obviously wrong" begin
   @test typeof(student_year) <: Number
   @test typeof(student_year) <: Integer
   @test 1 <= student_year <= 10
end;

end; # Exercise 1


text:
using KoLesky
using Plots
import Random.seed!
seed!(123)

uptos = [1, 100, 500]

x = rand(2, 5)

dofs = mat2points(x)

# @time P_r, ℓ_r = r_maximin(dofs; n_bins=100)
@time P, ℓ = maximin(dofs)

# pl = scatter(aspect_ratio=:equal, size=(500,500))
#
# for k = 2 : length(uptos)
#     scatter!(pl, vec(x[1, P[uptos[k-1] : uptos[k]]]), vec(x[2, P[uptos[k-1] : uptos[k]]]))
# end
#
# display(pl)
#

ρ = 6.0
λ = 1.5

@time 𝒩 = partition_into_supernodes(dofs, ℓ, P, ρ, λ)


text:
using Test
using LogicCircuits

include("../helper/little_circuits.jl")

@testset "Jlc write and save" begin

  sdd = readme_sdd()

  mktempdir() do tmp

    # write as a unstructured logic circuit
    jlc_path = "$tmp/example.jlc"
    write(jlc_path, sdd)

    # read as a unstructured logic circuit
    jlc2 = read(jlc_path, LogicCircuit, JlcFormat())

    @test jlc2 isa PlainLogicCircuit
    @test num_nodes(sdd) == num_nodes(jlc2)
    @test prob_equiv(sdd, jlc2, 10)

    # write/read as an instructured logic circuit, compressed
    write("$jlc_path.gz", sdd)
    jlc2 = read("$jlc_path.gz", LogicCircuit)

    @test jlc2 isa PlainLogicCircuit
    @test num_nodes(sdd) == num_nodes(jlc2)
    @test prob_equiv(sdd, jlc2, 10)

    # write with vtree
    vtree_path = "$tmp/example.vtree"
    paths = (jlc_path, vtree_path)
    write(paths, sdd)

    # read as a structured logic circuit
    formats = (JlcFormat(), VtreeFormat())
    jlc3 = read(paths, StructLogicCircuit, formats)

    @test jlc3 isa PlainStructLogicCircuit
    @test num_nodes(sdd) == num_nodes(jlc3)
    @test prob_equiv(sdd, jlc3, 10)
    @test Vtree(mgr(sdd)) == vtree(jlc3)

    # write/read with vtree, compressed
    pathsgz = ("$jlc_path.gz", vtree_path)
    write(pathsgz, sdd)
    jlc3 = read(pathsgz, StructLogicCircuit)

    @test jlc3 isa PlainStructLogicCircuit
    @test num_nodes(sdd) == num_nodes(jlc3)
    @test prob_equiv(sdd, jlc3, 10)
    @test Vtree(mgr(sdd)) == vtree(jlc3)

  end

end

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges
const MOIBC = MathOptInterface.Bridges.Constraint

include("../utilities.jl")

mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()))
config = MOIT.TestConfig()

@testset "SemiToBinary" begin
    bridged_mock = MOIBC.SemiToBinary{Float64}(mock)

    bridge_type = MOIBC.SemiToBinaryBridge{Float64,MOI.Semiinteger{Float64}}
    @test MOI.supports_constraint(
        bridge_type,
        MOI.SingleVariable,
        MOI.Semiinteger{Float64},
    )
    @test MOIBC.concrete_bridge_type(
        bridge_type,
        MOI.SingleVariable,
        MOI.Semiinteger{Float64},
    ) == bridge_type

    @test MOI.supports(bridged_mock, MOI.ConstraintPrimalStart(), bridge_type)

    MOIT.basic_constraint_tests(
        bridged_mock,
        config,
        include = [
            (F, S) for F in [MOI.SingleVariable] for
            S in [MOI.Semiinteger{Float64}, MOI.Semicontinuous{Float64}]
        ],
    )

    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> begin
            MOI.set(mock, MOI.ObjectiveBound(), 0.0)
            MOIU.mock_optimize!(mock, [0.0, 0.0, 0.0])
        end,
        (mock::MOIU.MockOptimizer) -> begin
            MOI.set(mock, MOI.ObjectiveBound(), 2.0)
            MOIU.mock_optimize!(mock, [2.0, 1.0, 1.0])
        end,
        (mock::MOIU.MockOptimizer) -> begin
            MOI.set(mock, MOI.ObjectiveBound(), 2.0)
            MOIU.mock_optimize!(mock, [2.0, 2.0, 1.0])
        end,
        (mock::MOIU.MockOptimizer) -> begin
            MOI.set(mock, MOI.ObjectiveBound(), 2.5)
            MOIU.mock_optimize!(mock, [2.5, 2.5, 1.0])
        end,
        (mock::MOIU.MockOptimizer) -> begin
            MOI.set(mock, MOI.ObjectiveBound(), 3.0)
            MOIU.mock_optimize!(mock, [3.0, 3.0, 1.0])
        end,
        (mock::MOIU.MockOptimizer) ->
            MOI.set(mock, MOI.TerminationStatus(), MOI.INFEASIBLE),
    )
    MOIT.semiconttest(bridged_mock, config)

    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.SingleVariable,
                MOI.Semicontinuous{Float64},
            }(),
        ),
    )

    test_delete_bridge(
        bridged_mock,
        ci,
        2,
        (
            (MOI.SingleVariable, MOI.EqualTo{Float64}, 1),
            (MOI.SingleVariable, MOI.ZeroOne, 0),
            (MOI.SingleVariable, MOI.Integer, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 1),
        ),
    )

    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> begin
            MOI.set(mock, MOI.ObjectiveBound(), 0.0)
            MOIU.mock_optimize!(mock, [0.0, 0.0, 0.0])
        end,
        (mock::MOIU.MockOptimizer) -> begin
            MOI.set(mock, MOI.ObjectiveBound(), 2.0)
            MOIU.mock_optimize!(mock, [2.0, 1.0, 1.0])
        end,
        (mock::MOIU.MockOptimizer) -> begin
            MOI.set(mock, MOI.ObjectiveBound(), 2.0)
            MOIU.mock_optimize!(mock, [2.0, 2.0, 1.0])
        end,
        (mock::MOIU.MockOptimizer) -> begin
            MOI.set(mock, MOI.ObjectiveBound(), 3.0)
            MOIU.mock_optimize!(mock, [3.0, 2.5, 1.0])
        end,
        (mock::MOIU.MockOptimizer) -> begin
            MOI.set(mock, MOI.ObjectiveBound(), 3.0)
            MOIU.mock_optimize!(mock, [3.0, 3.0, 1.0])
        end,
        (mock::MOIU.MockOptimizer) ->
            MOI.set(mock, MOI.TerminationStatus(), MOI.INFEASIBLE),
    )
    MOIT.semiinttest(bridged_mock, config)

    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.SingleVariable,
                MOI.Semiinteger{Float64},
            }(),
        ),
    )

    @test MOI.get(bridged_mock, MOI.ConstraintPrimal(), ci) == 3
    new_set = MOI.Semiinteger{Float64}(19.0, 20.0)
    MOI.set(bridged_mock, MOI.ConstraintSet(), ci, new_set)
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), ci) == new_set

    @testset "$attr" for attr in [MOI.ConstraintPrimalStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        value = 2.0
        MOI.set(bridged_mock, attr, ci, value)
        @test MOI.get(bridged_mock, attr, ci) ≈ value
    end

    test_delete_bridge(
        bridged_mock,
        ci,
        2,
        (
            (MOI.SingleVariable, MOI.EqualTo{Float64}, 1),
            (MOI.SingleVariable, MOI.ZeroOne, 0),
            (MOI.SingleVariable, MOI.Integer, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 1),
        ),
    )

    s = """
    variables: x, y
    cy: y == 4.0
    cx: x in Semiinteger(2.0, 3.0)
    minobjective: x
    """
    model = MOIU.Model{Float64}()
    MOIU.loadfromstring!(model, s)
    sb = """
    variables: x, y, z
    cy: y == 4.0
    bin: z in ZeroOne()
    int: x in Integer()
    clo: x + -2.0z >= 0.0
    cup: x + -3.0z <= 0.0
    minobjective: x
    """
    modelb = MOIU.Model{Float64}()
    MOIU.loadfromstring!(modelb, sb)

    MOI.empty!(bridged_mock)
    @test MOI.is_empty(bridged_mock)
    MOIU.loadfromstring!(bridged_mock, s)
    MOIU.test_models_equal(bridged_mock, model, ["x", "y"], ["cy", "cx"])

    # setting names on mock
    ci = first(
        MOI.get(
            mock,
            MOI.ListOfConstraintIndices{MOI.SingleVariable,MOI.ZeroOne}(),
        ),
    )
    MOI.set(mock, MOI.ConstraintName(), ci, "bin")
    z = MOI.VariableIndex(ci.value)
    MOI.set(mock, MOI.VariableName(), z, "z")
    ci = first(
        MOI.get(
            mock,
            MOI.ListOfConstraintIndices{MOI.SingleVariable,MOI.Integer}(),
        ),
    )
    MOI.set(mock, MOI.ConstraintName(), ci, "int")
    ci = first(
        MOI.get(
            mock,
            MOI.ListOfConstraintIndices{
                MOI.ScalarAffineFunction{Float64},
                MOI.GreaterThan{Float64},
            }(),
        ),
    )
    MOI.set(mock, MOI.ConstraintName(), ci, "clo")
    ci = first(
        MOI.get(
            mock,
            MOI.ListOfConstraintIndices{
                MOI.ScalarAffineFunction{Float64},
                MOI.LessThan{Float64},
            }(),
        ),
    )
    MOI.set(mock, MOI.ConstraintName(), ci, "cup")

    MOIU.test_models_equal(
        mock,
        modelb,
        ["x", "y", "z"],
        ["cy", "bin", "int", "clo", "cup"],
    )
end


text:
using StatsBase, JLD,  Hwloc,  DynNets, StatsFuns
initWorkers()
@everywhere using DynNets, StatsFuns
N_est = 4
NGWpoints = 3
Tvals = [50 250 1250]# [ 100 200  300  ]
N_Tind = length(Tvals)[1]
scoreRescType = ""#"FISHER-EWMA"#
useStartVal = false# true# false
Nmin = 10
Nspacing =  5
Nmax = 20
Nvals = Vector(Nmin:Nspacing:Nmax)
N_Nind = length(Nvals)

est_times = SharedArray{Float64}(zeros(N_Tind,N_Nind,NGWpoints,N_est))
est_conv_flag = SharedArray{Bool,4}(N_Tind,N_Nind,NGWpoints,N_est)
est_flag = SharedArray{Bool,4}(N_Tind,N_Nind,NGWpoints,N_est)
GWvalStore  = Array{Array{Int,1},1}(N_Nind)

@progress "N" for indN = N_Nind:-1:1
#indN = 1
N = Nvals[indN] ;
GWvals = unique(Vector{Int}(floor.(linspace(1,N,NGWpoints))))
GWvalStore[indN] = GWvals
N_GWind = length(GWvals)

@progress "GW"  for indGW =N_GWind:-1:1
#indGW = 2
GW = GWvals[indGW]
GBAdgp =1; GWdgp = GW ;
GBAest =1; GWest = GW;
#Distribuisci i nodi tra i gruppi in parti uguali

groupsindTsDgp = [[DynNets.distributeAinVecN(Array{Int}(1:GWdgp),N)];
                    [DynNets.distributeAinVecN(Array{Int}(1:GBAdgp),GWdgp)]]
groupsindTsEst = [[DynNets.distributeAinVecN(Array{Int}(1:GWest),N)];
                   [DynNets.distributeAinVecN(Array{Int}(1:GBAest),GWest)]]

#Initialize the variable to store
realUm = DynNets.linSpacedPar(DynNets.fooSdErgm1, N;NgroupsW = GW )[1]
sim_B = Array(linspace(0.4,0.9,GBAdgp+1)[2:end]) #
sim_A = Array(linspace(0.02,0.2,GBAdgp+1)[2:end]) #
sim_W = realUm.*(1-sim_B[groupsindTsDgp[2]])
[sim_W,sim_B,sim_A]

@progress "T" for indT = N_Tind:-1:1
#  indT=1
T = Tvals[indT]
@sync @parallel for n = 1:N_est
#n=1

        Model2Sample = DynNets.SdErgm1(ones(T,N),[sim_W,sim_B,sim_A],
                                                    groupsindTsDgp,scoreRescType)
        Model2Sample,Y_T,Fitness_T  =  DynNets.sampl(Model2Sample,T)
        sampleDegsT = squeeze(sum(Y_T,2),2)

        Model2Est =  DynNets.SdErgm1(sampleDegsT,[zeros(GWest),zeros(GBAest),zeros(GBAest)],
                                                    groupsindTsEst,scoreRescType)
        #estimate
        useStartVal ?   start_val = [sim_W, sim_B, sim_A] : start_val = zeros(3,3)
        est_times[indT,indN,indGW,n] = @elapsed  hat_pars , conv_flag =
                            DynNets.estimate(Model2Est,start_values = start_val )
        est_conv_flag[indT,indN,indGW,n] = conv_flag
        est_flag[indT,indN,indGW,n] = true

    end
end
end
    # save
    save_fold = "./data/estimatesTest/"
    file_name = "EstTest_varN_Nest_$(N_est)_$(Nmin)_$(Nspacing)_$(Nmax)_varGW_$(NGWpoints)_T_$(Tvals)_corectSpec"* scoreRescType *".jld"
    save_path = save_fold*file_name#
    @save(save_path,est_conv_flag,est_times,est_flag,Nvals,GWvalStore ,useStartVal)
    println([N,sum(est_times[:,indN,:,:]),sum(est_conv_flag[:,indN,:,:])/sum(est_flag[:,indN,:,:])])
end

##

 ##


text:
############################################################
# judiRHS ## #################################################
############################################################

# Authors: Philipp Witte (pwitte@eos.ubc.ca), Henryk Modzelewski (hmodzelewski@eos.ubc.ca)
# Date: January 2017

export judiRHS, judiRHSexception

############################################################

struct judiRHS{vDT<:Number} <: joAbstractLinearOperator{vDT,vDT}
    name::String
    m::Integer
    n::Integer
    info::Info
    geometry::Geometry
    data
end

mutable struct judiRHSexception <: Exception
    msg :: String
end

############################################################

## outer constructors

"""
    judiRHS
        name::String
        m::Integer
        n::Integer
        info::Info
        geometry::Geometry
        data

Abstract sparse vector for right-hand-sides of the modeling operators. The `judiRHS` vector has the\\
dimensions of the full time history of the wavefields, but contains only the data defined at the \\
source or receiver positions (i.e. wavelets or shot records).

Constructor
==========

    judiRHS(info, geometry, data)

Examples
========

Assuming `Pr` and `Ps` are projection operators of type `judiProjection` and `dobs` and `q` are\\
seismic vectors of type `judiVector`, then a `judiRHS` vector can be created as follows:

    rhs = Pr'*dobs    # right-hand-side with injected observed data
    rhs = Ps'*q    # right-hand-side with injected wavelet

"""
function judiRHS(info,geometry,data;vDT::DataType=Float32)
    vDT == Float32 || throw(judiRHSexception("Domain type not supported"))
    # length of vector
    m = info.n * sum(info.nt)
    n = 1
    return judiRHS{Float32}("judiRHS",m,n,info,geometry,data)
end

####################################################################
## overloaded Base functions

# conj(jo)
conj(A::judiRHS{vDT}) where vDT =
    judiRHS{vDT}("conj("*A.name*")",A.m,A.n,A.info,A.geometry,A.data)

# transpose(jo)
transpose(A::judiRHS{vDT}) where vDT =
    judiRHS{vDT}(""*A.name*".'",A.n,A.m,A.info,A.geometry,A.data)

# adjoint(jo)
adjoint(A::judiRHS{vDT}) where vDT =
    judiRHS{vDT}(""*A.name*"'",A.n,A.m,A.info,A.geometry,A.data)

####################################################################

# +(judiRHS,judiRHS)
function +(A::judiRHS{avDT}, B::judiRHS{bvDT}) where {avDT, bvDT}

    # Error checking
    size(A) == size(B) || throw(judiRHSexception("shape mismatch"))
    compareInfo(A.info, B.info) == true || throw(judiRHSexception("info mismatch"))
    isequal(A.geometry.nt,B.geometry.nt) == true || throw(judiRHSexception("sample number mismatch"))
    isequal(A.geometry.dt,B.geometry.dt) == true || throw(judiRHSexception("sample interval mismatch"))
    isequal(A.geometry.t,B.geometry.t) == true || throw(judiRHSexception("recording time mismatch"))

    # Size
    m = A.info.n * sum(A.info.nt)
    n = 1

    # merge geometries and data
    xloc = Array{Any}(undef, A.info.nsrc)
    yloc = Array{Any}(undef, A.info.nsrc)
    zloc = Array{Any}(undef, A.info.nsrc)
    dt = Array{Any}(undef, A.info.nsrc)
    nt = Array{Any}(undef, A.info.nsrc)
    t = Array{Any}(undef, A.info.nsrc)
    data = Array{Any}(undef, A.info.nsrc)

    for j=1:A.info.nsrc
        xloc[j] = [vec(A.geometry.xloc[j])' vec(B.geometry.xloc[j])']
        yloc[j] = [vec(A.geometry.yloc[j])' vec(B.geometry.yloc[j])']
        zloc[j] = [vec(A.geometry.zloc[j])' vec(B.geometry.zloc[j])']
        dt[j] = A.geometry.dt[j]
        nt[j] = A.geometry.nt[j]
        t[j] = A.geometry.t[j]
        data[j] = [A.data[j] B.data[j]]
    end
    geometry = Geometry(xloc,yloc,zloc,dt,nt,t)
    nvDT = promote_type(avDT,bvDT)

    return judiRHS{nvDT}("judiRHS",m,n,A.info,geometry,data)
end

# -(judiRHS,judiRHS)
function -(A::judiRHS{avDT}, B::judiRHS{bvDT}) where {avDT, bvDT}

    # Error checking
    size(A) == size(B) || throw(judiRHSexception("shape mismatch"))
    compareInfo(A.info, B.info) == true || throw(judiRHSexception("info mismatch"))
    isequal(A.geometry.nt,B.geometry.nt) == true || throw(judiRHSexception("sample number mismatch"))
    isequal(A.geometry.dt,B.geometry.dt) == true || throw(judiRHSexception("sample interval mismatch"))
    isequal(A.geometry.t,B.geometry.t) == true || throw(judiRHSexception("recording time mismatch"))

    # Size
    m = A.info.n * sum(A.info.nt)
    n = 1

    # merge geometries and data
    xloc = Array{Any}(undef, A.info.nsrc)
    yloc = Array{Any}(undef, A.info.nsrc)
    zloc = Array{Any}(undef, A.info.nsrc)
    dt = Array{Any}(undef, A.info.nsrc)
    nt = Array{Any}(undef, A.info.nsrc)
    t = Array{Any}(undef, A.info.nsrc)
    data = Array{Any}(undef, A.info.nsrc)

    for j=1:A.info.nsrc
        xloc[j] = [vec(A.geometry.xloc[j])' vec(B.geometry.xloc[j])']
        yloc[j] = [vec(A.geometry.yloc[j])' vec(B.geometry.yloc[j])']
        zloc[j] = [vec(A.geometry.zloc[j])' vec(B.geometry.zloc[j])']
        dt[j] = A.geometry.dt[j]
        nt[j] = A.geometry.nt[j]
        t[j] = A.geometry.t[j]
        data[j] = [A.data[j] -B.data[j]]
    end
    geometry = Geometry(xloc,yloc,zloc,dt,nt,t)
    nvDT = promote_type(avDT,bvDT)

    return judiRHS{nvDT}("judiRHS",m,n,A.info,geometry,data)
end

function subsample(a::judiRHS{avDT},srcnum) where avDT
    info = Info(a.info.n, length(srcnum), a.info.nt[srcnum])
    geometry = subsample(a.geometry,srcnum)     # Geometry of subsampled data container
    return judiRHS(info,geometry,a.data[srcnum];vDT=avDT)
end

getindex(x::judiRHS,a) = subsample(x,a)


text:
include("common.jl")

@testset "Optimize single objective problems in 5, 10, and 30 dimensions with DE" begin
    simple_problems = ["Sphere", "Schwefel2.22", "Schwefel2.22"]
    @testset "Simple problem $pr" for pr in simple_problems
        p = BlackBoxOptim.example_problems[pr]

        @test fitness_for_opt(p,  5, 20, 5e3, :de_rand_1_bin) < 0.01
        @test fitness_for_opt(p,  5, 20, 5e3, :de_rand_1_bin_radiuslimited) < 0.01

        @test fitness_for_opt(p, 10, 20, 1e4, :de_rand_1_bin) < 0.01
        @test fitness_for_opt(p, 10, 20, 1e4, :de_rand_1_bin_radiuslimited) < 0.01

        @test fitness_for_opt(p, 30, 25, 3e4, :de_rand_1_bin) < 0.01
        @test fitness_for_opt(p, 30, 25, 3e4, :de_rand_1_bin_radiuslimited) < 0.01
    end

    @testset "Schwefel1.2" begin
        problem = "Schwefel1.2"
        p = BlackBoxOptim.example_problems[problem]
        @test fitness_for_opt(p, 5,  20, 5e3, :de_rand_1_bin_radiuslimited) < 0.01
        @test fitness_for_opt(p, 10, 50, 5e4, :de_rand_1_bin_radiuslimited) < 0.01

        #DE/rand/1/bin seems to have troubles...
        @test fitness_for_opt(p, 30, 50, 2e5, :de_rand_1_bin_radiuslimited) < 10.0
        @test fitness_for_opt(p, 30, 50, 2e5, :adaptive_de_rand_1_bin) < 10.0
        @test fitness_for_opt(p, 30, 50, 2e5, :adaptive_de_rand_1_bin_radiuslimited) < 10.0
    end

    @testset "Rosenbrock" begin
        problem = "Rosenbrock"
        p = BlackBoxOptim.example_problems[problem]
        @test fitness_for_opt(p, 5,  20, 1e4, :de_rand_1_bin_radiuslimited) < 100.0
        @test fitness_for_opt(p, 10, 20, 5e4, :de_rand_1_bin_radiuslimited) < 100.0
        @test fitness_for_opt(p, 30, 40, 2e5, :de_rand_1_bin_radiuslimited) < 100.0

        @test fitness_for_opt(p, 30, 40, 2e5, :adaptive_de_rand_1_bin) < 100.0
        @test fitness_for_opt(p, 30, 40, 2e5, :adaptive_de_rand_1_bin_radiuslimited) < 100.0

        @test fitness_for_opt(p, 50, 40, 3e5, :adaptive_de_rand_1_bin_radiuslimited) < 100.0
    end
end


text:
using Test
using DataValues

@testset "DataValueArray: Show" begin

io = IOBuffer()
disp = TextDisplay(IOBuffer())

for typ in (Float64, Int, UInt, Char)
    for l in (10, 100, 1000)
        A1 = rand(typ, l)
        M1 = rand(Bool, l)
        X1 = DataValueArray(A1)
        Y1 = DataValueArray(A1, M1)

        show(io, X1)
        show(io, Y1)
        display(disp, X1)
        display(disp, Y1)
        typeof(X1)
        typeof(Y1)
    end

    for m in (10, 100), n in (10, 100)
        A2 = rand(typ, m, n)
        M2 = rand(Bool, m, n)
        X2 = DataValueArray(A2)
        Y2 = DataValueArray(A2, M2)

        show(io, X2)
        show(io, Y2)
        display(disp, X2)
        display(disp, Y2)
        typeof(X2)
        typeof(Y2)
    end

    nd = rand(3:5)
    sz = [ rand(3:5) for i in 1:nd ]
    A3 = rand(typ, sz...)
    M3 = rand(Bool, sz...)
    X3 = DataValueArray(A3)
    Y3 = DataValueArray(A3, M3)

    show(io, X3)
    show(io, Y3)
    display(disp, X3)
    display(disp, Y3)
    typeof(X3)
    typeof(Y3)

    X = DataValueArray{typ}()
    show(io, X)
    display(disp, X)
    typeof(X)
end

end


text:
function run_BMA(h, h_lag, J, J_cv, level, head_line, Y, Y_lag, X_lag, X_lag2, oos_i)
    name = "BMA"
    df_file = "../../data/results/rolling_cv/level$level-h$h-J$J/$name"
    ## Forecast using Bayesian Model Averaging
    θ_vec = 0.05:0.01:0.99
    pre_conv = 1_000
    post_conv = 2_000
    θ, ŷ, RMSE_cv = bma_oos_forc(J, h, h_lag, J_cv, Y, Y_lag, X_lag, θ_vec, pre_conv, post_conv, oos_i)
    writedlm("$df_file/forecast/$oos_i.csv", ŷ, ',')
    writedlm("$df_file/tuning/$oos_i.csv", θ, ',')
    writedlm("$df_file/RMSE_cv/$oos_i.csv", RMSE_cv, ',')
end

function run_BMA2(h, h_lag, J, J_cv, level, head_line, Y, Y_lag, X_lag, X_lag2, oos_i)
    name = "BMA2"
    df_file = "../../data/results/rolling_cv/level$level-h$h-J$J/$name"
    ## Forecast using Bayesian Model Averaging with interaction terms
    θ_vec = 0.05:0.01:0.99
    pre_conv = 1_000
    post_conv = 2_000
    X_lag = hcat(X_lag, X_lag2)
    θ, ŷ, RMSE_cv = bma_oos_forc(J, h, h_lag, J_cv, Y, Y_lag, X_lag, θ_vec, pre_conv, post_conv, oos_i)
    writedlm("$df_file/forecast/$oos_i.csv", ŷ, ',')
    writedlm("$df_file/tuning/$oos_i.csv", θ, ',')
    writedlm("$df_file/RMSE_cv/$oos_i.csv", RMSE_cv, ',')
end

function run_LAS(h, h_lag, J, J_cv, level, head_line, Y, Y_lag, X_lag, X_lag2, oos_i)
    name = "LAS"
    df_file = "../../data/results/rolling_cv/level$level-h$h-J$J/$name"
    #Forecast using a Ridge model with interaction terms
    α = 1.0
    λ_vec = 0.001:0.001:0.75
    β̂, λ, ŷ, RMSE_cv = glmneoos_i_forc(J, h, h_lag, J_cv, Y, Y_lag, X_lag, λ_vec, α, oos_i)
    writedlm("$df_file/forecast/$oos_i.csv", ŷ, ',')
    writedlm("$df_file/tuning/$oos_i.csv", λ, ',')
    writedlm("$df_file/RMSE_cv/$oos_i.csv", RMSE_cv, ',')
    writedlm("$df_file/beta/$oos_i.csv", β̂, ',')
end

function run_LAS2(h, h_lag, J, J_cv, level, head_line, Y, Y_lag, X_lag, X_lag2, oos_i)
    name = "LAS2"
    df_file = "../../data/results/rolling_cv/level$level-h$h-J$J/$name"
    #Forecast using a LASSO model with interaction terms
    α = 1.0
    λ_vec = 0.1:0.01:5.
    X_lag = hcat(X_lag, X_lag2)
    β̂, λ, ŷ, RMSE_cv = glmneoos_i_forc(J, h, h_lag, J_cv, Y, Y_lag, X_lag, λ_vec, α, oos_i)
    writedlm("$df_file/forecast/$oos_i.csv", ŷ, ',')
    writedlm("$df_file/tuning/$oos_i.csv", λ, ',')
    writedlm("$df_file/RMSE_cv/$oos_i.csv", RMSE_cv, ',')
    writedlm("$df_file/beta/$oos_i.csv", β̂, ',')
end

function run_RDG(h, h_lag, J, J_cv, level, head_line, Y, Y_lag, X_lag, X_lag2, oos_i)
    name = "RDG"
    df_file = "../../data/results/rolling_cv/level$level-h$h-J$J/$name"
    #Forecast using a Ridge model with interaction terms
    α = 0.0
    λ_vec = 1.:0.1:25.
    β̂, λ, ŷ, RMSE_cv = glmneoos_i_forc(J, h, h_lag, J_cv, Y, Y_lag, X_lag, λ_vec, α, oos_i)
    writedlm("$df_file/forecast/$oos_i.csv", ŷ, ',')
    writedlm("$df_file/tuning/$oos_i.csv", λ, ',')
    writedlm("$df_file/RMSE_cv/$oos_i.csv", RMSE_cv, ',')
end

function run_RDG2(h, h_lag, J, J_cv, level, head_line, Y, Y_lag, X_lag, X_lag2, oos_i)
    name = "RDG2"
    df_file = "../../data/results/rolling_cv/level$level-h$h-J$J/$name"
    #Forecast using a Ridge model with interaction terms
    α = 0.0
    λ_vec = 1.:0.1:50.
    X_lag = hcat(X_lag, X_lag2)
    β̂, λ, ŷ, RMSE_cv = glmneoos_i_forc(J, h, h_lag, J_cv, Y, Y_lag, X_lag, λ_vec, α, oos_i)
    writedlm("$df_file/forecast/$oos_i.csv", ŷ, ',')
    writedlm("$df_file/tuning/$oos_i.csv", λ, ',')
    writedlm("$df_file/RMSE_cv/$oos_i.csv", RMSE_cv, ',')
end

function run_RDF(h, h_lag, J, J_cv, level, head_line, Y, Y_lag, X_lag, X_lag2, oos_i)
    name = "RDF"
    df_file = "../../data/results/rolling_cv/level$level-h$h-J$J/$name"
    #Forecast using a Random Forrest Regression model
    m_vec = 1:1:25
    n_est = 50
    m, ŷ, RMSE_cv = rfr_oos_forc(J , h, h_lag, J_cv, Y, Y_lag, X_lag, m_vec, n_est, oos_i)
    writedlm("$df_file/forecast/$oos_i.csv", ŷ, ',')
    writedlm("$df_file/tuning/$oos_i.csv", m, ',')
    writedlm("$df_file/RMSE_cv/$oos_i.csv", RMSE_cv, ',')
end

function run_RDF2(h, h_lag, J, J_cv, level, head_line, Y, Y_lag, X_lag, X_lag2, oos_i)
    name = "RDF2"
    df_file = "../../data/results/rolling_cv/level$level-h$h-J$J/$name"
    #Forecast using a Random Forrest Regression model
    m_vec = 1:1:25
    n_est = 50
    X_lag = hcat(X_lag, X_lag2)
    m, ŷ, RMSE_cv = rfr_oos_forc(J , h, h_lag, J_cv, Y, Y_lag, X_lag, m_vec, n_est, oos_i)
    writedlm("$df_file/forecast/$oos_i.csv", ŷ, ',')
    writedlm("$df_file/tuning/$oos_i.csv", m, ',')
    writedlm("$df_file/RMSE_cv/$oos_i.csv", RMSE_cv, ',')
end

function dim_1to2(i)
    h_vec = [3, 6, 12]
    oos_i = mod(i, 408)
    h_indx = i ÷ 408
    if oos_i == 0
        oos_i = 408
    else
        h_indx += 1
    end
    h = h_vec[h_indx]
    return oos_i, h
end

## Set the number of cores to use
include("forc-helpers/prepare.jl")
level = 4
h_lag = 12
J = 301
## Different values for different cases
s = ArgParseSettings()
@add_arg_table s begin
    "i"
        arg_type = Int
        required = true
        help = "oos_i"
end
ps = parse_args(s)
i = ps["i"]
oos_i, h = dim_1to2(i)
head_line, Y, Y_lag, X_lag, X_lag2, J_cv = prepare(level, h, J, oos_i)
oos_vec = J - (h_lag + h):length(Y)
oos_num = length(oos_vec)
run_RDG(h, h_lag, J, J_cv, level, head_line, Y, Y_lag, X_lag, X_lag2, oos_i)

#=
h = 6
J = 301
name = "short_LAS"
tunes = readdlm("../../data/results/rolling_cv/level4-h$h-J$J/$name/tuning/combined.csv", ',')
ŷ_rolling = readdlm("../../data/results/rolling_cv/level4-h$h-J$J/$name/forecast/combined.csv", ',')
name = "LAS"
ŷ_baseline = convert(Vector{Float64}, readdlm("../../data/results/baseline/level4-h$h-J$J/$name.csv", ',')[2:end])
y_real = convert(Vector{Float64}, readdlm("../../data/results/baseline/level4-h$h-J$J/TRD.csv", ',')[2:end, 1])
sqrt(mean((ŷ_rolling - y_real).^2))
sqrt(mean((ŷ_baseline - y_real).^2))
plot(ŷ_rolling, label = "rolling")
plot!(y_real, lw = 2, label = "")
plot!(ŷ_baseline, label = "baseline")
plot(tunes)
=#

max_stars_repo_name: jungangl/forecast-machine-learning

max_issues_repo_name: jungangl/forecast-machine-learning

max_forks_repo_name: jungangl/forecast-machine-learning


text:
## Container widgets defined here

## methods
connect(object::Container, signal::Union(Symbol,String), obj, slot::Function) = connect(object.model, string(signal), obj, slot)
connect(object::Container, signal::Union(Symbol,String), slot::Function) = connect(object.model, string(signal), slot)
disconnect(o::Container, id::String) = disconnect(o.model, id)

## think of containers as Widget[] containers

@doc "Return all children of a container" ->
children(object::Container) = object.children

@doc "access children of a container by index" ->
getindex(object::Container, i::Union(Int, Range, Vector)) = children(object)[i]

## `push!(parent, child)` is the typical way a child is added to a
## container. This is a convenience as the child holds a reference to
## the parent.
push!(child::Widget) = push!(child.parent, child)

## Bin containers only allow one child
abstract BinContainer <: Container

## Children are added to a bin container by `push!`
function push!(parent::BinContainer, child::Widget)
    if parent != child.parent
        error("Window is not parent's child?")
    end

    if length(children(parent)) >= 1
        error("Can only add one child to a $(typeof(parent)) object")
    end
    ##
    push!(parent.children, child)

    set_child(parent.toolkit, parent, child)
end

@doc "Remove (alias for `pop!`) child of a bin container" ->
clear(parent::BinContainer) = pop!(parent)
length(parent::BinContainer) = length(children(parent))

##################################################
##
## Top-level window
##
type Window <: BinContainer
    o                           # frame
    block                       # window
    parent
    model
    toolkit
    children
    attrs::Dict
end

@doc """
Constructor for a toplevel window

A toplevel window is a `BinContainer` and can hold only one child.

## arguments

* `visible::Bool=true` if `false` suppress the initial drawing of the window

## Examples
```julia
w = window(title="title")
b = button(w, "delete")
connect(b, :clicked, () -> destroy(w))
raise(w)
```
""" ->
function window(;toolkit::MIME=default_toolkit, visible::Bool=true, kwargs...)
    widget, block = window(toolkit, visible=visible)
    obj = Window(widget, block, nothing, EventModel(), toolkit, {}, Dict())
    obj[:icontheme] = :default
    for (k, v) in kwargs
        obj[k] = v
    end
    obj
end

@doc "Raise window to top of viewable stack" ->
raise(o::Window) = raise(o.toolkit, o)
raise(o::Widget) = raise(o.parent)

@doc "lower window to bottom of visible stack" ->
lower(o::Widget) = lower(o.parent)
lower(o::Window) = lower(o.toolkit, o)

@doc "Destroy the window" ->
function destroy(o::Window)
    ## XXX call method that can intercept?
    destroy_window(o.toolkit, o)
end
destroy(o::Widget) = destroy(o.parent)

## properties
getTitle(o::Window) = getTitle(o.toolkit, o)
setTitle(o::Window, value::String) = setTitle(o.toolkit, o, value)
getModal(o::Window) = getModal(o.toolkit, o)
setModal(o::Window, value::Bool) = setModal(o.toolkit, o, value)
getPosition(o::Window) = getPosition(o.toolkit, o)
setPosition(o::Window, value::Vector{Int}) = setPosition(o.toolkit, o, value)
list_props(::@PROP("Window")) = {:title => "window title",
                                 :modal => "make window modal",
                                 :position => "move to position [x,y]"
                                 }

## XXX addToolBar, addMenuBar, addStatusBar...

## what name for this?
## not that useful as is, need to use Module.DataType, as these aren't exported
function lookup(o::Container, by::DataType)
    ## return all children of the given data type recursing over children
    function add_children(obj)
        for child in children(obj)
            isa(child, by) && push!(out, child)
            isa(child, Container) && add_children(child)
        end
    end
    out = {}
    add_children(o)
    out
end

## labelframe
type LabelFrame <: BinContainer
    o
    block
    parent
    model
    toolkit
    children
    attrs::Dict
end

@doc """
Label frame is a bin container to surround another container in a frame with a title.

## Arguments:

* `(parent, label::String; kwargs...)`
* `alignment::MaybeSymbol`, one of (nothing, :left, :center, :right)

## Examples
```julia
w = window(title="labelframe")
lf = labelframe(w, "A label")
push!(w, lf)
b = button(lf, "button")
push!(lf, b)
```

""" ->
function labelframe(parent::Container, label::String; alignment::Union(Nothing, Symbol)=nothing, kwargs...)
    (widget, block) = labelframe(parent.toolkit, parent, label, alignment=alignment)
    obj = LabelFrame(widget, block, parent, EventModel(), parent.toolkit, {}, Dict())
    for (k, v) in kwargs
        obj[k] = v
    end
    obj
end

##################################################
##
## Box containers

type BoxContainer <: Container
    o
    block
    parent
    toolkit
    children
    attrs::Dict
end

@doc """
A boxcontainer holds child objects either horizontally (`hbox`) or vertically (`vbox`)

Like Gtk, and not Qt or TclTk the container and the layout are coupled. Box containers are basically
frames with a box-like layout. By contrast, grid containers are frames with grid-like layouts.

Conceptually a box container is treated like a queue. Children are
pushed, appended, prepended, or inserted into it. They can be
popped, shifted, or spliced by index. Child components may be
referenced by integer indices.

Children are placed in box with size policies set by the
`sizepolicy` property. These are specified by direction `x` or `y`. The
value can be `nothing`, `:fixed`, or `:expand`. By default, widgets
will expand in the direction orthogonal to the packing direction
(horizontally for vertical boxes).

When a child does not expand to fill the allocated space, it can be
aligned within the cell. The `alignment` property sets attributes
that control this, again with `x` and `y` values in `(:left,
:right, :center, :justify)`, `(:top, :bottom, :center)`
respectively.

## Methods

The boxcontainer is viewed as a queue. These methods allow on to place children, as though into a vector

* `push!`
* `insert!`
* `unshift!`
* `append!`
* `prepend!`
* `splice!`
* `delete!`
* `pop!`
* `shift!`

""" ->
function boxcontainer(parent::Container; direction::Symbol=:horizontal, kwargs...)
    ## Toolkit
    widget, block = boxcontainer(parent.toolkit, parent, direction)

    ##
    obj = BoxContainer(widget, block, parent, parent.toolkit, {}, Dict())

    obj.attrs[:direction] = direction
    obj.attrs[:spacing] = [2,2]

    for (k, v) in kwargs
        obj[k] = v
    end
    obj
end

@doc "short form of `boxcontainer`" ->
box(parent::Container, direction=:h; kwargs...) = boxcontainer(parent, direction=direction==:h ? :horizontal : :vertical, kwargs...)

@doc "short form for a horizontal `boxcontainer`" ->
hbox(parent::Container; kwargs...) = boxcontainer(parent, direction=:horizontal, kwargs...)

@doc "short form for a vertical `boxcontainer`" ->
vbox(parent::Container; kwargs...) = boxcontainer(parent, direction=:vertical, kwargs...)

@doc "Returns number of children in box container" ->
length(object::BoxContainer) = length(children(object))

@doc "Remove all children in box container" ->
clear(object::BoxContainer)  = splice!(object, 1:length(object))

@doc "add strut (minimum height or width) to box container" ->
addstrut(parent::BoxContainer, px::Int) = addstruct(parent.toolkit, parent, px)

@doc "addStretch a stretching blank box" ->
addstretch(parent::BoxContainer, stretch::Int) = addstretch(parent.toolkit, parent, stretch)

@doc "addSpacing add blank space to layout" ->
addspacing(parent::BoxContainer, spacing::Int) = addspacing(parent.toolkit, parent, spacing)

## properties
## may be dynamic and adjust for just current chidren
## spacing around each child
getSpacing(object::BoxContainer) = object.attrs[:spacing]
setSpacing(object::BoxContainer, px::Int) = setSpacing(object, [px,px])
function setSpacing(object::BoxContainer, px::Vector{Int})
    object.attrs[:spacing] = px
    setSpacing(object.toolkit, object, px)
end
## margin around insider of container
setMargin(object::BoxContainer, px::Int) = setMargin(object::BoxContainer, [px,px])
setMargin(object::BoxContainer, px::Vector{Int}) = setMargin(object.toolkit, object, px)

list_props(::@PROP("BoxContainer")) = {:margin  => "[x,y] area in pixels around interior of box",
                                       :spacing => "[x,y] padding around each child widget"
                                       }

## findin
## return index of child in parent or 0 if not present
function findin(child::Widget, parent::Container)
    kids = children(parent)
    n = length(kids)
    out = [1:n][child .== kids]
    length(out) == 0 ? 0 : out[1]
end

## add child to parent

## add child widget to box container"
function push!(parent::BoxContainer, child::Widget)
    insert!(parent, length(parent) + 1, child)
end

## insert child into box container. index in `1:(n+1)`
function insert!(parent::BoxContainer, index::Int, child::Widget)
    insert!(parent.children, index, child)
    insert_child(parent.toolkit, parent, index, child)
    child
end

## insert child into first position of a box container
unshift!(parent::BoxContainer, child::Widget) = insert!(parent, 1, child)

## append several children to a box container
function append!{T <: Widget}(parent::BoxContainer, children::Vector{T})
    for child in children
        push!(parent, child)
    end
    children
end

## prepend severak children to a box container
function prepend!{T <: Widget}(parent::BoxContainer, children::Vector{T})
    for child in reverse(children)
        insert!(parent, 1, child)
    end
    children
end
prepend!(parent::BoxContainer, child) = prepend!(parent, [child])

@doc "Remove a child from a box container by index" ->
function splice!(parent::Container, index::Int)
    child = children(parent)[index]
    splice!(parent.children, index)
    remove_child(parent.toolkit, parent, child)
    child
end

## remove a child
function delete!(parent::Container, child::Widget)
    splice!(parent, findin(child, parent))
end

## Remove last child from a container
function pop!(parent::Container)
    splice!(parent, length(parent))
end

## Remove first child from a box container
shift!(parent::Container) = splice!(parent, 1)

##################################################
##
## Grid container
##

type GridContainer <: Container
    o
    block
    parent
    toolkit
    children
    attrs::Dict
end

@doc """
Grid

a grid-like container

Conceptually a grid is treated like a matrix.

One can add children by index:
- `parent[i,j] = child`
- `parent[:,:] = [child1 child2; nothing child3]`

Children can span multiple cells if added as follows
- `parent[i:(i+k), j:(j+l)] = child`

children can be referenced by index, as in `parent[i,j]`.

Children can be removed via `pop!(parent, child)`.

the size policy and alignment policies are used to position children within a cell

Rows and columns can be configured to have a minimum height, width or relative size through
`row_minimum_height`, `column_minimum_width`, `row_stretch`, `column_stretch`.

## Examples:
```julia
w = window(title="grid")
g = grid(w); push!(w, g)
g[1:2,1:2] = button(g, "1:2,1:2")
g[1:2, 3] = slider(g, 1:10, orientation=:vertical, size=[-1, 200])
g[3, 1:3] = slider(g, 1:10, size=[200, -1])
```

""" ->
function grid(parent::Container; kwargs...)

    widget, block = grid(parent.toolkit, parent)
    obj = GridContainer(widget, block, parent, parent.toolkit, {}, Dict())

    obj.attrs[:size] = [0,0]
    obj.attrs[:spacing] = [2,2]
    for (k, v) in kwargs
        obj[k] = v
    end
    obj
end

## XXX Document, but should these relate to addStretch, addSpacing for box containers (they play the same role)
function column_minimum_width(object::GridContainer, j::Int, width::Int)
    column_minimum_width(object.toolkit, object, j, width)
end

function row_minimum_height(object::GridContainer, i::Int, height::Int)
    row_minimum_width(object.toolkit, object, i, height)
end

function column_stretch(object::GridContainer, j::Int, weight::Int)
    column_stretch(object.toolkit, object, j, weight)
end
function row_stretch(object::GridContainer, i::Int, weight::Int)
    row_stretch(object.toolkit, object, i, weight)
end

## Properties
getSpacing(object::GridContainer) = object.attrs[:spacing]
setSpacing(object::GridContainer, px::Int) = setSpacing(object, [px, px])
function setSpacing(object::GridContainer, px::Vector{Int})
    object.attrs[:spacing] = px
    setSpacing(object.toolkit, object, px)
end

list_props(::@PROP("GridContainer")) = {:spacing => "[x,y] padding around each child widget"
                                       }

## return number of rows and columns
@doc "number of rows and columns of a grid container" ->
size(object::GridContainer) = grid_size(object.toolkit, object)
size(object::GridContainer, i::Int) = grid_size(object.toolkit, object)[i]
endof(object::GridContainer) = prod(size(object))
ndims(object::GridContainer) = 2

@doc "remove all children from a grid container" ->
function clear(object::GridContainer)
    for child in children(object)
        remove_child(child.toolkit, child.parent, child)
    end
    object.children = {}
end

@doc """
add children to grid container using array notation

* `parent[i,j] = child` places child at grid location i, j
* `parent[1:4, 2:3] = child` will place over multiple cells
""" ->
function setindex!(parent::GridContainer, child::Widget, i::Union(Int, Range1), j::Union(Int, Range1))
    ## must check that parents agree
    if parent != child.parent
        error("Parents are not the same")
    end

    ## How to organize children here? Could be a matrix, could carry ind info.... Useful for getindex..
    push!(parent.children, child) ## just a vector for now
    grid_add_child(parent.toolkit, parent, child, i, j)
end

## place children through grid notation
##
## * `parent[:,:] = [w1 w2; nothing w3]` will place the three children in cells (1,1), (1,2) and (2,2).
function setindex!(parent::GridContainer, children::Array, i::Union(Int, Range1), j::Union(Int, Range1))
    nr, nc = size(parent)
    if length(i) == nr && length(j) == nc
        clear(parent)
    end

    nr, nc = size(children)
    for i in 1:nr, j in 1:nc
        child = children[i, j]
        if isa(child, Widget)
            parent[i, j] = child
        end
    end
end

## retrieve children from a grid container through array notation
function getindex(object::GridContainer, i::Int, j::Int)
    grid_get_child_at(object.toolkit, object, i, j)
end

function getindex(object::GridContainer, i::Union(Vector, Range1), j::Union(Vector, Range1))
    [object[ii, jj] for j in jj, ii in i]
end

@doc "Remove a child from a grid container" ->
function delete!(parent::GridContainer, child::Widget)
    filter!(x -> !(x == child), parent.children)
    remove_child(parent.toolkit, parent, child)
end

###
type FormLayout <: Container
    o
    block
    parent
    toolkit
    children
    child_labels
    attrs::Dict
end

@doc """
formlayout

container to layout children in simple "label: control" form

children are pushed onto the end with an additional label (which may be null):
- `push!(parent, child, label)`

The `getValue` method will return a dictionary of control values keyed by the labels.

## Examples:
```julia
w = window(title = "formlayout")
fl = formlayout(w); push!(w, fl)
ed = lineedit(fl, "lineedit 1")
push!(fl, ed, "lineedit 1")
ed = lineedit(fl, "lineedit 2")
push!(fl, ed, "lineedit 2")
b = button(fl, "click")
push!(fl, b, "")
connect(b, :clicked, () -> println(fl[:value])) ## current values
```

""" ->
function formlayout(parent::Container;  kwargs...)

    widget, block = formlayout(parent.toolkit, parent)

    obj = FormLayout(widget, block, parent, parent.toolkit, {}, {}, Dict())

    obj.attrs[:nrows] = 0
    obj.attrs[:spacing] = [5,2]
    obj[:sizepolicy] = (:expand, :expand)
    for (k, v) in kwargs
        obj[k] = v
    end
    obj
end

@doc "add children to a `formlayout` is done via `push!(lty, child, label)`:" ->
function push!(parent::FormLayout, child::Widget, label::Union(Nothing, String))
    push!(parent.children, child)
    push!(parent.child_labels, label)
    formlayout_add_child(parent.toolkit, parent, child, label)
end

push!(parent::FormLayout, child::Widget; label::String="") = push!(parent, child, label)

## return dictionary of control values
function getValue(widget::FormLayout)
    kids = children(widget)
    labels = widget.child_labels
    d = Dict()
    for i in 1:length(kids)
        if !isa(kids[i], Control) next end
        label = labels[i]
        key = isa(label, Nothing) ? string("control i") : label
        d[key] = getValue(kids[i])
    end
    d
end

## Properties
getSpacing(object::FormLayout) = object.attrs[:spacing]
setSpacing(object::FormLayout, px::Int) = setSpacing(object, [px, px])
function setSpacing(object::FormLayout, px::Vector{Int})
    object.attrs[:spacing] = px
    setSpacing(object.toolkit, object, px)
end

list_props(::@PROP("FormLayout")) = {:spacing => "[x,y] padding around each child widget",
                                     :value => "read only. Returns dictionary of widget values"}

##################################################
##
## notebook

type NoteBook <: Container
    o
    block
    parent
    toolkit
    children
    model
    attrs::Dict
end

@doc """
notebook container

## Arugments

* `parent::Container` parent container

## Signals

* `valueChanged (value)` value is currently selected tab index

## Methods

* `length` number of tables
* `insert!` insert at `i` with label i in {1, 2, ..., length + 1}

## TODO: rename label

## Examples:
```julia
w = window(title="notebook")
nb = notebook(w); push!(w, nb)

push!(nb, button(nb, "one"), "one")
push!(nb, button(nb, "two"), "two")
b =  button(nb, "pop! last child")
push!(nb, b, "three")

connect(nb, :valueChanged, value -> println("On tab ", value))
connect(b, :clicked, () -> pop!(nb))
```

""" ->
function notebook(parent::Container, kwargs...)
    model = ItemModel(0)
    widget, block = notebook(parent.toolkit, parent, model)

    obj = NoteBook(widget, block, parent, parent.toolkit, {}, model, Dict())

    for (k, v) in kwargs
        obj[k] = v
    end
    obj
end

## interface
length(object::NoteBook) = length(object.children)
function setValue(object::NoteBook, i::Int; signal::Bool=true)
    if i < 1 || i > length(object)
        return
    end
    setValue(object.model, i; signal=signal)
end
getValue(object::NoteBook) = getValue(object.model)

function insert!(parent::NoteBook,  i::Int, child::Widget, label::String)
    if i < 0 || i > length(parent) + 1
        error("Index needs to be between 1 and n+1")
    end

    notebook_insert_child(parent.toolkit, parent, child, i, label)
    insert!(parent.children, i, child)

end
push!(parent::NoteBook, child::Widget, label::String) = insert!(parent,  length(parent) + 1, child, label)
unshift!(parent::NoteBook, child::Widget, label::String) = insert!(parent, 1, child, label)

## Remove a child by index
function splice!(parent::NoteBook, index::Int)
    child = children(parent)[index]
    notebook_remove_child(parent.toolkit, parent, child)
    splice!(parent.children, index)
    child
end
pop!(parent::NoteBook) = splice!(parent, length(parent))
shift!(parent::NoteBook) = splice!(parent, 1)

pop!(parent::NoteBook, child::Widget) = splice!(parent, findin(child, parent))

## panedgroup
## expandgroup...


text:
# violin

#Jl Plots

using Plots

y = rand(100, 4) # Four series of 100 points each

Plots.violin(["Series 1" "Series 2" "Series 3" "Series 4"], y, alpha=.7, leg = false)

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
"""
Produce an elliptic torus with circular tube cross
section
"""
module MajorEllipticTorus

# minimum distance between a point and ellipse
function min_distance(a, b, coordinates)
    x = abs(coordinates[1])
    y = abs(coordinates[2])
    z = coordinates[3]
    function aux(phi, accum)
        new_phi = atan(((a^2 - b^2)*sin(phi) + y*b)/(x*a))
        if abs(phi - new_phi) < 1e-4
            return new_phi
        else
            return aux(new_phi, accum+1)
        end
    end
    phi = aux(0, 0)
    return sqrt(z^2 + (x - a * cos(phi))^2 + (y - b*sin(phi))^2)
end

# check whether a point is inside the torus

function get_fn(argd::Dict)
    a = float(argd["major"]["a"])
    b = float(argd["major"]["b"])
    r = float(argd["minor"]["radius"])
    center = [float(x) for x in argd["center"]]
    force_constant = float(argd["k"])

    function fn_compute(raw_coordinates)
        coordinates = raw_coordinates - center
        if norm(coordinates[1:2]) < 1e-6
            d = sqrt(min(a,b)^2 + coordinates[3]^2)
        else
            d = min_distance(a, b, coordinates)
        end
        if d < r
            return 0
        else
            return force_constant * d
        end
    end
    return fn_compute
end

export get_fn

end

max_stars_repo_name: yuhangwang/SimShape

max_issues_repo_name: yuhangwang/SimShape

max_forks_repo_name: yuhangwang/SimShape


text:
module DistributionsCommonTest

using Test
using ReactiveMP
using Distributions
using Random

@testset "Distributions Common" begin

    # Categorical comes from Distributions.jl and most of the things should be covered there
    # Here we test some extra ReactiveMP.jl specific functionality

    @testset "Bernoulli × Categorical" begin
        @test prod(ProdAnalytical(), Bernoulli(0.3), Categorical([ 1/2, 1/2 ])) ≈ prod(ProdAnalytical(), Bernoulli(0.3), Bernoulli(0.5))
        @test prod(ProdAnalytical(), Bernoulli(0.1), Categorical([ 0.2, 0.8 ])) ≈ prod(ProdAnalytical(), Bernoulli(0.1), Bernoulli(0.2))
        @test prod(ProdAnalytical(), Bernoulli(0.9), Categorical([ 0.8, 0.2 ])) ≈ prod(ProdAnalytical(), Bernoulli(0.9), Bernoulli(0.8))
    end

end

end


text:
using CUDA
using FailureRepresentation
using FailureRepresentation.SystemUnderTest
using FailureRepresentation.Autoencoder
using Flux

sut = SystemUnderTest.load("models/sut.bson")
autoencoder = Autoencoder.load("models/autoencoder.bson")
encoder = Autoencoder.load("models/encoder.bson", :encoder)


text:
render(::Inline, m::Model) = m.data

view(x::AbstractString) = x
view(x::AbstractDict) = x

view(n::Node) =
  d(:type     => :dom,
    :tag      => tag(n),
    :attrs    => n.attrs,
    :contents => map(view, n.children))

view(n::Node{:code}) =
  d(:type  => :code,
    :attrs => n.attrs,
    :text  => join(n.children, '\n'))

view(n::Hiccup.Node{:latex}) =
  Dict(:type  => :latex,
       :attrs => n.attrs,
       :text  => join(n.children, ' '))

render(::Inline, n::Node) = view(n)

render(::Inline, x::HTML) = view(x)

render(::Console, x::Node) =
  @msg result(view(x))

render(::Editor, x::Node) = view(x)

render(i::Inline, t::Tree) =
  isempty(t.children) ?
    render(i, t.head) :
    d(:type     => :tree,
      :head     => render(i, t.head),
      :children => map(x->render(i, x),
                       t.children))

render(i::Inline, t::SubTree) =
  d(:type  => :subtree,
    :label => render(i, t.label),
    :child => render(i, t.child))

@render i::Inline t::Table begin
  xs = map(x -> render(i, x), t.xs)
  table("", vec(mapslices(xs -> tr("", map(x->td("", [x]), xs)), xs, 2)))
end

render(i::Inline, x::Copyable) =
  d(:type => :copy,
    :view => render(i, x.view),
    :text => x.text)

render(i::Inline, l::Link) =
  d(:type => :link,
    :file => l.file,
    :line => l.line-1,
    :contents => map(x->render(i, x), l.contents))

render(::Clipboard, l::Link) =
  "$(l.file):$(l.line)"


text:
#На прямоугольном поле произвольных размеров расставить маркеры в виде "шахматных" клеток,
#начиная с юго-западного угла поля, когда каждая отдельная "шахматная" клетка имеет размер n x n клеток поля
#(n - это параметр функции). Начальное положение Робота - произвольное, конечное - совпадает с начальным.
#Клетки на севере и востоке могут получаться "обрезанными" - зависит от соотношения размеров поля и "шахматных" клеток.
#(Подсказка: здесь могут быть полезными две глобальных переменных,
#в которых будут содержаться текущие декартовы координаты Робота относительно начала координат в левом нижнем углу поля, например)

function moves!(r,side) #двигает робота и возвращает количество шагов
    num = 0
    while isborder(r,side) == false
        move!(r,side)
        num+=1
    end
    return num
end

function reverse(side)
    for i=0:3
        if side == HorizonSide(i)
            return HorizonSide(mod(i+2, 4))
        end
    end
end

function mark_chek(r, x_coord, y_coord, size) #проверка координаты на наличие маркера(size-размер отступа)
    if ((mod(x_coord,2*size) in 0:size-1)&&(mod(y_coord,2*size) in 0:size-1))
        putmarker!(r)
    end
    if ((mod(x_coord,2*size) in size:2*size-1)&&(mod(y_coord,2*size) in size:2*size-1))
        putmarker!(r)
    end
end

function move_back!(r,side,num) #передвигает робота в сторону на число шагов
    for i = 1:num
        move!(r,side)
    end
end

function move_coord!(r,side) #двигает робота и отсчитывает координаты
    global x_coord, y_coord
    if side == Nord
        y_coord += 1
    end
    if side == Sud
        y_coord -= 1
    end
    if side == Ost
        x_coord += 1
    end
    if side == West
        x_coord -= 1
    end
    move!(r,side)
end

function mark_move!(r, side) #функция двигает робота и ставит метки
    global x_coord, y_coord
    mark_chek(r, x_coord, y_coord, size)
    while isborder(r,side) == false
        move_coord!(r, side)
        mark_chek(r, x_coord, y_coord, size)
    end
end

function marking(r)
    mark_move!(r, Ost)
    side = West
    while isborder(r, Nord) == false
        move_coord!(r, Nord)
        mark_move!(r, side)
        side = reverse(side)
    end
end

x_coord = 0
y_coord = 0
size = 0
function chest_n(r, n) # - главная функция
    global size, x_coord, y_coord
    size = n
    y = moves!(r, Sud)
    x = moves!(r, West)
    marking(r)
    move_back!(r, Sud, y_coord - y)
    move_back!(r, West, x_coord - x)
end

max_stars_repo_name: gardend/mirea-progs

max_issues_repo_name: gardend/mirea-progs

max_forks_repo_name: gardend/mirea-progs