text:
using DifferentialEquations
using Plots; gr()
using Formatting

function logistic(r,x)
    r*x*(1.0-x)
end

function logisticmap(;x0=0.21,r=4,n=10)
    results = [(x0,logistic(r,x0))];
    #N = 100; Note this line should not be here
    for i in 1:(n-1) # Note change of case N -> n.  Why?
      xold = results[end][2]
      xnew = logistic(r,xold)
      append!(results,[(xold,xold)])
      append!(results,[(xold,xnew)])
  end
  results
end

function vary(dname;x0_vals=[],r_vals=[], nvals=[])
 if !isdir(dname)
  mkdir(dname)
 end
 fno = 1
 for x0 in x0_vals
  for r in r_vals
      for n in nvals
          fname = format(dname*"/frame{:05d}.png",fno)
          fno = fno + 1
          simres = logisticmap(x0=x0, r=r, n=n)
          xvals = [x for (x, y) in simres]
          yvals = [y for (x, y) in simres]
          plot(xvals,yvals,legend=false,xaxis=false,yaxis=false,xlim=(0,1),ylim=(0,1));
          savefig(fname)
        end
    end
end
end

# Figure out how to make/deleate directories and remove their contents

r = 4.0
x = (r+sqrt(r-4)*sqrt(r))/(2*r)
vary("xn+2";x0_vals=(x+0.00001),r_vals=r, nvals=1:50)


text:
info("Testing ", KernelAffinity)
for kernelobj in (additive_kernels..., composition_kernels...)
    for T in FloatingPointTypes

        k1 = convert(Kernel{T}, (kernelobj)())
        a = convert(T,2)
        c = convert(T,3)

        k = KernelAffinity(a, c, k1)
        @test k.a.value == a
        @test k.c.value == c

        @test eltype(convert(KernelAffinity{Float16},  k)) == Float16
        @test eltype(convert(KernelAffinity{Float32},  k)) == Float32
        @test eltype(convert(KernelAffinity{Float64},  k)) == Float64
        @test eltype(convert(KernelAffinity{BigFloat}, k)) == BigFloat

        @test MOD.attainszero(k) == MOD.attainszero(k.kappa)
        @test MOD.attainspositive(k) == MOD.attainspositive(k.kappa)
        @test MOD.attainsnegative(k) == MOD.attainsnegative(k.kappa)
        @test ismercer(k) === ismercer(k.kappa)
        @test isnegdef(k) === isnegdef(k.kappa)

        @test_approx_eq MOD.phi(k, one(T)) (a*one(T) + c)

        k = k1 + c
        @test k.a == one(T)
        @test k.c == c
        @test k.kappa == k1

        k = c + ((c + k1) + c)
        @test k.a == one(T)
        @test k.c == 3c
        @test k.kappa == k1

        k = k1 * a
        @test k.a == a
        @test k.c == zero(T)
        @test k.kappa == k1

        k = a * ((a * k1) * a)
        @test k.a == a^3
        @test k.c == zero(T)
        @test k.kappa == k1

        k = (a * k1) + c
        @test k.a == a
        @test k.c == c
        @test k.kappa == k1

        k = c + (k1 * a)
        @test k.a == a
        @test k.c == c
        @test k.kappa == k1

        k = a * (a * k1 + c)
        @test k.a == a^2
        @test k.c == a*c
        @test k.kappa == k1

        k = (a * k1 + c) * a
        @test k.a == a^2
        @test k.c == a*c
        @test k.kappa == k1

       #= k = a * convert(Kernel{T}, SquaredDistanceKernel()) + c
        @test k^(one(T)/2) == KernelComposition(PowerClass(a, c, one(T)/2), k.kappa)

        k = a * convert(Kernel{T}, ScalarProductKernel()) + c
        @test k^3     == KernelComposition(PolynomialClass(a, c, 3one(T)), k.kappa)
        @test exp(k)  == KernelComposition(ExponentiatedClass(a, c), k.kappa)
        @test tanh(k) == KernelComposition(SigmoidClass(a, c), k.kappa)=#

        # Test that output does not create error
        @test show(DevNull, k) == nothing

    end
end

info("Testing ", KernelSum)
for kernelobj1 in (additive_kernels..., composition_kernels...)
    for kernelobj2 in (additive_kernels..., composition_kernels...)
        for T in FloatingPointTypes
            k1 = convert(Kernel{T}, (kernelobj1)())
            k2 = convert(Kernel{T}, (kernelobj2)())
            if ismercer(k1) && ismercer(k2) || isnegdef(k1) && isnegdef(k2)
                k = k1 + k2
                @test k.c.value == zero(T)
                @test k.kappa1 == k1
                @test k.kappa2 == k2

                @test eltype(convert(Kernel{Float32}, k))  == Float32
                @test eltype(convert(Kernel{Float64}, k))  == Float64
                @test eltype(convert(Kernel{BigFloat}, k)) == BigFloat

                @test ismercer(k) == (ismercer(k1) && ismercer(k2))
                @test isnegdef(k) == (isnegdef(k1) && isnegdef(k2))

                # Test that output does not create error
                @test show(DevNull, k) == nothing

                c = convert(T,2)

                k = (k1 + c) + (k2 + c)
                @test k.c.value == 2c
                @test k.kappa1 == k1
                @test k.kappa2 == k2

                k = (k1 + c) + k2
                @test k.c == c
                @test k.kappa1 == k1
                @test k.kappa2 == k2

                k = k1 + (k2 + c)
                @test k.c == c
                @test k.kappa1 == k1
                @test k.kappa2 == k2

                k = 2k1 + 2k2
                @test k.c.value == zero(T)
                @test k.kappa1 == 2k1
                @test k.kappa2 == 2k2

                k = k1 + 2k2
                @test k.c.value == zero(T)
                @test k.kappa1 == k1
                @test k.kappa2 == 2k2

                k = 2k1 + k2
                @test k.c.value == zero(T)
                @test k.kappa1 == 2k1
                @test k.kappa2 == k2

                k = (k1 + k2) + 1
                @test k.c == one(T)
                @test k.kappa1 == k1
                @test k.kappa2 == k2

            else
                @test_throws ErrorException k1 + k2
            end
        end
    end
end

info("Testing ", KernelProduct)
for kernelobj1 in (additive_kernels..., composition_kernels...)
    for kernelobj2 in (additive_kernels..., composition_kernels...)
        for T in FloatingPointTypes
            k1 = convert(Kernel{T}, (kernelobj1)())
            k2 = convert(Kernel{T}, (kernelobj2)())
            if ismercer(k1) && ismercer(k2)
                k = k1 * k2
                @test k.a.value == one(T)
                @test k.kappa1 == k1
                @test k.kappa2 == k2

                @test eltype(convert(Kernel{Float32}, k))  == Float32
                @test eltype(convert(Kernel{Float64}, k))  == Float64
                @test eltype(convert(Kernel{BigFloat}, k)) == BigFloat

                @test ismercer(k) == (ismercer(k1) && ismercer(k2))

                # Test that output does not create error
                @test show(DevNull, k) == nothing

                a = convert(T,2)

                k = (k1 * a) * (k2 * a)
                @test k.a.value == a^2
                @test k.kappa1 == k1
                @test k.kappa2 == k2

                k = (k1 * a) * k2
                @test k.a == a
                @test k.kappa1 == k1
                @test k.kappa2 == k2

                k = k1 * (k2 * a)
                @test k.a == a
                @test k.kappa1 == k1
                @test k.kappa2 == k2

                k = (k1 + 1) * (k2 + 1)
                @test k.a.value == one(T)
                @test k.kappa1 == (k1 + 1)
                @test k.kappa2 == (k2 + 1)

                k = k1 * (k2 + 1)
                @test k.a.value == one(T)
                @test k.kappa1 == k1
                @test k.kappa2 == (k2 + 1)

                k = (k1 + 1) * k2
                @test k.a.value == one(T)
                @test k.kappa1 == (k1 + 1)
                @test k.kappa2 == k2

                k = (k1 * k2) * 2
                @test k.a == convert(T,2)
                @test k.kappa1 == k1
                @test k.kappa2 == k2

            else
                @test_throws ErrorException k1 * k2
            end
        end
    end
end

ext: jl


text:
import Term: tprint, highlight, tprintln

print("\n\n"^10)

txt = """
{bold underline bright_blue}Term.jl{/bold underline bright_blue} highlights numbers: 123 1.5e10 1/2 and "strings too"

Operators are colored red: + - * % ^ / while symbols are orange: :x :y :z

Code looks like this:
```
    do α
    do β
```

And you can highlight expression like :(x + y) too.
"""
tprintln(txt)

# tprintln("1, 2, 3")

# tprintln(zeros(4); highlight=true)

# tprintln(Dict(:name=>"test", 1=>:v, :k=>[1, 2, 3]))

# x, y = :z, :(x+Y)
# tprintln("These are the results: x = $(highlight(x)), y = $(highlight(y))")src/layout.jl
# println()

# tprintln(highlight_syntax("""
# import Term

# tprintln("ok")
# """))


text:
Δt = 0.01
h = 0.5

ρ = 1
μ = 0.001
ν = μ / ρ

# ϕϕ = 1/9
# ϕϕx = 0.0
# ϕxϕx = 4/3

# # LHS * u = RHS
# # LHS = [ϕϕ 0.0;0.0 ϕϕ]/Δt +

# sol = - 2 * μ * ϕxϕx / (ϕϕ/Δt + 2 * μ * ϕxϕx)


text:
module FractionalCalculus

using QuadGK
using SpecialFunctions: gamma
using LinearAlgebra, InvertedIndices, SpecialMatrices

include("./Derivative/Derivative.jl")
include("./Integral/Integral.jl")

export FracDiffAlg, FracIntAlg
export fracdiff, fracint

# Export fractional derivative releating API
export Caputo, GL, RLDiff, Hadamard

# Caputo sense fractional derivative
export Caputo_Direct, Caputo_Direct_First_Diff_Known, Caputo_Direct_First_Second_Diff_Known
export Caputo_Piecewise
export Caputo_Diethelm
export Caputo_High_Precision

# Grünwald Letnikov sense fractional derivative
export GL_Direct, GL_Multiplicative_Additive, GL_Lagrange_Three_Point_Interp, GL_High_Precision

# Riemann Liouville sense fractional derivative
export RLDiff_Approx, RLDiff_Matrix, RL_Linear_Spline_Interp

# Hadamard sense fractional derivative
export Hadamard_LRect, Hadamard_RRect, Hadamard_Trap

# Riesz sense fractional derivative
export Riesz_Symmetric

# Macros for convenient computing
export @fracdiff, @semifracdiff

#Export fractional integral relating API
export RLInt
export RL_Direct, RL_Direct_First_Diff_Known, RL_Piecewise
export RLInt_Approx, RL_LinearInterp, RLInt_Matrix
export RLInt_Simpson, RLInt_Trapezoidal, RLInt_Rectangular, RLInt_Cubic_Spline_Interp

# Macros for convenient computing
export @fracint, @semifracint

end


text:
using Test, Random
include("tmp.jl")
include("src/types.jl")
include("src/params.jl")

using .Types
using .Params
using .Tmp

# TODO: Create keybinding "switch to console, execute last command, switch back"
# TODO: Create keybinding "paste selection in console, execute, switch back"

stockreturns = Func.assetreturns(stocks.value)

nstocks = size(stockreturns,1)
kurtoses = zeros(nstocks)
for stock in 1:size(stockreturns,1)
        kurtoses[stock] = StatsBase.kurtosis(stockreturns[stock,:])
end


text:
"""
    estimation_input_test(estim::VARSettings, benchmark_data::Tuple)

Check whether the settings in `estim` are equal to those in `benchmark_data`.

    estimation_input_test(estim_a::VARSettings, estim_b::VARSettings)

Check whether the settings in `estim_a` are equal to those in `estim_b`.
"""
function estimation_input_test(estim::VARSettings, benchmark_data::Tuple)

    benchmark_Y, benchmark_n, benchmark_T, benchmark_q, benchmark_nq, benchmark_m, benchmark_λ,
        benchmark_α, benchmark_β, benchmark_Γ, benchmark_ε, benchmark_tol, benchmark_max_iter,
            benchmark_prerun, benchmark_verb = benchmark_data;

    return ~(false in [estim.Y === benchmark_Y;
                      estim.n == benchmark_n;
                      estim.T == benchmark_T;
                      estim.q == benchmark_q;
                      estim.nq == benchmark_nq;
                      estim.m == benchmark_m;
                      estim.λ == benchmark_λ;
                      estim.α == benchmark_α;
                      estim.β == benchmark_β;
                      estim.Γ == benchmark_Γ;
                      estim.ε == benchmark_ε;
                      estim.tol == benchmark_tol;
                      estim.max_iter == benchmark_max_iter;
                      estim.prerun == benchmark_prerun;
                      estim.verb == benchmark_verb]);
