text:
using CSFML.LibCSFML

using ModernGL

mode = sfVideoMode(700, 720, 32)
window = sfRenderWindow_create(mode, "SFML window", sfResize | sfClose, C_NULL)
@assert window != C_NULL
sfWindow_setVerticalSyncEnabled(window,true)
sfWindow_setActive(window, sfTrue)

texture = sfTexture_createFromFile(joinpath("/julia/playgrounds/resources/ghost.jpg"), C_NULL)
@assert texture != C_NULL

sprite = sfSprite_create()
sfSprite_setTexture(sprite, texture, sfTrue)

font = sfFont_createFromFile(joinpath("/julia/playgrounds/resources/NotoSans-Black.ttf"))
@assert font != C_NULL

text = sfText_create()
sfText_setString(text, "The quick brown fox jumps over the lazy dog.")
sfText_setPosition(text, sfVector2f(100,5))
sfText_setFont(text, font)
sfText_setCharacterSize(text, 20)

event_ref = Ref{sfEvent}()

# complete list of event types :
#     sfEvtClosed = 0
#     sfEvtResized = 1
#     sfEvtLostFocus = 2
#     sfEvtGainedFocus = 3
#     sfEvtTextEntered = 4
#     sfEvtKeyPressed = 5
#     sfEvtKeyReleased = 6
#     sfEvtMouseWheelMoved = 7
#     sfEvtMouseWheelScrolled = 8
#     sfEvtMouseButtonPressed = 9
#     sfEvtMouseButtonReleased = 10
#     sfEvtMouseMoved = 11
#     sfEvtMouseEntered = 12
#     sfEvtMouseLeft = 13
#     sfEvtJoystickButtonPressed = 14
#     sfEvtJoystickButtonReleased = 15
#     sfEvtJoystickMoved = 16
#     sfEvtJoystickConnected = 17
#     sfEvtJoystickDisconnected = 18
#     sfEvtTouchBegan = 19
#     sfEvtTouchMoved = 20
#     sfEvtTouchEnded = 21
#     sfEvtSensorChanged = 22
#     sfEvtCount = 23
clock = sfClock_create()

# Let's express time in seconds.
function get_time(clock ::Ptr{Nothing})
    sfTime_asSeconds(sfClock_getElapsedTime(clock))
end
function restart(clock ::Ptr{Nothing})
     sfTime_asSeconds(sfClock_restart(clock))
end

target_frequency = 24.0
running = true
t0 = 0
while (running)

    while Bool(sfRenderWindow_isOpen(window)) && running
        frame_timestamp = get_time(clock)
        frame_time = frame_timestamp - t0
        t0 = frame_timestamp
        actual_frequency = 1 / frame_time

        # process events
        while Bool(sfRenderWindow_pollEvent(window, event_ref))
            # close window : exit
            event_ref[].type == sfEvtClosed && sfRenderWindow_close(window)
            event_ref[].type == sfEvtResized && println("Trigger sfEvtResized.")
            event_ref[].type == sfEvtLostFocus && println("Trigger sfEvtLostFocus.")
            event_ref[].type == sfEvtGainedFocus && println("Trigger sfEvtGainedFocus.")
            event_ref[].type == sfEvtTextEntered && println("Trigger sfEvtTextEntered: $(event_ref[].text.unicode)")
            event_ref[].type == sfEvtKeyPressed && println("Trigger sfEvtKeyPressed: $(event_ref[].key.code)")
            event_ref[].type == sfEvtKeyReleased && println("Trigger sfEvtKeyReleased: $(event_ref[].key.code)")
            event_ref[].type == sfEvtMouseWheelMoved && println("Trigger sfEvtMouseWheelMoved: $(event_ref[].mouseWheel.x), $(event_ref[].mouseWheel.y)")
            event_ref[].type == sfEvtMouseWheelScrolled && println("Trigger sfEvtMouseWheelScrolled: $(event_ref[].mouseWheelScroll.wheel)")
            event_ref[].type == sfEvtMouseButtonPressed && println("Trigger sfEvtMouseButtonPressed: $(event_ref[].mouseButton.button)")
            event_ref[].type == sfEvtMouseButtonReleased && println("Trigger sfEvtMouseButtonReleased: $(event_ref[].mouseButton.x), $(event_ref[].mouseButton.y)")
            event_ref[].type == sfEvtMouseMoved && println("Trigger sfEvtMouseMoved: $(event_ref[].mouseMove.x), $(event_ref[].mouseMove.y)")

            if event_ref[].type == sfEvtClosed
                running = false
                println("Render window closed.")
            elseif event_ref[].type == sfEvtResized
                glViewport(9, 9, event_ref[].size.width, event_ref[].size.height)
                println("Viewport updated.")
            elseif event_ref[].type == sfEvtKeyPressed
                if event_ref[].key.code == sfKeyQ
                    running = false
                    println("Quiting application.")
                end
            end
        end
        time_event_processing = get_time(clock) - t0
        #println("Event processing took $time_event_processing seconds.")

        # Render scene

        sfRenderWindow_clear(window, sfColor_fromRGBA(24,20,18,255))
        sfRenderWindow_drawSprite(window, sprite, C_NULL)
        sfRenderWindow_drawText(window, text, C_NULL)

        # circle = sfCircleShape_create()
        # sfCircleShape_setRadius(circle, 30)

        # sfCircleShape_setPosition(circle, sfVector2f(80+cos(frame_timestamp*22)*20, 240+sin(frame_timestamp*3)*200))
        # sfCircleShape_setFillColor(circle,sfColor_fromRGBA(255,155,105,255) )
        # sfCircleShape_setOutlineColor(circle , sfColor_fromRGBA(255,155,125,255) )

        square = sfRectangleShape_create()
        sfRectangleShape_setSize(square, sfVector2f(50,50) )
        t_warp = frame_timestamp + sin(frame_timestamp*2.2)*1.5
        sfRectangleShape_setPosition(square, sfVector2f(80+cos(t_warp*2.3)*20, 240+sin(t_warp*1.8)*200))
        sfRectangleShape_setFillColor(square,sfColor_fromRGBA(255,155,105,255) )

        println(sfRenderWindow_getSize(window))
        sfRenderWindow_drawRectangleShape(window, square, C_NULL)
        # --sfRenderStates
        # state.shader = shader;
        # state.blendMode = sfBlendAlpha;
        # state.transform = sfTransform_Identity;
        # state.texture = NULL;
        # struct sfRenderStates
        #     blendMode::sfBlendMode
        #     transform::sfTransform
        #     texture::Ptr{sfTexture}
        #     shader::Ptr{sfShader}
        # end

        # --shader
        # string frag = r"void main() { gl_FragColor = vec4(1,0,0,1); }";
        # sfShader* shader = sfShader_createFromMemory(null, frag);
        # sf::Shader shader; shader.loadFromMemory(frag, sf::Shader::Fragment);

        time_render = get_time(clock) - time_event_processing - t0

        sf_text = sfText_create()
        sfText_setPosition(sf_text, sfVector2f(4,700))
        sfText_setString(sf_text, "frame render time : $time_render seconds.")
        sfText_setFont(sf_text, font)
        sfText_setCharacterSize(sf_text, 14)
        sfRenderWindow_drawText(window, sf_text, C_NULL)

        sf_text2 = sfText_create()
        sfText_setPosition(sf_text2, sfVector2f(4,680))
        sfText_setString(sf_text2, "fps : $actual_frequency")
        sfText_setFont(sf_text2, font)
        sfText_setCharacterSize(sf_text2, 14)
        sfRenderWindow_drawText(window, sf_text2, C_NULL)

        sfRenderWindow_display(window)
        time_to_sleep = 1 / target_frequency - time_render - time_event_processing

        if time_to_sleep > 0.005
            sleep(time_to_sleep-0.005)
        end
    end

end

sfText_destroy(text)
sfFont_destroy(font)
sfSprite_destroy(sprite)
sfTexture_destroy(texture)
sfRenderWindow_destroy(window)

max_stars_repo_name: leroy-walton/playgrounds

max_issues_repo_name: leroy-walton/playgrounds

max_forks_repo_name: leroy-walton/playgrounds


text:
function beta(u::UGrid, v::VGrid, dt)
    umax = maximum(u)
    vmax = maximum(v)

    (abs(umax) + abs(vmax)) / min(u.grid.dx, u.grid.dy) * dt
end
beta(sim::Simulation) = beta(sim.u, sim.v, sim.dt)

function diagr(sim::Simulation)
    getnu(sim.params, sim.u.grid) * sim.dt / min(sim.u.grid.dx, sim.u.grid.dy)^2
end

function reynh(sim::Simulation)
    beta(sim) / diagr(sim)
end

function reynhvort(sim::Simulation)
    vort = VortGrid(sim.u, sim.v)
    maxvort = maximum(vort)
    maxvort * min(sim.u.grid.dx, sim.v.grid.dy)^2 / getnu(sim.params, sim.u.grid)
end

ext: jl


text:
# implementation of the GPUCompiler interfaces for generating SPIR-V code

# https://github.com/llvm/llvm-project/blob/master/clang/lib/Basic/Targets/SPIR.h
# https://github.com/KhronosGroup/LLVM-SPIRV-Backend/blob/master/llvm/docs/SPIR-V-Backend.rst
# https://github.com/KhronosGroup/SPIRV-LLVM-Translator/blob/master/docs/SPIRVRepresentationInLLVM.rst

const SPIRV_LLVM_Translator_jll = LazyModule("SPIRV_LLVM_Translator_jll", UUID("4a5d46fc-d8cf-5151-a261-86b458210efb"))
const SPIRV_Tools_jll = LazyModule("SPIRV_Tools_jll", UUID("6ac6d60f-d740-5983-97d7-a4482c0689f4"))

## target

export SPIRVCompilerTarget

Base.@kwdef struct SPIRVCompilerTarget <: AbstractCompilerTarget
end

llvm_triple(::SPIRVCompilerTarget) = Int===Int64 ? "spir64-unknown-unknown" : "spirv-unknown-unknown"

# SPIRV is not supported by our LLVM builds, so we can't get a target machine
llvm_machine(::SPIRVCompilerTarget) = nothing

llvm_datalayout(::SPIRVCompilerTarget) = Int===Int64 ?
    "e-i64:64-v16:16-v24:32-v32:32-v48:64-v96:128-v192:256-v256:256-v512:512-v1024:1024" :
    "e-p:32:32-i64:64-v16:16-v24:32-v32:32-v48:64-v96:128-v192:256-v256:256-v512:512-v1024:1024"

## job

# TODO: encode debug build or not in the compiler job
#       https://github.com/JuliaGPU/CUDAnative.jl/issues/368
runtime_slug(job::CompilerJob{SPIRVCompilerTarget}) = "spirv"

function process_module!(job::CompilerJob{SPIRVCompilerTarget}, mod::LLVM.Module)
    # calling convention
    for f in functions(mod)
        # JuliaGPU/GPUCompiler.jl#97
        #callconv!(f, LLVM.API.LLVMSPIRFUNCCallConv)
    end
end

function process_entry!(job::CompilerJob{SPIRVCompilerTarget}, mod::LLVM.Module, entry::LLVM.Function)
    if job.source.kernel
        # HACK: Intel's compute runtime doesn't properly support SPIR-V's byval attribute.
        #       they do support struct byval, for OpenCL, so wrap byval parameters in a struct.
        entry = wrap_byval(job, mod, entry)

        # calling convention
        callconv!(entry, LLVM.API.LLVMSPIRKERNELCallConv)
    end

    return entry
end

function finish_module!(job::CompilerJob{SPIRVCompilerTarget}, mod::LLVM.Module)
    # SPIR-V does not support trap, and has no mechanism to abort compute kernels
    # (OpKill is only available in fragment execution mode)
    ModulePassManager() do pm
        add!(pm, ModulePass("RemoveTrap", rm_trap!))
        run!(pm, mod)
    end
