text:
module TerminalRegressionTests
    using VT100
    import REPL

    function load_outputs(file)
        outputs = String[]
        decorators = String[]
        is_output = true
        is_continuation = false
        nlines = 0
        for line in readlines(file)
            if line[1] == '+'
                is_continuation = line[2] != '+'
                if is_continuation
                  until = findfirst(c -> c == '+', line[2:end])+1
                  nlines = parse(Int, line[2:until - 1])
                  push!(outputs,
                    join(split(outputs[end],'\n')[1:end-nlines],'\n'))
                else
                  push!(outputs,"")
                end
                is_output = true
                continue
            elseif line[1] == '-'
                if is_continuation
                  push!(decorators,
                    join(split(decorators[end],'\n')[1:end-nlines],'\n'))
                else
                  push!(decorators,"")
                end
                is_output = false
                continue
            elseif line[1] == '|'
                array = is_output ? outputs : decorators
                array[end] = string(array[end],isempty(array[end]) ? "" : "\n",line[2:end])
            else
                error("Unrecognized first character \"$(line[1])\"")
            end
        end
        outputs, decorators
    end

    mutable struct EmulatedTerminal <: REPL.Terminals.UnixTerminal
        input_buffer::IOBuffer
        out_stream::Base.TTY
        pty::VT100.PTY
        terminal::VT100.ScreenEmulator
        waiting::Bool
        step::Condition
        filled::Condition
        # Yield after every write, e.g. to test for flickering issues
        aggressive_yield::Bool
        function EmulatedTerminal()
            pty = VT100.create_pty(false)
            new(
            IOBuffer(UInt8[], true, true, true, true, typemax(Int)),
            Base.TTY(pty.slave), pty,
            pty.em, false, Condition(), Condition())
        end
    end
    function Base.wait(term::EmulatedTerminal)
        if !term.waiting || bytesavailable(term.input_buffer) != 0
            wait(term.step)
        end
    end
    for T in (Vector{UInt8}, Array, AbstractArray, String, Symbol, Any, Char, UInt8)
        function Base.write(term::EmulatedTerminal,a::T)
            write(term.out_stream, a)
            if term.aggressive_yield
                notify(term.step)
            end
        end
    end
    Base.eof(term::EmulatedTerminal) = false
    function Base.read(term::EmulatedTerminal, ::Type{Char})
        if bytesavailable(term.input_buffer) == 0
            term.waiting = true
            notify(term.step)
            wait(term.filled)
        end
        term.waiting = false
        read(term.input_buffer, Char)
    end
    function Base.readuntil(term::EmulatedTerminal, delim::UInt8; kwargs...)
        if bytesavailable(term.input_buffer) == 0
            term.waiting = true
            notify(term.step)
            wait(term.filled)
        end
        term.waiting = false
        readuntil(term.input_buffer, delim; kwargs...)
    end
    REPL.Terminals.raw!(t::EmulatedTerminal, raw::Bool) =
        ccall(:jl_tty_set_mode,
                 Int32, (Ptr{Cvoid},Int32),
                 t.out_stream.handle, raw) != -1
    REPL.Terminals.pipe_reader(t::EmulatedTerminal) = t.input_buffer
    REPL.Terminals.pipe_writer(t::EmulatedTerminal) = t.out_stream

    function _compare(output, outbuf)
        result = outbuf == output
        if !result
            println("Test failed. Expected result written to expected.out,
                actual result written to failed.out")
            open("failed.out","w") do f
                write(f,outbuf)
            end
            open("expected.out","w") do f
                write(f,output)
            end
            for (i,c) in enumerate(output)
                if c == Char(outbuf[i])
                    print(Char(c))
                elseif c == '\n'
                    println()
                else
                    printstyled(stdout, "█", color=:red)
                end
            end
            error()
        end
        return result
    end

    function compare(em, output, decorator = nothing)
        buf = IOBuffer()
        decoratorbuf = IOBuffer()
        VT100.dump(buf,decoratorbuf,em)
        outbuf = take!(buf)
        decoratorbuf = take!(decoratorbuf)
        _compare(Vector{UInt8}(codeunits(output)), outbuf) || return false
        decorator === nothing && return true
        _compare(Vector{UInt8}(codeunits(decorator)), decoratorbuf)
    end

    function process_all_buffered(emuterm)
        # Since writes to the tty are asynchronous, there's an
        # inherent race condition between them being sent to the
        # kernel and being available to epoll. We write a sentintel value
        # here and wait for it to be read back.
        sentinel = Ref{UInt32}(0xffffffff)
        ccall(:write, Cvoid, (Cint, Ptr{UInt32}, Csize_t), emuterm.pty.slave, sentinel, sizeof(UInt32))
        Base.process_events(false)
        # Read until we get our sentinel
        while bytesavailable(emuterm.pty.master) < sizeof(UInt32) ||
            reinterpret(UInt32, emuterm.pty.master.buffer.data[(emuterm.pty.master.buffer.size-3):emuterm.pty.master.buffer.size])[] != sentinel[]
            emuterm.aggressive_yield || yield()
            Base.process_events(false)
            sleep(0.01)
        end
        data = IOBuffer(readavailable(emuterm.pty.master)[1:(end-4)])
        while bytesavailable(data) > 0
            VT100.parse!(emuterm.terminal, data)
        end
    end

    function automated_test(f, outputpath, inputs; aggressive_yield = false)
        emuterm = EmulatedTerminal()
        emuterm.aggressive_yield = aggressive_yield
        emuterm.terminal.warn = true
        outputs, decorators = load_outputs(outputpath)
        c = Condition()
        @async Base.wait_readnb(emuterm.pty.master, typemax(Int64))
        yield()
        t1 = @async try
            f(emuterm)
            Base.notify(c)
        catch err
            Base.showerror(stderr, err, catch_backtrace())
            Base.notify_error(c, err)
        end
        t2 = @async try
            for input in inputs
                wait(emuterm);
                emuterm.aggressive_yield || @assert emuterm.waiting
                output = popfirst!(outputs)
                decorator = isempty(decorators) ? nothing : popfirst!(decorators)
                @assert !eof(emuterm.pty.master)
                process_all_buffered(emuterm)
                compare(emuterm.terminal, output, decorator)
                print(emuterm.input_buffer, input); notify(emuterm.filled)
            end
            Base.notify(c)
        catch err
            Base.showerror(stderr, err, catch_backtrace())
            Base.notify_error(c, err)
        end
        while !istaskdone(t1) || !istaskdone(t2)
            wait(c)
        end
    end

    function create_automated_test(f, outputpath, inputs; aggressive_yield=false)
        emuterm = EmulatedTerminal()
        emuterm.aggressive_yield = aggressive_yield
        emuterm.terminal.warn = true
        c = Condition()
        @async Base.wait_readnb(emuterm.pty.master, typemax(Int64))
        yield()
        t1 = @async try
            f(emuterm)
            Base.notify(c)
        catch err
            Base.showerror(stderr, err, catch_backtrace())
            Base.notify_error(c, err)
        end
        t2 = @async try
            outs = map(inputs) do input
                wait(emuterm);
                emuterm.aggressive_yield || @assert emuterm.waiting
                process_all_buffered(emuterm)
                out = IOBuffer()
                decorator = IOBuffer()
                VT100.dump(out, decorator, emuterm.terminal)
                print(emuterm.input_buffer, input); notify(emuterm.filled)
                out, decorator
            end
            open(outputpath, "w") do io
                print(io,"+"^50,"\n",
                    join(map(outs) do x
                        sprint() do io
                            out, dec = x
                            print(io, "|", replace(String(take!(out)),"\n" => "\n|"))
                            println(io, "\n", "-"^50)
                            print(io, "|", replace(String(take!(dec)),"\n" => "\n|"))
                        end
                    end, string('\n',"+"^50,'\n')))
            end
            Base.notify(c)
        catch err
            Base.showerror(stderr, err, catch_backtrace())
            Base.notify_error(c, err)
        end
        while !istaskdone(t1) || !istaskdone(t2)
            wait(c)
        end
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Imath_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Imath")
JLLWrappers.@generate_main_file("Imath", UUID("905a6f67-0a94-5f89-b386-d35d92009cd1"))
end  # module Imath_jll

ext: jl


text:
using SupremumInfimum
using InteractiveUtils
using Test

# write your own tests here
@test leaf_type(Int) == Int
@test leaf_type(UInt8) == UInt8
@test leaf_type(Real) == Float64
@test leaf_type(Integer) == Int
@test leaf_type(Unsigned) == UInt

@test supremum(Set{Int}[]) == Set{Int}()

@testset "Real type $T basics" for T in [subtypes(AbstractFloat); subtypes(Signed); subtypes(Unsigned)]

  let  E = one(T), Z = zero(T), IP = supremum(T), IN = infimum(T)

  @test supremum(T[]) == IN
  @test infimum(T[]) == IP # Inf - Note the names match!
  @test sup(E) == E
  @test inf(E) == E
  @test supremum([Z]) == Z
  @test infimum([Z]) == Z
  @test sup(IP) == IP
  @test inf(IP) == IP
  @test supremum([IN]) == IN
  @test infimum([IN]) == IN

  @test sup(E, Z) == E
  @test inf(E, Z) == Z
  @test sup(Z, E) == E
  @test inf(Z, E) == Z
  @test supremum([E, Z]) == E
  @test infimum([E, Z]) == Z
  @test supremum([Z, E]) == E
  @test infimum([Z, E]) == Z
  @test supremum([E, IP]) == IP
  @test infimum([E, IP]) == E
  @test supremum([IP, Z]) == IP
  @test infimum([IP, Z]) == Z
  @test supremum([E, IN]) == E
  @test infimum([E, IN]) == IN
  @test supremum([IN, Z]) == Z
  @test infimum([IN, Z]) == IN
  end
end

@testset "Set operations" begin
  let Z = Set{AbstractFloat}(), E = Set([Float16(1.0), 2]), A = Set([1.0, 3])

  @test inf(Z, E) == inf(E, Z) == Z
  @test sup(Z, E) == sup(Z, E) == E
  @test sup(E, A, Z) == Set([1.0, 2, 3])
  @test sup(Z, A, E) == Set([1.0, 2, 3])
  @test inf(E, A, Z) == Set{Float64}()
  @test inf(Z, A, E) == Set{Float64}()

  @test supremum([E, A]) == Set([1.0, 2, 3])
  @test supremum(maximum, [E, A]) === 3.0
  @test_throws ArgumentError supremum(maximum, [E, Z]) === 2.0
  @test supremum(supremum, [E, Z]) === 2.0
  @test supremum(length, [E, Z]) === 2
  @test infimum(length, [E, Z]) === 0

  end
end
@test infimum(Float64[]) == Inf # Note the names match!
@test supremum(Float64[]) == -Inf
@test sup(BitSet([1,42]), BitSet([2,42])) == BitSet([1, 2, 42])
@test inf(BitSet([1,42]), BitSet([2,42])) == BitSet([42])


text:
using Lorenz96
using Test

@testset "Bounded" begin
    X = L96(N=1000)
    @test maximum(X) < 100.0
    @test minimum(X) > -100.0
end

@testset "Chaotic" begin
    X = L96(N=1000,output=true)
    @test all(X[:,end-1] .!= X[:,end])
end

@testset "NoForcing" begin
    X = L96(N=1000,F=0.0,output=true)
    @test all(X[2:end,:] .== 0.0)
end

@testset "InitialCond" begin
    ini = randn(36)
    Xf64 = L96(Float64,X=ini,N=10,output=true)
    Xf32 = L96(Float32,X=ini,N=10,output=true)
    Xf16 = L96(Float16,X=ini,N=10,output=true)

    @test Xf64[:,1] == ini
    @test Xf32[:,1] == Float32.(ini)
    @test Xf16[:,1] == Float16.(ini)
end

@testset "Mixed Precision" begin
    X1 = L96(Float32,Float32,N=1000)
    X2 = L96(Float32,N=1000)

    @test X1 == X2

    X3 = L96(Float16,Float32,N=1000)
    @test X1 != X3
end


text:
# ## Data

import Random

Random.seed!(1)

n = 9
k = 6;

# ## Model

using Turing

@model function globe_toss(n, k)
    θ ~ Beta(1, 1)
    k ~ Binomial(n, θ)
    return k, θ
end;

# ## Output

using Random

Random.seed!(1)
chns = sample(globe_toss(n, k), NUTS(), 1000)

# \defaultoutput{}


text:
using GNSSTools
using JLD
using PyPlot
pygui(true)
include("ch4_data_plot_functions.jl")

directory = "/home/sjbilardi/Dropbox/Apps/Overleaf/MastersThesis/"
folder = "figures/l5_data/"

# file = "/home/sjbilardi/Documents/SeNSeLab/GNSSTools.jl/test/data/20200622_203230_000100_g5_ref.sc8"
# data = loaddata(:sc8, file, 25e6, L5_freq, L5_freq, 1, skip_to=0)
file = "/home/sjbilardi/Documents/SeNSeLab/GNSSTools.jl/test/data/20201211_211300_GPSL5_1176.45_25.0.sc4"
data = loaddata(:sc4, file, 25e6, L5_freq, L5_freq, 2, skip_to=2)

channel="Q"
sigtype = define_l5_code_type(; channel=channel)
# sigtype = definesignaltype(definecodetype([l5q_codes, nh20],
#                                           [L5_chipping_rate, nh_chipping_rate]),
#                                           L5_freq, "Q")
replica = definereplica(sigtype, 25e6, 20e-3)
detectable_prns = []
for i in 1:32
    fd_est, n0_est, SNR_est,
    P_d, above_threshold = courseacquisition(data, replica, i; M=1,
                                             Δfd=1/(1*replica.t_length))
    if above_threshold && (P_d >= 0.9)
        detectable = true
        push!(detectable_prns, i)
    else
        detectable = false
    end
    f_code_d, f_code_dd = GNSSTools.calc_doppler_code_rate(nh_chipping_rate,
                                                           replica.signal_type.sig_freq,
                                                           fd_est, 0)
    n0 = calcinitcodephase(length(nh20), f_code_d, f_code_dd, data.f_s,
                           n0_est)
    println("PRN $i\t$SNR_est\t$fd_est\t$n0_est\t$n0\t$(Float64(round(P_d, digits=2)))\t$detectable")
end
println("$(length(detectable_prns))) Detected PRNs: $detectable_prns")

fig = figure(figsize=(7.5, 4))
acquisition_T = 20e-3
fine_acq_T = 100e-3
tracking_T = 20e-3
M = 1
linewidth = 0.5
tracking_σᵩ = Array{Float64}(undef, length(detectable_prns))
tracking_snr = Array{Float64}(undef, length(detectable_prns))
tracking_code_phase_err = Array{Float64}(undef, length(detectable_prns))
for i in 1:length(detectable_prns)
    prn = detectable_prns[i]
    ax = fig.add_subplot(2, 3, i)
    tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
    title("PRN $prn")
    acqresults, trackresults = process(data, sigtype, prn, channel;
                                       show_plot=false,
                                       fine_acq_method=:fft, state_num=3,
                                       dynamickf=true,
                                    #    h₀=h_parms[3], h₋₂=h_parms[1],
                                       dll_b=1,
                                       acquisition_T=acquisition_T,
                                       fine_acq_T=fine_acq_T,
                                       tracking_T=tracking_T,
                                       use_fine_acq=true,
                                       σω=1,
                                       M=M)
    ax.plot(trackresults.t, real.(trackresults.ZP), linewidth=linewidth)
    ax.plot(trackresults.t, imag.(trackresults.ZP), linewidth=linewidth)
    N = floor(Int, length(trackresults.t)/2)
    tracking_σᵩ[i] = std(trackresults.dphi_meas[N:end]) * 180 / π
    tracking_snr[i] = mean(trackresults.SNR[N:end])
    tracking_code_phase_err[i] = std(trackresults.code_err_meas[N:end])
end
subplots_adjust(hspace=0.3, wspace=0.4,  top=0.86, left=0.02, right=0.98,
                bottom=0.02)
plot_title = string("Course Acq. T=$(floor(Int, acquisition_T*1000))ms    ",
                    "Fine Acq. T=$(floor(Int, fine_acq_T*1000))ms    ",
                    "Tracking T=$(floor(Int, tracking_T*1000))ms    ",
                    "M=$(M)")
suptitle(plot_title)
T = floor(Int, tracking_T*1e3)
savefig(string(directory, folder, "ch4_l5_T=$(T)ms_acquired_prns.svg"), dpi=300)

for prn in detectable_prns
    acqresults, trackresults, corr_result, SNR_est,
    P_d, above_threshold = process(data,
                                   sigtype,
                                   prn, channel;
                                   fine_acq_method=:fft,
                                   return_corrresult=true,
                                   return_Pd=true,
                                   state_num=3,
                                   #    h₀=h_parms[3],
                                   #    h₋₂=h_parms[1],
                                   q_a=1,
                                   dll_b=1,
                                   acquisition_T=20e-3,
                                   fine_acq_T=100e-3,
                                   fine_acq_sub_T=1e-3,
                                   tracking_T=20e-3,
                                   show_plot=false,
                                   σω=1,
                                   M=1,
                                   use_fine_acq=true);
    T = floor(Int, trackresults.T*1e3)
    name = "pdfs/ch4_l5$(lowercase(channel))_T=$(T)ms_prn_$(trackresults.prn).pdf"
    file_name = string(directory, folder, name)
    plot_example_track_results(trackresults; f_if=0, code_length=L5_code_length,
                               save_to_file=file_name,
                               bottom=0.075, top=0.95, left=0.125, right=0.875,
                               hspace=0.75, wspace=0.35, marker_size=1,
                               line_width=0.8)
end

# acqresults, trackresults, corr_result, SNR_est,
# P_d, above_threshold = process(data,
#                                sigtype,
#                                3, "Q";
#                                fine_acq_method=:fft,
#                                return_corrresult=true,
#                                return_Pd=true,
#                                state_num=3,
#                                dynamickf=true,
#                             #    h₀=h_parms[3],
#                             #    h₋₂=h_parms[1],
#                                q_a=1,
#                                dll_b=0.1,
#                                acquisition_T=20e-3,
#                                fine_acq_T=100e-3,
#                                tracking_T=20e-3,
#                                M=1,
#                                show_plot=true,
#                                use_fine_acq=true);


text:
# utility for softmax
struct OneHotVector <: AbstractVector{Bool}
  ix::UInt32
  of::UInt32
end

Base.size(xs::OneHotVector) = (Int64(xs.of),)
Base.getindex(xs::OneHotVector, i::Integer) = i == xs.ix
Base.getindex(xs::OneHotVector, ::Colon) = OneHotVector(xs.ix, xs.of)

function onehot(l, labels)
  i = something(findfirst(isequal(l), labels), 0)
  i > 0 || error("Value $l is not in labels")
  OneHotVector(i, length(labels))
end

# linear
function update_grads!(loss::Linear, α::T, pred::Vector{SVector{L,T}}, target::AbstractVector{T}, δ::Vector{SVector{L,T}}, δ²::Vector{SVector{L,T}}, 𝑤::Vector{SVector{1,T}}) where {T <: AbstractFloat, L}
    @inbounds for i in eachindex(δ)
        δ[i] = SVector(2 * (pred[i][1] - target[i]) * 𝑤[i][1])
        δ²[i] = SVector(2 * 𝑤[i][1])
    end
end

# logistic - on linear predictor
function update_grads!(loss::Logistic, α::T, pred::Vector{SVector{L,T}}, target::AbstractVector{T}, δ::Vector{SVector{L,T}}, δ²::Vector{SVector{L,T}}, 𝑤::Vector{SVector{1,T}}) where {T <: AbstractFloat, L}
    @inbounds for i in eachindex(δ)
        δ[i] = (sigmoid(pred[i][1]) * (1 - target[i]) - (1 - sigmoid(pred[i][1])) * target[i]) * 𝑤[i]
        δ²[i] = sigmoid(pred[i][1]) * (1 - sigmoid(pred[i][1])) * 𝑤[i]
    end
end

# Poisson
function update_grads!(loss::Poisson, α::T, pred::Vector{SVector{L,T}}, target::AbstractVector{T}, δ::Vector{SVector{L,T}}, δ²::Vector{SVector{L,T}}, 𝑤::Vector{SVector{1,T}}) where {T <: AbstractFloat, L}
    @inbounds for i in eachindex(δ)
        δ[i] = (exp(pred[i][1]) .- target[i]) * 𝑤[i]
        δ²[i] = exp(pred[i][1]) * 𝑤[i]
    end
end

# L1
function update_grads!(loss::L1, α::T, pred::Vector{SVector{L,T}}, target::AbstractArray{T, 1}, δ::Vector{SVector{L,T}}, δ²::Vector{SVector{L,T}}, 𝑤::Vector{SVector{1,T}}) where {T <: AbstractFloat, L}
    @inbounds for i in eachindex(δ)
        δ[i] =  (α * max(target[i] - pred[i][1], 0) - (1-α) * max(pred[i][1] - target[i], 0)) * 𝑤[i]
    end
end

# Softmax
function update_grads!(loss::Softmax, α::T, pred::Vector{SVector{L,T}}, target::AbstractVector{S}, δ::Vector{SVector{L,T}}, δ²::Vector{SVector{L,T}}, 𝑤::Vector{SVector{1,T}}) where {T <: AbstractFloat, L, S <: Integer}
    # pred = pred - maximum.(pred)
    # sums = sum(exp.(pred), dims=2)
    @inbounds for i in 1:size(pred,1)
        pred[i] = SVector{L,T}(pred[i] .- maximum(pred[i]))
        sums = sum(exp.(pred[i]))
        δ[i] = SVector{L,T}((exp.(pred[i]) / sums - (onehot(target[i], 1:L))) .* 𝑤[i][1])
        δ²[i] = SVector{L,T}(1 / sums .* (1 .- exp.(pred[i]) / sums) .* 𝑤[i][1])
    end
end

# Quantile
function update_grads!(loss::Quantile, α::T, pred::Vector{SVector{L,T}}, target::AbstractVector{T}, δ::Vector{SVector{L,T}}, δ²::Vector{SVector{L,T}}, 𝑤::Vector{SVector{1,T}}) where {T <: AbstractFloat, L}
    @inbounds for i in eachindex(δ)
        δ[i] = target[i] > pred[i][1] ? SVector(α * 𝑤[i][1]) : SVector((α - 1) * 𝑤[i][1])
        δ²[i] = SVector(target[i] - pred[i][1]) # δ² serves to calculate the quantile value - hence no weighting on δ²
    end
end

# Gaussian - http://jrmeyer.github.io/machinelearning/2017/08/18/mle.html
# pred[i][1] = μ
# pred[i][2] = log(σ²)
function update_grads!(loss::Gaussian, α, pred::Vector{SVector{L,T}}, target::AbstractArray{T, 1}, δ::Vector{SVector{L,T}}, δ²::Vector{SVector{L,T}}, 𝑤::Vector{SVector{1,T}}) where {T <: AbstractFloat, L}
    @inbounds @threads for i in eachindex(δ)
        δ[i] = SVector((pred[i][1] - target[i]) / max(1e-8, exp(pred[i][2])) * 𝑤[i][1], 𝑤[i][1] / 2 * (1 - (pred[i][1] - target[i])^2 / max(1e-8, exp(pred[i][2]))))
        δ²[i] = SVector(𝑤[i][1] / max(1e-8, exp(pred[i][2])), 𝑤[i][1] / max(1e-8, exp(pred[i][2])) * (pred[i][1] - target[i])^2)
    end
end

# utility functions
function logit(x::AbstractArray{T, 1}) where T <: AbstractFloat
    @. x = log(x / (1 - x))
    return x
end

function logit(x::T) where T <: AbstractFloat
    x = log(x / (1 - x))
    return x
end

function sigmoid(x::AbstractArray{T, 1}) where T <: AbstractFloat
    @. x = 1 / (1 + exp(-x))
    return x
end

function sigmoid(x::T) where T <: AbstractFloat
    x = 1 / (1 + exp(-x))
    return x
end

function softmax(x::AbstractVector{T}) where T <: AbstractFloat
    x .-= maximum(x)
    x = exp.(x) ./ sum(exp.(x))
    return x
end

##############################
# get the gain metric
##############################
# GradientRegression
function get_gain(loss::S, ∑δ::SVector{L,T}, ∑δ²::SVector{L,T}, ∑𝑤::SVector{1,T}, λ::T) where {S <: GradientRegression, T <: AbstractFloat, L}
    gain = sum((∑δ .^ 2 ./ (∑δ² .+ λ .* ∑𝑤)) ./ 2)
    return gain
end

# MultiClassRegression
function get_gain(loss::S, ∑δ::SVector{L,T}, ∑δ²::SVector{L,T}, ∑𝑤::SVector{1,T}, λ::T) where {S <: MultiClassRegression, T <: AbstractFloat, L}
    gain = sum((∑δ .^ 2 ./ (∑δ² .+ λ .* ∑𝑤)) ./ 2)
    return gain
end

# L1 Regression
function get_gain(loss::S, ∑δ::SVector{L,T}, ∑δ²::SVector{L,T}, ∑𝑤::SVector{1,T}, λ::T) where {S <: L1Regression, T <: AbstractFloat, L}
    gain = sum(abs.(∑δ))
    return gain
end

# QuantileRegression
function get_gain(loss::S, ∑δ::SVector{L,T}, ∑δ²::SVector{L,T}, ∑𝑤::SVector{1,T}, λ::T) where {S <: QuantileRegression, T <: AbstractFloat, L}
    gain = sum(abs.(∑δ) ./ (1 .+ λ))
    return gain
end

# GaussianRegression
function get_gain(loss::S, ∑δ::SVector{L,T}, ∑δ²::SVector{L,T}, ∑𝑤::SVector{1,T}, λ::T) where {S <: GaussianRegression, T <: AbstractFloat, L}
    gain = sum((∑δ .^ 2 ./ (∑δ² .+ λ .* ∑𝑤)) ./ 2)
    return gain
end


text:
using FdeSolver
using Plots

## inputs
tSpan = [0, 5]     # [intial time, final time]
y0 = 1             # intial value
β = 0.6            # order of the derivative

λ = -10
par = λ
F(t, y, par)= par * y
JF(t, y, par) = par

## Numerical solution
t,y = FDEsolver(F, tSpan, y0, β, par)
t_J, y_J = FDEsolver(F, tSpan, y0, β, par, JF = JF)

using MittagLeffler # for the exact solution

plot(t, y, linewidth = 5, title = "Solution to a linear fractional equation",
    xaxis = "Time (t)", yaxis = "y(t)", ls = :dot, label = "Approximation")
plot!(t, t -> mittleff(β, λ * t .^ β), linewidth = 5, label = "Exact solution")
plot!(t_J, y_J, linewidth = 5, ls = :dash, label = "Approximation with jacobian")


text:
module TestUtilities

using Test
using MLJBase

export @testset_accelerated, include_everywhere, @test_mach_sequence,
    @test_model_sequence

using ComputationalResources
using ComputationalResources: CPUProcesses

macro testset_accelerated(name::String, var, ex)
    testset_accelerated(name, var, ex)
end
macro testset_accelerated(name::String, var, opts::Expr, ex)
    testset_accelerated(name, var, ex; eval(opts)...)
end
function testset_accelerated(name::String, var, ex; exclude=[])
    final_ex = quote
       local $var = CPU1()
        @testset $name $ex
    end

    resources = AbstractResource[CPUProcesses(), CPUThreads()]

    for res in resources
        if any(x->typeof(res)<:x, exclude)
            push!(final_ex.args, quote
               local $var = $res
               @testset $(name*" ($(typeof(res).name))") begin
                   @test_broken false
               end
            end)
        else
            push!(final_ex.args, quote
               local $var = $res
               @testset $(name*" ($(typeof(res).name))") $ex
            end)
        end
    end
    # preserve outer location if possible
    if ex isa Expr && ex.head === :block && !isempty(ex.args) && ex.args[1] isa LineNumberNode
        final_ex = Expr(:block, ex.args[1], final_ex)
    end
    return esc(final_ex)
end

"""
    sedate!(fit_ex)

The input is a fit expression as `fit!(mach, kws...)`. This function
throws an error if the verbosity level is set and sets the verbosity level
to -5000 otherwise.
"""
function sedate!(fit_ex)
    kwarg_exs = filter(fit_ex.args) do arg
        arg isa Expr && arg.head == :kw
    end
    keys = map(kwarg_exs) do arg
        arg.args[1]
    end
    :verbosity in keys &&
        error("You cannot specify `verbosity` in @test_mach_sequence "*
              "or @test_model_sequence. ")
    push!(fit_ex.args, Expr(:kw, :verbosity, -5000))
    return fit_ex
end

macro test_mach_sequence(fit_ex, sequence_exs...)
    sedate!(fit_ex)
    seq = gensym(:sequence)
    esc(quote
        MLJBase.flush!(MLJBase.MACHINE_CHANNEL)
        $fit_ex
        local $seq = MLJBase.flush!(MLJBase.MACHINE_CHANNEL)
        # for s in $seq
        #     println(s)
        # end
        @test $seq in [$(sequence_exs...)]
    end)
end

# function weakly_in(object:Tuple{Symbol,Model}, itr)
#     for tup in itr
#         tup[1] === object[1] && tup[2] == tup

macro test_model_sequence(fit_ex, sequence_exs...)
    sedate!(fit_ex)
    seq = gensym(:sequence)
    esc(quote
        MLJBase.flush!(MLJBase.MACHINE_CHANNEL)
        $fit_ex
        local $seq = map(MLJBase.flush!(MLJBase.MACHINE_CHANNEL)) do tup
            (tup[1], tup[2].model)
        end
        # for s in $seq
        #     println(s)
        # end
        @test $seq in [$(sequence_exs...)]
    end)
end

###############################################################################
#####    THE FOLLOWINGS ARE USED TO TEST SERIALIZATION CAPACITIES         #####
###############################################################################

function test_args(mach)
    # Check source nodes are empty if any
    for arg in mach.args
        if arg isa Source
            @test arg == source()
        end
    end
end

function test_data(mach)
    @test !isdefined(mach, :old_rows)
    @test !isdefined(mach, :data)
    @test !isdefined(mach, :resampled_data)
    @test !isdefined(mach, :cache)
end

function generic_tests(mach₁, mach₂)
    test_args(mach₂)
    test_data(mach₂)
    @test mach₂.state == -1
    for field in (:frozen, :model, :old_model, :old_upstream_state, :fit_okay)
        @test getfield(mach₁, field) == getfield(mach₂, field)
    end
end

end


text:
if !@isdefined(AbstractModel)
    include(Pkg.dir("DifferentialDynamicProgramming","src","interfaces.jl"))
end

"""
    forward_kalman(y, C, R1, R2, P0)

#Arguments:
- `y`: DESCRIPTION
- `C`: DESCRIPTION
- `R1`: State noise
- `R2`: Meas noise
- `P0`: Initial state cov
"""
function forward_kalman(y,C,R1,R2, P0)
    na,n  = size(C)
    ma  = (n-na^2) ÷ na
    sa  = na+ma
    T   = size(y,2)
    xkk = zeros(n,T);
    Pkk = zeros(n,n,T)
    for i = 0:na-1 # TODO: This should be an input, maybe only initialize with first n_param datapoints?
        ran = (i*sa+1):((i+1)*sa)
        data_to_use = 1:min(10n, size(y,2))
        xkk[ran,1]    = C[i+1,ran,data_to_use]'\y[i+1,data_to_use]  # Initialize to semi-global ls solution
    end
    R2d = MvNormal(R2)
    Pkk[:,:,1] .= P0
    xk         = copy(xkk)
    Pk         = copy(Pkk)
    i          = 1
    Ck         = C[:,:,i]
    e          = y[:,i]-Ck*xk[:,i]
    ll         = logpdf(R2d,e)
    S          = Ck*Pk[:,:,i]*Ck' + R2
    K          = (Pk[:,:,i]*Ck')/S
    xkk[:,i]   = xk[:,i] + K*e
    Pkk[:,:,i] = (I - K*Ck)*Pk[:,:,i]
    @views for i = 2:T
        Ak         = 1 # This just assumes a random walk, no additional dynamics
        Ck         = C[:,:,i]
        xk[:,i]    = Ak*xkk[:,i-1]
        Pk[:,:,i]  = Ak*Pkk[:,:,i-1]*Ak' + R1
        e          = y[:,i]-Ck*xk[:,i]
        ll        += logpdf(R2d,e)
        S          = Ck*Pk[:,:,i]*Ck' + R2
        K          = (Pk[:,:,i]*Ck')/S
        xkk[:,i]   = xk[:,i] + K*e
        Pkk[:,:,i] = (I - K*Ck)*Pk[:,:,i]
    end

    return xkk,xk,Pkk,Pk,ll
end

function kalman_smoother(y, C, R1, R2, P0)
    T                = size(y,2)
    xkk,xk,Pkk,Pk,ll = forward_kalman(y,C,R1,R2, P0)
    xkn              = similar(xkk)
    Pkn              = similar(Pkk)
    xkn[:,end]       = xkk[:,end]
    Pkn[:,:,end]     = Pkk[:,:,end]
    @views for i = T-1:-1:1
        Ck          = Pkk[:,:,i]/Pk[:,:,i+1]
        xkn[:,i]    = xkk[:,i] + Ck*(xkn[:,i+1] - xk[:,i+1])
        Pkn[:,:,i]  = Pkk[:,:,i] + Ck*(Pkn[:,:,i+1] - Pk[:,:,i+1])*Ck'
    end
    return xkn, Pkn,ll
end

# TODO: Enable imposing of known structure with e.g. a boolean matrix and a coefficient matrix to tell the algorithm which entries are known to be ==1, ==0, ==h etc.
# Use this matrix to either set some values in C, xkn, Pkn, At,Bt to zero
eye(n) = Matrix{Float64}(I,n,n)

function KalmanModel(model::KalmanModel, d::InputOutputStateData, R1,R2, P0=100R1; extend=false, D=1)::KalmanModel
    T = length(d)-1
    n = nstates(d)
    m = ninputs(d)
    xi = oftype(Matrix, state(d))
    u = oftype(Matrix, input(d))

    x,u,xnew = xi[:,1:end-1],u[:,1:end-1],xi[:,2:end]
    @assert T > n "The calling convention for x and u is that time is the second dimention"
    Ta  = extend ? T+1 : T
    N   = n^2+n*m
    y   = copy(xnew)
    C   = zeros(n,D*N,T)
    if D == 2
        R1 = kron([1/4 1/2; 1/2 1]+0.001I, R1)
        P0 = kron([1/4 1/2; 1/2 1]+0.001I, P0)
    end
    @views for t = 1:T
        C[:,1:N,t] = kron(eye(n),[x[:,t]; u[:,t]]')
    end
    xkn, Pkn,ll    = kalman_smoother(y, C, R1, R2, P0)
    model.ll = ll
    @views for t = 1:T
        ABt      = reshape(xkn[1:N,t],n+m,n)'
        model.At[:,:,t] .= ABt[:,1:n]
        model.Bt[:,:,t] .= ABt[:,n+1:end]
    end
    @views if extend # Extend model one extra time step (primitive way)
        model.At[:,:,end] .= model.At[:,:,end-1]
        model.Bt[:,:,end] .= model.Bt[:,:,end-1]
        Pkn = cat(Pkn, Pkn[:,:,end], dims=3)
    end
    model.extended = extend
    model.Pt = Pkn

    return model
end

function KalmanModel(model::KalmanModel, prior::KalmanModel, d::AbstractIdData,args...; kwargs...)::KalmanModel
    model = KalmanModel(model, d,args...; kwargs...) # Fit model in the standard way without prior
    T = length(d)
    n = nstates(d)
    m = ninputs(d)
    # @views model2statevec(model,t) = [vec(model.At[:,:,t]);  vec(model.Bt[:,:,t])][:]
    # @views model2statevec(model,t) = [model.At[:,:,t]  model.Bt[:,:,t]] |> vec
    @views for t = 1:T     # Incorporate prior
        Pkk               = model.Pt[:,:,t]
        K̄                 = Pkk/(Pkk + prior.Pt[:,:,t])
        model.Pt[:,:,t] .-= K̄*Pkk
        x                 = model2statevec(model,t)
        xp                = model2statevec(prior,t)
        x               .+= K̄*(xp-x)
        ABt               = reshape(x,n+m,n)'
        model.At[:,:,t]  .= ABt[:,1:n]
        model.Bt[:,:,t]  .= ABt[:,n+1:end]

    end
    model
end

# function modelfit(model,d)
#     x = oftype(Matrix, state(d))
#     yhat = predict(model, d)
#     fit = nrmse(x[:,2:end],yhat)
# end

# """
#     KalmanModel(xi, R1, R2, P0=100R1; extend=false, D=1, printfit=true)::KalmanModel
#
# Estimate a Kalman model without control input
#
# #Arguments:
# - `xi`: states or iddata containing states
# - `R1`: Parameter transition noise
# - `R2`: State transition noise
# - `P0`: Initial state cov
# - `extend`: add one data point at the end to have length of parameter vector match length of input.
# - `D`: Order of integration of the noise.
# """
# function KalmanModel(model::KalmanModel, xi::AnyInput,R1,R2, P0=100R1; extend=false, D=1)::KalmanModel
#     xi = oftype(Matrix, state(xi))
#     x,xnew = xi[:,1:end-1],xi[:,2:end]
#     T = length(d)
#     n = nstates(d)
#     m = ninputs(d)
#     @assert T > n "The calling convention for x and u is that time is the second dimention"
#     Ta  = extend ? T+1 : T
#     N   = n^2
#     y   = copy(xnew)
#     C   = zeros(n,D*N,T)
#     if D == 2
#         R1 = kron([1/4 1/2; 1/2 1]+0.001I, R1)
#         P0 = kron([1/4 1/2; 1/2 1]+0.001I, P0)
#     end
#     @views for t = 1:T
#         C[:,1:N,t] = kron(eye(n),x[:,t]')
#     end
#     xkn, Pkn,ll    = kalman_smoother(y, C, R1, R2, P0)
#     model.ll = ll
#     @views for t = 1:T
#         model.At[:,:,t] .= reshape(xkn[1:N,t],n,n)'
#     end
#     @views if extend # Extend model one extra time step (primitive way)
#         model.At[:,:,end] .= model.At[:,:,end-1]
#         Pkn = cat(Pkn, Pkn[:,:,end], dims=3)
#     end
#     model.extended = extend
#     model.Pt = Pkn
#     if printfit
#         yhat = predict(model, x)
#         fit = nrmse(xnew,yhat)
#         println("Modelfit: ", round.(fit,digits=3))
#     end
#
#     return model
# end

function LTVAutoRegressive(model::LTVAutoRegressive, yi,R1,R2, P0=100R1; extend=false, D=1)
    yi = oftype(Matrix, output(yi))
    y,ynew = yi[1:end-1],yi[2:end]
    T = length(y)
    n = size(R1,1)
    ϕ   = zeros(1,D*n,T)
    T  = extend ? T : T-1
    if D == 2
        R1 = kron([1/4 1/2; 1/2 1]+0.001I, R1)
        P0 = kron([1/4 1/2; 1/2 1]+0.001I, P0)
    end
    inds = 1:n
    y0 = [zeros(n-1);y]
    size(y0)
    @views for t = 1:T
        ϕ[1,1:n,t] = y0[inds .+ (t-1)]
    end
    xkn, Pkn,ll = kalman_smoother(ynew', ϕ, R1, R2, P0)
    model.ll = ll
    @views for t = 1:T
        model.θ[:,t] .= xkn[1:n,t]
    end
    @views if extend # Extend model one extra time step (primitive way)
        model.θ[:,end] .= model.θ[:,end-1]
        Pkn = cat(Pkn, Pkn[:,:,end], dims=3)
    end
    model.extended = extend
    model.Pt = Pkn

    return model
end

function predict(model::LTVAutoRegressive, y)
    y = output(y)
    T  = length(y)
    n = size(model.θ, 1)
    @assert T<=length(model) "Can not predict further than the number of time steps in the model"
    y0 = [zeros(n-1);y]
    ynew = zeros(eltype(y), T)
    inds = 1:n
    @views for t = 1:T
        ynew[t] = model.θ[:,t]'y0[inds .+ (t-1)]
    end
    ynew
end

function rootspectrogram(model, fs)
    roots = map(1:length(model)) do i
        sys = tf(1,[1;-reverse(model.θ[:,i])],1)
        (sort(pole(sys), by=imag, rev=true)[1:end÷2])
    end
    S = reduce(hcat,roots)
    fs/(2pi) .* angle.(S)'
end


text:
using PCLCommon
using PCLIO
using Cxx
using Base.Test

import PCLCommon: StdVector

milk_cartoon_path = joinpath(dirname(@__FILE__), "data",
    "milk_cartoon_all_small_clorox.pcd")

@testset "PoindCloud with PCD File IO" begin
    milk_cloud = PointCloud{PointXYZRGBA}(milk_cartoon_path)
    @test (@cxx (milk_cloud.handle)->get()) != C_NULL

    @test width(milk_cloud) == 640
    @test height(milk_cloud) == 480
    @test !is_dense(milk_cloud)

    ps = points(milk_cloud)
    @test length(ps) == 640 * 480
    @test isa(ps, StdVector)
    @test isa(ps[1], PointXYZRGBARef)
    @test isa(ps[1], PointXYZRGBAValOrRef)

    cloud = PointCloud{PointXYZ}(2,3)
    @test width(cloud) == 2
    @test height(cloud) == 3

    @test typeof(PointCloud{eltype(cloud)}()) == typeof(cloud)
    @test typeof(similar(cloud)) == typeof(cloud)
end

@testset "PointCloud conversions" begin
    milk_cloud = PointCloud{PointXYZRGBA}(milk_cartoon_path)
    non_nan_idx = 100000
    @assert non_nan_idx < length(milk_cloud)

    milk_cloud_xyzrgb = convert(PointCloud{PointXYZRGB}, milk_cloud)
    @test isa(milk_cloud_xyzrgb, typeof(PointCloud{PointXYZRGB}()))
    for (s, T) in [(:x, Float32), (:y, Float32), (:z, Float32),
            (:r, UInt8), (:g, UInt8), (:b, UInt8)]
        @test T(milk_cloud_xyzrgb[non_nan_idx,s]) == T(milk_cloud[non_nan_idx,s])
    end

    milk_cloud_xyz = convert(PointCloud{PointXYZ}, milk_cloud)
    @test isa(milk_cloud_xyz, typeof(PointCloud{PointXYZ}()))
    for (s, T) in [(:x, Float32), (:y, Float32), (:z, Float32)]
        @test T(milk_cloud_xyzrgb[non_nan_idx,s]) == T(milk_cloud[non_nan_idx,s])
    end
end

@testset "loadPCDFile" begin
    cloudxyz = PointCloud{PointXYZ}()
    @test loadPCDFile(milk_cartoon_path, cloudxyz) == 0
    cloud_xyzrgb = PointCloud{PointXYZRGB}()
    @test loadPCDFile(milk_cartoon_path, cloud_xyzrgb) == 0
end


text: # Place any deprecations prior to a 2.0 release here


text:
"""
    OlsSureFn(Y,X,m=0)

LS of Y on X; for one n dependent variables, same regressors

# Usage
(b,res,Yhat,Covb,R2) = OlsSureFn(Y,X,m)

# Input
- `Y::Array`:     Txn, the n dependent variables
- `X::Array`:     Txk matrix of regressors (including deterministic ones)
- `m::Int`:       scalar, bandwidth in Newey-West

# Output
- `b::Array`:     n*kx1, regression coefficients
- `u::Array`:     Txn, residuals Y - Yhat
- `Yhat::Array`:  Txn, fitted values X*b
- `V::Array`:     matrix, covariance matrix of vec(b)
- `R2::Number`:  n vector, R2 value

"""
function OlsSureFn(Y,X,m=0)

    (T,n) = (size(Y,1),size(Y,2))
    k     = size(X,2)

    b     = X\Y
    Yhat  = X*b
    u     = Y - Yhat

    g     = zeros(T,n*k)
    for i = 1:n
      vv      = (1+(i-1)*k):(i*k)   #1:k,(1+k):2k,...
      g[:,vv] = X.*u[:,i]           #moment conditions for Y[:,i] regression
    end

    S0    = NWFn(g,m)            #Newey-West covariance matrix
    Sxx   = X'X/T
    Sxx_1 = kron(I(n),inv(Sxx))
    V     = Sxx_1 * S0 * Sxx_1/T

    R2   = 1 .- var(u,dims=1)./var(Y,dims=1)

    return b, u, Yhat, V, R2

end

max_issues_repo_name: PaulSoderlind/EmpiricalFinancePhD

max_forks_repo_name: PaulSoderlind/EmpiricalFinancePhD


text:
using OpenCV

const cv = OpenCV

# chess1.png is at https://raw.githubusercontent.com/opencv/opencv_extra/master/testdata/cv/cameracalibration/chess1.png
img = cv.imread("chess1.png",cv.IMREAD_GRAYSCALE)
climg = cv.cvtColor(img, cv.COLOR_GRAY2BGR)

# Find the chess board corners
ret, corners = cv.findChessboardCorners(img, cv.Size{Int32}(7,5))

# If found, add object points, image points (after refining them)
if ret
    climg = cv.drawChessboardCorners(climg, cv.Size{Int32}(7,5), corners,ret)
    cv.imshow("img",climg)
    cv.waitKey(Int32(0))

    cv.destroyAllWindows()
end


text:
module Behavior

export i_dt, i_dt_wrapper, identify_fixatedTargets, plot_eye_events

using DataFrames
using MAT
using JLD
using Dates
using LinearAlgebra
using DSP
using Plots
using Statistics
using KernelDensity
using StatsBase
using Images

using ..MyUtils

### functions

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

function i_dt_wrapper(
    df,
    sf,
    fix_duration_threshold,
    fix_dispersion_threshold,
    saccade_velocity_threshold,
    saccade_duration_threshold;
    plot_saccades=false,
    plot_fixations=false,
    which_trials="all",
    how_ = "mean",
    which="Fixations_clean"
    )

    ####
    # functions takes a dataframe containg trial data of an experimental session as input
    # and returns the dataframe with filled :Fixsation and :Saccades columns
    ####

    for trl in range(1, size(df, 1), step=1)

        xy_, fix_idx, sac_idx, blink_idx = i_dt( df[!, "EyeXY"][trl], sf, fix_duration_threshold, fix_dispersion_threshold, saccade_velocity_threshold, saccade_duration_threshold);#, blink_velocity_threshold );
        df.Fixations[trl] = fix_idx
        df.Saccades[trl] = sac_idx

    end

end

function plot_eye_events(
    df;
    what="both", # fixations / saccades
    how_="mean",
    which_fixations="all_fixations",
    p=plot(legend=false))

    if what == "both" || what == "saccades"
        ### plot saccades
        for (itrl,trls_saccades) in enumerate(df[:Saccades])
            for sac_idxs in trls_saccades
                xx = df[itrl,:EyeXY][1][sac_idxs]
                yy = df[itrl,:EyeXY][2][sac_idxs]
                plot!(xx,yy)
            end
        end
    end

    if which_fixations == "all_fixations"
        clean_fixations_bool = [[fix in [-1,0,1,2,3,4,5,6,7,8] for fix in el] for el in df.LookedAtTarget]
    elseif which_fixations == "clean_fixations"
        clean_fixations_bool = [[fix in [0,1,2,3,4,5,6,7,8] for fix in el] for el in df.LookedAtTarget]
    end

    if what == "both" || what == "fixations"
        for (itrl,trls_fix) in enumerate(df.Fixations)
            for (ifix, fix) in enumerate(trls_fix)
                if clean_fixations_bool[itrl][ifix]
                    scatter!( [mean( df.EyeXY[itrl][1][fix] )], [mean( df.EyeXY[itrl][2][fix] )] )
                end
            end
        end

        # for (itrl,trls_fixations) in enumerate(df[:Fixations])
        #     for fix_idxs in trls_fixations
        #         # print(df[itrl,:EyeXY][1][sacd])
        #         xx = df[itrl,:EyeXY][1][fix_idxs]
        #         yy = df[itrl,:EyeXY][2][fix_idxs]
        #         if how_ == "mean"
        #             scatter!([mean(xx)],[mean(yy)])
        #         elseif how_ == "raw"
        #             scatter!(xx,yy, markersize=0.5)
        #         end
        #     end
        # end
    end

    xlims!(-15,15)
    ylims!(-15,15)
    return p
end

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

function identify_fixatedTargets(
    df;
    nbins_length=150,
    nbins_angles=150,
    kde_bandwidth_length=1/3,
    kde_bandwidth_angle=1/15,
    sf=2000,
    plot_distributions=false,
    plot_only=false,
    target_rotation=22.5)

    #
    # also identifies Invalid trials, i.e. those that are inititated and Fritz went for the WRONG target
    #
    #

    if !isempty(df)

        ### rotate targets to prevent splitting a cluster between 0 and 360 degrees
        rot_angl = target_rotation
        rot_rad = 2*pi/360 * rot_angl ### rotate the data points by 22.5 degrees to get target3 as one cluster
        rot_matr = [cos(rot_rad) -sin(rot_rad); sin(rot_rad) cos(rot_rad) ]
        ref_dist = (2pi)/360 * 45 ### theoretical distance of targets in rad

        ### determin theoretical angles of targets
        angle_of_first_target = (2*pi)/360 * rot_angl
        angle_of_last_target = 7*(2*pi)/360*45 + (2*pi)/360*rot_angl
        theoretical_x_vals_targets_va = collect(range(angle_of_first_target, angle_of_last_target; step=ref_dist))

        ### get vector length and angles of fixations
        vec_length = Array{Any}(undef, 0)
        vec_angl = Array{Any}(undef, 0)

        for (itrl,trls_fix) in enumerate(df.Fixations)
            push!(vec_length, Array{Float64,1}(undef,0))
            push!(vec_angl, Array{Float64,1}(undef,0))

            for (ifix,fix) in enumerate(trls_fix)

                # vec = [mean(xy_nb[itrl][1][fix]), mean(xy_nb[itrl][2][fix])]
                vec = [mean(df.EyeXY[itrl][1][fix]), mean(df.EyeXY[itrl][2][fix])]

                ### rotate fixations to prevent the cluster at 0deg to be splitted between >0 and <360
                vec_rotated = rot_matr*vec

                ### length
                push!(vec_length[itrl], round(norm(vec_rotated), digits=1))

                ### angle relative to x axis
                if vec_rotated[2] >= 0
                    push!(vec_angl[itrl], acos(dot([1,0],vec_rotated)/(norm([1,0])*norm(vec_rotated))))
                else
                    push!(vec_angl[itrl], acos(dot([1,0],vec_rotated)/(norm([1,0])*norm(vec_rotated)))+pi)
                end

            end
        end

        ### fit histograms and kernel densitys to length
        vec_length_ = MyUtils.flatten(vec_length)
        h_vl = fit(Histogram, vec_length_, nbins=nbins_length)
        h_norm_vl = normalize(h_vl, mode=:pdf)
        vec_len_dens = kde(vec_length_, bandwidth=kde_bandwidth_length)
        ### find modes of fitted distributions
        x_idx_all_modes_vl = findlocalmaxima(vec_len_dens.density)
        x_idx_two_largest_modes_vl = x_idx_all_modes_vl[sortperm(vec_len_dens.density[x_idx_all_modes_vl], rev=true)][1:2]
        x_val_two_largest_modes_vl = vec_len_dens.x[x_idx_two_largest_modes_vl]
        y_val_two_largest_modes_vl = vec_len_dens.density[x_idx_two_largest_modes_vl]
        ### sort saccades according fixation vector length
        # min and max vector lengths of central fixations and target fixations
        c_fix = [x_val_two_largest_modes_vl[1]-1, x_val_two_largest_modes_vl[1]+1]
        margin_t_fix_i = 3
        margin_t_fix_o = 3
        t_fix = [x_val_two_largest_modes_vl[2]-margin_t_fix_i, x_val_two_largest_modes_vl[2]+margin_t_fix_o]

        fixation_types = Array{Any,1}(undef, 0)

        for (itrl,trls_fix) in enumerate(df.Fixations)
            clean_fixations = [] ### Fixations that fall into the defined areas of central fixation point and targets

            push!(fixation_types, Array{Float64,1}(undef,0))

            for (ifix,fix) in enumerate(trls_fix)

                ### catch central fixations
                if c_fix[1] < vec_length[itrl][ifix] < c_fix[2]
                    push!(fixation_types[itrl], 0)
                elseif t_fix[1] < vec_length[itrl][ifix] < t_fix[2]
                    push!(fixation_types[itrl], 10)
                else
                    push!(fixation_types[itrl], -1)
                end

            end
        end

        bool_list_central_fixations = [ el.==10.0 for el in fixation_types ]
        vec_angl_targets_only = [vec_angl[ii][el] for (ii,el) in enumerate(bool_list_central_fixations)]

        ### fit histograms and kernel densitys to angles distributions
        vec_angl_ = MyUtils.flatten(vec_angl_targets_only)
        h_va = fit(Histogram, vec_angl_, nbins=nbins_angles)
        h_norm_va = normalize(h_va, mode=:pdf)
        vec_angl_dens = kde(vec_angl_, bandwidth=kde_bandwidth_angle)
        ### find modes of fitted distributions
        x_idx_all_modes_va = findlocalmaxima(vec_angl_dens.density)
        x_val_all_modes_va = vec_angl_dens.x[x_idx_all_modes_va]
        which_modes = [argmin(abs.(tx .- x_val_all_modes_va)) for tx in theoretical_x_vals_targets_va] ### find modes closest to theoretical location of modes
        x_idx_eight_largest_modes_va = x_idx_all_modes_va[which_modes]
        x_val_eight_largest_modes_va = x_val_all_modes_va[which_modes]
        y_val_eight_largest_modes_va = vec_angl_dens.density[x_idx_eight_largest_modes_va]
        ### correct for potentially wrong modes: replace modes that deviate more than 0.5 from theoretical mode
        false_modes = x_val_eight_largest_modes_va .- theoretical_x_vals_targets_va .> 0.5
        x_val_eight_largest_modes_va[false_modes] .=  theoretical_x_vals_targets_va[false_modes]

        ### angular range of target fixations
        margin_ts = 0.2
        t1 = [x_val_eight_largest_modes_va[1]-margin_ts, x_val_eight_largest_modes_va[1]+margin_ts]
        t2 = [x_val_eight_largest_modes_va[2]-margin_ts, x_val_eight_largest_modes_va[2]+margin_ts]
        t3 = [x_val_eight_largest_modes_va[3]-margin_ts, x_val_eight_largest_modes_va[3]+margin_ts]
        t4 = [x_val_eight_largest_modes_va[4]-margin_ts, x_val_eight_largest_modes_va[4]+margin_ts]
        t5 = [x_val_eight_largest_modes_va[5]-margin_ts, x_val_eight_largest_modes_va[5]+margin_ts]
        t6 = [x_val_eight_largest_modes_va[6]-margin_ts, x_val_eight_largest_modes_va[6]+margin_ts]
        t7 = [x_val_eight_largest_modes_va[7]-margin_ts, x_val_eight_largest_modes_va[7]+margin_ts]
        t8 = [x_val_eight_largest_modes_va[8]-margin_ts, x_val_eight_largest_modes_va[8]+margin_ts]

        for (itrl,trls_fix) in enumerate(fixation_types)
            clean_fixations = Array{Array{Float64,1},1}(undef,0)

            for (ifix,fix) in enumerate(trls_fix)

                ### catch target fixations, could either be 0 for central fixations or -1 for noise fixations or 10 for target fixations
                if fix == 10.0

                    ### sort into traget IDs
                    if t1[1] < vec_angl[itrl][ifix] < t1[2]
                        fixation_types[itrl][ifix] = 3
                    elseif t2[1] < vec_angl[itrl][ifix] < t2[2]
                        fixation_types[itrl][ifix] = 2
                    elseif t3[1] < vec_angl[itrl][ifix] < t3[2]
                        fixation_types[itrl][ifix] = 1
                    elseif t4[1] < vec_angl[itrl][ifix] < t4[2]
                        fixation_types[itrl][ifix] = 8
                    elseif t5[1] < vec_angl[itrl][ifix] < t5[2]
                        fixation_types[itrl][ifix] = 4
                    elseif t6[1] < vec_angl[itrl][ifix] < t6[2]
                        fixation_types[itrl][ifix] = 5
                    elseif t7[1] < vec_angl[itrl][ifix] < t7[2]
                        fixation_types[itrl][ifix] = 6
                    elseif t8[1] < vec_angl[itrl][ifix] < t8[2]
                        fixation_types[itrl][ifix] = 7
                    else
                        fixation_types[itrl][ifix] = -1
                    end

                end

                ### store fixations belonging to target or central fixation dot
                if !(fixation_types[itrl][end] == -1)
                    push!( clean_fixations, df.Fixations[itrl][ifix] )
                end

            end # end iterate trials fixations

            if !plot_only
                df.LookedAtTarget[itrl] = fixation_types[itrl]

                ### add columns Values
                df.Valid[itrl] = (!df.WasMissed[itrl]) & (any(unique(df[itrl, :LookedAtTarget]) .> 0))
                df.Invalid[itrl] = (!df.WasMissed[itrl]) & (!any(unique(df[itrl, :LookedAtTarget]) .> 0))
            #     # df.Hit[itrl] = (df.Valid[itrl]) & (df[itrl, :LookedAtTarget][end] == df[itrl, :ValidTargetIndex])
                if !isempty(df[itrl, :LookedAtTarget])
                    df.FA[itrl] = ((df.Valid[itrl]) & (df[itrl, :LookedAtTarget][end] != df[itrl, :ValidTargetIndex])) & !(df.Hit[itrl])
                end

                ### Get reaction times
                targetFixation_idxs = findall([ !(el in [-1,0]) for el in df[itrl,:LookedAtTarget]]) ### list of all target fixations
                if !isempty(targetFixation_idxs)
                    firstTargetFixation_id = minimum(targetFixation_idxs) ### the first target fixation
                    firstTargetFixation_time = df[itrl,:Fixations][firstTargetFixation_id] / sf
                    RT_firstTarget = minimum(firstTargetFixation_time .- df[itrl, :LuminanceChangeTime])

                    if df[itrl,:LookedAtTarget][firstTargetFixation_id] == df[itrl,:ValidTargetIndex]
                        RT_CorrectTarget = RT_firstTarget
                    else
                        correctTargetFixation_idxs = findall(df[itrl,:LookedAtTarget] .== df[itrl,:ValidTargetIndex]) ### list of all fixations on correct target
                        if !isempty(correctTargetFixation_idxs)
                            correctTargetFixation_idx = minimum(correctTargetFixation_idxs) ### the first fixation on correct target
                            correctTargetFixation_time = df[itrl,:Fixations][correctTargetFixation_idx] / sf
                            RT_CorrectTarget = minimum(correctTargetFixation_time .- df[itrl, :LuminanceChangeTime])
                        else ### if correct target was not fixated
                            RT_CorrectTarget = NaN
                        end
                    end

                else ### if no target was fixated
                    RT_firstTarget = NaN
                    RT_CorrectTarget = NaN
                end
                df.RTfirstFixation[itrl] = RT_firstTarget
                df.RTcorrectFixation[itrl] = RT_CorrectTarget

            end

        end # end iterate trials

        if plot_distributions | plot_only
            p = plot(layout=2, size=(1000,500), xlabel="saccade length / AU")
            plot!(h_norm_vl, subplot=1, xlims=[0, 50])
            plot!(vec_len_dens.x, vec_len_dens.density, linewidth=3, subplot=1)
            scatter!(x_val_two_largest_modes_vl, y_val_two_largest_modes_vl, subplot=1)
            #
            plot!(h_norm_va, subplot=2, xlabel="saccade angles / rad")
            plot!(vec_angl_dens.x, vec_angl_dens.density, linewidth=3, subplot=2)
            scatter!(x_val_eight_largest_modes_va, y_val_eight_largest_modes_va, subplot=2)
            return p
        end

    end

end

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

function i_dt(
    xy_o::Array{Array{Float64,1},1},
    sf,
    fix_duration_threshold_sec,
    fix_dispersion_threshold,
    saccade_velocity_threshold,
    saccade_duration_threshold
    )

    ###########
    #
    # Implements I-DT algorithm to detect fixations, sacades
    # extended to identify and interpolate blinks
    #
    # INPUT
    # fix_duration_threshold_sec: minimum duration of a fixation to be counted in sec
    # fix_dispersion_threshold: can be left empty then it is inferred from central fixation period in beginning of trials
    # saccade_velocity_threshold: not normalized to sampling rate!
    # saccade_duration_threshold: in sec
    #
    # OUTPUT
    #   xy_: eye trace with blinks interpolated
    #   sac_idx:
    #
    #
    #
    ###########

    if size(xy_o[1],1) > 1000 # if number of samples is lower empty arrays are returned

        xy = deepcopy(xy_o)
        [deleteat!(xy[ii], findall(isnan.(xy[ii]))) for ii in 1:length(xy)] # remove NaNs
        xy_ = deepcopy(xy) # to be returned with blinks interpolated
        xy_idxs = [collect(1:length(xy[1])), collect(1:length(xy[2]))] # array with corresponding indices
        fix_window_length = Int(fix_duration_threshold_sec * sf)

        if fix_dispersion_threshold == []

            fix_dispersion_threshold = maximum(xy[1][1:499])-minimum(xy[1][1:499]) + maximum(xy[2][1:499])-minimum(xy[2][1:499]) # deg

        end

        #
        sac_idx = Array{Array{Int,1}, 1}(undef,0)#[]
        fix_idx = Array{Array{Int,1}, 1}(undef,0)#[]
        blink_idx = Array{Array{Int,1}, 1}(undef,0)#[]

        ### I-DT loop:
        while !isempty(xy_idxs[1])

            if length(xy_idxs[1]) > fix_window_length # check if length(xy) > fix_window_length (lower duration threshold of saccades)
                current_window_idxs = collect(1:fix_window_length)
            else
                # @show length(xy_idxs[1])
                current_window_idxs = collect(1:length(xy_idxs[1])) # if length(xy) > fix_window_length
            end

            # println(current_window_idxs)

            current_window_vals = [xy[ii][current_window_idxs] for ii in 1:2]
            wd = (maximum(current_window_vals[1])-minimum(current_window_vals[1])) + (maximum(current_window_vals[2])-minimum(current_window_vals[2]))

            ### if dispersion in current window is smaller than threshold
            if wd <= fix_dispersion_threshold

                # while this is the case and window does not exceed amount of data points add next point to window
                while wd<=fix_dispersion_threshold && length(current_window_idxs)<length(xy[1])

                    ### append to indices of current window
                    current_window_idxs = append!( current_window_idxs, maximum(current_window_idxs)+1 )
                    ### append to values of current widnow
                    current_window_vals[1] = xy[1][current_window_idxs]
                    current_window_vals[2] = xy[2][current_window_idxs]
                    wd = (maximum(current_window_vals[1])-minimum(current_window_vals[1])) + (maximum(current_window_vals[2])-minimum(current_window_vals[2]))

                end

                ### add fixation to array and delete window from data
                # push!( fix_idx, [xy_idxs[ii][current_window_idxs] for ii in 1:2] )
                push!( fix_idx, xy_idxs[1][current_window_idxs] )
                [deleteat!( xy[ii], current_window_idxs ) for ii in 1:length(xy)]
                [deleteat!( xy_idxs[ii], current_window_idxs ) for ii in 1:length(xy_idxs)]

            else # if dispersion in current window is greater than threshold delete first data point

                [deleteat!( xy[ii], 1 ) for ii in 1:length(xy)]
                [deleteat!( xy_idxs[ii], 1 ) for ii in 1:length(xy_idxs)]

            end

        end

        ### get saccade-like indices, i.e. those not identified as fixations
        for (idx,fixs) in enumerate(fix_idx[1:end-1])
            #if (fix_idx[idx+1][1][1]-fixs[1][end]) > blink_velocity_threshold # blink velocity threshold
                # push!( sac_idx, fixs[1][end]:fix_idx[idx+1][1][1] )
                push!( sac_idx, fixs[end]:fix_idx[idx+1][1] )

            #end
        end

        ### filter blinks: saccades with velocities > threshold
        saccade_velocities = [diff(xy_[1][sac]) for sac in sac_idx]
        blinks = findall(maximum.(saccade_velocities) .> saccade_velocity_threshold)

        for bl in blinks

            push!( blink_idx, sac_idx[bl] )

        end
        # delete blinks from saccades collection
        deleteat!( sac_idx, blinks )
        deleteat!(sac_idx, findall(length.(sac_idx) .< saccade_duration_threshold*sf)) # remove implausible short saccades

        ### inerpolate blinks (linear)
        interpolated_blinks = []

        for bl in 1:length(blink_idx)

            push!(interpolated_blinks, range( xy_[1][minimum(blink_idx[bl])], xy_[1][maximum(blink_idx[bl])], length=length(blink_idx[bl]) ))

        end

        [xy_[dim][blink_idx[bl]] .= interpolated_blinks[bl] for bl in 1:length(blink_idx), dim=1:2 ] # dim is x,y dimensions of eyedata -> bot needed

        return xy_, fix_idx, sac_idx, blink_idx

    else
        return [], [], [], []
    end

end # i_dt

end # end Behavior.jl

max_stars_repo_name: maalaria/Newbe

max_issues_repo_name: maalaria/Newbe

max_forks_repo_name: maalaria/Newbe


text:
immutable Poisson <: DiscreteUnivariateDistribution
    lambda::Float64
    function Poisson(l::Real)
    	l > zero(l) || error("lambda must be positive")
        new(float64(l))
    end
    Poisson() = new(1.0)
end

@_jl_dist_1p Poisson pois

isupperbounded(::Union(Poisson, Type{Poisson})) = false
islowerbounded(::Union(Poisson, Type{Poisson})) = true
isbounded(::Union(Poisson, Type{Poisson})) = false

minimum(::Union(Poisson, Type{Poisson})) = 0
maximum(::Union(Poisson, Type{Poisson})) = Inf

insupport(::Poisson, x::Real) = isinteger(x) && zero(x) <= x
insupport(::Type{Poisson}, x::Real) = isinteger(x) && zero(x) <= x

function probs(d::Poisson, rgn::UnitRange)
    λ = d.lambda
    f, l = rgn[1], rgn[end]
    0 <= f <= l || throw(BoundsError())
    r = Array(Float64, l - f + 1)
    v = r[1] = pdf(d, f)
    if l > f
        b = f - 1
        for x = f+1:l
            c = λ / x
            r[x - b] = (v *= c)
        end
    end
    return r
end

mean(d::Poisson) = d.lambda

median(d::Poisson) = quantile(d, 0.5)

mode(d::Poisson) = ifloor(d.lambda)
modes(d::Poisson) = isinteger(d.lambda) ? [int(d.lambda)-1,int(d.lambda)] : [ifloor(d.lambda)]

var(d::Poisson) = d.lambda

skewness(d::Poisson) = 1.0 / sqrt(d.lambda)

kurtosis(d::Poisson) = 1.0 / d.lambda

function entropy(d::Poisson)
    λ = d.lambda
    if λ < 50.0
        s = 0.0
        for k in 1:100
            s += λ^k * lgamma(k + 1.0) / gamma(k + 1.0)
        end
        return λ * (1.0 - log(λ)) + exp(-λ) * s
    else
        return 0.5 * log(2 * pi * e * λ) -
               (1 / (12 * λ)) -
               (1 / (24 * λ * λ)) -
               (19 / (360 * λ * λ * λ))
    end
end

function mgf(d::Poisson, t::Real)
    l = d.lambda
    return exp(l * (exp(t) - 1.0))
end

function cf(d::Poisson, t::Real)
    l = d.lambda
    return exp(l * (exp(im * t) - 1.0))
end

# model fitting

immutable PoissonStats <: SufficientStats
    sx::Float64   # (weighted) sum of x
    tw::Float64   # total sample weight
end

suffstats(::Type{Poisson}, x::Array) = PoissonStats(float64(sum(x)), float64(length(x)))

function suffstats(::Type{Poisson}, x::Array, w::Array{Float64})
    n = length(x)
    n == length(w) || throw(ArgumentError("Inconsistent array lengths."))
    sx = 0.
    tw = 0.
    for i = 1 : n
        @inbounds wi = w[i]
        @inbounds sx += x[i] * wi
        tw += wi
    end
    PoissonStats(sx, tw)
end

fit_mle(::Type{Poisson}, ss::PoissonStats) = Poisson(ss.sx / ss.tw)


text:
# sample implementations of the Abstract Types
# used for testing, not to be used for real code

mutable struct ConcreteDGMesh{Tmsh} <: AbstractDGMesh{Tmsh}
  # counts
  numVert::Int
  numEl::Int
  numNodes::Int
  numDof::Int
  numDofPerNode::Int
  numNodesPerElement::Int
  order::Int
  numNodesPerFace::Int

  # parallel counts
  npeers::Int
  numGlobalEl::Int
  numSharedEl::Int
  peer_face_counts::Array{Int, 1}
  local_element_counts::Array{Int, 1}
  remote_element_counts::Array{Int, 1}

  # MPI info
  comm::Int  # this should really be an MPI.Comm
  myrank::Int
  commsize::Int
  peer_parts::Array{Int, 1}

  # Discretization type
  isDG::Bool
  isInterpolated::Bool

  # Mesh data
  coords::Array{Tmsh, 3}
  dxidx::Array{Tmsh, 4}
  jac::Array{Tmsh, 2}

  # interpolated data
  coords_bndry::Array{Tmsh, 3}
  dxidx_bndry::Array{Tmsh, 4}
  jac_bndry::Array{Tmsh, 2}
  dxidx_face::Array{Tmsh, 4}
  jac_face::Array{Tmsh, 2}

  # parallel data
  coords_sharedface::Array{Array{Tmsh, 3}, 1}
  dxidx_sharedface::Array{Array{Tmsh, 4}, 1}
  jac_sharedface::Array{Array{Tmsh, 2}, 1}

  # boundary condition data
  numBC::Int
  numBoundaryEdges::Int
  bndryfaces::Array{Boundary, 1}
  bndry_offsets::Array{Int, 1}
  bndry_funcs::Array{BCType, 1}

  # interior edge data
  numInterfaces::Int
  interfaces::Array{Interface, 1}

  # dof number data
  dofs::Array{Int, 3}
  dof_offset::Int
  sparsity_bnds::Array{Int, 2}
  sparsity_nodebnds::Array{Int, 2}

  # mesh coloring data
  numColors::Int
  maxColors::Int
  color_masks::Array{BitArray{1}, 1}
  shared_element_colormasks::Array{Array{BitArray{1}, 1}, 1}
  pertNeighborEls::Array{Int, 2}

  # parallel bookkeeping
  bndries_local::Array{Array{Boundary, 1}, 1}
  bndries_remote::Array{Array{Boundary, 1}, 1}
  shared_interfaces::Array{Array{Interface, 1}, 1}
  shared_element_offsets::Array{Int, 1}
  local_element_lists::Array{Array{Int, 1}, 1}

  function ConcreteDGMesh{Tmsh}(opts) where {Tmsh}

    mesh = new{Tmsh}()  # incomplete initialization

    # make up values
    mesh.numVert = 5
    mesh.numEl = 7
    mesh.numNodesPerElement = 4
    mesh.numNodes = mesh.numEl*mesh.numNodesPerElement
    mesh.numDofPerNode = 5
    mesh.numDof = mesh.numNodes*mesh.numDofPerNode
    mesh.order = 2
    mesh.numNodesPerFace = 3

    mesh.npeers = 2
    mesh.numGlobalEl = 42
    mesh.numSharedEl = 4
    mesh.peer_face_counts = [2, 2]
    mesh.local_element_counts = [2, 2]
    mesh.remote_element_counts = [1, 2]

    mesh.comm = 666
    mesh.myrank = 0
    mesh.commsize = 4
    mesh.peer_parts = [1, 3]

    mesh.isDG = true
    mesh.isInterpolated = true

    mesh.coords = rand(2, mesh.numNodesPerElement, mesh.numEl)
    mesh.dxidx = rand(2, 2, mesh.numNodesPerElement, mesh.numEl)
    mesh.jac = rand(mesh.numNodesPerElement, mesh.numEl)

    mesh.numBoundaryEdges = 2
    mesh.coords_bndry = rand(2, mesh.numNodesPerFace, mesh.numBoundaryEdges)
    mesh.dxidx_bndry = rand(2, 2, mesh.numNodesPerFace, mesh.numBoundaryEdges)
    mesh.jac_bndry = rand(mesh.numNodesPerFace, mesh.numBoundaryEdges)

    mesh.numInterfaces = 2
    mesh.dxidx_face = rand(2, 2, mesh.numNodesPerFace, mesh.numInterfaces)
    mesh.jac_face = rand(mesh.numNodesPerFace, mesh.numInterfaces)

    mesh.numBC = 2
    mesh.bndryfaces = [Boundary(2, 1), Boundary(3, 2)]
    mesh.bndry_offsets = [1, 2, 3]
    mesh.bndry_funcs = Array{BCType}(mesh.numBC)

    mesh.interfaces = [Interface(2, 3, 2, 3, 0), Interface(4, 7, 1, 3, 2)]

    mesh.dofs = Array{Int}(mesh.numDofPerNode, mesh.numNodesPerElement, mesh.numEl)
    for i=1:mesh.numDof
      mesh.dofs[i] = i
    end

    mesh.dof_offset = 0
    mesh.sparsity_bnds = Array{Int}(2, mesh.numDof)
    # make the matrix dense
    for i=1:mesh.numDof
      mesh.sparsity_bnds[1, i] = 1
      mesh.sparsity_bnds[2, i] = mesh.numDof
    end

    mesh.sparsity_nodebnds = Array{Int}(2, mesh.numNodes)
    for i=1:mesh.numNodes
      mesh.sparsity_nodebnds[1, i] = 1
      mesh.sparsity_nodebnds[2, i] = mesh.numNodes
    end

    mesh.numColors = 4  # like that's really possible...
    mesh.maxColors = 5
    mesh.color_masks = Array{BitArray{1}}(mesh.numColors)
    mesh.shared_element_colormasks = Array{Array{BitArray{1}, 1}}(mesh.npeers)
    mesh.pertNeighborEls = Array{Int}(4, mesh.numEl)

    mesh.bndries_local = Array{Array{Boundary, 1}}(mesh.npeers)
    mesh.bndries_remote = Array{Array{Boundary, 1}}(mesh.npeers)
    mesh.shared_interfaces = Array{Array{Interface, 1}}(mesh.npeers)
    mesh.shared_element_offsets = Array{Int}(2)
    mesh.local_element_lists = Array{Array{Int, 1}}(mesh.npeers)

    return mesh
  end

end

mutable struct ConcreteParamType{Tdim}
  t::Float64
  order::Int
  time  # this really should have a concrete type

  function ConcreteParamType{Tdim}(mesh, opts) where Tdim

    t = 0.0
    order = mesh.order
    time = 0

    return new(t, order, time)
  end
end

mutable struct ConcreteSolutionData{Tsol, Tres} <: AbstractSolutionData{Tsol, Tres}
  q::Array{Tsol, 3}
  q_vec::Array{Tsol, 1}
  # AbstractSharedFaceData should be a concrete type
  shared_data::Array{AbstractSharedFaceData{Tsol}, 1}
  res::Array{Tres, 3}
  res_vec::Array{Tres, 1}
  M::Array{Float64, 1}
  Minv::Array{Float64, 1}
  disassembleSolution::Function
  assembleSolution::Function
  multiplyA0inv::Function
  majorIterationCallback::Function
  params::ConcreteParamType{2}

  function ConcreteSolutionData{Tsol, Tres}(mesh, sbp, opts)  where {Tsol, Tres}

    numDofPerNode = mesh.numDofPerNode
    numNodesPerElement = mesh.numNodesPerElement
    numEl = mesh.numEl
    numNodesPerFace = mesh.numNodesPerFace
    numDof = mesh.numDof

    q = Array{Tsol}(numDofPerNode, numNodesPerElement, numEl)
    q_vec = Array{Tsol}(numDof)
    # the elements of this array should be populated as well
    shared_data = Array{AbstractSharedFaceData{Tsol}}(mesh.npeers)
    res = Array{Tres}(mesh.numDofPerNode, mesh.numNodesPerElement, mesh.numEl)
    res_vec = Array{Tres}(mesh.numDof)
    M = ones(numDof)
    Minv = ones(numDof)

    # create empty functions
    disassembleSolution = (mesh, sbp, eqn, opts) -> begin end
    assembleSolution = (mesh, sbp, eqn, opts) -> begin end
    multiplyA0inv = (mesh, sbp, eqn, opts, res_arr) -> begin end
    majorIterationCallback = (itr, mesh, sp, eqn, opts) -> begin end

    params = ConcreteParamType{2}(mesh, opts)

    return new(q, q_vec, shared_data, res, res_vec, M, Minv,
               disassembleSolution, assembleSolution, multiplyA0inv,
               majorIterationCallback, params)
  end
end

mutable struct ConcreteSBP
  Q::Array{Float64, 2}

  function ConcreteSBP()
    Q = rand(3, 3)
    return new(Q)
  end
end


text:
__precompile__()

module Polyhedra

using FixedSizeArrays, GeometryTypes

export PolyhedraLibrary, Polyhedron, getlibrary, getlibraryfor

abstract PolyhedraLibrary
abstract Polyhedron{N,T} <: GeometryPrimitive{N,T}

# Libraries can implement this function to be the default for a certain type T
# getlibraryfor{T<:Real}(::Type{T})
# This function should be implemented by libraries
# or it would be like not implementing similar for AbstractArray
getlibraryfor{T}(p::Polyhedron, ::Type{T}) = getlibraryfor(T)
getlibrary(p::Polyhedron) = getlibraryfor(p, eltype(p))

import Base.intersect, Base.+, Base.*, Base.isempty, Base.copy, Base.push!, Base.length, Base.eltype, Base.start, Base.done, Base.next

include("mycomp.jl")
include("elements.jl")
include("representation.jl")
include("repop.jl")
include("operations.jl")
include("opt.jl")
include("simplerep.jl")
include("simplepolyhedron.jl")
include("decompose.jl")

end # module


text:
@testset "basic" begin
    s = "ABC"
    parts = FP.md_partition(s)
    @test length(parts) == 1
    @test parts[1].name == :TEXT
    @test parts[1].ss == s
    s = "ABC<!--DEF-->"
    parts = FP.md_partition(s)
    @test length(parts) == 2
    @test parts[1].name == :TEXT
    @test parts[2].name == :COMMENT
    s = "ABC<!--DEF-->GHI"
    parts = FP.md_partition(s)
    @test length(parts) == 3
    @test all(parts[i].name == :TEXT for i in (1, 3))
    s = "ABC<!--DEF-->GHI<!--JKL-->MNO"
    parts = FP.md_partition(s)
    @test all(parts[i].name == :COMMENT for i in (2, 4))
end

@testset "cov" begin
    s = "abc"
    p = FP.md_partition(s)
    @test length(p) == 1
    @test p[1].name == :TEXT

    p = "abc @@d ef *AA* g @@ end" |> FP.md_partition
    @test p[2] isa FP.Block
    c = p[2] |> FP.md_partition
    @test c[1] // "ef"
    @test c[2] // "*AA*"
    @test c[3] // "g"
end

@testset "group" begin
    g = """
        abc *def* ghi
        | a | b | c |""" |> grouper
    @test g[1].role == :PARAGRAPH
    @test g[2].role == :TABLE
    g = """
        abc *def* ghi
        @@d klm @@""" |> grouper
    @test g[1].role == :PARAGRAPH
    @test g[2].role == :DIV
    g = """
        @@d klm @@
        abc *def* ghi""" |> grouper
    @test g[2].role == :PARAGRAPH
    @test g[1].role == :DIV
end

@testset "environment" begin
    g = raw"abc \begin{foo} bar @@d baz @@ and \end{foo} hello" |> grouper
    @test g[1] // "abc"
    @test g[2] // raw"\begin{foo} bar @@d baz @@ and \end{foo}"
    @test g[3] // "hello"

    g = raw"abc \begin{foo} bar @@d baz @@ \begin{foo} a \end{foo} and \end{foo} hello" |> grouper
    @test g[1] // "abc"
    @test g[3] // "hello"

    @test_throws FP.FranklinParserException raw"abc \begin{foo}" |> grouper
end

@testset "raw" begin
    s = """
        abc ~~~ def ~~~ ghi %%% klm %%% and ??? %%% foo ???.
        """
    g = s |> grouper
    @test g[1] // s
    b = g[1].blocks
    @test b[1] // "abc"
    @test b[2] // "~~~ def ~~~"
    @test b[3] // "ghi"
    @test b[4] // "%%% klm %%%"
    @test b[5] // "and"
    @test b[6] // "??? %%% foo ???"
    @test b[7] // "."
end

@testset "nop-par" begin
    s = """

        {{foo}}

        """
    g = s |> grouper
    @test g[1] // s
    @test g[1].role == :PARAGRAPH_NOP
    s = """
        ABC ~~~foo~~~

        ~~~foo~~~

        ABC <!--bar-->

        <!--bar-->
        """
    g = s |> grouper
    @test g[1].role == :PARAGRAPH
    @test g[2].role == :PARAGRAPH_NOP
    @test g[3].role == :PARAGRAPH
    @test g[4].role == :PARAGRAPH_NOP

    # from franklin docs

    s = raw"""
        +++
        a = 5
        +++

        \newcommand{\foo}{bar}

        ~~~
        baz
        ~~~

        ABC
        """
    g = s |> grouper
    @test g[1].role == :MD_DEF_BLOCK
    @test g[2].role == :PARAGRAPH_NOP
    @test g[3].role == :PARAGRAPH_NOP
    @test g[4].role == :PARAGRAPH
    @test g[4] // "ABC"

    s = raw"""
        +++
        a = 5
        +++

        \newcommand{\foo}{bar}
        ~~~
        baz
        ~~~

        ABC
        """
    g = s |> grouper
    @test g[1].role == :MD_DEF_BLOCK
    @test g[2].role == :PARAGRAPH_NOP
    @test g[3].role == :PARAGRAPH
    @test g[3] // "ABC"
end


text:
using ForwardDiff, OrdinaryDiffEq

include("utils.jl")

struct TestEuclidean{N} <: Manifold end
struct TestEuclideanMetric <: Metric end

@testset "scaled Euclidean metric" begin
    Manifolds.manifold_dimension(::TestEuclidean{N}) where {N} = N
    function Manifolds.local_metric(M::MetricManifold{<:TestEuclidean,<:TestEuclideanMetric}, x)
        return Diagonal(1.0:manifold_dimension(M))
    end

    n = 3
    E = TestEuclidean{n}()
    g = TestEuclideanMetric()
    M = MetricManifold(E, g)
    G = Diagonal(1.0:n)
    invG = inv(G)
    @test manifold_dimension(M) == n
    @test base_manifold(M) === E
    @test metric(M) === g

    for vtype in (Vector, SVector{n}, MVector{n})
        x, v, w = vtype(randn(n)), vtype(randn(n)), vtype(randn(n))

        @test local_metric(M, x) ≈ G
        @test inverse_local_metric(M, x) ≈ invG
        @test det_local_metric(M, x) ≈ *(1.0:n...)
        @test log_local_metric_density(M, x) ≈ sum(log.(1.0:n)) / 2
        @test inner(M, x, v, w) ≈ dot(v, G * w)
        @test norm(M, x, v) ≈ sqrt(dot(v, G * v))

        T = 0:.5:10
        @test exp(M, x, v, T) ≈ [x + t * v for t in T]
        @test geodesic(M, x, v, T) ≈ [x + t * v for t in T]

        @test christoffel_symbols_first(M, x) ≈ zeros(n, n, n)
        @test christoffel_symbols_second(M, x) ≈ zeros(n, n, n)
        @test riemann_tensor(M, x) ≈ zeros(n, n, n, n)
        @test ricci_tensor(M, x) ≈ zeros(n, n)
        @test ricci_curvature(M, x) ≈ 0
        @test gaussian_curvature(M, x) ≈ 0
        @test einstein_tensor(M, x) ≈ zeros(n, n)
    end
end

struct TestSphere{N,T} <: Manifold
    r::T
end

struct TestSphericalMetric <: Metric end

@testset "scaled Sphere metric" begin
    Manifolds.manifold_dimension(::TestSphere{N}) where {N} = N
    function Manifolds.local_metric(M::MetricManifold{<:TestSphere,<:TestSphericalMetric}, x)
        r = base_manifold(M).r
        d = similar(x)
        d[1] = r^2
        d[2] = d[1] * sin(x[1])^2
        return Diagonal(d)
    end
    sph_to_cart(θ, ϕ) = [cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)]

    n = 2
    r = 10 * rand()
    θ, ϕ = π * rand(), 2π * rand()
    Sr = TestSphere{n,Float64}(r)
    S = Manifolds.Sphere(n)
    g = TestSphericalMetric()
    M = MetricManifold(Sr, g)

    @test manifold_dimension(M) == n
    @test base_manifold(M) === Sr
    @test metric(M) === g

    for vtype in (Vector, SVector{n}, MVector{n})
        x = vtype([θ, ϕ])
        G = Diagonal(vtype([1, sin(θ)^2])) .* r^2
        invG = Diagonal(vtype([1, 1 / sin(θ)^2])) ./ r^2
        v, w = normalize(randn(n)), normalize(randn(n))

        @test local_metric(M, x) ≈ G
        @test inverse_local_metric(M, x) ≈ invG
        @test det_local_metric(M, x) ≈ r^4 * sin(θ)^2
        @test log_local_metric_density(M, x) ≈ 2log(r) + log(sin(θ))
        @test inner(M, x, v, w) ≈ dot(v, G * w)
        @test norm(M, x, v) ≈ sqrt(dot(v, G * v))

        xcart = sph_to_cart(θ, ϕ)
        vcart = [cos(ϕ)*cos(θ) -sin(ϕ)*sin(θ);
                 sin(ϕ)*cos(θ)  cos(ϕ)*sin(θ);
                 -sin(θ) 0] * v

        if !Sys.iswindows() || Sys.ARCH == :x86_64
            @testset "numerically integrated geodesics for $vtype" begin
                T = 0:.1:1
                @test isapprox([sph_to_cart(yi...) for yi in exp(M, x, v, T)],
                               exp(S, xcart, vcart, T); atol=1e-3, rtol=1e-3)
                @test isapprox([sph_to_cart(yi...) for yi in geodesic(M, x, v, T)],
                               geodesic(S, xcart, vcart, T); atol=1e-3, rtol=1e-3)
            end
        end

        Γ₁ = christoffel_symbols_first(M, x)
        for i=1:n, j=1:n, k=1:n
            if (i,j,k) == (1,2,2) || (i,j,k) == (2,1,2)
                @test Γ₁[i,j,k] ≈ r^2*cos(θ)*sin(θ)
            elseif (i,j,k) == (2,2,1)
                @test Γ₁[i,j,k] ≈ -r^2*cos(θ)*sin(θ)
            else
                @test Γ₁[i,j,k] ≈ 0
            end
        end

        Γ₂ = christoffel_symbols_second(M, x)
        for l=1:n, i=1:n, j=1:n
            if (l,i,j) == (1,2,2)
                @test Γ₂[l,i,j] ≈ -cos(θ)*sin(θ)
            elseif (l,i,j) == (2,1,2) || (l,i,j) == (2,2,1)
                @test Γ₂[l,i,j] ≈ cot(θ)
            else
                @test Γ₂[l,i,j] ≈ 0
            end
        end

        R = riemann_tensor(M, x)
        for l=1:n, i=1:n, j=1:n, k=1:n
            if (l,i,j,k) == (2,1,1,2)
                @test R[l,i,j,k] ≈ -1
            elseif (l,i,j,k) == (2,1,2,1)
                @test R[l,i,j,k] ≈ 1
            elseif (l,i,j,k) == (1,2,1,2)
                @test R[l,i,j,k] ≈ sin(θ)^2
            elseif (l,i,j,k) == (1,2,2,1)
                @test R[l,i,j,k] ≈ -sin(θ)^2
            else
                @test R[l,i,j,k] ≈ 0
            end
        end

        @test ricci_tensor(M, x) ≈ G ./ r^2
        @test ricci_curvature(M, x) ≈ 2 / r^2
        @test gaussian_curvature(M, x) ≈ 1 / r^2
        @test einstein_tensor(M, x) ≈ ricci_tensor(M, x) - gaussian_curvature(M, x)  .* G
    end
end

struct BaseManifold{N} <: Manifold end
struct BaseManifoldMetric{M} <: Metric end

@testset "Metric decorator" begin
    Manifolds.manifold_dimension(::BaseManifold{N}) where {N} = N
    Manifolds.inner(::BaseManifold, x, v, w) = 2 * dot(v,w)
    Manifolds.exp!(::BaseManifold, y, x, v) = y .= x + 2 * v
    Manifolds.log!(::BaseManifold, v, x, y) = v .= (y - x) / 2
    Manifolds.project_tangent!(::BaseManifold, w, x, v) = w .= 2 .* v
    Manifolds.local_metric(::MetricManifold{BaseManifold{N},BaseManifoldMetric{N}},x) where N = 2*one(x*x')
    Manifolds.exp!(::MetricManifold{BaseManifold{N},BaseManifoldMetric{N}}, y, x, v) where N = exp!(base_manifold(M), y, x, v)
    function Manifolds.flat!(::BaseManifold, v::FVector{Manifolds.CotangentSpaceType}, x, w::FVector{Manifolds.TangentSpaceType})
        v.data .= 2 .* w.data
        return v
    end
    function Manifolds.sharp!(::BaseManifold, v::FVector{Manifolds.TangentSpaceType}, x, w::FVector{Manifolds.CotangentSpaceType})
        v.data .= w.data ./ 2
        return v
    end

    M = BaseManifold{3}()
    g = BaseManifoldMetric{3}()
    MM = MetricManifold(M, g)

    x = [0.1 0.2 0.4]
    v = [0.5 0.7 0.11]
    w = [0.13 0.17 0.19]
    y = similar(x)

    @test inner(M, x, v, w) == 2 * dot(v,w)
    @test inner(MM, x, v, w) === inner(M, x, v, w)
    @test norm(MM, x, v) === norm(M, x, v)
    @test distance(MM, x, y) === distance(M, x, y)
    @test exp(M, x, v) == x + 2 * v
    @test exp!(MM, y, x, v) === exp!(M, y, x, v)
    @test log(M, x, y) == (y - x) / 2
    @test log!(MM, v, x, y) === log!(M, v, x, y)
    @test retract!(MM, y, x, v) === retract!(M, y, x, v)
    @test retract!(MM, y, x, v, 1) === retract!(M, y, x, v, 1)
    @test project_tangent!(MM, w, x, v) === project_tangent!(M, w, x, v)
    @test zero_tangent_vector!(MM, v, x) === zero_tangent_vector!(M, v, x)
    @test injectivity_radius(MM, x) === injectivity_radius(M, x)
    @test injectivity_radius(MM) === injectivity_radius(M)
    @test is_manifold_point(MM, x) === is_manifold_point(M, x)
    @test is_tangent_vector(MM, x, v) === is_tangent_vector(M, x, v)

    cov = flat(M, x, FVector(TangentSpace, v))
    cow = flat(M, x, FVector(TangentSpace, w))
    @test cov.data ≈ flat(MM, x, FVector(TangentSpace, v)).data
    cotspace = CotangentBundleFibers(M)
    @test cov.data ≈ 2 * v
    @test inner(M, x, v, w) ≈ inner(cotspace, x, cov.data, cow.data)
    @test inner(MM, x, v, w) ≈ inner(cotspace, x, cov.data, cow.data)
    @test sharp(M, x, cov).data ≈ v
end


text:
A1d3 = readtimeaxisarray("1dtest.csv", date(), header=false)
A1d = TimeAxisArray(randn(10), collect(1.0:10.0))
A1d[5] = NaN
A2d = readtimeaxisarray("2dtest.csv", date())
A4d1 = readtimeaxisarray("4dtest.csv", date(), headlines=3)

# split
@test split(A4d1, week) == split(week, A4d1)
@test map(size, split(A1d3, week)) == [(3,), (7,), (4,)]
@test map(size, split(A2d, week)) == [(3,3), (7,3), (4,3)]

# collapse
@test collapse(A1d1, first, mean) ≈ TimeAxisArray([mean(A1d1.data)], timestamps(A1d1)[1:1])
@test collapse(A2d1, last, std) ≈ TimeAxisArray(std(A2d1.data, 1), [timestamps(A2d1)[end]], A2d1.axes[2])
@test collapse(A2d1, last, std) ≈ collapse(A2d1, last, x->std(x,1), lift=false)

# downsample
@test downsample(A1d3, week, first) == TimeAxisArray(A1d3.data[[1,4,11]], timestamps(A1d3)[[1,4,11]])
@test downsample(A2d, week, last, maximum) == TimeAxisArray([A2d.data[[2,4,13],1] A2d.data[[2,4,11],2] A2d.data[[1,7,13],3]], timestamps(A2d)[[3,10,14]], A2d.axes[2])
@test downsample(A2d, week, last, maximum) == downsample(A2d, week, last, x->maximum(x,1), lift=false)

# moving
@test moving(A1d2, mean, 3) ≈ TimeAxisArray(map(mean, TimeAxisArray[A1d2[1:3], A1d2[2:4], A1d2[3:5]]), timestamps(A1d2)[3:5])
@test moving(A2d, x->maximum(x,1), 10, lift=false) == moving(A2d, maximum, 10)
@test moving(A2d, maximum, 10) == TimeAxisArray(
                                                               [[A2d[2,1] A2d[4,2] A2d[7,3]];
                                                                [A2d[2,1] A2d[4,2] A2d[7,3]];
                                                                [A2d[3,1] A2d[4,2] A2d[7,3]];
                                                                [A2d[4,1] A2d[4,2] A2d[13,3]];
                                                                [A2d[13,1] A2d[5,2] A2d[13,3]]
                                                                ], timestamps(A2d)[10:end], A2d.axes[2:end]...)

# lead
T = size(A4d1, 1)
@test lead(A1d1) == TimeAxisArray(A1d1.data[2:end], timestamps(A1d1)[1:end-1])
@test lead(A1d1) == lead(A1d1, 1)
@test lead(A1d1, 3) == TimeAxisArray(A1d1.data[4:end], timestamps(A1d1)[1:end-3])
@test lead(A4d1) == TimeAxisArray(A4d1[Axis{:Timestamp}(2:T)].data, timestamps(A4d1)[1:end-1], A4d1.axes[2:end]...)
@test lead(A4d1) == lead(A4d1, 1)
@test lead(A4d1, 2) == TimeAxisArray(A4d1[Axis{:Timestamp}(3:T)].data, timestamps(A4d1)[1:end-2], A4d1.axes[2:end]...)

# lag
@test lag(A1d1) == TimeAxisArray(A1d1.data[1:end-1], timestamps(A1d1)[2:end])
@test lag(A1d1) == lag(A1d1, 1)
@test lag(A1d1, 3) == TimeAxisArray(A1d1.data[1:end-3], timestamps(A1d1)[4:end])
@test lag(A4d1) == TimeAxisArray(A4d1[Axis{:Timestamp}(1:T-1)].data, timestamps(A4d1)[2:end], A4d1.axes[2:end]...)
@test lag(A4d1) == lag(A4d1, 1)
@test lag(A4d1, 2) == TimeAxisArray(A4d1[Axis{:Timestamp}(1:T-2)].data, timestamps(A4d1)[3:end], A4d1.axes[2:end]...)

# diff
@test isequal(diff(A1d), TimeAxisArray(A1d[2:end] .- A1d[1:end-1], timestamps(A1d)[2:end]))
@test isequal(diff(A1d, 2), diff(diff(A1d)))
@test diff(A4d1) == TimeAxisArray(A4d1[2:end,:,:,:] .- A4d1[1:end-1,:,:,:], timestamps(A4d1)[2:end], A4d1.axes[2:end]...)
@test diff(A4d1, 2) == diff(diff(A4d1))

# percentchange
T = size(A2d, 1)
@test isapprox(percentchange(TimeAxisArray([1,2,3,4,5], 1:5)), TimeAxisArray([1,1/2,1/3,1/4], 2:5))
@test isapprox(percentchange(TimeAxisArray([1,2,3,4,5], 1:5), logdiff=true), TimeAxisArray(log1p([1,1/2,1/3,1/4]), 2:5))
@test isapproxequal(percentchange(A2d), TimeAxisArray((A2d[2:end,:] .- A2d[1:end-1,:]) ./ A2d[1:end-1,:], timestamps(A2d)[2:end], A2d.axes[2:end]...))
@test isapproxequal(percentchange(A2d, logdiff=true), map(log1p, percentchange(A2d)))

# dropif
@test isequal(dropif(any, x->x<=0, A1d), A1d[find((A1d .> 0) | isnan(A1d))])
@test dropif(any, x->x>0, A2d) == TimeAxisArray(zeros(0,3), Date[], A2d.axes[2])
@test isequal(dropif(any, x->x<0, A2d), A2d)

# dropnan
most(f::Function, x::AbstractArray) = sum(f(x)) > length(x)/2
@test dropnan(A1d) == A1d[[1:4;6:10]]
@test isequal(dropnan(A2d), A2d)
@test isequal(dropnan(A2d, all), A2d)
@test isequal(dropnan(A2d, most), A2d[[1:10;12:14], :])
@test dropnan(A2d, any) == A2d[[2:3;5:6;8:10;12;14], :]

max_stars_repo_name: GordStephen/TimeAxisArrays

max_issues_repo_name: GordStephen/TimeAxisArrays

max_forks_repo_name: GordStephen/TimeAxisArrays


text:
function AUSM(Qbase,Qcon,cellxmax,cellymax,vecAx,vecAy,specific_heat_ratio)
    E_adv_hat = zeros(cellxmax+1,cellymax)
    F_adv_hat = zeros(cellxmax,cellymax+1)

    for i in 2:cellxmax+1 -1
        for j in 2:cellymax -1
            # i-1セル
            cell_vecAx = zeros(2)
            cell_vecAx[1] = 0.5*(vecAx[i,j,1]+vecAx[i+1,j,1])
            cell_vecAx[2] = 0.5*(vecAx[i,j,2]+vecAx[i+1,j,2])

            rho = Qbase[i-1,j,1]
            U = Qbase[i-1,j,2]*cell_vecAx[1] + Qbase[i-1,j,3]*cell_vecAx[2]
            p = Qbase[i-1,j,4]
            a = (specific_heat_ratio * rho / p)^0.5
            M = u/a

            if M <= 1
                u_p = 0.25 * M*a * (M+1)^2
                p_p = 0.25 *p * (1+M)^2 *(2-M)
            elseif M > 1
                u_p = 0.5 * a * (M+abs(M))^2
                p_p = 0.5 *p * (M+abs(M))/M
            end

            # iセル
            cell_vecAx = zeros(2)
            cell_vecAx[1] = 0.5*(vecAx[i+1,j,1]+vecAx[i+2,j,1])
            cell_vecAx[2] = 0.5*(vecAx[i+1,j,2]+vecAx[i+2,j,2])

            rho = Qbase[i,j,1]
            U = Qbase[i,j,2]*cell_vecAx[1] + Qbase[i,j,3]*cell_vecAx[2]
            p = Qbase[i,j,4]
            a = (specific_heat_ratio * rho / p)^0.5
            M = u/a

            if M <= 1
                u_m = -0.25 * M*a * (M-1)^2
                p_m = 0.25 *p * (1-M)^2 *(2+M)
            elseif M > 1
                u_m = 0.5 * a * (M-abs(M))^2
                p_m = 0.5 *p * (M-abs(M))/M
            end

            # bound
            u_half = u_p + u_m
            p_half = p_p + p_m
            temp_vecX = zeros(4)
            temp_vecX[2] = vecAx[i,j,1]
            temp_vecX[3] = vecAx[i,j,2]

            for k in 1:4
                #----------------------------------------
                E_adv_hat[i,j,k] = 0.5 * ((u_half+abs(u_half))*Qcon_hat[i-1,j,k]+(u_half-abs(u_half))*Qcon_hat[i,j,k]) + p_half * temp_vecX[k]
            end
        end
    end

    for i in 2:cellxmax -1
        for j in 2:cellymax+1 -1
            # j-1セル
            cell_vecAy = zeros(2)
            cell_vecAy[1] = 0.5*(vecAy[i,j,1]+vecAy[i,j+1,1])
            cell_vecAy[2] = 0.5*(vecAy[i,j,2]+vecAy[i,j+1,2])

            rho = Qbase[i,j-1,1]
            #-------------
            V = Qbase[i,j-1,2]*cell_vecAy[1] + Qbase[i,j-1,3]*cell_vecAy[2]
            p = Qbase[i,j-1,4]
            a = (specific_heat_ratio * rho / p)^0.5
            M = u/a

            if M <= 1
                v_p = 0.25 * M*a * (M+1)^2
                p_p = 0.25 *p * (1+M)^2 *(2-M)
            elseif M > 1
                v_p = 0.5 * a * (M+abs(M))^2
                p_p = 0.5 *p * (M+abs(M))/M
            end

            # iセル
            cell_vecAy = zeros(2)
            cell_vecAy[1] = 0.5*(vecAy[i,j+1,1]+vecAy[i,j+2,1])
            cell_vecAy[2] = 0.5*(vecAy[i,j+1,2]+vecAy[i,j+2,2])

            rho = Qbase[i,j,1]
            V = Qbase[i,j,2]*cell_vecAy[1] + Qbase[i,j,3]*cell_vecAy[2]
            p = Qbase[i,j,4]
            a = (specific_heat_ratio * rho / p)^0.5
            M = u/a

            if M <= 1
                v_m = -0.25 * M*a * (M-1)^2
                p_m = 0.25 *p * (1-M)^2 *(2+M)
            elseif M > 1
                v_m = 0.5 * a * (M-abs(M))^2
                p_m = 0.5 *p * (M-abs(M))/M
            end

            # bound
            v_half = v_p + v_m
            p_half = p_p + p_m
            temp_vecX = zeros(4)
            temp_vecX[2] = vecAy[i,j,1]
            temp_vecX[3] = vecAy[i,j,2]

            for k in 1:4
                #----------------------------------------
                F_adv_hat[i,j,k] = 0.5 * ((v_half+abs(v_half))*Qcon_hat[i,j-1,k]+(v_half-abs(v_half))*Qcon_hat[i,j,k]) + p_half * temp_vecX[k]
            end
        end
    end

    return E_adv_hat, F_adv_hat
end

function setup_cell_flux_hat(Qcon, cellxmax, cellymax, cellzmax, cell_Ahat_plas, cell_Ahat_minus, cell_Bhat_plas, cell_Bhat_minus, cell_Chat_plas, cell_Chat_minus)
    cell_E_hat_plas  = zeros(cellxmax, cellymax, cellzmax, 5)
    cell_E_hat_minus = zeros(cellxmax, cellymax, cellzmax, 5)
    cell_F_hat_plas  = zeros(cellxmax, cellymax, cellzmax, 5)
    cell_F_hat_minus = zeros(cellxmax, cellymax, cellzmax, 5)
    cell_G_hat_plas  = zeros(cellxmax, cellymax, cellzmax, 5)
    cell_G_hat_minus = zeros(cellxmax, cellymax, cellzmax, 5)

    for i in 1:cellxmax
        for j in 1:cellymax
            for k in 1:cellzmax
                for l in 1:5
                    for m in 1:5
                        # Ehat = Ahat * Qhat = Ahat * Q/J = 1/J*Ahat * Q
                        # 1/J*Ahatの計算をvecAを使ったため，ここではQconを掛けてる
                        cell_E_hat_plas[i,j,k,l]  += cell_Ahat_plas[i,j,k,l,m]*Qcon[i,j,k,m]
                        cell_E_hat_minus[i,j,k,l] += cell_Ahat_minus[i,j,k,l,m]*Qcon[i,j,k,m]
                        cell_F_hat_plas[i,j,k,l]  += cell_Bhat_plas[i,j,k,l,m]*Qcon[i,j,k,m]
                        cell_F_hat_minus[i,j,k,l] += cell_Bhat_minus[i,j,k,l,m]*Qcon[i,j,k,m]
                        cell_G_hat_plas[i,j,k,l]  += cell_Chat_plas[i,j,k,l,m]*Qcon[i,j,k,m]
                        cell_G_hat_minus[i,j,k,l] += cell_Chat_minus[i,j,k,l,m]*Qcon[i,j,k,m]
                    end
                end
            end
        end
    end

    return cell_E_hat_plas, cell_E_hat_minus, cell_F_hat_plas, cell_F_hat_minus, cell_G_hat_plas, cell_G_hat_minus
end

function cal_jacobi(Qbase, Qcon, cellxmax, cellymax, cellzmax, specific_heat_ratio, vecAx, vecAy, vecAz, volume)
    """
    Qbase=[rho,u,v,p]
    Qcon=[rho,rhou,rhov,e]

    q=u^2+v^2
    c=(rp/rho)^0.5
    theta=Z=U= A^xi_x * u   +   A^xi_y * v
    H=(e+p)/rho
    b1=q^2/2 * (r-1)/(c^2)
    b2=(r-1)/(c^2)

    xi_x_bar= xi_x / (xi_x^2+xi_y^2)^0.5
    xi_y_bar= xi_y / (xi_x^2+xi_y^2)^0.5
    U_bar=Z_bar=theta_bar=

    A = Rleft * Lambda * Rright
    """

    cell_Ahat_plas  = zeros(cellxmax, cellymax, cellzmax, 5, 5)
    cell_Ahat_minus = zeros(cellxmax, cellymax, cellzmax, 5, 5)
    cell_Bhat_plas  = zeros(cellxmax, cellymax, cellzmax, 5, 5)
    cell_Bhat_minus = zeros(cellxmax, cellymax, cellzmax, 5, 5)
    cell_Chat_plas  = zeros(cellxmax, cellymax, cellzmax, 5, 5)
    cell_Chat_minus = zeros(cellxmax, cellymax, cellzmax, 5, 5)
    r = specific_heat_ratio

    for i in 1:cellxmax
        for j in 1:cellymax
            for k in 1:cellzmax
                q = (Qbase[i,j,k,2]^2+Qbase[i,j,k,3]^2+Qbase[i,j,k,4]^2)^0.5
                c = 0
                try c = (r*Qbase[i,j,k,5]/Qbase[i,j,k,1])^0.5
                catch
                    println("\n"*string(i)*","*string(j)*","*string(k)*" Qbase error")
                    println(Qbase[i,j,k,:])
                    println("\n")
                    throw(UndefVarError(:x))
                end

                H  = (Qcon[i,j,k,5]+Qbase[i,j,k,5])/Qbase[i,j,k,1]
                b1 = q^2/2 * (r-1)/(c^2)
                b2 = (r-1)/(c^2)
                u  = Qbase[i,j,k,2]
                v  = Qbase[i,j,k,3]
                w  = Qbase[i,j,k,4]

                # cell_Ahat
                k_x = 0.5*(vecAx[i,j,k,1]+vecAx[i+1,j,k,1])
                k_y = 0.5*(vecAx[i,j,k,2]+vecAx[i+1,j,k,2])
                k_z = 0.5*(vecAx[i,j,k,3]+vecAx[i+1,j,k,3])

                Rleft_E, Lambda_E, Rright_E = Eigenvalue_vector(q,c,H,b1,b2,u,v,w,k_x,k_y,k_z)
                A_plas, A_minus = RLmbdaR(Rleft_E,Lambda_E,Rright_E)

                for l in 1:5
                    for m in 1:5
                        cell_Ahat_plas[i,j,k,l,m]  = A_plas[l,m]
                        cell_Ahat_minus[i,j,k,l,m] = A_minus[l,m]
                    end
                end

                # cell_Bhat
                k_x = 0.5*(vecAy[i,j,k,1]+vecAy[i,j+1,k,1])
                k_y = 0.5*(vecAy[i,j,k,2]+vecAy[i,j+1,k,2])
                k_z = 0.5*(vecAy[i,j,k,3]+vecAy[i,j+1,k,3])

                Rleft_E, Lambda_E, Rright_E = Eigenvalue_vector(q,c,H,b1,b2,u,v,w,k_x,k_y,k_z)
                A_plas, A_minus = RLmbdaR(Rleft_E,Lambda_E,Rright_E)

                for l in 1:5
                    for m in 1:5
                        cell_Bhat_plas[i,j,k,l,m]  = A_plas[l,m]
                        cell_Bhat_minus[i,j,k,l,m] = A_minus[l,m]
                    end
                end

                # cell_Chat
                k_x = 0.5*(vecAz[i,j,k,1]+vecAz[i,j,k+1,1])
                k_y = 0.5*(vecAz[i,j,k,2]+vecAz[i,j,k+1,2])
                k_z = 0.5*(vecAz[i,j,k,3]+vecAz[i,j,k+1,3])

                Rleft_E, Lambda_E, Rright_E = Eigenvalue_vector(q,c,H,b1,b2,u,v,w,k_x,k_y,k_z)
                A_plas, A_minus = RLmbdaR(Rleft_E,Lambda_E,Rright_E)

                for l in 1:5
                    for m in 1:5
                        cell_Chat_plas[i,j,k,l,m]  = A_plas[l,m]
                        cell_Chat_minus[i,j,k,l,m] = A_minus[l,m]
                    end
                end
            end
        end
    end

    return cell_Ahat_plas, cell_Ahat_minus, cell_Bhat_plas, cell_Bhat_minus, cell_Chat_plas, cell_Chat_minus
end

function Eigenvalue_vector(q,c,H,b1,b2,u,v,w,k_x,k_y,k_z)
    Z       = k_x*u + k_y*v + k_z*w
    k_x_bar = k_x / (k_x^2 + k_y^2 + k_z^2)^0.5
    k_y_bar = k_y / (k_x^2 + k_y^2 + k_z^2)^0.5
    k_z_bar = k_z / (k_x^2 + k_y^2 + k_z^2)^0.5
    Z_bar   =   Z / (k_x^2 + k_y^2 + k_z^2)^0.5

    Lambda = [Z - c*(k_x^2+k_y^2+k_z^2)^0.5
              Z
              Z
              Z
              Z + c*(k_x^2+k_y^2+k_z^2)^0.5]

    temp1 = k_y_bar*q^2/2+k_z_bar*u*c-k_x_bar*w*c
    temp2 = k_z_bar*q^2/2-k_y_bar*u*c+k_x_bar*v*c
    temp3 = k_x_bar*q^2/2-k_z_bar*v*c+k_y_bar*w*c

    Rleft = [[1           k_y_bar             k_z_bar             k_x_bar             1          ]
             [u-k_x_bar*c k_y_bar*u+k_z_bar*c k_z_bar*u-k_y_bar*c k_x_bar*u           u+k_x_bar*c]
             [v-k_y_bar*c k_y_bar*v           k_z_bar*v+k_x_bar*c k_x_bar*v-k_z_bar*c v+k_y_bar*c]
             [w-k_z_bar*c k_y_bar*w-k_x_bar*c k_z_bar*w           k_x_bar*w+k_y_bar*c w+k_z_bar*c]
             [H-c*Z_bar   temp1               temp2               temp3               H+c*Z_bar  ]]

    temp1 = k_y_bar*(1-b1)+(k_x_bar*w-k_z_bar*u)/c
    temp2 = k_z_bar*(1-b1)+(k_y_bar*u-k_x_bar*v)/c
    temp3 = k_x_bar*(1-b1)+(k_z_bar*v-k_y_bar*w)/c

    Rright = [[(b1+Z_bar/c)/2 -(k_x_bar/c+b2*u)/2    -(k_y_bar/c+b2*v)/2    -(k_z_bar/c+b2*w)/2    b2/2       ]
              [temp1          k_y_bar*b2*u+k_z_bar/c k_y_bar*b2*v           k_y_bar*b2*w-k_x_bar/c -k_y_bar*b2]
              [temp2          k_z_bar*b2*u-k_y_bar/c k_z_bar*b2*v+k_x_bar/c k_z_bar*b2*w           -k_z_bar*b2]
              [temp3          k_x_bar*b2*u           k_x_bar*b2*v-k_z_bar/c k_x_bar*b2*w+k_z_bar/c -k_x_bar*b2]
              [(b1-Z_bar/c)/2 (k_x_bar/c-b2*u)/2     (k_y_bar/c-b2*v)/2     (k_z_bar/c-b2*w)/2     b2/2       ]]

    return Rleft, Lambda, Rright
end

function RLmbdaR(R, Lam, Rm)
    Lamp = zeros(5,5)
    Lamm = zeros(5,5)
    for i in 1:length(Lam)
        Lamp[i,i] = 0.5*(Lam[i]+abs(Lam[i]))
        Lamm[i,i] = 0.5*(Lam[i]-abs(Lam[i]))
    end
    Ap = nn_inner_product(nn_inner_product(R,Lamp),Rm)
    Am = nn_inner_product(nn_inner_product(R,Lamm),Rm)
    return Ap, Am
end

function nn_inner_product(a, b)
    temp = zeros(size(a)[1],size(a)[1])
    for i in 1:size(a)[1]
        for j in 1:size(a)[1]
            for k in 1:size(a)[1]
                temp[i,j] += a[i,k]*b[k,j]
            end
        end
    end
    return temp
end

function FVS(cell_E_hat_plas, cell_E_hat_minus, cell_F_hat_plas, cell_F_hat_minus, cell_G_hat_plas, cell_G_hat_minus, cellxmax, cellymax, cellzmax)
    E_hat = zeros(cellxmax+1,   cellymax,   cellzmax, 5)
    F_hat = zeros(  cellxmax, cellymax+1,   cellzmax, 5)
    G_hat = zeros(  cellxmax,   cellymax, cellzmax+1, 5)

    Threads.@threads for i in 2:cellxmax+1-1
        for j in 2:cellymax-1
            for k in 2:cellzmax-1
                for l in 1:5
                    E_hat[i,j,k,l] = cell_E_hat_plas[i-1,j,k,l] + cell_E_hat_minus[i,j,k,l]
                end
            end
        end
    end
    Threads.@threads for i in 2:cellxmax-1
        for j in 2:cellymax+1-1
            for k in 2:cellzmax-1
                for l in 1:5
                    F_hat[i,j,k,l] = cell_F_hat_plas[i,j-1,k,l] + cell_F_hat_minus[i,j,k,l]
                end
            end
        end
    end
    Threads.@threads for i in 2:cellxmax-1
        for j in 2:cellymax-1
            for k in 2:cellzmax+1-1
                for l in 1:5
                    G_hat[i,j,k,l] = cell_G_hat_plas[i,j,k-1,l] + cell_G_hat_minus[i,j,k,l]
                end
            end
        end
    end

    return E_hat, F_hat, G_hat
end

max_stars_repo_name: hide-dog/test

max_issues_repo_name: hide-dog/test

max_forks_repo_name: hide-dog/test


text:
using FdeSolver
using Test
using SpecialFunctions
using Statistics

@testset "FdeSolver.jl" begin

    tSpan = [0, 1]
    y0 = 0
    β = 0.9

    par = β
    F(t, y, par) = (40320 ./ gamma(9 - par) .* t .^ (8 - par) .- 3 .* gamma(5 + par / 2)
               ./ gamma(5 - par / 2) .* t .^ (4 - par / 2) .+ 9/4 * gamma(par + 1) .+
               (3 / 2 .* t .^ (par / 2) .- t .^ 4) .^ 3 .- y .^ (3 / 2))

    JacobF(t, y, β) = -(3 / 2) .* y .^ (1 / 2)

    t, Yapp = FDEsolver(F, tSpan, y0, β, par, StopIt = "Convergence", tol = 10e-8, itmax = 15)
    t1, Yapp1 = FDEsolver(F, tSpan, y0, β, par, JF = JacobF, StopIt = "Convergence", tol = 10e-8, itmax = 15)

    @test @isdefined(t)
    @test @isdefined(Yapp)

    @test @isdefined(t1)
    @test @isdefined(Yapp1)

    # @test abs(mean(Yapp) - 0.99971) < 0.05
    # @test abs(mean(Yapp1) - 0.99971) < 0.05

end