end

function estimation_input_test(estim_a::VARSettings, estim_b::VARSettings)

    return ~(false in [estim_a.Y === estim_b.Y;
                      estim_a.n == estim_b.n;
                      estim_a.T == estim_b.T;
                      estim_a.q == estim_b.q;
                      estim_a.nq == estim_b.nq;
                      estim_a.m == estim_b.m;
                      estim_a.λ == estim_b.λ;
                      estim_a.α == estim_b.α;
                      estim_a.β == estim_b.β;
                      estim_a.Γ == estim_b.Γ;
                      estim_a.ε == estim_b.ε;
                      estim_a.tol == estim_b.tol;
                      estim_a.max_iter == estim_b.max_iter;
                      estim_a.prerun == estim_b.prerun;
                      estim_a.verb == estim_b.verb]);
end

"""
    var_estimation_test(Y::JArray, q::Int64, λ::Number, α::Number, β::Number, benchmark_data::Tuple)

Run basic tests on VARSettings.
"""
function var_estimation_test(Y::JArray, q::Int64, λ::Number, α::Number, β::Number, benchmark_data::Tuple)

    #=
    Tests on VARSettings
    =#

    # First group of settings

    estim1 = VARSettings(Y, q, λ, α, β);
    estim2 = VARSettings(Y, q, λ, α, β, ε=1e-4);
    estim3 = VARSettings(Y, q, λ, α, β, ε=1e-4, tol=1e-4);
    estim4 = VARSettings(Y, q, λ, α, β, ε=1e-4, tol=1e-4, max_iter=1000);
    estim5 = VARSettings(Y, q, λ, α, β, ε=1e-4, tol=1e-4, max_iter=1000, prerun=2);
    estim6 = VARSettings(Y, q, λ, α, β, ε=1e-4, tol=1e-4, max_iter=1000, prerun=2, verb=true);

    @test estimation_input_test(estim1, benchmark_data);
    @test estimation_input_test(estim1, estim2);
    @test estimation_input_test(estim1, estim3);
    @test estimation_input_test(estim1, estim4);
    @test estimation_input_test(estim1, estim5);
    @test estimation_input_test(estim1, estim6);

    # Second group of settings

    estim1b = VARSettings(Y, q, λ, α, β, ε=1e-8);
    estim2b = VARSettings(Y, q, λ, α, β, ε=1e-8, tol=1e-4);
    estim3b = VARSettings(Y, q, λ, α, β, ε=1e-8, tol=1e-4, max_iter=1000);
    estim4b = VARSettings(Y, q, λ, α, β, ε=1e-8, tol=1e-4, max_iter=1000, prerun=2);
    estim5b = VARSettings(Y, q, λ, α, β, ε=1e-8, tol=1e-4, max_iter=1000, prerun=2, verb=true);

    @test estimation_input_test(estim1b, estim2b);
    @test estimation_input_test(estim1b, estim3b);
    @test estimation_input_test(estim1b, estim4b);
    @test estimation_input_test(estim1b, estim5b);

    # Third group of settings

    estim1c = VARSettings(Y, q, λ, α, β, ε=1e-8, tol=1e-16);
    estim2c = VARSettings(Y, q, λ, α, β, ε=1e-8, tol=1e-16, max_iter=1000);
    estim3c = VARSettings(Y, q, λ, α, β, ε=1e-8, tol=1e-16, max_iter=1000, prerun=2);
    estim4c = VARSettings(Y, q, λ, α, β, ε=1e-8, tol=1e-16, max_iter=1000, prerun=2, verb=true);

    @test estimation_input_test(estim1c, estim2c);
    @test estimation_input_test(estim1c, estim3c);
    @test estimation_input_test(estim1c, estim4c);

    # Fourth group of settings

    estim1d = VARSettings(Y, q, λ, α, β, ε=1e-8, tol=1e-16, max_iter=10000);
    estim2d = VARSettings(Y, q, λ, α, β, ε=1e-8, tol=1e-16, max_iter=10000, prerun=2);
    estim3d = VARSettings(Y, q, λ, α, β, ε=1e-8, tol=1e-16, max_iter=10000, prerun=2, verb=true);

    @test estimation_input_test(estim1d, estim2d);
    @test estimation_input_test(estim1d, estim3d);

    # Fifth group of settings

    estim1e = VARSettings(Y, q, λ, α, β, ε=1e-8, tol=1e-16, max_iter=10000, prerun=3);
    estim2e = VARSettings(Y, q, λ, α, β, ε=1e-8, tol=1e-16, max_iter=10000, prerun=3, verb=true);

    @test estimation_input_test(estim1e, estim2e);

    # Sixth group of settings

    estim1f = VARSettings(Y, q, λ, α, β, ε=1e-8, tol=1e-16, max_iter=10000, prerun=3, verb=false);

    @test estim1f.verb == false;

    #=
    Set default settings
    =#

    estim = estim1;

    #=
    Tests on initialisation
    =#

    sspace, B_star, C_star, Q_view, coordinates_measurement_states, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs, coordinates_transition_P0, coordinates_free_params_B, coordinates_free_params_C = MessyTimeSeriesOptim.initialise(estim);

    @test sspace.B[:, 1:estim.n] == Matrix(1.0I, estim.n, estim.n);
    @test sspace.B[:, estim.n+1:end] == zeros(estim.n, estim.m);
    @test sspace.C[1:estim.n, 1:estim.m] == C_star;
    @test sspace.C[1:estim.n, estim.m+1:end] == zeros(estim.n, estim.n);

    #=
    Tests on ksmoother_ecm!
    =#

    # Run Kalman filter
    status = SizedKalmanStatus(sspace.Y.T);
    kfilter_full_sample!(sspace, status);

    # Compute Kalman stats from the ksmoother function
    X_sm, P_sm, X0_sm, P0_sm = ksmoother(sspace, status);

    # Memory pre-allocation
    smoother_arrays_benchmark = MessyTimeSeriesOptim.SmootherArrays(estim, sspace, coordinates_measurement_states, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs);

    # Compute F_benchmark, G_benchmark and H_benchmark
    for t=length(X_sm)-1:-1:1
        copyto!(smoother_arrays_benchmark.Xs_leading, X_sm[t+1]);
        copyto!(smoother_arrays_benchmark.Ps_leading, P_sm[t+1]);
        Xs = X_sm[t];
        Ps = P_sm[t];
        MessyTimeSeriesOptim.update_ecm_stats_transition!(smoother_arrays_benchmark.F, smoother_arrays_benchmark.G, estim, smoother_arrays_benchmark, Xs, Ps, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs);
    end

    copyto!(smoother_arrays_benchmark.Xs_leading, X_sm[1]);
    copyto!(smoother_arrays_benchmark.Ps_leading, P_sm[1]);
    MessyTimeSeriesOptim.update_ecm_stats_transition!(smoother_arrays_benchmark.F, smoother_arrays_benchmark.G, estim, smoother_arrays_benchmark, X0_sm, P0_sm, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs);

    # Compute Kalman stats from the ksmoother_ecm! function
    smoother_arrays = MessyTimeSeriesOptim.SmootherArrays(estim, sspace, coordinates_measurement_states, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs);
    MessyTimeSeriesOptim.ksmoother_ecm!(estim, sspace, status, smoother_arrays, coordinates_measurement_states, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs, coordinates_transition_P0);

    # Tests
    @test smoother_arrays.F == smoother_arrays_benchmark.F;
    @test smoother_arrays.G == smoother_arrays_benchmark.G;
    @test smoother_arrays.H == smoother_arrays_benchmark.H;
    @test sspace.X0 == X0_sm;
    @test sspace.P0 == P0_sm;
end