end

@unlocked function mcgen(job::CompilerJob{SPIRVCompilerTarget}, mod::LLVM.Module,
                         format=LLVM.API.LLVMAssemblyFile)
    # write the bitcode to a temporary file (the SPIRV Translator library doesn't have a C API)
    mktemp() do input, input_io
        write(input_io, mod)
        flush(input_io)

        # compile to SPIR-V
        mktemp() do output, output_io
            SPIRV_LLVM_Translator_jll.llvm_spirv() do translator
                cmd = `$translator`
                if format == LLVM.API.LLVMAssemblyFile
                    cmd = `$cmd -spirv-text`
                end
                cmd = `$cmd --spirv-debug-info-version=ocl-100 -o $output $input`
                run(cmd)
            end

            # read back the file
            if format == LLVM.API.LLVMAssemblyFile
                read(output_io, String)
            else
                read(output_io)
            end
        end
    end
end

# reimplementation that uses `spirv-dis`, giving much more pleasant output
function code_native(io::IO, job::CompilerJob{SPIRVCompilerTarget}; raw::Bool=false, dump_module::Bool=false)
    if raw
        # The SPIRV Tools don't handle Julia's debug info, rejecting DW_LANG_Julia...
        # so just return what LLVM gives us in that case (which is also more faithful).
        asm, _ = codegen(:asm, job; strip=false, only_entry=!dump_module, validate=false)
        print(io, asm)
        return
    end

    obj, _ = codegen(:obj, job; strip=!raw, only_entry=!dump_module, validate=false)
    mktemp() do input_path, input_io
        write(input_io, obj)
        flush(input_io)

        SPIRV_Tools_jll.spirv_dis() do disassembler
            if io == stdout
                run(`$disassembler $input_path`)
            else
                mktemp() do output_path, output_io
                    run(`$disassembler $input_path -o $output_path`)
                    asm = read(output_io, String)
                    print(io, asm)
                end
            end
        end
    end
end

## LLVM passes

# remove llvm.trap and its uses from a module
function rm_trap!(mod::LLVM.Module)
    job = current_job::CompilerJob
    changed = false
    @timeit_debug to "remove trap" begin

    if haskey(functions(mod), "llvm.trap")
        trap = functions(mod)["llvm.trap"]

        for use in uses(trap)
            val = user(use)
            if isa(val, LLVM.CallInst)
                unsafe_delete!(LLVM.parent(val), val)
                changed = true
            end
        end

        @compiler_assert isempty(uses(trap)) job
        unsafe_delete!(mod, trap)
    end

    end
    return changed
end

# wrap byval pointers in a single-value struct
function wrap_byval(@nospecialize(job::CompilerJob), mod::LLVM.Module, entry_f::LLVM.Function)
    ctx = context(mod)
    entry_ft = eltype(llvmtype(entry_f)::LLVM.PointerType)::LLVM.FunctionType
    @compiler_assert return_type(entry_ft) == LLVM.VoidType(ctx) job

    args = classify_arguments(job, entry_f)
    filter!(args) do arg
        arg.cc != GHOST
    end

    # generate the wrapper function type & definition
    wrapper_types = LLVM.LLVMType[]
    for arg in args
        typ = if arg.cc == BITS_REF
            st = LLVM.StructType([eltype(arg.codegen.typ)]; ctx)
            LLVM.PointerType(st, addrspace(arg.codegen.typ))
        else
            convert(LLVMType, arg.typ; ctx)
        end
        push!(wrapper_types, typ)
    end
    wrapper_fn = LLVM.name(entry_f)
    LLVM.name!(entry_f, wrapper_fn * ".inner")
    wrapper_ft = LLVM.FunctionType(LLVM.VoidType(ctx), wrapper_types)
    wrapper_f = LLVM.Function(mod, wrapper_fn, wrapper_ft)

    # emit IR performing the "conversions"
    let builder = Builder(ctx)
        entry = BasicBlock(wrapper_f, "entry"; ctx)
        position!(builder, entry)

        wrapper_args = Vector{LLVM.Value}()

        # perform argument conversions
        for arg in args
            param = parameters(wrapper_f)[arg.codegen.i]
            attrs = parameter_attributes(wrapper_f, arg.codegen.i)
            if arg.cc == BITS_REF
                push!(attrs, EnumAttribute("byval", 0; ctx))
                ptr = struct_gep!(builder, param, 0)
                push!(wrapper_args, ptr)
            else
                push!(wrapper_args, param)
                for attr in collect(attrs)
                    push!(parameter_attributes(wrapper_f, arg.codegen.i), attr)
                end
            end
        end

        call!(builder, entry_f, wrapper_args)

        ret!(builder)

        dispose(builder)
    end

    # early-inline the original entry function into the wrapper
    delete!(function_attributes(entry_f), EnumAttribute("noinline", 0; ctx))
    push!(function_attributes(entry_f), EnumAttribute("alwaysinline", 0; ctx))
    linkage!(entry_f, LLVM.API.LLVMInternalLinkage)

    ModulePassManager() do pm
        always_inliner!(pm)
        run!(pm, mod)
    end

    return wrapper_f
end


text:
# ML parameter estimation with filter result log-likelihood (via Nelder-Mead)

function fit{T}(y::Matrix{T}, build::Function, theta0::Vector{T};
          u::Array{T}=zeros(size(y,1), build(theta0).nu))
    objective(theta) = kalman_filter(y, build(theta), u=u).loglik
    kfit = Optim.optimize(objective, theta0)
    return (kfit.minimum, build(kfit.minimum), filter(y, build(kfit.minimum)))
end #fit

# Expectation-Maximization (EM) parameter estimation