"""
    naive_var_simulation(true_coefficient::FloatMatrix, n::Int64, T::Int64, q::Int64, λ::Number, α::Number, β::Number; no_sim::Int64=500, verb::Bool=false)

Run a VAR(q) naive simulation.
"""
function naive_var_simulation(true_coefficient::FloatMatrix, n::Int64, T::Int64, q::Int64, λ::Number, α::Number, β::Number; no_sim::Int64=500, verb::Bool=false)

    # Run simulation
    rng = StableRNG(1);
    C_sim = zeros(n, n*q, no_sim);
    params_accuracy = zeros(no_sim);
    oos_fc_error = zeros(no_sim);

    for sim=1:no_sim

        # Simulate data with mean zero and unit standard deviation
        v = randn(rng, n, T+q);
        Y = zeros(n, T+q);
        for t=1:T+q
            if t <= q
                Y[:,t] .= 0.0;
            else
                Y[:,t] = true_coefficient*vcat([Y[:,t-lag][:] for lag=1:q]...) + v[:,t];
            end
        end
        Y = standardise(Y[:,q+1:end]);

        # Estimate
        t0 = floor(T/2) |> Int64;
        estim = VARSettings(Y[:,1:t0], q, λ, α, β, verb=false);
        sspace = ecm(estim, output_sspace_data=Y);

        # Run Kalman filter
        status = kfilter_full_sample(sspace);

        # Forecast
        forecast = mapreduce(Xt -> sspace.B*Xt, hcat, status.history_X_prior);

        # Store estimated parameters
        C_sim[:,:,sim] = sspace.C[1:n, 1:n*q];

        # Store params_accuracy
        difference_C_truth = @views C_sim[:,:,sim] - true_coefficient;
        params_accuracy[sim] = sqrt(tr(difference_C_truth'*difference_C_truth));

        # Store oos_fc_error
        oos_fc_error[sim] = @views mean(mean(abs.(Y[:,t0+1:end] - forecast[:,t0+1:end]), dims=2)./mean(abs.(Y[:,t0+1:end]), dims=2));
    end

    if verb == true
        println("=============================================================================================================================================")
        println("True parameters:")
        println(true_coefficient);
        println("=============================================================================================================================================")
        println("Estimated parameters with λ=$(λ), α=$(α), and β=$(β):")
        println(round.(median(C_sim, dims=3)[:,:], digits=5));
        println("=============================================================================================================================================")
        println("");
    end

    return C_sim, params_accuracy, oos_fc_error;
end

"""
    var_simulation_tests(model_id::String, λ::Number, α::Number, β::Number)

Default tests based on VAR simulations.
"""
function var_simulation_tests(model_id::String, λ::Number, α::Number, β::Number)

    # Run simulations
    _, bivariate_var1_full, _ = naive_var_simulation([0.3 -0.2; -0.2 0.5], 2, 200, 1, λ, α, β, no_sim=10);
    _, bivariate_var2_full, _ = naive_var_simulation([0.5 0.4 0.2 0.1; 0.3 -0.2 0.1 0.1], 2, 200, 2, λ, α, β, no_sim=10);
    _, bivariate_var2_sparse, _ = naive_var_simulation([0.3 0.0 0.0 0.0; 0.0 -0.1 0.0 0.0], 2, 200, 2, λ, α, β, no_sim=10);
    _, trivariate_var2_groups, _ = naive_var_simulation([0.5 0.4 0.0 0.2 0.1 0.0; 0.3 -0.2 0.0 0.1 0.1 0.0; 0.0 0.0 0.9 0.0 0.0 0.0], 3, 200, 2, λ, α, β, no_sim=10);

    # Load benchmarks
    benchmark_bivariate_var1_full = read_test_input("./input/var/$(model_id)_bivariate_var1_full");
    benchmark_bivariate_var2_full = read_test_input("./input/var/$(model_id)_bivariate_var2_full");
    benchmark_bivariate_var2_sparse = read_test_input("./input/var/$(model_id)_bivariate_var2_sparse");
    benchmark_trivariate_var2_groups = read_test_input("./input/var/$(model_id)_trivariate_var2_groups");

    # Perform tests
    @test norm(bivariate_var1_full-benchmark_bivariate_var1_full, Inf) <= 1e-8;
    @test norm(bivariate_var2_full-benchmark_bivariate_var2_full, Inf) <= 1e-8;
    @test norm(bivariate_var2_sparse-benchmark_bivariate_var2_sparse, Inf) <= 1e-8;
    @test norm(trivariate_var2_groups-benchmark_trivariate_var2_groups, Inf) <= 1e-8;
end

@testset "VAR basic functionalities" begin

    # Initialise data and state-space parameters
    Y = [0.72 missing 1.86 missing missing 2.52 2.98 3.81 missing 4.36;
         0.95 0.70 missing missing missing missing 2.84 3.88 3.84 4.63];

    q = 2;
    nq = 4;
    m = 4;
    λ = 0.0;
    α = 0.5;
    β = 1.0;
    ε = 1e-4;
    tol = 1e-4;
    max_iter = 1000;
    prerun = 2;
    verb = true;

    Γ_1 = build_Γ(2, q, λ, β);
    Γ_2 = build_Γ(2, q, λ+1, β);
    Γ_3 = build_Γ(2, q, λ+1, β+1);

    # Benchmark data
    benchmark_data_1 = (Y, size(Y,1), size(Y,2), q, nq, m, λ, α, β, Γ_1, ε, tol, max_iter, prerun, verb);
    benchmark_data_2 = (Y, size(Y,1), size(Y,2), q, nq, m, λ+1, α, β, Γ_2, ε, tol, max_iter, prerun, verb);
    benchmark_data_3 = (Y, size(Y,1), size(Y,2), q, nq, m, λ+1, α, β+1, Γ_3, ε, tol, max_iter, prerun, verb);

    # Run tests
    var_estimation_test(Y, q, λ, α, β, benchmark_data_1);
    var_estimation_test(Y, q, λ+1, α, β, benchmark_data_2);
    var_estimation_test(Y, q, λ+1, α, β+1, benchmark_data_3);
end

@testset "VAR simulations: MLE" begin
    var_simulation_tests("mle", 0.0, 0.0, 1.0);
end

@testset "VAR simulations: penalised MLE" begin

    models_triplets = Dict("ridge" => (10.0, 0.0, 1.0), "lasso" => (10.0, 1.0, 1.0), "en" => (10.0, 0.5, 1.0));

    for (key, value) in models_triplets
        var_simulation_tests(key, value...);
    end
end


text:
# http://stronglyconvex.com/blog/accelerated-gradient-descent.html
# TODO: Need to specify alphamax on each iteration
# Flip notation relative to Duckworth
# Start with x_{0}
# y_{t} = x_{t - 1} - alpha g(x_{t - 1})
# If converged, return y_{t}
# x_{t} = y_{t} + (t - 1.0) / (t + 2.0) * (y_{t} - y_{t - 1})

immutable AcceleratedGradientDescent{L<:Function} <: Optimizer
    linesearch!::L
end

#= uncomment for v0.8.0
AcceleratedGradientDescent(; linesearch = LineSearches.hagerzhang!) =
  AcceleratedGradientDescent(linesearch)
=#
function AcceleratedGradientDescent(; linesearch! = nothing,
                                      linesearch = LineSearches.hagerzhang!)
    linesearch = get_linesearch(linesearch!, linesearch)
    AcceleratedGradientDescent(linesearch)
end

type AcceleratedGradientDescentState{T}
    @add_generic_fields()
    x_previous::Array{T}
    g::Array{T}
    f_x_previous::T
    iteration::Int64
    y::Array{T}
    y_previous::Array{T}
    s::Array{T}
    @add_linesearch_fields()
end

function initial_state{T}(method::AcceleratedGradientDescent, options, d, initial_x::Array{T})
    g = similar(initial_x)
    f_x = d.fg!(initial_x, g)

    AcceleratedGradientDescentState("Accelerated Gradient Descent",
                         length(initial_x),
                         copy(initial_x), # Maintain current state in state.x
                         f_x, # Store current f in state.f_x
                         1, # Track f calls in state.f_calls
                         1, # Track g calls in state.g_calls
                         0, # Track h calls in state.h_calls
                         copy(initial_x), # Maintain current state in state.x_previous
                         g, # Store current gradient in state.g
                         T(NaN), # Store previous f in state.f_x_previous
                         0, # Iteration
                         copy(initial_x), # Maintain intermediary current state in state.y
                         copy(initial_x), # Maintain intermediary state in state.y_previous
                         similar(initial_x), # Maintain current search direction in state.s
                         @initial_linesearch()...) # Maintain a cache for line search results in state.lsr
end

function update_state!{T}(d, state::AcceleratedGradientDescentState{T}, method::AcceleratedGradientDescent)
    lssuccess = true
    state.iteration += 1
    # Search direction is always the negative gradient
    @simd for i in 1:state.n
        @inbounds state.s[i] = -state.g[i]
    end

    # Refresh the line search cache
    dphi0 = vecdot(state.g, state.s)
    LineSearches.clear!(state.lsr)
    push!(state.lsr, zero(T), state.f_x, dphi0)

    # Determine the distance of movement along the search line
    try
        state.alpha, f_update, g_update =
        method.linesearch!(d, state.x, state.s, state.x_ls, state.g_ls, state.lsr,
                           state.alpha, state.mayterminate)
        state.f_calls, state.g_calls = state.f_calls + f_update, state.g_calls + g_update
    catch ex
        if isa(ex, LineSearches.LineSearchException)
            lssuccess = false
            state.f_calls, state.g_calls = state.f_calls + ex.f_update, state.g_calls + ex.g_update
            state.alpha = ex.alpha
            Base.warn("Linesearch failed, using alpha = $(state.alpha) and exiting optimization.")
        else
            rethrow(ex)
        end
    end

    # Make one move in the direction of the gradient
    copy!(state.y_previous, state.y)
    @simd for i in 1:state.n
        @inbounds state.y[i] = state.x_previous[i] + state.alpha * state.s[i]
    end

    # Record previous state
    copy!(state.x_previous, state.x)

    # Update current position with Nesterov correction
    scaling = (state.iteration - 1) / (state.iteration + 2)
    @simd for i in 1:state.n
        @inbounds state.x[i] = state.y[i] + scaling * (state.y[i] - state.y_previous[i])
    end

    # Update the function value and gradient
    state.f_x_previous, state.f_x = state.f_x, d.fg!(state.x, state.g)
    state.f_calls, state.g_calls = state.f_calls + 1, state.g_calls + 1

    (lssuccess == false) # break on linesearch error
end


text: const Float = Float64

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


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

using Markdown
using InteractiveUtils

# ╔═╡ 0d696b90-14b9-11eb-19c6-7110e9eded9e
using PotentialFlow

# ╔═╡ 40e53300-14b9-11eb-28ca-c39ecf8c05ef
using Plots

# ╔═╡ 59e4bc40-14b9-11eb-3670-cdc8fa8f0ded
begin
	X = range(-1, 1, length=101)
	Y = range(-1, 1, length=101)
	Z = [x + im*y for y in Y, x in X]
end

# ╔═╡ 27ebaaf0-14be-11eb-3b42-d3c031377d61
begin

	v2 = Vortex.Point(0.0, 2)

	ψ = streamfunction(Z,v2)

	contour(X, Y, ψ,
		color = cgrad([:red, :blue]),
		grid = false,
		legend = false,
		ratio = 1,
		xlim = (-1, 1),
		ylim = (-1, 1),)

end

# ╔═╡ 2e843420-1538-11eb-14b2-d9ab38fa091b
begin

	vs = Vortex.Point.(rand(ComplexF64, 10), π*rand(10))
    fs = Freestreams.Freestream(rand(ComplexF64))

	ψs = streamfunction(Z,(vs,fs))

	contour(X, Y, ψs,
		color = cgrad([:red, :blue]),
		grid = false,
		legend = false,
		ratio = 1,
		xlim = (0, 1),
		ylim = (0, 1),)

end

# ╔═╡ e0cf3580-1538-11eb-1ba8-1baff87ba9e1
md"""
### Flow over a cylinder
"""

# ╔═╡ Cell order:
# ╠═0d696b90-14b9-11eb-19c6-7110e9eded9e
# ╠═40e53300-14b9-11eb-28ca-c39ecf8c05ef
# ╠═59e4bc40-14b9-11eb-3670-cdc8fa8f0ded
# ╠═27ebaaf0-14be-11eb-3b42-d3c031377d61
# ╠═2e843420-1538-11eb-14b2-d9ab38fa091b
# ╟─e0cf3580-1538-11eb-1ba8-1baff87ba9e1


text:
# Codec Interfaces
# ================

"""
An abstract codec type.

Any codec supporting the transcoding protocol must be a subtype of this type.

Transcoding protocol
--------------------

Transcoding proceeds by calling some functions in a specific way. We call this
"transcoding protocol" and any codec must implement it as described below.

There are six functions for a codec to implement:
- `expectedsize`: return the expected size of transcoded data
- `minoutsize`: return the minimum output size of `process`
- `initialize`: initialize the codec
- `finalize`: finalize the codec
- `startproc`: start processing with the codec
- `process`: process data with the codec.

These are defined in the `TranscodingStreams` and a new codec type must extend
these methods if necessary.  Implementing a `process` method is mandatory but
others are optional.  `expectedsize`, `minoutsize`, `initialize`, `finalize`,
and `startproc` have a default implementation.

Your codec type is denoted by `C` and its object by `codec`.

Errors that occur in these methods are supposed to be unrecoverable and the
stream will go to the panic mode. Only `Base.isopen` and `Base.close` are
available in that mode.

### `expectedsize`

The `expectedsize(codec::C, input::Memory)::Int` method takes `codec` and
`input`, and returns the expected size of transcoded data. This method will be
used as a hint to determine the size of a data buffer when `transcode` is
called. A good hint will reduce the number of buffer resizing and hence result
in better performance.

### `minoutsize`

The `minoutsize(codec::C, input::Memory)::Int` method takes `codec` and `input`,
and returns the minimum required size of the output memory when `process` is
called.  For example, an encoder of base64 will write at least four bytes to the
output and hence it is reasonable to return 4 with this method.

### `initialize`

The `initialize(codec::C)::Void` method takes `codec` and returns `nothing`.
This is called once and only once before starting any data processing.
Therefore, you may initialize `codec` (e.g. allocating memory needed to process
data) with this method. If initialization fails for some reason, it may throw an
exception and no other methods (including `finalize`) will be called. Therefore,
you need to release the memory before throwing an exception.

### `finalize`

The `finalize(codec::C)::Void` method takes `codec` and returns `nothing`.  This
is called once and only only once just before the transcoding stream goes to the
close mode (i.e. when `Base.close` is called) or just after `startproc` or
`process` throws an exception. Other errors that happen inside the stream (e.g.
`EOFError`) will not call this method. Therefore, you may finalize `codec` (e.g.
freeing memory) with this method. If finalization fails for some reason, it may
throw an exception. You should release the allocated memory in codec before
returning or throwing an exception in `finalize` because otherwise nobody cannot
release the memory. Even when an exception is thrown while finalizing a stream,
the stream will become the close mode for safety.

### `startproc`

The `startproc(codec::C, mode::Symbol, error::Error)::Symbol` method takes
`codec`, `mode` and `error`, and returns a status code. This is called just
before the stream starts reading or writing data. `mode` is either `:read` or
`:write` and then the stream starts reading or writing, respectively.  The
return code must be `:ok` if `codec` is ready to read or write data.  Otherwise,
it must be `:error` and the `error` argument must be set to an exception object.

### `process`

The `process(codec::C, input::Memory, output::Memory,
error::Error)::Tuple{Int,Int,Symbol}` method takes `codec`, `input`, `output`
and `error`, and returns a consumed data size, a produced data size and a status
code. This is called repeatedly while processing data. The input (`input`) and
output (`output`) data are a `Memory` object, which is a pointer to a contiguous
memory region with size. You must read input data from `input`, transcode the
bytes, and then write the output data to `output`.  Finally you need to return
the size of read data, the size of written data, and `:ok` status code so that
the caller can know how many bytes are consumed and produced in the method.
When transcoding reaches the end of a data stream, it is notified to this method
by empty input. In that case, the method need to write the buffered data (if
any) to `output`. If there is no data to write, the status code must be set to
`:end`. The `process` method will be called repeatedly until it returns `:end`
status code. If an error happens while processing data, the `error` argument
must be set to an exception object and the return code must be `:error`.
"""
abstract type Codec end

# Methods
# -------

"""
    expectedsize(codec::Codec, input::Memory)::Int

Return the expected size of the transcoded `input` with `codec`.

The default method returns `input.size`.
"""
function expectedsize(codec::Codec, input::Memory)::Int
    return input.size
end

"""
    minoutsize(codec::Codec, input::Memory)::Int

Return the minimum output size to be ensured when calling `process`.

The default method returns `max(1, div(input.size, 4))`.
"""
function minoutsize(codec::Codec, input::Memory)::Int
    return max(1, div(input.size, 4))
end

"""
    initialize(codec::Codec)::Void

Initialize `codec`.

The default method does nothing.
"""
function initialize(codec::Codec)
    return nothing
end

"""
    finalize(codec::Codec)::Void

Finalize `codec`.

The default method does nothing.
"""
function finalize(codec::Codec)::Nothing
    return nothing
end

"""
    startproc(codec::Codec, mode::Symbol, error::Error)::Symbol

Start data processing with `codec` of `mode`.

The default method does nothing and returns `:ok`.
"""
function startproc(codec::Codec, mode::Symbol, error::Error)::Symbol
    return :ok
end

"""
    process(codec::Codec, input::Memory, output::Memory, error::Error)::Tuple{Int,Int,Symbol}

Do data processing with `codec`.

There is no default method.
"""
function process(codec::Codec, input::Memory, output::Memory, error::Error)::Tuple{Int,Int,Symbol}
    # no default method
    throw(MethodError(process, (codec, input, output, error)))
end


text:
# Файл roblib.jl

"""
    moves!(r::Robot, side::HorizonSide)

-- перемещает Робота до упора в заданном направлении
"""
moves!(r::Robot, side::HorizonSide) =
while isborder(r,side)
    move!(r,side)
end

"""
    moves!(r::Robot, side::HorizonSide, num_steps::Int)

-- перемещает Робота в заданном направлении на заданное число шагов (если на пути - перегородка, то - ошибка)
"""
moves!(r::Robot, side::HorizonSide, num_steps::Int) =
for _ in 1:num_steps
    move!(r,side)
end

"""
    find_border!(r::Robot,side_to_border::HorizonSide, side_of_movement::HorizonSide)

-- останавливает робота у перегородки, которая ожидается с направления side_to_border, при движении робота "змейкой" в сторону перегородки (от упора до упора в поперечном этому напавлении).

-- side_of_movement - начальное "поперечное" направление
"""
find_border!(r::Robot,side_to_border::HorizonSide, side_of_movement::HorizonSide) =
while isborder(r,side_to_border)==false
    if isborder(r,side_of_movement)==false
        move!(r,side_of_movement)
    else
        move!(r,side_to_border)
        side_of_movement=inverse(side_of_movement)
    end
end

"""
    inverse(side::HorizonSide)

-- возвращает сторону горизонта, противоположную заданной
"""
inverse(side::HorizonSide) = HorizonSide(mod(Int(side)+2, 4))

"""
    putmarkers!(r::Robot, side::HorizonSide)

-- ставит маркеры, пермещая Робота до упора в заданном направлении (в начальной клетке маркер не ставится)
"""
putmarkers!(r::Robot, side::HorizonSide) =
while isborder(r,side)==false
    move!(r,side)
    putmarker!(r)
end

"""
    putmarkers!(r::Robot,side_of_movement::HorizonSide,side_to_border::HorizonSide)

-- Ставит маркеры и перемещает Робота в направлении side_of_movement пока рядом с ним в направлении side_to_border имеется перегородка (эти два направления должны быть взаимно перпендикулярными)
"""
putmarkers!(r::Robot,side_of_movement::HorizonSide,side_to_border::HorizonSide) =
while isborder(r,side_to_border)==true
    move!(r,side_of_movement)
end

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

# Файл с наиболее часто использумыми функциями

# Используем более привычные английские названия сторон света вместо немецких
"""
Север
"""
North = Nord
"""
Юг
"""
South = Sud
#Запад - уже определён как West
"""
Восток
"""
East = Ost

"""
Возвращает инвертированное направление
"""
function inverse_side(side::HorizonSide)::HorizonSide
    return HorizonSide(mod(Int(side)+2, 4))
end

"""
Возвращает следующее по часовой стрелке направление
"""
function clockwise_side(side::HorizonSide)::HorizonSide
    return HorizonSide(mod(Int(side)-1,4))
end

"""
Возвращает следующее против часовой стрелки направление
"""
function counterclockwise_side(side::HorizonSide)::HorizonSide
    return HorizonSide(mod(Int(side)+1,4))
end

"""
go!(r::Robot, side::HorizonSide; steps::Int = 1, go_around_barriers::Bool = false, markers = false)::Int

Перемещает робота в направлении. Возвращает 0, если невозможно пройти в данном направлении, отрицательное число (Если робот обходил препятствия, но в нужной клетке был прямоугольник. Тогда это отрицательное число, равное числу шагов, сделанных роботом, возвращающимся в исходную точку),либо положительное число (Количество реально сделанных шагов в данном направлении).

Именованные параметры:
steps - количество шагов в направлении;

go_around_barriers - робот будет по пути обходить перегородки, если такие имеются;

markers - робот будет ставить маркеры по пути (не ставит их при обходе перегородок).
"""
function go!(r::Robot, side::HorizonSide; steps::Int = 1, go_around_barriers::Bool = false, markers = false)::Int
    my_ans = 0
    if markers
        putmarker!(r)
    end
    if (go_around_barriers)
        path = around_move_return_path!(r, side; steps, markers)
        my_ans = get_path_length_in_direction(path, side)
    else
        for i ∈ 1:steps

            if (markers)
                putmarker!(r)
            end

            if !isborder(r, side)
                move!(r, side)
                my_ans += 1
            else
                for i ∈ 1:my_ans
                    move!(r, inverse_side(side))
                end
                my_ans = 0
                break
            end
        end
        if (markers)
            putmarker!(r)
        end
    end

    return my_ans
end

"""
around_move_return_path!(r::Robot, side::HorizonSide; steps::Int = 1, markers = false)::Array{Tuple{HorizonSide,Int64},1}

Вспомогательная функция. Возвращает путь после перемещения робота
"""
function around_move_return_path!(r::Robot, side::HorizonSide; steps::Int = 1, markers = false)::Array{Tuple{HorizonSide,Int64},1}
    path = [ (Nord, 0) ]
    steps_to_do = steps

    while steps_to_do > 0

        if markers
            putmarker!(r)
        end

        path_now = go_around_barrier_and_return_path!(r, side)

        for i in path_now
            push!(path, i)
        end

        steps_to_do -= get_path_length_in_direction(path_now, side)

        if !isborder(r, side) && steps_to_do > 0
            push!(path, ( inverse_side(side), 1))
            move!(r, side)
            steps_to_do -= 1

            if markers
                putmarker!(r)
            end
        elseif get_path_length_in_direction(path_now, side) == 0
            steps_to_do = -1
            break
        end
        if markers && steps_to_do >= 0
            putmarker!(r)
        end
    end
    if steps_to_do < 0
        my_ans = 0
        # возвращаемся назад
        go_by_path!(r, path)
        path = [ (North, 0) ]
    end

    return path
end

"""
go_to_border_come_back_and_return_distance!(r::Robot, side::HorizonSide; go_around_barriers::Bool = false, markers = false)::Int

Перемещает робота до границы и обратно, возвращая расстояние до границы (ортогональные шаги при обходе перегородок не учитываются)

Именованные параметры:
go_around_barriers - робот будет по пути обходить перегородки, если такие имеются;

markers - робот будет ставить маркеры по пути (не ставит их при обходе перегородок).
"""
function go_to_border_come_back_and_return_distance!(r::Robot, side::HorizonSide; go_around_barriers::Bool = false, markers = false)::Int
    my_ans = 0
    if go_around_barriers
        # здесь и далее интерпритатор понимает, что я хочу записать в именованный параметр markers вызываемой функции значение из параметра markers текущей функции
        if markers
            putmarker!(r)
        end
        if !isborder(r, side)
            move!(r, side)
            steps = 1
        else
            steps = get_path_length_in_direction(go_around_barrier_and_return_path!(r, side), side)
        end
        my_ans += steps
        if markers
            putmarker!(r)
        end
        while steps > 0
            if !isborder(r, side)
                move!(r, side)
                steps = 1
            else
                steps = get_path_length_in_direction(go_around_barrier_and_return_path!(r, side), side)
            end
            my_ans += steps
        end
        if markers
            putmarker!(r)
        end
        go!(r, inverse_side(side); steps = my_ans, go_around_barriers = true)
    else
        while go!(r, side; markers) > 0
            my_ans += 1
        end
        go!(r, inverse_side(side); steps = my_ans)
    end
    return my_ans
end

"""
go_to_border_and_return_path!(r::Robot, side::HorizonSide; go_around_barriers::Bool = false, markers = false)::Array{Tuple{HorizonSide,Int64},1}

Перемещает робота до границы и возвращает путь для обраного следования в виде массива пар типа (направление, количество шагов)

Именованные параметры:
go_around_barriers - робот будет по пути обходить перегородки, если такие имеются;

markers - робот будет ставить маркеры по пути (не ставит их при обходе перегородок).
"""
function go_to_border_and_return_path!(r::Robot, side::HorizonSide; go_around_barriers::Bool = false, markers = false)::Array{Tuple{HorizonSide,Int64},1}
    my_ans = [ (North, 0) ]
    if go_around_barriers
        steps = 0
        if markers
            putmarker!(r)
        end
        if !isborder(r, side)
            move!(r, side)
            steps = 1
            push!(my_ans, (inverse_side(side), 1) )
        else
            path = go_around_barrier_and_return_path!(r, side)
            steps = get_path_length_in_direction(path, side)
            for i in path
                push!(my_ans, i)
            end
        end
        if markers
            putmarker!(r)
        end
        while steps > 0
            if !isborder(r, side)
                move!(r, side)
                steps = 1
                push!(my_ans, (inverse_side(side), 1) )
                if markers
                    putmarker!(r)
                end
            else
                path = go_around_barrier_and_return_path!(r, side)
                steps = get_path_length_in_direction(path, side)
                for i in path
                    push!(my_ans, i)
                end
                if markers
                    putmarker!(r)
                end
            end
        end

    else
        steps=0
        steps_now = go!(r,side; markers)
        while steps_now > 0
            steps += steps_now
            steps_now = go!(r,side; markers)
        end
        push!(my_ans, (inverse_side(side), steps) )
    end
    return my_ans
end

"""
go_to_west_south_corner!(r::Robot; go_around_barriers::Bool = false, markers = false)::Array{Tuple{HorizonSide,Int64},1}

Перемещает робота в левый нижний угол и возвращает путь для обраного следования в виде массива пар типа (направление, количество шагов)

Именованные параметры:
go_around_barriers - робот будет по пути обходить перегородки, если такие имеются;

markers - робот будет ставить маркеры по пути (не ставит их при обходе перегородок).
"""
function go_to_west_south_corner_and_return_path!(r::Robot; go_around_barriers::Bool = false, markers = false)::Array{Tuple{HorizonSide,Int64},1}
    my_ans = []
    a = go_to_border_and_return_path!(r, West; go_around_barriers, markers)
    b = go_to_border_and_return_path!(r, South; go_around_barriers, markers)

    for i in a
        push!(my_ans, i)
    end
    for i in b
        push!(my_ans, i)
    end
    return my_ans
end

function go_to_some_corner_and_return_path!(r::Robot, side_1::HorizonSide, side_2::HorizonSide; go_around_barriers::Bool = false, markers = false)::Array{Tuple{HorizonSide,Int64},1}
    my_ans = []
    a = go_to_border_and_return_path!(r, side_1; go_around_barriers, markers)
    b = go_to_border_and_return_path!(r, side_2; go_around_barriers, markers)

    for i in a
        push!(my_ans, i)
    end
    for i in b
        push!(my_ans, i)
    end
    return my_ans
end

"""
go_by_path!(r::Robot, path::Array{Tuple{HorizonSide,Int64},1})

Перемещает робота по пути. Путь должен быть в формате массива пар (направление, количество шагов).
Путь на корректность не проверяется, что может привести к ошибкам в ходе выполнения.
"""
function go_by_path!(r::Robot, path::Array{Tuple{HorizonSide,Int64},1})
    new_path = reverse(path)
    for i in new_path
        go!(r, i[1]; steps = i[2])
    end
end

"""
go_around_barrier_and_return_path!(r::Robot, direct_side::HorizonSide)::Array{Tuple{HorizonSide,Int64},1}

Обходит барьер, если таковой имеется перед роботом. Возвращает путь для возвращения в формате массива пар (направление, количество шагов).
"""
function go_around_barrier_and_return_path!(r::Robot, direct_side::HorizonSide)::Array{Tuple{HorizonSide,Int64},1}
    my_ans = []
    orthogonal_side = clockwise_side(direct_side)
    reverse_side = inverse_side(orthogonal_side)
    num_of_orthohonal_steps = 0
    num_of_direct_steps = 0

    if !isborder(r, direct_side)
        my_ans = [ (North, 0) ]
    else
        while isborder(r,direct_side) == true
            if isborder(r, orthogonal_side) == false
                move!(r, orthogonal_side)
                num_of_orthohonal_steps += 1
            else
                break
            end
        end

        if isborder(r,direct_side) == false
            move!(r,direct_side)
            num_of_direct_steps += 1
            while isborder(r,reverse_side) == true
                num_of_direct_steps += 1
                move!(r,direct_side)
            end
            push!(my_ans, (inverse_side(orthogonal_side), num_of_orthohonal_steps) )
            push!(my_ans, (inverse_side(direct_side), num_of_direct_steps) )
            push!(my_ans, (inverse_side(reverse_side), num_of_orthohonal_steps) )
        else
            my_ans = [ (North, 0) ]
        end

        while num_of_orthohonal_steps > 0
            num_of_orthohonal_steps=num_of_orthohonal_steps-1
            move!(r,reverse_side)
        end

    end
    return my_ans
end

"""
get_path_length(path::Array{Tuple{HorizonSide,Int64},1})::Int

Возвращает общую количество шагов в массиве. Массив должен быть в формате пар (направление, количество шагов).
"""
function get_path_length(path::Array{Tuple{HorizonSide,Int64},1})::Int
    my_ans = 0
    for i in path
        my_ans += i[2]
    end
    return my_ans
end

"""
get_path_length_in_direction(path::Array{Tuple{HorizonSide,Int64},1}, direction::HorizonSide)::Int

Возвращает количество шагов в массиве в данном направлении, или в противоположном ему. Массив должен быть в формате пар (направление, количество шагов).
"""
function get_path_length_in_direction(path::Array{Tuple{HorizonSide,Int64},1}, direction::HorizonSide)::Int
    my_ans = 0
    for i in path
        (i[1] == direction || i[1] == inverse_side(direction)) ? my_ans += i[2] : Nothing
    end
    return my_ans
end

"""
move_if_possible!(r::Robot, direct_side::HorizonSide)::Bool

Робот двигается в направлении, если это возможно. Если препятствия нет, робот пройдёт 1 клетку. Иначе, обойдёт препятствие
"""
function move_if_possible!(r::Robot, direct_side::HorizonSide)::Bool
    orthogonal_side = counterclockwise_side(direct_side)
    reverse_side = inverse_side(orthogonal_side)
    num_steps=0
    if isborder(r,direct_side)==false
        move!(r,direct_side)
        result=true
    else
        while isborder(r,direct_side) == true
            if isborder(r, orthogonal_side) == false
                move!(r, orthogonal_side)
                num_steps += 1
            else
                break
            end
        end
        if isborder(r,direct_side) == false
            move!(r,direct_side)
            while isborder(r,reverse_side) == true
                move!(r,direct_side)
            end
            result = true
        else
            result = false
        end
        while num_steps>0
            num_steps=num_steps-1
            move!(r,reverse_side)
        end
    end
    return result
end

max_stars_repo_name: Savko-fokin/mirea-progs

max_issues_repo_name: Savko-fokin/mirea-progs

max_forks_repo_name: Savko-fokin/mirea-progs


text:
"""
Clamped cylinder, free vibration.

The reference gives the frequencies and mode shape indicators in the table
![](Table-freq-1.gif). Unfortunately, no indication of the accuracy of these
numbers is provided: we don't know how many elements were used, or what type
of elements.

The mode shapes are described with M = number of circumferential waves, N =
number of longitudinal half waves, using the tuple (M, N).

For instance, for the thickness = 0.01, the fundamental frequency should be
approximately 1701 Hz (2, 1), mode 13 frequency should be approximately 3926
Hz (torsion). The torsional mode does not occur for thinner shells.

Reference
http://www.adina.com/newsgH53.shtml

"""
module clamp_cyl_vibration_examples

using Arpack
using LinearAlgebra
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json
using FinEtools.MeshExportModule.VTKWrite: vtkwrite

E = 200e9;
nu = 0.3;
rho = 7800.0
R = 0.1
L = 0.4

cylindrical!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) = begin
    r = vec(XYZ); r[2] = 0.0;
    csmatout[:, 3] .= vec(r)/norm(vec(r))
    csmatout[:, 2] .= (0.0, 1.0, 0.0) #  this is along the axis
    cross3!(view(csmatout, :, 1), view(csmatout, :, 2), view(csmatout, :, 3))
    return csmatout
end

function _execute(n = 2, thickness = 0.01, visualize = true)
    formul = FEMMShellT3FFModule

    n = Int(round(n/2)*2)
    @info "Mesh: $n elements per side"

    # Mesh
    tolerance = min(R, L) / n  / 100

    tolerance = R/n/1000
    fens, fes = T3block(360.0, L, n, n);
    fens.xyz = xyz3(fens)
    for i in 1:count(fens)
        a=fens.xyz[i, 1]; y=fens.xyz[i, 2];
        fens.xyz[i, :] .= (R*sin(a/180*pi), y-L/2, R*cos(a/180*pi))
    end
    fens, fes = mergenodes(fens, fes, tolerance)
    bfes = meshboundary(fes)

    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)
    ocsys = CSys(3, 3, cylindrical!)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), ocsys, mater)

    associategeometry! = formul.associategeometry!
    stiffness = formul.stiffness
    mass = formul.mass

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Clamp the cross sections
    nl = connectednodes(bfes)
    for i in 1:6
        setebc!(dchi, nl, true, i)
    end

    numberdofs!(dchi);

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);
    M = mass(femm, geom0, dchi);

    # Solve
    OmegaShift = 0.0*2*pi
    neigvs = 20
    evals, evecs, nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM, explicittransform=:none)
    evals[:] = evals .- OmegaShift;
    fs = real(sqrt.(complex(evals)))/(2*pi)
    @show fs

    for ev in 1:length(fs)
        U = evecs[:, ev]
        scattersysvec!(dchi, 1.0/maximum(abs.(U)).*U)
        vtkwrite("clamped_cylinder_vibration-$thickness-mode-$(ev).vtu", fens, fes; vectors = [("u", dchi.values[:, 1:3]), ("ur", dchi.values[:, 4:6])])
    end

    return true
end

function test_convergence()
    @info "Clamped cylinder vibration"
    for n in [64]
        _execute(n, 0.01, !false)
    end
    return true
end

function allrun()
    println("#####################################################")
    println("# test_convergence ")
    test_convergence()
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing


text:
# # Comparison to Neumann analytic solution

# Before reading this tutorial,
# we recommend that you look over the coupled energy
# and water [tutorial](../Coupled/equilibrium_test.md)
# and the freezing front [tutorial](freezing_front.md).
# The former shows how to solve the heat equation for soil volumetric
# internal energy `ρe_int` simultaneously
# with Richards equation for volumetric liquid water fraction `ϑ_l`, assuming zero
# volumetric ice fraction `θ_i` for all time, everywhere in the domain[^a].
# The latter shows how to include freezing and thawing, and explains the freeze thaw model employed
# by CliMA Land. This tutorial compares a simulated temperature profile
# from a freezing front with an analytic solution.

# The analytic solution applies to a freezing front propagating under certain assumptions. It assumes
# that there is no water movement, which we mimic by setting `K_sat=0`. It also assumes a semi-infinite
# domain, which we approximate by making the domain larger than the extent to which the front propagates.
# The solution also assumes that all water freezes. Our model does not satisfy that assumption, as discussed
# [here](freezing_front.md#Discussion-and-Model-Explanation), but as we will see, the model still matches the
# analytic expectation well in the frozen region.