function fit{T}(y::Array{T}, pmodel::ParametrizedSSM, params::SSMParameters;
          u::Array{T}=zeros(size(y,1), pmodel.nu), eps::Float64=1e-6, niter::Int=typemax(Int))

    n = size(y, 1)
    @assert n == size(u, 1)

    y_orig = copy(y)
    y = y'
    y_notnan = !isnan(y)
    y = y .* y_notnan

    u_orig = copy(u)
    u = u'

    I_nx = eye(pmodel.nx)
    I0_nx = 0I_nx
    I_ny = eye(pmodel.ny)
    I0_ny = 0I_ny
    I_nu = eye(pmodel.nu)

    # Requirement: zeros-rows in G and H remain consistent
    x_deterministic = all(pmodel.G(1) .== 0, 2)
    all_x_deterministic = all(x_deterministic)
    y_deterministic = all(pmodel.H(1) .== 0, 2)
    all_y_deterministic = all(y_deterministic)

    if any(x_deterministic)
        OQ0, OQp = I_nx[find(x_deterministic), :], I_nx[find(!x_deterministic), :]
        Id(M_t::Array{Int,2}) = diagm(OQ0' * all((OQ0 * M_t * OQp') .== 0, 2)[:])
    else
        Id(_::Array{Int,2}) = I0_nx
    end

    phi(t)  = (pmodel.G(t)' * pmodel.G(t)) \ pmodel.G(t)'
    xi(t) = (pmodel.H(t)' * pmodel.H(t)) \ pmodel.H(t)'

    estimate_A   = length(params.A) > 0
    estimate_B   = length(params.B) > 0
    estimate_Q   = length(params.Q) > 0
    estimate_C   = length(params.C) > 0
    estimate_D   = length(params.D) > 0
    estimate_R   = length(params.R) > 0
    estimate_x1  = length(params.x1) > 0
    estimate_P1  = length(params.P1) > 0

    function em_kernel!{T}(params::SSMParameters{T})

        m = pmodel(params)

        #print("Expectations (smoothing)... ")
        #tic()

        if estimate_A | estimate_Q
            exs, P, Plag1, loglik = lag1_smooth(y_orig, u_orig, m)
        else
            smoothed = kalman_smooth(y_orig, m, u=u_orig)
            exs, P, loglik = smoothed.smoothed', smoothed.error_cov, smoothed.loglik
        end

        #toc()
        #println("Negative log-likelihood: ", loglik)

        #print("Maximizations... ")
        #tic()

        if estimate_A | estimate_B | estimate_Q
            Qinv    = inv(pmodel.Q(params.Q))
            Q_      = all_x_deterministic ? I0_nx : phi(1)' * Qinv * phi(1)
        end #if ABQ

        if estimate_B | estimate_C | estimate_D | estimate_R
            Rinv  = inv(pmodel.R(params.R))
            R_    = all_y_deterministic ? I0_ny : xi(1)' * Rinv * xi(1)
        end #if

        if estimate_x1
            Linv = inv(m.P1) #TODO: Degenerate accomodations?
        end #if

        HRH(t) = pmodel.H(t) * pmodel.R(params.R) * pmodel.H(t)'
        HRH_nonzero_rows(t) = diag(HRH(t)) .!= 0

        function N(t::Int)
            HRHt = HRH(t)
            O = I_ny[find(y_notnan[:, t] & HRH_nonzero_rows(t)), :]
            return I - HRHt * O' * inv(O * HRHt * O') * O
        end #N

        ex  = exs[:, 1]
        exx = zeros(m.nx, m.nx)

        yt = y[:, 1]
        Nt = N(1)
        ey = yt - Nt * (yt - m.C(1) * exs[:, 1] - m.D(1) * u[:, 1])

        But = m.B(1) * u[:, 1]
        Dut = m.D(1) * u[:, 1]

        if estimate_A # A setup
            A_S1 = zeros(length(params.A), length(params.A))
            A_S2 = zeros(params.A)
        end #if A

        if estimate_B # B setup
            pkp = kron(I_nu, pmodel.B1(1))
            f_Bt = pkp * pmodel.B2.f
            D_Bt = pkp * pmodel.B2.D

            B_kp = kron(u[:, 1]', I_nx)
            fB = B_kp * f_Bt
            DB = B_kp * D_Bt

            Idt = I_nx
            A_  = I_nx
            f_  = B_kp * f_Bt * 0
            D_  = B_kp * D_Bt * 0
            M   = 1 * (m.A(1) .!= 0) #Assumption: constant position of nonzero values in A
            M_t = copy(M)
            EX0 = exs[:, 1]

            Dt1 = ey - m.C(1) * (I - Idt) * ex - m.C(1) * Idt * (A_ * EX0 + f_) - Dut
            Dt2 = m.C(1) * Idt * D_
            Dt3 = zeros(pmodel.nx)
            Dt4 = DB * 0

            Idt1 = Idt
            At1_ = A_
            ft1_ = f_
            Dt1_ = D_

            B_S1 = Dt2' * R_ * Dt2
            B_S2 = Dt2' * R_ * Dt1
        end #if B

        if estimate_Q # Q setup
            Q_S1 = zeros(length(params.Q), length(params.Q))
            Q_S2 = zeros(params.Q)
        end #if

        if estimate_C # C setup
            C_S1 = zeros(length(params.C), length(params.C))
            C_S2 = zeros(params.C)
        end #if

        if estimate_D # D setup
            D_S1 = zeros(length(params.D), length(params.D))
            D_S2 = zeros(params.D)
        end #if

        if estimate_R # R setup
            R_S1 = zeros(length(params.R), length(params.R))
            R_S2 = zeros(params.R)
        end #if

        for t in 1:n

            At      = m.A(t)
            Pt      = P[:, :, t]
            ut      = u[:, t]

            ex_prev = ex
            ex      = exs[:, t]

            exx_prev  = exx
            exx       = Pt + ex * ex'

            !all_x_deterministic && (estimate_A || estimate_B || estimate_Q) && (Q_ = phi(t)' * Qinv * phi(t))
            !all_y_deterministic && (estimate_B || estimate_C || estimate_D || estimate_R) && (R_ = xi(t)' * Rinv * xi(t))

            if estimate_A || estimate_Q || estimate_C

                if (estimate_A || estimate_Q) && t > 1

                    exx1 = Plag1[:, :, t] + ex * ex_prev'
                    But_prev = But
                    But = m.B(t) * ut

                    if estimate_A
                        pkp = kron(pmodel.A3(t)', pmodel.A1(t))
                        D_At = pkp * pmodel.A2.D
                        f_At = pkp * pmodel.A2.f
                        A_kp = kron(exx_prev, Q_)
                        A_S1 += D_At' * A_kp * D_At
                        A_S2 += D_At' * (vec(Q_ * exx1) -
                                      A_kp * f_At - vec(Q_ * But_prev * ex_prev'))
                    end #if B

                    if estimate_Q
                        Q_S1 += pmodel.Q.D' * pmodel.Q.D
                        Q_S2 += pmodel.Q.D' * vec(phi(t) * (exx - exx1 * At' -
                                    At * exx1' - ex * But_prev' - But_prev * ex' + At * exx_prev * At' +
                                    At * ex_prev * But_prev' + But_prev * ex_prev' * At' + But_prev * But_prev') * phi(t)')
                    end #if Q

                end #if BQ

            end #if BQZ

            if estimate_C | estimate_R | estimate_B | estimate_D
                yt = y[:, t]
                Ct = m.C(t)
                Nt = N(t)
                Dut = m.D(t) * ut
                ey = yt - Nt * (yt - Ct * ex - Dut)

                if estimate_C | estimate_R

                    eyx = Nt * Ct * Pt + ey * ex'

                    if estimate_C
                        pkp = kron(pmodel.C3(t)', pmodel.C1(t))
                        f_Ct = pkp * pmodel.C2.f
                        D_Ct = pkp * pmodel.C2.D
                        C_kp = kron(exx, R_)
                        C_S1 += D_Ct' * C_kp * D_Ct
                        C_S2 += D_Ct' * (vec(R_ * eyx) - C_kp * f_Ct - vec(R_ * Dut * ex'))
                    end #if Z

                    if estimate_R
                        I2 = diagm(!y_notnan[:, t])
                        eyy = I2 * (Nt * HRH(t) + Nt * Ct * Pt * Ct' * Nt') * I2 + ey * ey'
                        R_S1 += pmodel.R.D' * pmodel.R.D
                        R_S2 += pmodel.R.D' * vec(xi(t) * (eyy - eyx * Ct' -
                                    Ct * eyx' - ey * Dut' - Dut * ey' + Ct * exx * Ct' +
                                    Ct * ex * Dut' + Dut * ex' * Ct' + Dut * Dut') * xi(t)')
                    end #if R

                end #if ZR

                if estimate_B

                    Idt = Id(M_t)
                    B_kp = kron(ut', I_nx)
                    A_ = At * At1_
                    fB = B_kp * f_Bt
                    DB = B_kp * D_Bt
                    f_ = At * ft1_ + fB
                    D_ = At * Dt1_ + DB

                    Dt1 = ey - Ct * (I - Idt) * ex - Ct * Idt * (A_ * EX0 + f_) - Dut
                    Dt2 = Ct * Idt * D_
                    Dt3 = ex - At * (I - Idt1) * ex_prev -
                                At * Idt1 * (At1_ * EX0 + ft1_) - fB
                    Dt4 = DB + At * Idt1 * Dt1_

                    B_S1 += Dt4' * Q_ * Dt4 + Dt2' * R_ * Dt2
                    B_S2 += Dt4' * Q_ * Dt3 + Dt2' * R_ * Dt1

                    t <= (pmodel.nx +1) && (M_t *= M)
                    Idt1 = Idt
                    At1_ = A_
                    ft1_ = f_
                    Dt1_ = D_

                end #if B

                if estimate_D
                    pkp = kron(I_nu, pmodel.D1(t))
                    f_Dt = pkp * pmodel.D2.f
                    D_Dt = pkp * pmodel.D2.D
                    D_kp = kron(ut', I_ny)
                    D_S1 += D_Dt' * D_kp' * R_ * D_kp * D_Dt
                    D_S2 += D_Dt' * D_kp' * R_ * (ey - Ct * ex - D_kp * f_Dt)
                end #if D

            end #if CRBD

        end #for

        params.A[:]  = estimate_A ? A_S1 \ A_S2 : T[]
        params.B[:]  = estimate_B ? B_S1 \ B_S2 : T[]
        params.Q[:]  = estimate_Q ? Q_S1 \ Q_S2 : T[]

        params.C[:]  = estimate_C ? C_S1 \ C_S2 : T[]
        params.D[:]  = estimate_D ? D_S1 \ D_S2 : T[]
        params.R[:]  = estimate_R ? R_S1 \ R_S2 : T[]

        params.x1[:] = estimate_x1 ? (pmodel.x1.D' * Linv * pmodel.x1.D) \
                        pmodel.x1.D' * Linv * (exs[:, 1] - pmodel.x1.f) : T[]
        params.P1[:] = estimate_P1 ? (pmodel.P1.D' * pmodel.P1.D) \ pmodel.P1.D' *
                        vec(V[1] + exs[:, 1] * exs[:, 1]' -
                        exs[:, 1]*m.x1' - m.x1*exs[:, 1]' + m.x1*m.x1') : T[]

        #toc()

        return loglik

    end #em_kernel

    fraction_change(ll_prev, ll_curr) = isinf(ll_prev) ?
        1 : 2 * (ll_prev - ll_curr) / (ll_prev + ll_curr + 1e-6)

    ll, ll_prev = Inf, Inf

    while (niter > 0) && (fraction_change(ll_prev, ll) > eps)
        ll_prev = ll
        ll      = em_kernel!(params)
        niter -= 1
    end #while

    (niter == 0) && warn("Parameter estimation timed out - results may have failed to converge")

    return params, pmodel(params)
end #fit

function fit{T}(y::Array{T}, model::StateSpaceModel{T}; u::Array{T}=zeros(size(y,1), model.nu),
              eps::Float64=1e-6, niter::Int=typemax(Int))

    # B, Z, x1 default to parametrized as fully unconstrained
    A, A_params = parametrize_full(model.A(1))
    C, C_params = parametrize_full(model.C(1))
    x1, x1_params = parametrize_full(model.x1)

    # U, A default to fixed
    B, B_params = parametrize_none(model.B(1))
    D, D_params = parametrize_none(model.D(1))

    # Q, R, P1 default to parametrized as diagonal with independent elements - any other values
    #   are ignored / set to zero
    Q, Q_params = parametrize_diag(diag(model.V(1)))
    R, R_params = parametrize_diag(diag(model.W(1)))
    P1, P1_params = parametrize_diag(diag(model.P1))

    pmodel = ParametrizedSSM(A, Q, C, R, x1, P1, B=B, D=D)
    params = SSMParameters(A=A_params, B=B_params, Q=Q_params,
                            C=C_params, D=D_params, R=R_params,
                            x1=x1_params, P1_params)
    fit(y, pmodel, params, u=u, eps=eps, niter=niter)

end #fit


text:
mutable struct DEStats
  nf::Int
  nf2::Int
  nw::Int
  nsolve::Int
  njacs::Int
  nnonliniter::Int
  nnonlinconvfail::Int
  ncondition::Int
  naccept::Int
  nreject::Int
  maxeig::Float64
end

DEStats(x::Int = -1) = DEStats(x, x, x, x, x, x, x, x, x, x, 0.0)

function Base.show(io::IO, s::DEStats)
  println(io, summary(s))
  @printf io "%-50s %-d\n" "Number of function 1 evaluations:" s.nf
  @printf io "%-50s %-d\n" "Number of function 2 evaluations:" s.nf2
  @printf io "%-50s %-d\n" "Number of W matrix evaluations:" s.nw
  @printf io "%-50s %-d\n" "Number of linear solves:" s.nsolve
  @printf io "%-50s %-d\n" "Number of Jacobians created:" s.njacs
  @printf io "%-50s %-d\n" "Number of nonlinear solver iterations:" s.nnonliniter
  @printf io "%-50s %-d\n" "Number of nonlinear solver convergence failures:" s.nnonlinconvfail
  @printf io "%-50s %-d\n" "Number of rootfind condition calls:" s.ncondition
  @printf io "%-50s %-d\n" "Number of accepted steps:"  s.naccept
  @printf io "%-50s %-d" "Number of rejected steps:" s.nreject
  iszero(s.maxeig) || @printf io "\n%-50s %-d" "Maximum eigenvalue recorded:" s.maxeig
end


text:
@testset "DWD climate raw" begin
	temp_out = joinpath(dirname(@__DIR__),"test","output","temp_dwdclimateraw.txt");
	temp_in = joinpath(dirname(@__DIR__),"test","input")
	dwd = dwdclimateraw(DateTime(2012,04),67890,"P",downto=temp_out,
						url=temp_in);
	@test dwd[:datetime][end] == DateTime(2012,04)
	@test dwd[:P][1] == 63
	@test size(dwd) == (1,2)

	timein = collect(DateTime(2012,04):Dates.Month(1):DateTime(2012,05));
	dwdvec = dwdclimateraw(timein,12345,"Po",downto=temp_out,
						url=temp_in);
	@test dwdvec[:datetime] == timein
	@test dwdvec[:Po] == [0,22222]
end

@testset "CDO daily data" begin
	file_in = joinpath(dirname(@__DIR__),"test","input","nndc_climate_cdo.txt");
	df = readgssd(file_in);
	@test names(df) == [:STN,:WBAN,:YEARMODA,:TEMP,:Count_TEMP,:DEWP,:Count_DEWP,
						:SLP,:Count_SLP,:STP,:Count_STP,:VISIB,:Count_VISIB,
						:WDSP,:Count_WDSP,:MXSPD,:GUST,:MAX,:Flag_MAX,:MIN,
						:Flag_MIN,:PRCP,:Flag_PRCP,:SNDP,:FRSHTT];
	x = [875760 99999 20160101 74.3 24 67.9 24 1011.3 6 1009.0 6 5.5 24 4.9 24 12.0 999.9 86.9 " " 66.2 "*" 0.20 "G" 999.9 010010;
	 	 875760 99999 20160102 78.6 24 66.9 24 1013.3 8 1011.0 8 6.2 18 7.9 24 18.1 999.9 87.8 "*" 67.1 " " 0.59 "G" 999.9 000000;
	 	 875760 99999 20160103 78.4 24 68.4 24 1013.0 8 1010.7 8 6.2 18 8.0 24 13.0 999.9 89.4 " " 74.5 " " 0.00 "G" 999.9 000000];
    for i in 1:ncol(df)
		@test df[i] == x[:,i]
	end

	df[:PRCP][2] = 99.99;
	dfc = convertgssd(df);
	@test names(dfc) == [:datetime,:TEMP,:DEWP,:MAX,:MIN,:WDSP,:MXSPD,:PRCP,:STP]
	@test dfc[:datetime] == [DateTime(2016,1,1),DateTime(2016,1,2),DateTime(2016,1,3)]
	s = DataFrame(TEMP = [23.5,25.88889,25.77778], DEWP = (x[:,6].-32)./1.8,
					WDSP = x[:,14]./10.0.*0.514444444444, STP = df[:STP],
					PRCP = x[:,end-3]*2.54*10.0);
	for i in names(s)
		if i != :PRCP
			@test isapprox(dfc[i],s[i], atol=1e-3);
		else
			@test isapprox(dfc[i][1],s[i][1],atol=1e-2)
			@test isapprox(dfc[i][3],s[i][3],atol=1e-2)
			@test isnan(dfc[i][2])
		end
	end
end


text:
#!/usr/bin/env julia

using Logging
@Logging.configure(level=DEBUG)

@debug("test")

max_stars_repo_name: abogushov/lplib

max_issues_repo_name: abogushov/lplib

max_forks_repo_name: abogushov/lplib


text:
"""
Block and braile rendering of julia arrays, for terminal graphics.
"""
module UnicodeGraphics

export blockize, brailize, blockize!, brailize!

"""
    brailize(a, cutoff=0)

Convert an array to a block unicode string, filling values above the cutoff point.
"""
blockize(a, cutoff=0) = blockize!(initblock(size(a)), a, cutoff)

# x and y are inverted: repl rows are columns.
initblock((y, x)) = initblock(y, x)
initblock(y, x) = Array{Char,2}(undef, x + 1, (y - 1) ÷ 2 + 1)

"""
    blockize!(out, a, cutoff=0)

Convert an array to a braile unicode string, filling the `out` array.
Calculation of array dims is a little complicated:
"""
blockize!(out, a, cutoff=0) = join(block_array!(out, a, cutoff))

function block_array!(out, a, cutoff)
    yrange, xrange = axes(a)
    for y in first(yrange):2:last(yrange)
        for x in xrange
            top = checkval(a, y, x, yrange, xrange, cutoff)
            bottom = checkval(a, y + 1, x, yrange, xrange, cutoff)
            if top
                ch = bottom ? '█' : '▀'
            else
                ch = bottom ? '▄' : ' '
            end
            out[x-first(xrange) + 1, (y-first(yrange)) ÷ 2 + 1] = Char(ch)
        end
        # Return after every column
        out[end, (y-first(yrange)) ÷ 2 + 1] = Char('\n')
    end
    # The last character is null
    out[end, end] = 0x00
    out
end

const braile_hex = ((0x01, 0x08), (0x02, 0x10), (0x04, 0x20), (0x40, 0x80))

"""
    brailize(a, cutoff=0)

Convert an array to a braile unicode string, filling values above the cutoff point.
"""
brailize(a, cutoff=0) = brailize!(initbraile(size(a)), a, cutoff)

# x and y are inverted: repl rows are columns.
initbraile((y, x)) = initbraile(y, x)
initbraile(y, x) = Array{Char,2}(undef, (x - 1) ÷ 2 + 2, (y - 1) ÷ 4 + 1)

"""
    brailize!(out, a, cutoff=0)

Convert an array to a braile unicode string, filling the `out` array.
"""
brailize!(out, a, cutoff=0) = join(braile_array!(out, a, cutoff))

function braile_array!(out, a, cutoff)
    yrange, xrange = axes(a)
    for y in first(yrange):4:last(yrange)
        for x in first(xrange):2:last(xrange)
            ch = 0x2800
            for j = 0:3, i = 0:1
                if checkval(a, y+j, x+i, yrange, xrange, cutoff)
                    ch += braile_hex[j % 4 + 1][i % 2 + 1]
                end
            end
            out[(x - first(xrange)) ÷ 2 + 1, (y-first(yrange)) ÷ 4 + 1] = ch
        end
        # Return after every column
        out[end, (y-first(yrange)) ÷ 4 + 1] = Char('\n')
    end
    # The last character is null
    out[end, end] = 0x00
    out
end

checkval(a, y, x, yrange, xrange, cutoff) = begin
    if x <= last(xrange) && y <= last(yrange)
        a[y, x] > cutoff
    else
        false
    end
end

end # module


text:
using AbstractActuators
using Test

@testset "AbstractActuators.jl" begin
    # Write your tests here.
end


text:
module CelesteModSampleTrigger

using ..Ahorn, Maple

@mapdef Trigger "CelesteMod/SampleTrigger" SampleTrigger(
    x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight,
    sampleProperty::Integer=0
)

const placements = Ahorn.PlacementDict(
    "Sample Trigger (CelesteMod)" => Ahorn.EntityPlacement(
        SampleTrigger,
        "rectangle",
    )
)

end

max_stars_repo_name: swoolcock/CelesteModTemplate

max_issues_repo_name: swoolcock/CelesteModTemplate

max_forks_repo_name: swoolcock/CelesteModTemplate


text:
import Base: filter, map, reduce

"""
    map(f, d::DTable) -> DTable

Applies `f` to each row of `d`.
The applied function needs to return a `Tables.Row` compatible object (e.g. `NamedTuple`).

# Examples
```julia
julia> d = DTable((a = [1, 2, 3], b = [1, 1, 1]), 2);

julia> m = map(x -> (r = x.a + x.b,), d)
DTable with 2 partitions
Tabletype: NamedTuple

julia> fetch(m)
(r = [2, 3, 4],)

julia> m = map(x -> (r1 = x.a + x.b, r2 = x.a - x.b), d)
DTable with 2 partitions
Tabletype: NamedTuple

julia> fetch(m)
(r1 = [2, 3, 4], r2 = [0, 1, 2])
```
"""
function map(f, d::DTable)
    chunk_wrap = (_chunk, _f) -> begin
        if isnonempty(_chunk)
            m = TableOperations.map(_f, _chunk)
            Tables.materializer(_chunk)(m)
        else
            _chunk
        end
    end
    chunks = map(c -> Dagger.spawn(chunk_wrap, c, f), d.chunks)
    DTable(chunks, d.tabletype)
end

"""
    map(f, gd::GDTable) -> GDTable

Applies `f` to each row of `gd`.
The applied function needs to return a `Tables.Row` compatible object (e.g. `NamedTuple`).

# Examples
```julia
julia> g = Dagger.groupby(DTable((a=repeat('a':'c', inner=2),b=1:6), 2), :a)
GDTable with 3 partitions and 3 keys
Tabletype: NamedTuple
Grouped by: [:a]

julia> m = map(r -> (a = r.a, b = r.b, c = r.a + r.b), g)
GDTable with 3 partitions and 3 keys
Tabletype: NamedTuple
Grouped by: [:a]

julia> fetch(m)
(a = ['a', 'a', 'c', 'c', 'b', 'b'], b = [1, 2, 5, 6, 3, 4], c = ['b', 'c', 'h', 'i', 'e', 'f'])
```
"""
function map(f, gd::GDTable)
    d = map(f, gd.dtable)
    GDTable(d, gd.cols, gd.index)
end

"""
    reduce(f, d::DTable; cols=nothing, [init]) -> NamedTuple

Reduces `d` using function `f` applied on all columns of the DTable.

By providing the kwarg `cols` as a `Vector{Symbol}` object it's possible
to restrict the reduction to the specified columns.
The reduced values are provided in a NamedTuple under names of reduced columns.

For the `init` kwarg please refer to `Base.reduce` documentation,
as it follows the same principles.

# Examples
```julia
julia> d = DTable((a = [1, 2, 3], b = [1, 1, 1]), 2);

julia> r1 = reduce(+, d)
EagerThunk (running)

julia> fetch(r1)
(a = 6, b = 3)

julia> r2 = reduce(*, d, cols=[:a])
EagerThunk (running)

julia> fetch(r2)
(a = 6,)
```
"""
function reduce(f, d::DTable; cols=nothing::Union{Nothing, Vector{Symbol}}, init=Base._InitialValue())
    # handle empty dtables
    nchunks(d) == 0 && return Dagger.@spawn NamedTuple()

    columns = cols === nothing ? _columnnames_svector(d) : cols

    chunk_reduce_results = _reduce_chunks(f, d.chunks, columns, init=init)

    construct_single_column = (_col, _chunk_results) -> getindex.(fetch.(_chunk_results), _col)
    result_columns = [Dagger.@spawn construct_single_column(c, chunk_reduce_results) for c in columns]

    reduce_result_column = (_f, _c, _init) -> reduce(_f, _c; init=_init)
    reduce_chunks = [Dagger.@spawn reduce_result_column(f, c, deepcopy(init)) for c in result_columns]

    construct_result = (_cols, _vals) -> (; zip(_cols, fetch.(_vals))...)
    Dagger.@spawn construct_result(columns, reduce_chunks)
end

function _reduce_chunks(f, chunks::Vector{Union{EagerThunk, Chunk}}, columns::Vector{Symbol}; init=Base._InitialValue())
    col_in_chunk_reduce = (_f, _c, _init, _chunk) -> reduce(_f, Tables.getcolumn(_chunk, _c); init=_init)

    chunk_reduce = (_f, _chunk, _cols, _init) -> begin
        # TODO: potential speedup enabled by commented code below by reducing the columns in parallel
        v = [col_in_chunk_reduce(_f, c, deepcopy(_init), _chunk) for c in _cols]
        (; zip(_cols, v)...)

        # TODO: uncomment and define a good threshold for parallelization when this get's resolved
        # https://github.com/JuliaParallel/Dagger.jl/issues/267
        # This piece of code (else option) below is causing the issue above
        # when reduce is repeatedly executed or @btime is used.
        # if length(_cols) <= 1
        #     v = [col_in_chunk_reduce(_f, c, _init, _chunk) for c in _cols]
        # else
        #     values = [Dagger.spawn(col_in_chunk_reduce, _f, c, _init, _chunk) for c in _cols]
        #     v = fetch.(values)
        # end
        # (; zip(_cols, v)...)
    end
    chunk_reduce_spawner = (_d, _f, _columns, _init) -> [Dagger.@spawn chunk_reduce(_f, c, _columns, _init) for c in _d]
    Dagger.@spawn chunk_reduce_spawner(chunks, f, columns, init)
end

"""
    reduce(f, gd::GDTable; cols=nothing, prefix="result_", [init]) -> EagerThunk -> NamedTuple

Reduces `gd` using function `f` applied on all columns of the DTable.
Returns results per group in columns with names prefixed with the `prefix` kwarg.
For more information on kwargs see `reduce(f, d::DTable)`

# Examples
```julia
julia> g = Dagger.groupby(DTable((a=repeat('a':'d', inner=2),b=1:8), 2), :a)
GDTable with 4 partitions and 4 keys
Tabletype: NamedTuple
Grouped by: [:a]

julia> fetch(reduce(*, g))
(a = ['a', 'c', 'd', 'b'], result_a = ["aa", "cc", "dd", "bb"], result_b = [2, 30, 56, 12])
```
"""
function reduce(
    f,
    gd::GDTable;
    cols=nothing::Union{Nothing, Vector{Symbol}},
    prefix::String="result_",
    init=Base._InitialValue())

    # handle empty dtables
    nchunks(gd) == 0 && return Dagger.@spawn NamedTuple()

    columns = cols === nothing ? _columnnames_svector(gd) : cols

    chunk_reduce_results = _reduce_chunks(f, gd.dtable.chunks, columns, init=init)

    reduce_result_column = (_f, _c, _init) -> reduce(_f, _c; init=deepcopy(_init))
    construct_single_column = (_col, _chunk_results, _index, _f, _init) -> begin
        r = getindex.(fetch.(_chunk_results), _col)
        [reduce_result_column(_f, getindex.(Ref(r), chunk_indices), _init) for (_, chunk_indices) in _index]
    end
    result_columns = [Dagger.@spawn construct_single_column(c, chunk_reduce_results, gd.index, f, init) for c in columns]

    construct_result = (_keys, _gcols, _columns, _results, _prefix) -> begin
        ks = [col => getindex.(_keys, i) for (i, col) in enumerate(_gcols)]
        rs = [Symbol(_prefix * string(r)) => fetch(_results[i]) for (i, r) in enumerate(_columns)]
        (;ks...,rs...)
    end

    Dagger.@spawn construct_result(keys(gd), grouped_cols(gd), columns, result_columns, prefix)
end

"""
    filter(f, d::DTable) -> DTable

Filter `d` using `f`.
Returns a filtered `DTable` that can be processed further.

# Examples
```julia
julia> d = DTable((a = [1, 2, 3], b = [1, 1, 1]), 2);

julia> f = filter(x -> x.a < 3, d)
DTable with 2 partitions
Tabletype: NamedTuple

julia> fetch(f)
(a = [1, 2], b = [1, 1])

julia> f = filter(x -> (x.a < 3) .& (x.b > 0), d)
DTable with 2 partitions
Tabletype: NamedTuple

julia> fetch(f)
(a = [1, 2], b = [1, 1])
```
"""
function filter(f, d::DTable)
    chunk_wrap = (_chunk, _f) -> begin
        m = TableOperations.filter(_f, _chunk)
        Tables.materializer(_chunk)(m)
    end
    DTable(map(c -> Dagger.spawn(chunk_wrap, c, f), d.chunks), d.tabletype)
end

"""
    filter(f, gd::GDTable) -> GDTable

Filter 'gd' using 'f', returning a filtered `GDTable`.
Calling `trim!` on a filtered `GDTable` will clean up the empty keys and partitions.

# Examples
```julia
julia> g = Dagger.groupby(DTable((a=repeat('a':'d', inner=2),b=1:8), 2), :a)
GDTable with 4 partitions and 4 keys
Tabletype: NamedTuple
Grouped by: [:a]

julia> f = filter(x -> x.a ∈ ['a', 'b'], g)
GDTable with 4 partitions and 4 keys
Tabletype: NamedTuple
Grouped by: [:a]

julia> fetch(f)
(a = ['a', 'a', 'b', 'b'], b = [1, 2, 3, 4])

julia> trim!(f)
GDTable with 2 partitions and 2 keys
Tabletype: NamedTuple
Grouped by: [:a]
```
"""
function filter(f, gd::GDTable)
    d = filter(f, gd.dtable)
    GDTable(d, gd.cols, gd.index)
end


text:
# Importing packages
using Plots
pyplot() #add pyplot from python
using Oceananigans
using Oceananigans.Units

L = 20meters #Length in meters (thus x1)
H = 4meters #Total depth

Ny = Int(2*L/0.1meters) # number of points in y
Nz = Int(H/0.1meters) # number of points in z

grid = RectilinearGrid( #creating grid
    size=(Ny, Nz),
    y=(-L*meters, L*meters),
    z=(-H*meters,0),
    topology=(Flat, Bounded, Bounded) # Solving for y and z where we define as "Bounded". "Flat" will not be solved.
)

closure = SmagorinskyLilly() #Choosing turbulence closure

model = NonhydrostaticModel(grid=grid, closure=closure, buoyancy=SeawaterBuoyancy(), tracers=(:T, :S)) # Setting the model configuration

initial_temperature(x, y, z) = 20 + 0.5*y/abs(y)
set!(model, T=initial_temperature, S=35) # Adding initial conditions to the model setting

simulation = Simulation(model, Δt = 0.025seconds, stop_time = 10minutes) #Creating the simulation

#Setting the format, frequency and path of the output file
simulation.output_writers[:temperature] = JLD2OutputWriter(
                    model, model.tracers, prefix = "../data/lockexchange",
                    schedule=TimeInterval(30second), force = true
)

run!(simulation)

####################################################################################################
T = FieldTimeSeries("../data/lockexchange.jld2", "T") # Reading the output file
x,y,z = nodes(T) #Extracting the grid

#Creating the animation
anim = @animate for (i, t) in enumerate(T.times)
    contour(y,z,T[1,:,:,i]'; xlabel="y [m]", ylabel="z [m]",clim=(19,20.5),c=cgrad(:dense,rev=true), levels=18.9:0.1:20.5, fill=true, linewidth=0.0)
end

#Saving the animation as a GIF
gif(anim, "../img/animation.gif", fps = 9)


text:
using MathProgBase

"
--- PostSolving Utilities ---
add_to_stack!           : function that will add one type of a presolve stack element to the stack.
postsolve_stack_element!: function that will post solve one presolve stack element using the corresponding function.
postsolve_lindep        : function that will post solve one Linear_Dependency element.
return_postsolved       : function that will take in the solution from solver for reduced problem and returns solution for original problem
"

# EMPTY ROW
function add_to_stack!(r::Empty_Row,  independentvar::BitArray{1}, active_constr::BitArray{1}, pstack::Array{Presolve_Stack,1})
    active_constr[r.i] = false
    push!(pstack,r)
end

function postsolve_stack_element!(r::Empty_Row, var_primalsol::Array{Float64,1}, var_dualsol::Array{Float64,1}, constr_primalsol::Array{Float64,1}, constr_dualsol::Array{Float64,1})
    constr_primalsol[r.i] = 0.0
    constr_dualsol[r.i] = 0.0
end

# FREE ROW
function add_to_stack!(r::Free_Row,  independentvar::BitArray{1}, active_constr::BitArray{1}, pstack::Array{Presolve_Stack,1})
    active_constr[r.i] = false
    push!(pstack,r)
end

function postsolve_stack_element!(r::Free_Row, var_primalsol::Array{Float64,1}, var_dualsol::Array{Float64,1}, constr_primalsol::Array{Float64,1}, constr_dualsol::Array{Float64,1})
    postsolve_lindep!(constr_primalsol,l)
    constr_dualsol[r.i] = 0.0
end

# SINGLETON ROW EQUALITY
function add_to_stack!(r::Singleton_Row_Equality,  independentvar::BitArray{1}, active_constr::BitArray{1}, pstack::Array{Presolve_Stack,1})
    active_constr[r.i] = false
    push!(pstack,r)
end

function postsolve_stack_element!(r::Singleton_Row_Equality, var_primalsol::Array{Float64,1}, var_dualsol::Array{Float64,1}, constr_primalsol::Array{Float64,1}, constr_dualsol::Array{Float64,1})
    var_primalsol[r.j] = r.x_j
    constr_primalsol[r.i] = r.y_i
    constr_dualsol[r.row] = (var_dualsol[r.j] * r.x_j)/ r.y_i
    var_dualsol[r.j] = 0.0
end

# SINGLETON ROW INEQUALITY
function add_to_stack!(r::Singleton_Row_Inequality,  independentvar::BitArray{1}, active_constr::BitArray{1}, pstack::Array{Presolve_Stack,1})
    #active_constr[r.i] = false
    push!(pstack,r)
end

function postsolve_stack_element!(r::Singleton_Row_Inequality, var_primalsol::Array{Float64,1}, var_dualsol::Array{Float64,1}, constr_primalsol::Array{Float64,1}, constr_dualsol::Array{Float64,1})
    l = r.l
    constr_primalsol[r.i] = l.value
    for j in 1:length(l.vec1)
        constr_primalsol[r.i] += l.vec2[j]*constr_dualsol[l.vec1[j]]
    end

    # col bound
    if(var_primalsol[r.j] = r.l_new)
        col_bound = 1
    elseif (l_new <= var_primalsol[r.j] && var_primalsol[r.j] <= u_new)
        col_bound = 2
    elseif (var_primalsol[r.j] == r.u_new)
        col_bound = 3
    end

    row_bound = r.row_bound

    if (col_bound == 2)
        constr_dualsol[r.i] = 0.0
    end

    if(matval > 0)
        if (col_bound == 1 && (row_bound == 1 || row_bound == 2))
            constr_dualsol[r.i] = 0.0
        elseif (col_bound == 1 && (row_bound == 3 || row_bound == 4))
            constr_dualsol[r.i] = var_dualsol[r.j] / r.matval
        elseif (col_bound == 3 && (row_bound == 1 || row_bound == 3))
            constr_dualsol[r.i] = 0.0
        elseif (col_bound == 3 && (row_bound == 2 || row_bound == 4))
            constr_dualsol[r.i] = var_dualsol[r.j] / r.matval
        end
    else
        if (col_bound == 1 && (row_bound == 1 || row_bound == 2))
            constr_dualsol[r.i] = 0.0
        elseif (col_bound == 1 && (row_bound == 3 || row_bound == 4))
            constr_dualsol[r.i] = -var_dualsol[r.j] / r.matval
        elseif (col_bound == 3 && (row_bound == 1 || row_bound == 3))
            constr_dualsol[r.i] = 0.0
        elseif (col_bound == 3 && (row_bound == 2 || row_bound == 4))
            constr_dualsol[r.i] = -var_dualsol[r.j] / r.matval
        end
    end
end

# Forcing ROW
function add_to_stack!(r::Forcing_Row,  independentvar::BitArray{1}, active_constr::BitArray{1}, pstack::Array{Presolve_Stack,1})
    active_constr[r.i] = false
    push!(pstack,r)
end

function postsolve_stack_element!(r::Forcing_Row, var_primalsol::Array{Float64,1}, var_dualsol::Array{Float64,1}, constr_primalsol::Array{Float64,1}, constr_dualsol::Array{Float64,1})
    # lastcol , val
    lastcol = -1
    val = 0
    for(j in 1:length(r.col_ind))
        dual = var_dualsol[r.col_ind[j]]
        temp = dual/r.mat_val[r.col_ind[j]]

        if((r.col_bound == -1) && (dual < 0 && val < temp))
            lastcol = r.col_ind[j]
            val = temp
        elseif((r.col_bound == 1) && (dual > 0 && val < temp))
            lastcol = r.col_ind[j]
            val = temp
        end
    end

    constr_primalsol[r.i] = r.fixed_val

    if(lastcol == -1)
        constr_dualsol[r.i] = 0.0
    else
        #constr_dualsol[r.i] = var_dualsol[lastcol] /r.mat_val[lastcol]
        constr_dualsol[r.i] = val
        for(j in 1:length(r.col_ind))
            if(j!=lastcol)
                var_dualsol[j] -= (constr_dualsol[r.i])*r.mat_val[j]
            end
        end
        var_dualsol[lastcol] = 0.0
    end

end

# EMPTY COL
function add_to_stack!(c::Empty_Col,  independentvar::BitArray{1}, active_constr::BitArray{1}, pstack::Array{Presolve_Stack,1})
    independentvar[c.j] = false
    push!(pstack,c)
end

function postsolve_stack_element!(c::Empty_Col, var_primalsol::Array{Float64,1}, var_dualsol::Array{Float64,1}, constr_primalsol::Array{Float64,1}, constr_dualsol::Array{Float64,1})
    var_primalsol[c.j] = x_j
    var_dualsol[c.j] = c_j
end

# FIXED COL
function add_to_stack!(c::Fixed_Col,  independentvar::BitArray{1}, active_constr::BitArray{1}, pstack::Array{Presolve_Stack,1})
    independentvar[c.j] = false
    push!(pstack,c)
end

function postsolve_stack_element!(c::Fixed_Col, var_primalsol::Array{Float64,1}, var_dualsol::Array{Float64,1}, constr_primalsol::Array{Float64,1}, constr_dualsol::Array{Float64,1})
    var_primalsol[c.j] = c.x_j

    l = c.l
    var_dualsol[c.j] = l.value
    for i in 1:length(l.vec1)
        var_dualsol[c.j] += l.vec2[i]*constr_dualsol[l.vec1[i]]
        constr_primalsol[l.vec1[i]] += l.vec2[i]*c.x_j
    end
end

"
# SINGLETON COL
function add_to_stack!(c::Fixed_Col,  independentvar::BitArray{1}, active_constr::BitArray{1}, pstack::Array{Presolve_Stack,1})
    independentvar[c.j] = false
    push!(pstack,c)
end

function postsolve_stack_element!(c::Fixed_Col, var_primalsol::Array{Float64,1}, var_dualsol::Array{Float64,1}, constr_primalsol::Array{Float64,1}, constr_dualsol::Array{Float64,1})
    var_primalsol[c.j] = c.x_j

    l = c.l
    var_dualsol[c.j] = l.value
    for i in 1:length(l.vec1)
        var_dualsol[c.j] += l.vec2[i]*constr_dualsol[l.vec1[i]]
        constr_primalsol[l.vec1[i]] += l.vec2[i]*c.x_j
    end
end
"

# General Linear Dependency postsolve
function postsolve_lindep!(post_solvedX::Array{Float64,1}, l::Linear_Dependency)
    post_solvedX[l.index] = l.value

    for i in 1:length(l.vec1)
        post_solvedX[l.index] += l.vec2[i]*post_solvedX[l.vec1[i]]
        #println("made postsolved at $(l.index) to value $(post_solvedX[l.index])")
    end
end

function return_postsolved(vpsol::Array{Float64,1}, vdsol::Array{Float64,1}, cpsol::Array{Float64,1}, cdsol::Array{Float64,1}, independent_var::BitArray{1}, active_constr::BitArray{1}, pstack :: Array{Presolve_Stack,1})
    n = length(independent_var)
    m = length(active_constr)

    var_primalsol = zeros(n)
    constr_primalsol = zeros(m)
    var_dualsol = zeros(n)
    constr_dualsol = zeros(m)

    newcols = find(independent_var)
    newrows = find(active_constr)

    for i in 1:length(newcols)
        var_primalsol[newcols[i]] = vpsol[i]
        var_dualsol[newcols[i]] = vdsol[i]
    end

    for i in 1:length(newrows)
        constr_primalsol[newrows[i]] = cpsol[i]
        constr_dualsol[newrows[i]] = cdsol[i]
    end

    sol = [var_primalsol,var_dualsol,constr_primalsol,constr_dualsol]

    #NOTE. What is the mapping? from different sols to flags. comment it.
    # 1 - var_primalsol
    # 2 - var_dualsol
    # 3 - constr_primalsol
    # 4 - constr_dualsol

    for i in reverse(collect(1:length(pstack)))
        postsolve_stack_element!(pstack[i],sol[1],sol[2],sol[3],sol[4])
    end
    return sol
end


text:
module Models

export
    IncompressibleModel, NonDimensionalIncompressibleModel, ShallowWaterModel,
    HydrostaticFreeSurfaceModel, VectorInvariant,
    ExplicitFreeSurface, ImplicitFreeSurface,
    HydrostaticSphericalCoriolis, VectorInvariantEnstrophyConserving,
    PrescribedVelocityFields

using Oceananigans: AbstractModel

abstract type AbstractIncompressibleModel{TS} <: AbstractModel{TS} end

include("IncompressibleModels/IncompressibleModels.jl")
include("HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl")
include("ShallowWaterModels/ShallowWaterModels.jl")

using .IncompressibleModels: IncompressibleModel, NonDimensionalIncompressibleModel
using .ShallowWaterModels: ShallowWaterModel

using .HydrostaticFreeSurfaceModels:
    HydrostaticFreeSurfaceModel, VectorInvariant,
    ExplicitFreeSurface, ImplicitFreeSurface,
    HydrostaticSphericalCoriolis, VectorInvariantEnstrophyConserving,
    PrescribedVelocityFields

end


text:
using DrWatson
@quickactivate
using ArgParse
using GenerativeAD
import StatsBase: fit!, predict
using StatsBase
using BSON
using Flux
using GenerativeModels

s = ArgParseSettings()
@add_arg_table! s begin
	"max_seed"
		default = 1
		arg_type = Int
		help = "max_seed"
	"category"
		default = "wood"
		arg_type = String
		help = "category"
	"contamination"
		arg_type = Float64
		help = "contamination rate of training data"
		default = 0.0
end
parsed_args = parse_args(ARGS, s)
@unpack category, max_seed, contamination = parsed_args

modelname = "fAnoGAN"

"""
	sample_params()

Should return a named tuple that contains a sample of model parameters.
"""
function sample_params()
	# first sample the number of layers
	nlayers = rand(2:4)
	kernelsizes = reverse((3,5,7,9)[1:nlayers])
	channels = reverse((16,32,64,128)[1:nlayers])
	scalings = reverse((1,2,2,2)[1:nlayers])

	par_vec = (
		2 .^(3:8),
		10f0 .^(-4:-3),
		10f0 .^(-4:-3),
		2 .^ (3:5),
		["relu", "swish", "tanh"],
		1:Int(1e8),
		[true, false],
		[0.001,0.005,0.01,0.05,0.1]
	)
	argnames = (
		:zdim,
		:lr_gan,
		:lr_enc,
		:batchsize,
		:activation,
		:init_seed,
		:batchnorm,
		:weight_clip,
	)
	parameters = (;zip(argnames, map(x->sample(x, 1)[1], par_vec))...)
	return merge(parameters, (nlayers=nlayers, kernelsizes=kernelsizes,
		channels=channels, scalings=scalings))
end

"""
	fit(data, parameters)

This is the most important function - returns `training_info` and a tuple or a vector of tuples `(score_fun, final_parameters)`.
`training_info` contains additional information on the training process that should be saved, the same for all anomaly score functions.
Each element of the return vector contains a specific anomaly score function - there can be multiple for each trained model.
Final parameters is a named tuple of names and parameter values that are used for creation of the savefile name.
"""
function fit(data, parameters)
	# construct model - constructor should only accept kwargs
	all_parameters = merge(
		parameters,
		(
			idim = size(data[1][1])[1:3],
			kappa = 1f0,
			weight_clip=0.1,
			lr_gan = 0.00005,
			lr_enc = 0.001,
			batch_size=32,
			patience=10,
			check_every=30,
			iters=10000,
			mtt=82800,
			n_critic=5,
			usegpu=true
		)
	)

	# construct model - constructor should only accept kwargs
	model = GenerativeAD.Models.fanogan_constructor(;idim = size(data[1][1])[1:3], parameters...) |> gpu

	max_iter = 50 # this should be enough

	# fit train data

	try
		global info, fit_t, _, _, _ = @timed fit!(model, data, all_parameters)
	catch e
		# return an empty array if fit fails so nothing is computed
		@info "Failed training due to \n$e"
		return (fit_t = NaN, history=nothing, npars=nothing, model=nothing), []
	end
	model = info[1]

	# construct return information - put e.g. the model structure here for generative models
	training_info = (
		fit_t = fit_t,
		history = info[2],
		npars = info[3],
		model = model |> cpu
		)

	# now return the different scoring functions
	training_info, [(x -> GenerativeAD.Models.anomaly_score_gpu(model, x), parameters)]
end

#___________________________________________________________________________________________________________________
if abspath(PROGRAM_FILE) == @__FILE__
	try_counter = 0
	max_tries = 10*max_seed
	cont_string = (contamination == 0.0) ? "" : "_contamination-$contamination"
	while try_counter < max_tries
		parameters = sample_params()

		for seed in 1:max_seed
			savepath = datadir("experiments/images_mvtec$cont_string/$(modelname)/$(category)/ac=1/seed=$(seed)")
			mkpath(savepath)

			# get data
			data = GenerativeAD.load_data("MVTec-AD", seed=seed, category=category,
					contamination=contamination, img_size=128)

			@info "Trying to fit $modelname on $category with parameters $(parameters)..."
			@info "Train/validation/test splits: $(size(data[1][1], 4)) | $(size(data[2][1], 4)) | $(size(data[3][1], 4))"
			@info "Number of features: $(size(data[1][1])[1:3])"

			# here, check if a model with the same parameters was already tested
			if GenerativeAD.check_params(savepath, parameters) & GenerativeAD.Models.check_scaling(size(data[1][1])[1:3], parameters.scalings)
				# fit model
				training_info, results = fit(data, parameters)
				# saving model separately
				if training_info.model !== nothing
					tagsave(joinpath(savepath, savename("model", parameters, "bson", digits=5)),
						Dict("model"=>training_info.model),
						safe = true)
					training_info = merge(training_info, (model = nothing,))
				end
				save_entries = merge(training_info, (modelname = modelname, seed = seed,
					category = category, dataset = "MVTec-AD_$category",
					contamination=contamination))
				# now loop over all anomaly score funs
				for result in results
					GenerativeAD.experiment(result..., data, savepath; save_entries...)
				end
				global try_counter = max_tries + 1
			else
				@info "Model already present, sampling new hyperparameters..."
				global try_counter += 1
			end
		end
	end
	(try_counter == max_tries) ? (@info "Reached $(max_tries) tries, giving up.") : nothing
end


text:
@testset "Event Probabilities" begin

    @testset "Basic Test" begin

        bg = 0.5
        kappa = 0.5
        kern = Distributions.Exponential(1/0.5)
        events = [1,2,3]
        startT = 4
        endT = 5
        prob = HawkesProcesses.event_probability(startT, endT, events, bg, kappa, kern)
        @test isa(prob, Number)
    end

    @testset "Value Test" begin
        bg = 0.5
        kappa = 0.5
        kern = Distributions.Exponential(1/0.5)
        events = [1,2,3]
        startT = 400
        endT = 500
        prob = HawkesProcesses.event_probability(startT, endT, events, bg, kappa, kern)
        @test prob == 50
    end

end


text:
### ============== ### ============== ###
##    3D Simple Vicksek Model          ##
##    Grid version (particle search)   ##
##    Martin Zumaya Hernandez          ##
##    30 / 06 / 2017                   ##
### ============== ### ============== ###

### ============== ### ============== ###

using CollectiveDynamics.SVM3D

### =============== ### =============== ###
###   DEFINITION OF INITIAL PARAMETERS  ###
### =============== ### =============== ###

N   = parse(Int, ARGS[1]) # number of particles
ρ   = parse(Float64, ARGS[2]) # density
η   = parse(Float64, ARGS[3]) # noise intensity
T   = parse(Int, ARGS[4]) # integration time steps
rep = parse(Int, ARGS[5]) # ensemble index

### =============== ### =============== ###

v0 = 1.0 # particle's speed
dt = 1.0 # integration time step

### =============== ### =============== ###

M = 10 # number of boxes per dimension

l = 0.5 # interaction range is double the distance each particle moves in one time step
# l = 0.1 # interaction range is ten times the distance each particle moves in one time step

r0 = (v0 * dt) / l # local interaction range

L  = cbrt(N / ρ) # size of box

### =============== ### =============== ###

cell_size = step(linspace(0., L , M))

### ============== ### ============== ###

box = Box(L, M)
flock = Flock(N, L, dt, v0, r0, η)

### ============== ### ============== ###

output_path = set_output_data_structure_vsk("SVM_3D", N, ρ)

pos_file = open(output_path * "/pos_$(rep).dat", "w+")
vel_file = open(output_path * "/vel_$(rep).dat", "w+")

times = [convert(Int, exp10(i)) for i in 0:T]

for i in 1:(length(times) - 1)

    if i > 1

        for t in (times[i]+1):times[i+1]

            evolve_system(flock, box, cell_size)

            if t % times[i] == 0 || t % times[i-1] == 0
                println("//////// ", t)
                write(pos_file, vcat(flock.pos...))
                write(vel_file, vcat(flock.vel...))
            end
        end

    else

        for t in (times[i]+1):times[i+1]

            evolve_system(flock, box, cell_size)

            if t % times[i] == 0
                println("//////// ", t)
                write(pos_file, vcat(flock.pos...))
                write(vel_file, vcat(flock.vel...))
            end
        end

    end

end

close(pos_file)
close(vel_file)

println("Done all")


text:
"Create a 2D square using Quad4 elements."
function buildSquare(xmin::Real, xmax::Real, ymin::Real, ymax::Real, n_elems_x::Integer, n_elems_y::Integer)
    startLog(main_perf_log, "buildSquare()")

    # Create the Mesh
    mesh = Mesh()

    # Set the sideset / nodeset numberings:
    boundary_info = mesh.boundary_info

    boundary_info.sidesets = [1,2,3,4]
    boundary_info.nodesets = [1,2,3,4]

    for bid in [1,2,3,4]
        boundary_info.side_list[bid] = Array{ElemSidePair}(0)
        boundary_info.node_list[bid] = Array{Node}(0)
    end

    # First: create all of the Nodes
    n_nodes_x = n_elems_x + 1
    n_nodes_y = n_elems_y + 1

    nodes = Array{Node}(n_nodes_x * n_nodes_y)

    x_range = xmax - xmin
    y_range = ymax - ymin

    x_increment = x_range / n_elems_x
    y_increment = y_range / n_elems_y

    node_id = 1
    for y_idx = 1:n_nodes_y
        # y position
        y = ymin + ((y_idx-1) * y_increment)

        for x_idx = 1:n_nodes_x
            # x position
            x = xmin + ((x_idx-1) * x_increment)

            # Create the node
            node = Node{2}(node_id, Vec{2}((x,y)), [], invalid_processor_id)

            # Save it off
            nodes[node_id] = node

            # Add this Node to any relevant nodesets...

            # Bottom
            if y_idx == 1
                push!(boundary_info.node_list[1], node)
            end

            # Right
            if x_idx == n_nodes_x
                push!(boundary_info.node_list[2], node)
            end

            # Top
            if y_idx == n_nodes_y
                push!(boundary_info.node_list[3], node)
            end

            # Left Side
            if x_idx == 1
                push!(boundary_info.node_list[4], node)
            end

            node_id += 1
        end
    end

    # Now create Elements
    elements = Array{Element}(n_elems_x * n_elems_y)

    element_id = 1

    # Idea taken from libMesh mesh_generation.C
    idx(x_idx, y_idx) = ( (x_idx-1) + (y_idx-1) * (n_nodes_x) ) + 1

    for y_idx = 1:n_elems_y
        for x_idx = 1:n_elems_x
            # Create the Element
            element = Element(element_id,
                              [nodes[idx(x_idx, y_idx)],
                               nodes[idx(x_idx+1, y_idx)],
                               nodes[idx(x_idx+1, y_idx+1)],
                               nodes[idx(x_idx, y_idx+1)]],
                              [],
                              invalid_processor_id)

            # Save it off
            elements[element_id] = element

            # Add it to any relevant sidesets...

            # Bottom
            if y_idx == 1
                push!(boundary_info.side_list[1], ElemSidePair(element, 1))
            end

            # Right
            if x_idx == n_elems_x
                push!(boundary_info.side_list[2], ElemSidePair(element, 2))
            end

            # Top
            if y_idx == n_elems_y
                push!(boundary_info.side_list[3], ElemSidePair(element, 3))
            end

            # Left Side
            if x_idx == 1
                push!(boundary_info.side_list[4], ElemSidePair(element, 4))
            end

            element_id += 1
        end
    end

    mesh.elements = elements
    mesh.nodes = nodes

    initialize!(mesh)

    stopLog(main_perf_log, "buildSquare()")

    return mesh
end


text:
module LatticeDiracOperators
    using Wilsonloop
    using Gaugefields
    using Requires
    using AlgRemez_jll

    import Gaugefields:add_U!
    import Gaugefields:Abstractfields,clear_U!
    import Gaugefields:println_verbose_level1,println_verbose_level2,println_verbose_level3,Verbose_print

    # Write your package code here.
    include("./rhmc/AlgRemez.jl")
    include("./rhmc/rhmc.jl")
    #include("./cgmethods.jl")

    include("Diracoperators.jl")
        include("./SakuraiSugiura/SSmodule.jl")

    import .Dirac_operators:Initialize_pseudofermion_fields,Dirac_operator,gauss_distribution_fermion!,
            Initialize_WilsonFermion,Initialize_4DWilsonFermion,
            DdagD_operator,solve_DinvX!,FermiAction,
            shift_fermion,bicg,bicgstab,
            sample_pseudofermions!,gauss_sampling_in_action!,evaluate_FermiAction,calc_UdSfdU,calc_UdSfdU!,bicgstab,Z2_distribution_fermion!,
            Wilson_Dirac_operator_evenodd,calc_p_UdSfdU!,Wilson_GeneralDirac_operator,set_wing_fermion!,eigensystem,eigensystem_old,
            construct_sparsematrix,initialize_Adjoint_fermion,calc_dSfdU!,
            Wilson_Dirac_operator_faster,Dx!,Ddagx!,setvalue_fermion!,setindex_global!,
            uniform_distribution_fermion!,γ5D,convert_to_normalvector,save_fermionfield,load_fermionfield!

    export Initialize_pseudofermion_fields,Dirac_operator,gauss_distribution_fermion!,cg,bicg
    export Initialize_WilsonFermion,Initialize_4DWilsonFermion
    export DdagD_operator,solve_DinvX!,FermiAction
    export shift_fermion
    export WilsonFermion_4D_wing
    export sample_pseudofermions!,gauss_sampling_in_action!,evaluate_FermiAction,calc_UdSfdU,calc_UdSfdU!,bicgstab,calc_p_UdSfdU!
    export Wilson_Dirac_operator_evenodd,
            Wilson_GeneralDirac_operator,set_wing_fermion!
    export println_verbose_level1,println_verbose_level2,println_verbose_level3,Verbose_print
    export bicg,bicgstab
    export eigensystem,calc_dSfdU!,eigensystem_old
    export construct_sparsematrix,initialize_Adjoint_fermion,Z2_distribution_fermion!
    export Wilson_Dirac_operator_faster
    export setvalue_fermion!
    export setindex_global!
    export uniform_distribution_fermion!,γ5D
    export convert_to_normalvector
    export save_fermionfield,load_fermionfield!
end


text:
# To run this script:
# pkg> build GoogleSheets

println("Instalation of GoogleSheets python package dependencies")

using PyCall
@pyimport pip
pip.main(["install","google-api-python-client","google-auth-httplib2","google-auth-oauthlib"])


text:
struct BranchIterator
    bs::Vector{sideal}
    vars::Vector{String}
    auxvars::Vector{String}
    length::Int
end

function BranchIterator(bs::Vector{Vector{MPolyElem}}, vars::Vector{Symbol}, auxvars::Vector{Symbol})
    vs = [map(Recurrences.initvar, vars); vars; auxvars]
    R, _ = Singular.PolynomialRing(Nemo.QQ, map(string, vs))
    σ = IdMap(R)
    ideals = [Singular.Ideal(R, map(σ, b)) for b in bs]
    @debug "Ideals in branch iterator" ideals
    len = length(ideals) == 1 ? 1 : length(ideals) * (length(vars) + 1)
    BranchIterator(ideals, map(string, vars), map(string, auxvars), len)
end

nbranches(iter::BranchIterator) = length(iter.bs)

splitgens(R, lss...) = let sf = Iterators.Stateful(gens(R))
    (collect(Iterators.take(sf, length(ls))) for ls in lss)
end

function stateideal(iter::BranchIterator, I::sideal, state::Int)
    is = [Recurrences.initvar(v, 0) for v in iter.vars]
    ms = [Recurrences.initvar(v, state + 1) for v in iter.vars]
    fs = [Recurrences.initvar(v, state + 2) for v in iter.vars]
    as = iter.auxvars
    R, _ = Singular.PolynomialRing(Nemo.QQ, map(string, [ms; fs; as; is]))
    mss, _, ass, _ = splitgens(R, ms, fs, as, is)
    σ = IndexMap(R)
    Singular.Ideal(R, [σ(g; force=true) for g in gens(I)]), [mss; ass]
end

function Base.iterate(iter::BranchIterator)
    if length(iter) == 0
        return nothing
    # elseif length(iter) == 1
    #     I = iter.bs[1]
    #     _, _, as = splitgens(base_ring(I), iter.vars, iter.vars, iter.auxvars)
    #     return (iter.bs[1], as), 1
    end
    stateideal(iter, iter.bs[1], 0), 1
end

function Base.iterate(iter::BranchIterator, state)
    if state >= length(iter)
        return nothing
    end
    i = state % length(iter.bs) + 1
    stateideal(iter, iter.bs[i], state), state + 1
end

Base.IteratorSize(::Type{BranchIterator}) = Base.HasLength()
Base.length(iter::BranchIterator) = iter.length

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

function initial_ideal(vars::Vector{String}, init::ValueMap)
    is = [Recurrences.initvar(v, 0) for v in vars]
    fs = [Recurrences.initvar(v, 1) for v in vars]

    R, svars = Singular.PolynomialRing(Nemo.QQ, [is; fs])
    iss, fss = splitgens(R, is, fs)

    iss = [haskey(init, Symbol(v)) ? Nemo.QQ(init[Symbol(v)]) : iss[i] for (i, v) in enumerate(vars)]

    basis = [x - y for (x, y) in zip(iss, fss)]
    Singular.Ideal(R, basis)
end

function invariants(bs::Vector{Vector{MPolyElem}}, vars::Vector{Symbol}, auxvars::Vector{Symbol}, init::ValueMap)
    biter = BranchIterator(bs, vars, auxvars)
    # if length(biter) == 1
    #     I, elim = first(biter)
    #     I = Singular.eliminate(I, elim...)
    #     return I
    # end
    fixedpoint(biter, map(string, vars), init)
end

function fixedpoint(biter::BranchIterator, vars::Vector{String}, init::ValueMap)
    is = [Recurrences.initvar(v, 0) for v in vars]
    vs = [vars; is]
    R, _ = Singular.PolynomialRing(Nemo.QQ, vs)
    bcount = nbranches(biter)
    I = Singular.Ideal(R)
    T = initial_ideal(vars, init)
    for (i, (J, elim)) in enumerate(biter)
        S = base_ring(J)
        T = imap(T, S) # map old ideal to new ring via variable name
        T = Singular.eliminate(T + J, elim...)

        if i % bcount == 0
            fs = [Recurrences.initvar(v, i + 1) for v in vars]
            RR, _ = Singular.PolynomialRing(Nemo.QQ, [fs; is])
            II = std(fetch(imap(T, RR), R)) # map current ideal to final ring
            if bcount == 1
                return II
            end
            @debug "Check if ideals are equal" T II I std(I) R isequal(I, II)
            if isequal(I, II)
                return I
            end
            I = II
        end
    end
    error("Fixed-point computation failed. Should not happen.")
end

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

function invariants(cs::Vector{Vector{T}}, _vars::Vector{Symbol}, init::ValueMap) where {T <: Recurrences.ClosedForm}
    __vars = [Set(first(c) for c in branch) for branch in cs]
    vars = [x for x in intersect(__vars...)]
    removed = setdiff(Set(_vars), Set(vars))
    if !isempty(removed)
        @debug "No closed form found for $removed"
    end
    ps = Vector{MPolyElem}[]
    as = Symbol[]
    for c in cs
        p, auxvars = polys(c, vars)
        push!(ps, p)
        push!(as, auxvars...)
    end
    invariants(ps, vars, as, init)
end

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

AbstractAlgebra.isconstant(x::FracElem) = AbstractAlgebra.isconstant(numerator(x)) && AbstractAlgebra.isconstant(denominator(x))
AbstractAlgebra.isconstant(x::fmpq) = true

get_constant(x::PolyElem) = get_constant(AbstractAlgebra.coeff(x, 0))
get_constant(x::MPolyElem) = iszero(x) ? zero(base_ring(x)) : get_constant(AbstractAlgebra.coeff(x, 1))
get_constant(x::FracElem) = get_constant(numerator(x)) // get_constant(denominator(x))
get_constant(x::RingElem) = x
get_constant(x::FieldElem) = x

AbstractAlgebra.show_minus_one(::Type{T}) where T <: FracElem = false

function polys(cs::Vector{ClosedForm}, _all_vars::Vector{Symbol})

    all_vars = copy(_all_vars)

    geometrics = Dict{fmpq,Symbol}()
    factorials = Dict{fmpq,Symbol}()

    function map_geom(x::fmpq)
        isone(x) && return x
        haskey(geometrics, x) && return geometrics[x]
        s = Recurrences.gensym_unhashed(:g)
        push!(geometrics, x=>s)
        return s
    end

    function map_fact(x::fmpq_poly)
        R = parent(x)
        v = gen(parent(x))
        c = x - v
        @assert isconstant(c) c
        c = Nemo.coeff(c, 0)
        for (k, s) in factorials
            d = k - c
            iszero(d) && return s, one(R)//one(R)
            if isone(denominator(d))
                n = numerator(d)
                p = prod(v + k + i for i in 1:abs(n))
                q = R(prod(c + i for i in 0:abs(n)-1))
                @debug "map_fact" n p q factorials
                if n < 0
                    return s, p // q
                else
                    return s, q // p
                end
            end
        end
        s = Recurrences.gensym_unhashed(:f)
        @assert isconstant(c)
        push!(factorials, c=>s)
        return s, one(R)//one(R)
    end

    function _poly(t::Recurrences.HyperTerm)
        c = Recurrences.coeff(t)
        g = Recurrences.geom(t)
        f = Recurrences.fact(t)
        n, d = numerator(f), denominator(f)
        # make fmpq_poly, needed for factorization
        nn = map_coeffs(get_constant, n)
        dd = map_coeffs(get_constant, d)

        function __factor(__x::fmpq_poly)
            __f = Nemo.factor(__x)
            c *= get_constant(unit(__f))
            __vs = Pair{Symbol,fmpz}[]
            for (__fac, __mul) in __f
                # poly should be monic
                lcoeff = Nemo.coeff(__fac, 1)
                __fac = __fac * (1//lcoeff)
                c *= lcoeff
                __var, __coeff = map_fact(__fac)
                __coeff = change_base_ring(base_ring(base_ring(c)), __coeff)
                c *= __coeff
                push!(__vs, __var=>__mul)
            end
            __vs
        end
        @assert AbstractAlgebra.isconstant(g)
        facts = __factor(nn)=>__factor(dd)
        c, map_geom(get_constant(g)), facts
    end

    function _poly(c::ClosedForm)
        first(c), [_poly(t) for t in Recurrences.terms(last(c))]
    end

    res = map(_poly, cs)

    final_vars = map(string, all_vars)
    init_vars = map(Recurrences.initvar, final_vars)
    rec_var = string(var(base_ring(last(first(cs)))))

    mvars = String[rec_var]
    append!(mvars, final_vars)
    append!(mvars, init_vars)
    append!(mvars, string(v) for v in values(geometrics))
    append!(mvars, string(v) for v in values(factorials))
    R, _ = Nemo.PolynomialRing(Nemo.QQ, collect(mvars))
    σ = IdMap(R)
    lcm_global = one(R)

    function _to_mpoly(v::Symbol, terms)
        filter!(x->x!=v, all_vars)
        lcm = one(R)
        res = zero(R)
        for (coeff, geom, (fact_num, fact_den)) in terms
            c, g = σ(coeff), σ(geom)
            fn = reduce(*, σ(first(x)) for x in fact_num; init=one(R))
            fd = reduce(*, σ(first(x)) for x in fact_den; init=one(R))
            AbstractAlgebra.lcm(lcm, denominator(c)*fd)
            res += c*g*(fn//fd)
        end
        AbstractAlgebra.lcm(lcm_global, lcm)
        rhs = lcm*res
        # TODO: compute saturated ideal
        denominator(rhs)*lcm*σ(v) - numerator(rhs)
    end

    polys = map(x->_to_mpoly(x...), res)
    for v in all_vars
        push!(polys, σ(v) - σ(Recurrences.initvar(v)))
    end

    if length(geometrics) > 0
        geoms = [replace(string(g), "//"=>"/") for g in keys(geometrics)]
        gvars = collect(values(geometrics))
        I = dependencies(sideal, geoms; variables = gvars)
        if !isnothing(I)
            aRing, _ = Singular.AsEquivalentAbstractAlgebraPolynomialRing(base_ring(I))
            for p in gens(I)
                push!(polys, σ(aRing(p); force=true))
            end
        end
    end
    @debug "Closed form polynomials" polys
    polys, Iterators.flatten(((v for v in values(geometrics)), (v for v in values(factorials)), [Symbol(rec_var)]))
end


text:
using DataStructures

function kmer_count(infile::String, k::Int)
    sequences = read_sequences(infile)

    counts = DefaultDict{String, Int8}(0)
    for seq in sequences
        for i = 1:length(seq)-k+1
            kmer = seq[i : i+k-1]
            counts[kmer] += 1
        end
    end

    for kmer in eachindex(counts)
        println(kmer, ' ', counts[kmer])
    end
end

function read_sequences(input_filename::String)
    sequences = String[]
    inseq = false
    currentseq = [""]
    input = open(input_filename, "r")
    for line in eachline(input)
        if line[1] == '>'
            if inseq
                push!(sequences, join(currentseq))
                currentseq = [""]
            end
            inseq = false
        else
            inseq = true
            push!(currentseq, line)
        end
    end
    if inseq
        push!(sequences, join(currentseq))
    end
    return sequences
end

input_filename = "input.fa"
k = 11
nargs = length(ARGS)
if nargs > 0
    input_filename = ARGS[1]
end
if nargs > 1
    k = parse(Int8, ARGS[2])
end
kmer_count(input_filename, k)


text:
using StatsBase
using Base.Test

fnecdf = ecdf(randn(10000000))
@test isapprox(fnecdf([-1.96, -1.644854, -1.281552, -0.6744898, 0, 0.6744898, 1.281552, 1.644854, 1.96]),
    [0.025, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.975], atol=1e-3)
@test isapprox(fnecdf(1.96), 0.975, atol=1e-3)
@test fnecdf([-1.96, -1.644854, -1.281552, -0.6744898, 0, 0.6744898, 1.281552, 1.644854, 1.96]) ≈
    map(fnecdf, [-1.96, -1.644854, -1.281552, -0.6744898, 0, 0.6744898, 1.281552, 1.644854, 1.96])

fnecdf = ecdf([0.5])
@test fnecdf([zeros(5000); ones(5000)]) == [zeros(5000); ones(5000)]


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.1
#   kernelspec:
#     display_name: Julia 1.1.1
#     language: julia
#     name: julia-1.1
# ---

# +
using ModelingToolkit, DifferentialEquations, Plots

# Define some variables
@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t
# -

eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

de = ODESystem(eqs)

generate_function(de, [x,y,z], [σ,ρ,β])

f = ODEFunction(de, [x,y,z], [σ,ρ,β])

u0 = [1.0,0.0,0.0]
tspan = (0.0,100.0)
p = [10.0,28.0,8/3]
prob = ODEProblem(f,u0,tspan,p)

sol = solve(prob)

plot(sol,vars=(1,2,3))


text:
#
# This file is part of the Guards.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

"""
    Guard{T}

A parametric struct returned by [`guard`](@ref).

# Parameter/Field

- `T`: the type of the guarded variable,
- `link::Union{Link,Symbol}`: `Link` or name of guard actor.
"""
struct Guard{T}
    link::Union{Link,Symbol}
end

id(x) = x

guardtype(::Guard{T}) where T = T

#
# return a Guard variable with a remote link
#
Actors._rlink(gd::Guard{T}) where T =
    Guard{guardtype(gd)}(Actors._rlink(gd.link))

"""
```
guard(var;  name=nothing, pid=myid(), thrd=false,
            sticky=false, taskref=nothing)
```

Start a `:guard` actor for the variable `var` and return
a [`Guard`](@ref) link to it.

# Parameters
- `var`: variable to guard for,
- `name=nothing`: if a `name::Symbol` is provided the server
    is registered and the name is returned,
- `remote=false`: if `remote=true` a guard with an
    remote link is returned,
- `pid=myid()`: worker pid to create the actor on,
- `thrd=false`: thread to create the actor on,
- `sticky=false`: if `true`, the actor is created in
    the same thread,
- `taskref=nothing`: if a `Ref{Task}` variable is
    provided, it gets the created `Task`.
"""
function guard(var; name=nothing, remote=false,
               pid=myid(), thrd=false,
               sticky=false, taskref=nothing)
    s = spawn(Bhv(id, var); mode=:guard, pid=pid,
            thrd=thrd, sticky=sticky, taskref=taskref)
    isnothing(name) || register(name, s)
    init!(s, id, var)
    return isnothing(name) ?
        Guard{typeof(var)}(remote ? Actors._rlink(s) : s) :
        Guard{typeof(var)}(name)
end


text:
using Base.Test
using CpuHints
using CpuId

@testset "Result Types" begin

    a = [1,2,3,4]
    p = pointer(a)

    @test isa( lfence(), Void )
    @test isa( mfence(), Void )
    @test isa( sfence(), Void )

    @test isa( prefetch(p),    Void )
    @test isa( prefetcht0(p),  Void )
    @test isa( prefetcht1(p),  Void )
    @test isa( prefetcht2(p),  Void )
    @test isa( prefetchnta(p), Void )
    @test isa( prefetchw(p),   Void )

    @test isa( clflush(p),     Void )
    @test isa( clflushopt(p),  Void )
    @test isa( clwb(p),        Void )

    @test isa( reorder_barrier(),      Void )
    @test isa( elimination_barrier(p), Void )

end


text:
using Test

@testset "mini-FL Tests" begin
   @testset "comp" begin
		@test Plasm.comp([acos,atan,tan,sin])(pi/4) * 4 - pi < 0.1e-14
		@test Plasm.comp([acos,atan,tan,sin])(1)==(acos∘atan∘tan∘sin)(1)
		@test Plasm.comp([atan,tan])(pi/4)*4 > 3.14159265358979
   end

   @testset "cons Tests" begin
		@test ((Plasm.cons([acos,atan,tan,sin])(pi/4) .<=
		[0.667458,0.665774,1.0,0.707107]))==[true,true,true,true]
		@test (Plasm.cons([sin,cos])(pi/2).>=[1.0,0.0])==[true,true]
		@test Plasm.cons([+,-])([2,3])==[[2,3],[-2,-3]]
		@test typeof(Plasm.cons([+,-])([2,3]))==Array{Array{Int64,1},1}
   end

   @testset "k Tests" begin
		@test Plasm.k(10)(100)==10
		@test Plasm.k(10)(sin)==10
		@test Plasm.k(sin)(cos)==sin
		@test Plasm.k(sin)(1000)==sin
		@test Plasm.k([1,2,3])(100)==[1,2,3]
		@test Plasm.k([1,2,3])(sin)==[1,2,3]
   end

   @testset "aa Tests" begin
		@test (Plasm.aa(sin)([0.,pi/2])==[0.0,1.0])
		@test typeof(Plasm.aa(sin)([0.,pi/6,pi/2]))==Array{Float64,1}
		@test Plasm.aa(sum)([[1,2],[3,4],[5,6],[7,8]])==[3,7,11,15]
   end
end


text:
# Social costs dynamic
using Mimi
using MimiFUND

m = MimiFUND.get_model()
update_param!(m, :climatedynamics, :climatesensitivity, 4.5)

list_scc_year = [i for i in 2000:2020]

scc_dynamic = []

for i in eachindex(list_scc_year)
    push!(scc_dynamic, MimiFUND.compute_scc(m, year = list_scc_year[i], eta = 1., prtp = 0.) * 1.68)
end

using Plots

plot(list_scc_year, scc_dynamic, label = nothing, xlabel = "Year", ylabel = "Social Cost of Carbon")

max_stars_repo_name: tlorans/SustainableFinance

max_issues_repo_name: tlorans/SustainableFinance

max_forks_repo_name: tlorans/SustainableFinance