# As such, our set of equations is

# ``
# \frac{∂ ρe_{int}}{∂ t} =  ∇ ⋅ κ(θ_l, θ_i; ν, ...) ∇T
# ``

# ``
# \frac{ ∂ ϑ_l}{∂ t} = -\frac{F_T}{ρ_l}
# ``

# ``
# \frac{ ∂ θ_i}{∂ t} = \frac{F_T}{ρ_i}
# ``

# Here

# ``t`` is the time (s),

# ``z`` is the location in the vertical (m),

# ``ρe_{int}`` is the volumetric internal energy of the soil (J/m^3),

# ``T`` is the temperature of the soil (K),

# ``κ`` is the thermal conductivity (W/m/K),

# ``ϑ_l`` is the augmented volumetric liquid water fraction,

# ``θ_i`` is the volumetric ice fraction,

# ``ν, ...`` denotes parameters relating to soil type, such as porosity, and

# ``F_T`` is the freeze-thaw term.

# Our domain is effectively 1-d, with ``z ∈ [-3,0]``,
# and with the following boundary and initial conditions:

# `` T(t, z=0) = 263.15 K``

# ``- κ ∇T(t, z= -3) =  0 ẑ ``

# `` T(t = 0, z) = 275.15 K``

# ``- K ∇h(t, z = 0) = 0 ẑ ``

# `` -K ∇h(t, z = -3) = 0 ẑ``

# `` ϑ(t = 0, z) = 0.33 ``.

# # Import necessary modules
# External (non - CliMA) modules
using MPI
using OrderedCollections
using StaticArrays
using Statistics
using Test
using DelimitedFiles
using Plots

# CliMA Parameters
using CLIMAParameters
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()
using CLIMAParameters.Planet: ρ_cloud_liq
using CLIMAParameters.Planet: ρ_cloud_ice
using CLIMAParameters.Planet: LH_f0

# ClimateMachine modules
using ClimateMachine
using ClimateMachine.Land
using ClimateMachine.Land.SoilWaterParameterizations
using ClimateMachine.Land.SoilHeatParameterizations
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.DGMethods: BalanceLaw, LocalGeometry
using ClimateMachine.MPIStateArrays
using ClimateMachine.GenericCallbacks
using ClimateMachine.SystemSolvers
using ClimateMachine.ODESolvers
using ClimateMachine.VariableTemplates
using ClimateMachine.SingleStackUtils
using ClimateMachine.BalanceLaws:
    BalanceLaw, Prognostic, Auxiliary, Gradient, GradientFlux, vars_state
using SpecialFunctions
using ArtifactWrappers

# # Preliminary set-up

# Get the parameter set, which holds constants used across CliMA models:
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet();
# Initialize and pick a floating point precision:
ClimateMachine.init()
const FT = Float64;

# # Simulation specific parameters
N_poly = 1
nelem_vert = 40
zmax = FT(0)
zmin = FT(-3)
t0 = FT(0)
dt = FT(50)
timeend = FT(3600 * 24 * 20)
n_outputs = 540
every_x_simulation_time = ceil(Int, timeend / n_outputs)
Δ = abs(zmin - zmax) / FT(nelem_vert);

# # Soil properties
# All units are mks.
porosity = FT(0.535)
vg_α = 1.11
vg_n = 1.48
Ksat = 0.0
ν_ss_quartz = FT(0.2)
ν_ss_minerals = FT(0.6)
ν_ss_om = FT(0.2)
ν_ss_gravel = FT(0.0);
κ_quartz = FT(7.7)
κ_minerals = FT(2.5)
κ_om = FT(0.25)
κ_liq = FT(0.57)
κ_ice = FT(2.29)
ρp = FT(2700)
κ_solid = k_solid(ν_ss_om, ν_ss_quartz, κ_quartz, κ_minerals, κ_om)
κ_sat_frozen = ksat_frozen(κ_solid, porosity, κ_ice)
κ_sat_unfrozen = ksat_unfrozen(κ_solid, porosity, κ_liq)
ρc_ds = FT((1 - porosity) * 2.3e6)
soil_param_functions = SoilParamFunctions(
    FT;
    porosity = porosity,
    ν_ss_gravel = ν_ss_gravel,
    ν_ss_om = ν_ss_om,
    ν_ss_quartz = ν_ss_quartz,
    ρc_ds = ρc_ds,
    ρp = ρp,
    κ_solid = κ_solid,
    κ_sat_unfrozen = κ_sat_unfrozen,
    κ_sat_frozen = κ_sat_frozen,
    water = WaterParamFunctions(FT; Ksat = Ksat, S_s = 1e-3),
);

# # Build the model
# Initial and Boundary conditions. The default initial condition for
# `θ_i` is zero everywhere, so we don't modify that. Furthermore, since
# the equation for `θ_i` does not involve spatial derivatives, we don't need
# to supply boundary conditions for it.
zero_flux = (aux, t) -> eltype(aux)(0.0)
ϑ_l0 = (aux) -> eltype(aux)(0.33)
surface_state = (aux, t) -> eltype(aux)(273.15 - 10.0)
T_init = aux -> eltype(aux)(275.15)
bc = LandDomainBC(
    bottom_bc = LandComponentBC(
        soil_heat = Neumann(zero_flux),
        soil_water = Neumann(zero_flux),
    ),
    surface_bc = LandComponentBC(
        soil_heat = Dirichlet(surface_state),
        soil_water = Neumann(zero_flux),
    ),
);

# Create the [`SoilWaterModel`](@ref ClimateMachine.Land.SoilWaterModel),
# [`SoilHeatModel`](@ref ClimateMachine.Land.SoilHeatModel),
# and the [`SoilModel`](@ref ClimateMachine.Land.SoilModel) instances.
# Note that we are still specifying a hydraulics model, because the matric potential
# and hydraulic conductivity functions are still evaluated (though they don't affect
# the outcome). Setting `Ksat =0` is just a
# hack for turning off water flow.
soil_water_model = SoilWaterModel(
    FT;
    hydraulics = vanGenuchten(FT; α = vg_α, n = vg_n),
    initialϑ_l = ϑ_l0,
);

soil_heat_model = SoilHeatModel(FT; initialT = T_init);

m_soil = SoilModel(soil_param_functions, soil_water_model, soil_heat_model);
# Create the source term instance. Our phase change model requires
# knowledge of the vertical spacing, so we pass
# that information in via an attribute of the
# [`PhaseChange`](@ref ClimateMachine.Land.PhaseChange) structure.
freeze_thaw_source = PhaseChange{FT}(Δz = Δ);

# Sources are added as elements of a list of sources. Here we just add freezing
# and thawing.
sources = (freeze_thaw_source,);

# Next, we define the required `init_soil!` function, which takes the user
# specified functions of space for `T_init` and `ϑ_l0` and initializes the state
# variables of volumetric internal energy and augmented liquid fraction. This requires
# a conversion from `T` to `ρe_int`.
function init_soil!(land, state, aux, localgeo, time)
    myFT = eltype(state)
    ϑ_l = myFT(land.soil.water.initialϑ_l(aux))
    θ_i = myFT(land.soil.water.initialθ_i(aux))
    state.soil.water.ϑ_l = ϑ_l
    state.soil.water.θ_i = θ_i
    param_set = land.param_set

    θ_l = volumetric_liquid_fraction(ϑ_l, land.soil.param_functions.porosity)
    ρc_ds = land.soil.param_functions.ρc_ds
    ρc_s = volumetric_heat_capacity(θ_l, θ_i, ρc_ds, param_set)

    state.soil.heat.ρe_int = volumetric_internal_energy(
        θ_i,
        ρc_s,
        land.soil.heat.initialT(aux),
        param_set,
    )
end;

# Lastly, package it all up in the `LandModel`:
m = LandModel(
    param_set,
    m_soil;
    boundary_conditions = bc,
    source = sources,
    init_state_prognostic = init_soil!,
);

# # Set up and run the simulation

driver_config = ClimateMachine.SingleStackConfiguration(
    "LandModel",
    N_poly,
    nelem_vert,
    zmax,
    param_set,
    m;
    zmin = zmin,
    numerical_flux_first_order = CentralNumericalFluxFirstOrder(),
);

solver_config =
    ClimateMachine.SolverConfiguration(t0, timeend, driver_config, ode_dt = dt);

state_types = (Prognostic(), Auxiliary(), GradientFlux())
all_data = Dict[dict_of_nodal_states(solver_config, state_types; interp = true)]
time_data = FT[0]
callback = GenericCallbacks.EveryXSimulationTime(every_x_simulation_time) do
    dons = dict_of_nodal_states(solver_config, state_types; interp = true)
    push!(all_data, dons)
    push!(time_data, gettime(solver_config.solver))
    nothing
end;

ClimateMachine.invoke!(solver_config; user_callbacks = (callback,));
z = get_z(solver_config.dg.grid; rm_dupes = true);

# # Analytic Solution of Neumann
# All details here are taken from [DallAmico2011](@cite) (see also [CarslawJaeger](@cite)), and the reader is referred to that
# for further information on the solution. It takes the form of a function
# for T(z) on each side of the freezing front interface, which depends on
# the thermal properties in that region, and which is also parameterized by
# a parameter (ζ), which we show how to solve for below. In computing the
# thermal properties, we evaluate the conductivity and heat capacity
# assuming that all of the water is either in liquid or frozen form, with total
# mass proportional to ``θ_{l,0}ρ_l`` (as we have no water flow).

# Compute the thermal conductivity and heat capacity in the frozen region - subscript 1.
θ_l0 = FT(0.33)
kdry = k_dry(param_set, soil_param_functions)
ksat = saturated_thermal_conductivity(
    FT(0.0),
    θ_l0 * ρ_cloud_liq(param_set) / ρ_cloud_ice(param_set),
    κ_sat_unfrozen,
    κ_sat_frozen,
)
kersten = kersten_number(
    θ_l0 * ρ_cloud_liq(param_set) / ρ_cloud_ice(param_set),
    θ_l0 * ρ_cloud_liq(param_set) / ρ_cloud_ice(param_set) / porosity,
    soil_param_functions,
)
λ1 = thermal_conductivity(kdry, kersten, ksat)
c1 = volumetric_heat_capacity(
    FT(0.0),
    θ_l0 * ρ_cloud_liq(param_set) / ρ_cloud_ice(param_set),
    ρc_ds,
    param_set,
)
d1 = λ1 / c1;

# Compute the thermal conductivity and heat capacity in the region
# with liquid water - subscript 2.
ksat =
    saturated_thermal_conductivity(θ_l0, FT(0.0), κ_sat_unfrozen, κ_sat_frozen)
kersten = kersten_number(FT(0.0), θ_l0 / porosity, soil_param_functions)
λ2 = thermal_conductivity(kdry, kersten, ksat)
c2 = volumetric_heat_capacity(θ_l0, FT(0.0), ρc_ds, param_set)
d2 = λ2 / c2;

# Initial T and surface T, in Celsius
Ti = FT(2)
Ts = FT(-10.0);

# The solution requires the root of the implicit equation below
# (you'll need to install `Roots` via the package manager).
# ``` julia
# using Roots
# function implicit(ζ)
#    term1 = exp(-ζ^2) / ζ / erf(ζ)
#    term2 =
#        -λ2 * sqrt(d1) * (Ti - 0) /
#        (λ1 * sqrt(d2) * (0 - Ts) * ζ * erfc(ζ * sqrt(d1 / d2))) *
#        exp(-d1 / d2 * ζ^2)
#    term3 =
#        -LH_f0(param_set) * ρ_cloud_liq(param_set) * θ_l0 * sqrt(π) / c1 /
#        (0 - Ts)
#    return (term1 + term2 + term3)
# end
# find_zero(implicit, (0.25,0.27), Bisection())
# ```

# The root is
ζ = 0.26447353269809687;

# This function plots the analytic solution
# and the simulated result, at the output time index
# `k`.
function f(k; ζ = ζ)
    T = all_data[k]["soil.heat.T"][:] .- 273.15
    plot(
        T,
        z[:],
        xlim = [-10.5, 3],
        ylim = [zmin, zmax],
        xlabel = "T",
        label = "simulation",
    )
    t = time_data[k]
    zf = 2.0 * ζ * sqrt(d1 * t)
    myz = zmin:0.001:0
    spatially_varying =
        (erfc.(abs.(myz) ./ (zf / ζ / (d1 / d2)^0.5))) ./
        erfc(ζ * (d1 / d2)^0.5)
    mask = abs.(myz) .>= zf
    plot!(
        Ti .- (Ti - 0.0) .* spatially_varying[mask],
        myz[mask],
        label = "analytic",
        color = "green",
    )
    spatially_varying = ((erf.(abs.(myz) ./ (zf / ζ)))) ./ erf(ζ)
    mask = abs.(myz) .< zf
    plot!(
        Ts .+ (0.0 - Ts) .* spatially_varying[mask],
        myz[mask],
        label = "",
        color = "green",
    )
end

# Now we will plot this and compare to other methods for modeling phase change without water movement.
# These solutions were produced by modifying the Supplemental Program 5:3 from [Bonan19a](@cite).

# Excess heat solution:
eh_dataset = ArtifactWrapper(
    @__DIR__,
    isempty(get(ENV, "CI", "")),
    "eh",
    ArtifactFile[ArtifactFile(
        url = "https://caltech.box.com/shared/static/6xs1r98wk7u1b0xjhpkdvt80q4sh16un.csv",
        filename = "bonan_data.csv",
    ),],
);
eh_dataset_path = get_data_folder(eh_dataset);
eh_data = joinpath(eh_dataset_path, "bonan_data.csv")
ds = readdlm(eh_data, ',');
# Apparent heat capacity solution:
ahc_dataset = ArtifactWrapper(
    @__DIR__,
    isempty(get(ENV, "CI", "")),
    "ahc",
    ArtifactFile[ArtifactFile(
        url = "https://caltech.box.com/shared/static/d6xciskzl2djwi03xxfo8wu4obrptjmy.csv",
        filename = "bonan_data_ahc.csv",
    ),],
);
ahc_dataset_path = get_data_folder(ahc_dataset);
ahc_data = joinpath(ahc_dataset_path, "bonan_data_ahc.csv")
ds2 = readdlm(ahc_data, ',');

k = n_outputs
f(k)
plot!(ds[:, 2], ds[:, 1], label = "Excess Heat")
plot!(ds2[:, 2], ds2[:, 1], label = "Apparent heat capacity")
plot!(legend = :bottomleft)
savefig("analytic_comparison.png")
# ![](analytic_comparison.png)

# # References
# - [DallAmico2011](@cite)
# - [CarslawJaeger](@cite)
# - [Bonan19a](@cite)

# [^a]:
#       Note that `θ_i` is always treated as a prognostic variable
#       in the `SoilWaterModel`, but with
#       zero terms on the RHS unless freezing and thawing is turn on, as demonstrated in this
#       tutorial. That means that the user could, in principle, set the initial condition to be nonzero
#       (`θ_i(x, y, z ,t=0) = 0` is the default), which in turn would allow a nonzero `θ_i`
#       profile to affect things like thermal conductivity, etc,
#       in a consistent way. However, it would not be enforced that ``θ_l+θ_i \leq ν``, because there would
#       be no physics linking the liquid and water content to each other, and they are independent
#       variables in our model. We don't envision this being a common use case.


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

using Markdown
using InteractiveUtils

# ╔═╡ bd3a666c-51c4-11eb-0ea5-3beaf9b624c2
begin
  import Pkg
  # Adding packages is required for deployment on Binder
  Pkg.add("CSV")
  Pkg.add("DataFrames")
  Pkg.add("PlutoUI")
  Pkg.add("StatsPlots")
  Pkg.add("RemoteFiles")
  using PlutoUI, CSV, DataFrames, StatsPlots, RemoteFiles
  plotly()
end

# ╔═╡ 3a965f9c-510e-11eb-08dc-8d1cdc52e218
with_terminal() do
	Pkg.status()
end

# ╔═╡ 30e459ae-526c-11eb-19b3-8f67825515c2
begin
	@RemoteFile data_csv "https://raw.githubusercontent.com/StatisticalMice/superforecasters/main/data2020_b.csv"
	download(data_csv, verbose=true)
	csv_path = path(data_csv)
end

# ╔═╡ 0cdb58ae-510f-11eb-12d4-5393fb3dec0c
df_raw = CSV.read(csv_path, DataFrame)[!,3:end];

# ╔═╡ 3493c430-51b9-11eb-2b59-6bd9bde30ed5

# ╔═╡ 9ebdd936-5187-11eb-2ed8-3fa4ec45d80c

# ╔═╡ c9f9c422-5185-11eb-3e8d-3f96592baac1
function split_column_name(s)
    rxa = r"^\d+\.\s*"
    n = match(rxa, s).match
    text = s[length(n)+1:end]
    n, text
end

# ╔═╡ ef460ef8-5186-11eb-360f-cfb255b92d0f
split_column_name("25. Ensimmainen laillisesti crispr-tekniikalla geenimuunneltu ihminen syntyy vuoden 2028 loppuun mennessa.")

# ╔═╡ 0bce4e76-5187-11eb-2fbd-a5a81560ac90
function clean_column_name(s)
    rxb = r"^\d+"
    "Q" * match(rxb, s).match
end

# ╔═╡ 93b434d4-5189-11eb-32be-e97d67e6f555
begin
	column_names = Dict{Symbol,String}()
	df = DataFrame()

	let cnames, nimi, tyyppi
		cnames = names(df_raw)
		nimi = popfirst!(cnames)
		tyyppi = popfirst!(cnames)
		column_names[:nimi] = nimi
		df[!, :nimi] = df_raw[!, nimi]
		column_names[:tyyppi] = tyyppi
		df[!, :tyyppi] = df_raw[!, tyyppi]
		for cn in cnames
			n,t = split_column_name(cn)
			n = Symbol(clean_column_name(n))
			df[!, n] = df_raw[!, cn]
			column_names[n] = t
		end
	end
end;

# ╔═╡ 93552228-5189-11eb-0a42-5776e305433f
column_names

# ╔═╡ 824e0940-526f-11eb-383a-915bba34cc29
df

# ╔═╡ 3006417c-5187-11eb-0d0c-1bcaaf659578
clean_column_name("27.  ")

# ╔═╡ 9c4bc926-518f-11eb-3bd9-312fdacf8168
Symbol(clean_column_name("27.  "))

# ╔═╡ 8a5eca56-5194-11eb-2e18-9d8b842aba0c
function string_with_na_to_float(s)
    s=="NA" && return missing
    parse(Float32, s)
end

# ╔═╡ cddd9434-526f-11eb-3510-bd2b0159f846
df2 = select(df, :nimi, :tyyppi,
	intersect(names(df, r"^Q"), names(df, Float64)),
	intersect(names(df, r"^Q"), names(df, String)) .=> ByRow(string_with_na_to_float),
	renamecols=false
)

# ╔═╡ 35013218-51b9-11eb-275e-d5a3685e2ea0
begin
	@df df2 histogram(:Q1, ticks=:native)
	#yaxis!("count")
	xaxis!(column_names[:Q1])
end

# ╔═╡ e936638a-51dd-11eb-01b5-672579b31497
describe(df2, :first, :median, :mean, :eltype, :nmissing)

# ╔═╡ 8e371412-5194-11eb-221f-152be903cc6e
string_with_na_to_float("20")

# ╔═╡ 8ded747e-5194-11eb-122d-a90e49de9b31
string_with_na_to_float("NA")

# ╔═╡ Cell order:
# ╠═35013218-51b9-11eb-275e-d5a3685e2ea0
# ╠═bd3a666c-51c4-11eb-0ea5-3beaf9b624c2
# ╠═3a965f9c-510e-11eb-08dc-8d1cdc52e218
# ╠═30e459ae-526c-11eb-19b3-8f67825515c2
# ╠═0cdb58ae-510f-11eb-12d4-5393fb3dec0c
# ╠═93b434d4-5189-11eb-32be-e97d67e6f555
# ╠═93552228-5189-11eb-0a42-5776e305433f
# ╠═824e0940-526f-11eb-383a-915bba34cc29
# ╠═cddd9434-526f-11eb-3510-bd2b0159f846
# ╠═e936638a-51dd-11eb-01b5-672579b31497
# ╠═3493c430-51b9-11eb-2b59-6bd9bde30ed5
# ╠═9ebdd936-5187-11eb-2ed8-3fa4ec45d80c
# ╠═c9f9c422-5185-11eb-3e8d-3f96592baac1
# ╠═ef460ef8-5186-11eb-360f-cfb255b92d0f
# ╠═0bce4e76-5187-11eb-2fbd-a5a81560ac90
# ╠═3006417c-5187-11eb-0d0c-1bcaaf659578
# ╠═9c4bc926-518f-11eb-3bd9-312fdacf8168
# ╠═8a5eca56-5194-11eb-2e18-9d8b842aba0c
# ╠═8e371412-5194-11eb-221f-152be903cc6e
# ╠═8ded747e-5194-11eb-122d-a90e49de9b31

max_stars_repo_name: StatisticalMice/superforecasters

max_issues_repo_name: StatisticalMice/superforecasters

max_forks_repo_name: StatisticalMice/superforecasters


text:
module TestTransducers

using Transducers
using Folds.Testing: parse_tests, test_with_sequential

length_partition_by(f) = ReducePartitionBy(f, Map(_ -> 1)'(+), 0)
length_partition_by(f, xs) = length_partition_by(f)(xs)

divby(y) = x -> x ÷ y

rawdata = """
maximum(length_partition_by(divby(8), 1:10))
maximum(length_partition_by(divby(8), 5:19))
extrema(length_partition_by(divby(8), 1:10))
extrema(length_partition_by(divby(8), 5:19))
reduce(TeeRF(min, max), 1:10)
reduce(TeeRF(min, max), (2x for x in 1:10))
reduce(TeeRF(min, max), (2x for x in 1:10 if isodd(x)))
reduce(+, 1:10 |> Map(tuple) |> Broadcasting())
reduce(+, ((x, 2x) for x in 1:10) |> Broadcasting())
sum(1:10 |> Map(tuple) |> Broadcasting())
sum(((x, 2x) for x in 1:10) |> Broadcasting())
"""

tests = parse_tests(rawdata, @__MODULE__)

test_with_sequential(tests, [
    SequentialEx(),
    ThreadedEx(),
    ThreadedEx(basesize = 1),
])

end  # module


text:
@doc """
Takes the singles sided FFT. returns (frequencies, complex 1-sided FT)
"""-> function FFT_1s(TimeDomain::Vector; fs = 1, Window = "Hanning", PlotOn = false)
    ## This function takes in a time domain vector, optional sampling rate "fs," and optional Hanning window
    L = length(TimeDomain)
    if Window == "Hanning"
        CorFac = 2
        Window = WindowHanning(L)
    elseif (Window == "none")|(Window == "None")
        CorFac = 1
        Window = ones(L)
    else
        error("Window not recognized")

    end
    WindowedData = TimeDomain .* Window
    if L % 2 == 1
        EndIndex = Int(round(L / 2))
    else
        EndIndex = Int(round(L / 2))
    end
    FTData_1s = fft(WindowedData)[1:EndIndex]
    FTData_1s = FTData_1s ./ L #Corrects for vector length
    freqs = fftfreq(L, fs)[1:EndIndex]

    if PlotOn
        plot(fftfreq(L), abs.(fft(WindowedData)) ./ L)
        plot(freqs, abs.(FTData_1s), "ro")
    end

    return freqs, FTData_1s
end

@doc """
Creates a Hanning window of length N
"""-> function WindowHanning(N)
    return 0.5 .- 0.5 .* cos.(2 .* pi .* collect(0:(N-1)) ./ (N - 1))
end

@doc """
This is a documentation test2
"""-> function known_aliasing(
    Data::Vector,
    fs,
    f_unwrap::Union{Float64,Vector};
    FreqBuf::Int = 1,
    oversamp::Int = 1,
    PlotOn=false
)

    freqs, ftData = FFT_1s(Data, fs = fs)
    L = length(Data)
    L_1s = length(ftData)
    T = 1 / fs
    T_vec = 0:T:((length(Data)-1)*T)
    df = fs / L

    MaxFreq = maximum(f_unwrap)
    ScaleFac = Int(ceil(MaxFreq / (fs / 2))) * oversamp
    new_fs = ScaleFac * fs
    new_L = ScaleFac * L
    new_L_1s = Int(round(new_L / 2))

    newFreqs = fftfreq(ScaleFac * L, new_fs)[1:new_L_1s]
    NewSpectra = Complex.(zeros(length(newFreqs)))
    NewSpectra[1:L_1s] = ftData[:]

    MovedIndex = [0] #freqs which have been moved to a new location due to aliasing already

for i = 1:length(f_unwrap)
    if f_unwrap[i] > fs / 2
        AliasedFreq = calcAlias(fs, f_unwrap[i]) #the freq. the frequency has aliased to
        AliasedIndex = Int(round(AliasedFreq / df)) #the index of the aliased frequency
        if AliasedIndex == 0
            AliasedIndex = 1
        end
        NewIndex = Int(round(f_unwrap[i] / df)) #The new (corrected) index

        if findfirst(x -> x == NewIndex, MovedIndex) == nothing
            push!(MovedIndex, NewIndex)

            IndexDistFromEdge = [
                AliasedIndex,
                length(AliasedFreq) - AliasedIndex,
                NewIndex,
                new_L_1s - NewIndex,
            ]
            if minimum(abs.(IndexDistFromEdge)) < FreqBuf
                FreqBuf = minimum(abs.(IndexDistFromEdge))
            end

            StartInd_Alias = AliasedIndex - FreqBuf
            StopInd_Alias = AliasedIndex + FreqBuf
            StartInd_New = NewIndex - FreqBuf
            StopInd_New = NewIndex + FreqBuf
            if StopInd_Alias>length(ftData)
                FreqBuf = FreqBuf-(StopInd_Alias-length(ftData))
            end
            if StartInd_Alias<0
                FreqBuf = 0
            end
            if StopInd_New>length(NewSpectra)
                FreqBuf = FreqBuf-(StopInd_New-length(NewSpectra))
            end
            if StartInd_New<0
                FreqBuf = 0
            end

            StartInd_Alias = AliasedIndex - FreqBuf
            StopInd_Alias = AliasedIndex + FreqBuf
            StartInd_New = NewIndex - FreqBuf
            StopInd_New = NewIndex + FreqBuf

            println(length(ftData))
            println(StartInd_Alias)
            println(StopInd_Alias)
            println(StartInd_New)
            println(StopInd_New)

            NewSpectra[StartInd_New:StopInd_New] =
                ftData[StartInd_Alias:StopInd_Alias]
            NewSpectra[StartInd_Alias:StopInd_Alias] =
                zeros(length(StartInd_Alias:StopInd_Alias))

            println("Actual freq: $(f_unwrap[i]). Aliased to: $(AliasedFreq)")
        end
    end
end

FullSpectraNew = FFT_1s_to2s(NewSpectra)
NewTimeDomain = ifft(FullSpectraNew .* length(FullSpectraNew))
T_vec_new = 0:T/ScaleFac:((length(NewTimeDomain)-1)*T/ScaleFac)
if PlotOn
    figure(1)
    subplot(211)
    plot(freqs, abs.(ftData), "g",label="Orig data")
    stem(newFreqs, abs.(NewSpectra),label="New data")
    legend()

    subplot(212)

    # figure(2)
    plot(T_vec, Data, "r",label="Orig data")
    plot(T_vec_new, NewTimeDomain, "g",label="New data")
    legend()
end
    return FullSpectraNew, NewTimeDomain, T_vec_new

end

@doc """
Calculates the frequency some other freq. will alias to given a sampling rate, fs
calcAlias(fs, f)
fs is sampling rate
f is the original freq
"""-> function calcAlias(fs, f)
    #fs is the sampling rate
    #df is the df in the DFT
    #f is the actual waveform frequency

    Nyqu = fs / 2

    Wraps = Int(floor(f / Nyqu)) #numer of full wraparounds
    WrapDist = f % Nyqu # how far past the last wrap it will be

    if (Wraps % 2) == 0#if it is even, the wrap will be "backwards"
        return AliasedFreq = WrapDist
    else
        return AliasedFreq = Nyqu - WrapDist
    end
end

function FFT_1s_to2s(ftData)
    return cat(ftData[:], conj.(reverse(ftData[2:end]))[:], dims = 1)
end


text:
using PuppyBook
using Test
using Logging

function test_running_proportion_length(ntests=10)
    n_trials = rand(1:100,ntests)
    has_correct_length = Bool[]
    for n in n_trials
        is_correct = length(compute_running_proportion(0.1,n)) == n
        is_correct || @error "Test failed with N = $n"
        append!(has_correct_length, is_correct)
    end
    return all(has_correct_length)
end

function test_running_proportion_convergence(;tol=1e-1, ntrials=1000)
    p_heads = rand(0.01:0.01:0.99, 10)
    converged = []
    for p in p_heads
        convd = isapprox(compute_running_proportion(p, ntrials)[end],p,atol=tol)
        convd || @error "Test failed with p_heads=$p and N=$ntrials"
        append!(converged, convd)
    end
    return all(converged)
end

@testset "PuppyBook.jl" begin
    @testset "03_running_proportion.jl" begin
        @test_throws AssertionError compute_running_proportion(-2.0, 10)
        @test_throws AssertionError compute_running_proportion(1.01, 10)
        @test test_running_proportion_length()
        @test test_running_proportion_convergence()
    end
end


text:
module SnaggedNecklace

using LinearAlgebra: Symmetric, dot, eigen, normalize
using Statistics: mean, std
using StatsBase: Histogram, fit

export
    NmmSystem,
    NmmParameters,
    nmm,

    PimcSystem,
    StepSize1D1,
    StepSize1Dj,
    State1D,
    StepSize3D1,
    StepSize3Dj,
    State3D,
    PimcParameters,
    pimc

include("units.jl")
include("stats.jl")

include("nmm_system.jl")
include("nmm.jl")

include("pimc_system.jl")
include("pimc_1d.jl")
include("pimc_3d.jl")
include("pimc_integrate.jl")

end


text:
# Structs and constants related to the codes_handle type

export Message, getbytes, writemessage, missingvalue, clone, data, maskedvalues

mutable struct codes_handle
end

mutable struct Message
    ptr::Ptr{codes_handle}
end

"""
    Message(index::Index)
    Message(file::GribFile)

Retrieve the next message.

Returns `nothing` if there are no more messages.
"""
function Message(index::Index)::Union{Message, Nothing}
    eref = Ref(Int32(0))
    ptr = ccall((:codes_handle_new_from_index, eccodes), Ptr{codes_handle},
                (Ptr{codes_index}, Ref{Cint}), index.ptr, eref)
    if eref[] != 0 && eref[] != CODES_END_OF_INDEX
        throw(ErrorException("Error: $(errors[-eref[]+1])"))
    end
    if ptr == C_NULL || eref[] == CODES_END_OF_INDEX
        nothing
    else
        finalizer(destroy, Message(ptr))
    end
end

function Message(file::GribFile)::Union{Message, Nothing}
    eref = Ref(Int32(0))
    ptr = ccall((:codes_grib_handle_new_from_file, eccodes), Ptr{codes_handle},
                (Ptr{codes_context}, Ptr{File}, Ref{Cint}),
                C_NULL, file.ptr, eref)
    file.pos += 1
    errorcheck(eref[])
    if ptr == C_NULL
        nothing
    else
        finalizer(destroy, Message(ptr))
    end
end

"""
    Message(raw::Vector{UInt8})
    Message(messages::Vector{Vector{UInt8}})

Create a Message from a bytes array or array of bytes arrays.

The Message will share its underlying data with the inputed data.
"""
function Message(raw::Vector{UInt8})
    ptr = ccall((:codes_handle_new_from_message, eccodes), Ptr{codes_handle},
                (Ptr{codes_context}, Ptr{UInt8}, Csize_t),
                C_NULL, raw, length(raw))
    if ptr == C_NULL
        throw(ErrorException("Invalid message, could not create Message"))
    end
    finalizer(destroy, Message(ptr))
end

function Message(messages::Vector{Vector{UInt8}})
    sizes = map(length, messages)
    eref = Ref(Int32(0))
    ptr = ccall((:codes_grib_handle_new_from_multi_message, eccodes), Ptr{codes_handle},
                (Ptr{codes_context}, Ptr{Ptr{UInt8}}, Ptr{Csize_t}, Ptr{Cint}),
                C_NULL, messages, sizes, eref)
    errorcheck(eref[])
    finalizer(destroy, Message(ptr))
end

"""
    getbytes(message::Message)

Get a coded representation of the Message.
"""
function getbytes(message::Message)
    lenref = Ref(Csize_t(0))
    vp = Ref{Ptr{UInt8}}(C_NULL)
    err = ccall((:codes_get_message, eccodes), Cint,
                (Ptr{codes_handle}, Ptr{Ptr{UInt8}}, Ref{Csize_t}),
                message.ptr, vp, lenref)
    errorcheck(err)
    unsafe_wrap(Vector{UInt8}, vp[], lenref[])
end

"""
    writemessage(handle::Message, filename::AbstractString; mode="w")

Write the message respresented by `handle` to the file at `filename`.

`mode` is a mode as described by `Base.open`.
"""
function writemessage(handle::Message, filename::AbstractString; mode="w")
    err = ccall((:codes_write_message, eccodes), Cint, (Ptr{codes_handle}, Cstring, Cstring),
                handle.ptr, filename, mode)
    errorcheck(err)
end

function Base.haskey(msg::Message, key::AbstractString)
    try
        msg[key]
        true
    catch e
        false
    end
end

"""
    missingvalue(msg::Message)

Return the missing value of the message. If one isn't included returns 1e30.
"""
function missingvalue(msg::Message)
    if haskey(msg, "missingValue")
        msg["missingValue"]
    else
        1e30 # Why do we do this? Could we use nothing?
    end
end

"""
    maskedvalues(msg::Message)

Return the values of the message masked so the missing value is `missing`.
"""
function maskedvalues(msg::Message)
    missingval = missingvalue(msg)
    vals = msg["values"]
    masked = Array{Union{Float64, Missing}, 2}(undef, size(vals)...)
    for i in eachindex(vals)
        @inbounds masked[i] = vals[i] == missingval ? missing : vals[i]
    end
    masked
end

""" Get the native type of a key. """
function getnativetype(handle::Message, key::AbstractString)::DataType
    typeref = Ref(Int32(0))
    err = ccall((:codes_get_native_type, eccodes), Cint, (Ptr{codes_handle}, Cstring, Ref{Cint}),
                handle.ptr, key, typeref)
    errorcheck(err)
    typeint = typeref[]

    if typeint == CODES_TYPE_LONG
        Clong
    elseif typeint == CODES_TYPE_STRING
        String
    elseif typeint == CODES_TYPE_DOUBLE
        Float64
    elseif typeint == CODES_TYPE_BYTES
        Vector{UInt8}
    elseif typeint == CODES_TYPE_MISSING
        Missing
    else
        Tuple{}
    end
end

""" Get the number of values in a key. """
function getsize(message::Message, key::AbstractString)
    lenref = Ref(Csize_t(0))
    err = ccall((:codes_get_size, eccodes), Cint, (Ptr{codes_handle}, Cstring, Ref{Csize_t}),
                 message.ptr, key, lenref)
    errorcheck(err)
    lenref[]
end

function getvalues(::Type{String}, handle, key)
    bufr = UInt(1024)
    bufref = Ref(bufr)
    mesg = zeros(Cuchar, bufr)
    err = ccall((:codes_get_string, eccodes), Cint,
                (Ptr{codes_handle}, Cstring, Ptr{Cuchar}, Ref{Csize_t}),
                handle.ptr, key, mesg, bufref)

    mesgstr = rstrip(transcode(String, mesg), '\0')
    errorcheck(err)
    mesgstr
end

function getvalues(::Type{Float64}, handle, key)
    len = getsize(handle, key)
    if len == 1
        valref = Ref(Float64(0))
        err = ccall((:codes_get_double, eccodes), Cint,
                    (Ptr{codes_handle}, Cstring, Ref{Cdouble}),
                    handle.ptr, key, valref)
        errorcheck(err)
        valref[]
    else
        lenref = Ref(len)
        vals = Vector{Float64}(undef, len)
        err = ccall((:codes_get_double_array, eccodes), Cint,
                   (Ptr{codes_handle}, Cstring, Ref{Cdouble}, Ref{Csize_t}),
                    handle.ptr, key, vals, lenref)
        errorcheck(err)
        if key == "values" && handle["Ni"] != (2^31 - 1)
            ni = handle["Ni"]
            nj = handle["Nj"]
            columnmajor = handle["jPointsAreConsecutive"] != 0
            if columnmajor
                reshape(vals, nj, ni)
            else
                reshape(vals, ni, nj)
            end
        else
            vals
        end
    end
end

function getvalues(::Type{Clong}, handle, key)
    len = getsize(handle, key)
    if len == 1
        valref = Ref(Clong(0))
        err = ccall((:codes_get_long, eccodes), Cint, (Ptr{codes_handle}, Cstring, Ref{Clong}),
                    handle.ptr, key, valref)
        errorcheck(err)
        valref[]
    else
        lenref = Ref(len)
        vals = Vector{Clong}(undef, len)
        err = ccall((:codes_get_long_array, eccodes), Cint,
                   (Ptr{codes_handle}, Cstring, Ref{Clong}, Ref{Csize_t}),
                    handle.ptr, key, vals, lenref)
        errorcheck(err)
        vals
    end
end

function getvalues(::Type{Vector{UInt8}}, handle, key)
    bufr = UInt(1024)
    bufref = Ref(bufr)
    mesg = Vector{UInt8}(undef, bufr)
    err = ccall((:codes_get_bytes, eccodes), Cint,
                (Ptr{codes_handle}, Cstring, Ptr{Cuchar}, Ref{Csize_t}),
                handle.ptr, key, mesg, bufref)
    errorcheck(err)
    mesg[1:bufref[]]
end

# Throw an error for an unspecialized type
getvalues(typ::Type{Any}, handle, key) = throw(ErrorException("Could not get key: unsupported type $typ"))

""" Get the value of a key. """
Base.getindex(handle::Message, key::AbstractString) = getvalues(getnativetype(handle, key), handle, key)

"""
    setindex!(message::Message, value, key::AbstractString)

Set the value of a key.
"""
function Base.setindex!(handle::Message, value::Integer, key::AbstractString)
    err = ccall((:codes_set_long, eccodes), Cint, (Ptr{codes_handle}, Cstring, Clong),
                 handle.ptr, key, value)
    errorcheck(err)
end

function Base.setindex!(handle::Message, value::AbstractFloat, key::AbstractString)
    err = ccall((:codes_set_double, eccodes), Cint, (Ptr{codes_handle}, Cstring, Cdouble),
                 handle.ptr, key, value)
    errorcheck(err)
end

function Base.setindex!(handle::Message, value::AbstractString, key::AbstractString)
    lenref = Ref(Csize_t(length(value)))
    err = ccall((:codes_set_string, eccodes), Cint,
                (Ptr{codes_handle}, Cstring, Cstring, Ref{Csize_t}),
                 handle.ptr, key, value, lenref)
    errorcheck(err)
end

function Base.setindex!(handle::Message, value::Vector{UInt8}, key::AbstractString)
    lenref = Ref(Csize_t(length(value)))
    err = ccall((:codes_set_bytes, eccodes), Cint,
                (Ptr{codes_handle}, Cstring, Ref{UInt8}, Ref{Csize_t}),
                 handle.ptr, key, value, lenref)
    errorcheck(err)
end

function Base.setindex!(handle::Message, value::Array{Float64}, key::AbstractString)
    len = Csize_t(length(value))
    err = ccall((:codes_set_double_array, eccodes), Cint,
                (Ptr{codes_handle}, Cstring, Ref{Cdouble}, Csize_t),
                handle.ptr, key, value, len)
    errorcheck(err)
end

function Base.setindex!(handle::Message, value::Array{Clong}, key::AbstractString)
    len = Csize_t(length(value))
    err = ccall((:codes_set_long_array, eccodes), Cint,
                (Ptr{codes_handle}, Cstring, Ref{Clong}, Csize_t),
                 handle.ptr, key, value, len)
    errorcheck(err)
end

function Base.setindex!(handle::Message, value::Array{String}, key::AbstractString)
    len = Csize_t(length(value))
    err = ccall((:codes_set_string_array, eccodes), Cint,
                (Ptr{codes_handle}, Cstring, Ref{Cstring}, Csize_t),
                 handle.ptr, key, value, len)
    errorcheck(err)
end

""" Duplicate a message. """
function clone(handle::Message)::Union{Message, Nothing}
    ptr = ccall((:codes_handle_clone, eccodes), Ptr{codes_handle}, (Ptr{codes_handle},), handle.ptr)
    if ptr == C_NULL
        nothing
    else
        finalizer(destroy, Message(ptr))
    end
end

# Aliases for clone for more Julian style
Message(handle::Message) = clone(handle)
Base.deepcopy(handle::Message) = clone(handle)

"""
    data(handle::Message)

Retrieve the longitudes, latitudes, and values from the message.

# Example
```
GribFile(filename) do f
    msg = Message(f)
    lons, lats, values = data(msg)
end
```
"""
function data(handle::Message)
    # Allocate arrays
    npoints = handle["numberOfPoints"]
    lons = Vector{Float64}(undef, npoints)
    lats = Vector{Float64}(undef, npoints)
    values = Vector{Float64}(undef, npoints)
    err = ccall((:codes_grib_get_data, eccodes), Cint,
                (Ptr{codes_handle}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
                handle.ptr, lats, lons, values)
    errorcheck(err)

    ni = handle["Ni"]
    if ni != (2^31 - 1)
        nj = handle["Nj"]
        columnmajor = handle["jPointsAreConsecutive"] != 0
        if columnmajor
            return reshape(lons, nj, ni), reshape(lats, nj, ni), reshape(values, nj, ni)
        else
            return reshape(lons, ni, nj), reshape(lats, ni, nj), reshape(values, ni, nj)
        end
    else
        lons, lats, values
    end
end

""" Safely destroy a message handle. """
function destroy(handle::Message)
    if handle.ptr != C_NULL
        err = ccall((:codes_handle_delete, eccodes), Cint, (Ptr{codes_handle},), handle.ptr)
        #errorcheck(err) Calling this throws, which is not allowed in finalizers
        handle.ptr == C_NULL
    end
end

# Functions related to printing
function Base.show(io::IO, mime::MIME"text/plain", message::Message)
    dispkeys = ["date", "gridType", "stepRange", "typeOfLevel", "level",
                 "shortName", "name"]
    dispkeys = [k for k in dispkeys if haskey(message, k)]
    offsets = [9, 15, 10, 18, 6, 10, 4]
    line1 = strip(join([rpad(k, offsets[i]) for (i, k) in enumerate(dispkeys)])) * "\n"
    line2 = strip(join([rpad(message[k], offsets[i]) for (i, k) in enumerate(dispkeys)]))
    write(io, line1 * line2)
end

# So it pretty-prints for arrays
Base.show(io::IO, message::Message) = show(io, MIME("text/plain"), message)

function Base.print(io::IO, message::Message)
    bytes = getbytes(message)
    str = bytes2hex(bytes)
    print(io, str)
end


text:
using Test, GetGene, DataFrames

snps = ["rs113980419", "rs17367504",
"rs113980419", "rs2392929", "rs11014171",
"78067132", "1465537", "Affx-4150211"]

loci = ["C1orf167", "MTHFR", "C1orf167",
"No gene listed", "CACNB2", "No gene listed",
"No gene listed", "snpid not in database"]

genedict = Dict("annotation_release" => "Homo sapiens Annotation Release 109",
"gene_name"          => "solute carrier family 39 member 8",
"gene_locus"         => "SLC39A8",
"gene_id"            => "64116",
"seq_id"             => "NC_000004.12",
"gene_is_pseudo"     => "0",
"gene_orientation"   => "1")

emptydict = Dict("Empty" => "No gene info listed for ref SNP ID")

df = DataFrame(rsid = snps)

@testset "snpid match" begin
    @test getgenes(snps) == loci
    @test getgenes(df; idvar = "rsid") == loci
    @test getgenes(2270993) == "PTPRJ"
end

@testset "Argument Error" begin
    @test_throws ArgumentError getgenes(df)
    @test_throws ArgumentError getgeneinfo("Affx-4150211")
end

@testset "get gene info" begin
    @test getgeneinfo("rs13107325") == genedict
    @test getgeneinfo("rs1173771") == emptydict
end