text:
using NKFtool
using Test

try
    success(`nkf -v`)
    global got_it = true
catch
    @test 1 == 1
    global got_it = false
end

if (got_it)
    @testset "NKFtool.jl" begin

        @testset "Installation" begin
            @test !isempty(nkf_version())
            @test !isempty(nkf_help())
        end

        text_iso2022jp = raw"ＩＳＯー２０２２ーJ Ｐコード"
        code_iso2022jp = raw"GyRCI0kjUyNPITwjMiMwIzIjMiE8GyhCSiAbJEIjUCUzITwlSRsoQg=="
        #
        text_sjis = raw"シフトJ I S コード"
        code_sjis = raw"g1aDdINnSiBJIFMgg1KBW4No"
        #
        text_eucjp = raw"ＥＵＣーＪＰコード"
        code_eucjp = raw"o8Wj1aPDobyjyqPQpbOhvKXJ"
        #
        text_utf8 = raw"ＵＴＦー８コード"
        code_utf8 = raw"77y177y077ym44O877yY44Kz44O844OJ"

        using Base64
        @testset "Guess" begin
            @test nkf_guess(String(base64decode(code_iso2022jp))) ==
                  "ISO-2022-JP"
            @test nkf_guess(String(base64decode(code_sjis))) == "Shift_JIS"
            @test nkf_guess(String(base64decode(code_eucjp))) == "EUC-JP"
            @test nkf_guess(String(base64decode(code_utf8))) == "UTF-8"
        end

        @testset "Convert to UTF-8" begin
            @test nkf_convert(String(base64decode(code_iso2022jp)), "-w -m0") ==
                  text_iso2022jp
            @test nkf_convert(String(base64decode(code_sjis)), "-w -m0") ==
                  text_sjis
            @test nkf_convert(String(base64decode(code_eucjp)), "-w -m0") ==
                  text_eucjp
            @test nkf_convert(String(base64decode(code_utf8)), "-w -m0") ==
                  text_utf8
        end

        @testset "Convert from UTF-8" begin
            @test base64encode(nkf_convert(text_iso2022jp, "-j m0")) ==
                  code_iso2022jp
            @test base64encode(nkf_convert(text_sjis, "-s -m0")) == code_sjis
            @test base64encode(nkf_convert(text_eucjp, "-e -m0")) == code_eucjp
            @test base64encode(nkf_convert(text_utf8, "-w -m0")) == code_utf8
        end
    end

end


text:
# ----------------------------------------------------------------------------------------------
function imCreateButtonImages()
	new_bits = Uint8[
		2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2
		,3,2,2,1,3,2,2,1,2,2,2,2,2,2,2,2
		,1,3,2,1,4,2,1,3,2,2,2,2,2,2,2,2
		,2,1,3,1,3,1,3,0,0,0,0,0,0,0,0,2
		,2,3,1,3,4,1,1,1,1,4,4,4,4,4,0,2
		,1,1,1,4,3,3,4,4,4,4,4,4,4,4,0,2
		,2,2,4,1,4,4,1,4,4,4,4,4,4,4,0,2
		,1,2,2,1,3,4,4,4,4,4,4,4,4,4,0,2
		,2,2,2,1,4,4,4,4,4,4,4,4,4,4,0,2
		,2,2,2,0,4,4,4,4,4,4,4,4,4,4,0,2
		,2,2,2,0,4,4,4,4,4,4,4,4,4,4,0,2
		,2,2,2,0,4,4,4,4,4,4,4,4,4,4,0,2
		,2,2,2,0,4,4,4,4,4,4,4,4,4,4,0,2
		,2,2,2,0,4,4,4,4,4,4,4,4,4,4,0,2
		,2,2,2,0,4,4,4,4,4,4,4,4,4,4,0,2
		,2,2,2,0,0,0,0,0,0,0,0,0,0,0,0,2
	]

	new_colors = [
		"0 0 0",
		"132 132 132",
		"BGCOLOR",
		"255 255 0",
		"255 255 255",
		""
	]

	open_bits = Uint8[
		1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
		1,1,1,1,1,1,1,1,1,1,0,0,0,1,1,1,
		1,1,1,1,1,1,1,1,1,0,1,1,1,0,1,0,
		1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,
		1,1,0,0,0,1,1,1,1,1,1,1,1,0,0,0,
		1,0,1,3,1,0,0,0,0,0,0,0,1,1,1,1,
		1,0,3,1,3,1,3,1,3,1,3,0,1,1,1,1,
		1,0,1,3,1,3,1,3,1,3,1,0,1,1,1,1,
		1,0,3,1,3,1,0,0,0,0,0,0,0,0,0,0,
		1,0,1,3,1,0,2,2,2,2,2,2,2,2,2,0,
		1,0,3,1,0,2,2,2,2,2,2,2,2,2,0,1,
		1,0,1,0,2,2,2,2,2,2,2,2,2,0,1,1,
		1,0,0,2,2,2,2,2,2,2,2,2,0,1,1,1,
		1,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,
		1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
		1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
	]

	open_colors = [
		"0 0 0",
		"BGCOLOR",
		"128 128 0",
		"255 255 0",
		""
	]

	line_bits = Uint8[
		1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
		1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
		1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,
		1,1,1,1,1,1,1,1,1,1,1,0,0,0,1,1,
		1,1,1,1,1,1,1,1,1,1,0,0,0,1,1,1,
		1,1,1,1,1,1,1,1,1,0,0,0,1,1,1,1,
		1,1,1,1,1,1,1,1,0,0,0,1,1,1,1,1,
		1,1,1,1,1,1,1,0,0,0,1,1,1,1,1,1,
		1,1,1,1,1,1,0,0,0,1,1,1,1,1,1,1,
		1,1,1,1,1,0,0,0,1,1,1,1,1,1,1,1,
		1,1,1,1,0,0,0,1,1,1,1,1,1,1,1,1,
		1,1,1,0,0,0,1,1,1,1,1,1,1,1,1,1,
		1,1,0,0,0,1,1,1,1,1,1,1,1,1,1,1,
		1,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,
		1,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,
		1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
	]

	CreateButtonImage(16,16, new_bits, new_colors, "imImgWinNewButton");
	CreateButtonImage(16,16, open_bits, open_colors, "imImgWinOpenButton");
	CreateButtonImage(16,16, line_bits, open_colors, "imImgWinLineButton");
end

# ----------------------------------------------------------------------------------------------
function CreateButtonImage(w::Int, h::Int, bits, colors, name::String)
	global total_images, images

	color_str = ["0","1","2","3","4","5","6","7","8","9","10","11","12","13","14","15"]
	i = 1

	iup_image = IupImage(w, h, bits)
	IupSetHandle(name, iup_image)

	while (colors[i] != "")
		if (i > 16)
			aux_color_str = @sprintf("%d", i)
			IupStoreAttribute(iup_image, aux_color_str, colors[i]);
		else
			IupStoreAttribute(iup_image, color_str[i], colors[i]);
		end
		i += 1
	end

end

# ---------------------------------------------------------------------------------------------
function CreateMainMenu()

	# Creates items of menu File
	item_new   = IupItem ("New","");
	item_open  = IupItem ("Open Grid/Image", "");
	item_close = IupItem ("Close", "");
	item_exit  = IupItem ("Exit", "item_exit_act");

	# Creates items of menu Image
	item_copy  = IupItem ("Copy", "");
	item_paste = IupItem ("Paste", "")

	# Creates items for menu triangle
	item_scalenus = IupItem ("Scalenus","")
	item_isoceles = IupItem ("Isoceles", "")
	item_equilateral = IupItem ("Equilateral", "")

	# Create menu triangle
	menu_triangle = IupMenu(item_equilateral, item_isoceles, item_scalenus)

	# Creates submenu triangle
	submenu_triangle = IupSubmenu("Triangle", menu_triangle)

	# Creates items for menu create
	item_line = IupItem ("Line", "")
	item_circle = IupItem ("Circle", "")

	# Creates menu create
	menu_create = IupMenu(item_line, item_circle, submenu_triangle)

	# Creates submenu create
	submenu_create = IupSubmenu ("Create", menu_create)

	# Creates items of menu help
	item_help = IupItem ("Help", "item_help_act")

	# Creates menus
	menu_file  = IupMenu (item_new, item_open, item_close, IupSeparator(), item_exit)
	menu_image = IupMenu (item_copy, item_paste, IupSeparator(), submenu_create)
	menu_tools = IupMenu (IupItem("Tools",""))
	menu_help  = IupMenu (item_help);

	# Creates submenus
	submenu_file  = IupSubmenu ("File", menu_file)
	submenu_image = IupSubmenu ("Image", menu_image)
	submenu_tools = IupSubmenu ("Tools", menu_tools)
	submenu_help  = IupSubmenu ("Help", menu_help)

	# The menu Draw and its descendents
	menu_draw = IupMenu(IupItem("Draw",""))
	submenu_draw = IupSubmenu ("Draw", menu_draw)

	# The menu Geography and its descendents
	menu_geog = IupMenu(IupItem("Plot-coastline",""))
	submenu_geog = IupSubmenu ("Geography", menu_geog)

	# The menu Plates and its descendents
	menu_plates = IupMenu(IupItem("Plate calculator",""))
	submenu_plates = IupSubmenu ("Plates", menu_plates)

	# The menu Mag/Grav and its descendents
	menu_mag = IupMenu(IupItem("IGRF calculator",""))
	submenu_mag = IupSubmenu ("Mag/Grav", menu_mag)

	# The menu Seismology and its descendents
	menu_seismic = IupMenu(IupItem("Seismicity",""))
	submenu_seismic = IupSubmenu ("Seismology", menu_seismic)

	# The menu Tsunamis and its descendents
	menu_tsu = IupMenu(IupItem("Tsunami Travel Time",""))
	submenu_tsu = IupSubmenu("Tsunamis", menu_tsu)

	# The menu GMT and its descendents
	menu_gmt = IupMenu(IupItem("grdfilter",""))
	submenu_gmt = IupSubmenu("GMT", menu_gmt)

	# The menu Grid Tools and its descendents
	menu_grdtool = IupMenu(IupItem("grid calculator",""))
	submenu_grdtool = IupSubmenu("Grid Tools", menu_grdtool)

	# The menu Projections and its descendents
	menu_proj = IupMenu(IupItem("None",""))
	submenu_proj = IupSubmenu("Projections", menu_proj)

	# Creates main menu with file menu
	mainMenu = IupMenu(submenu_file, submenu_image, submenu_tools, submenu_draw, submenu_geog, submenu_plates,
	                   submenu_mag, submenu_seismic, submenu_tsu, submenu_gmt, submenu_grdtool, submenu_proj,
	                   submenu_help)

	return mainMenu
end

ext: jl


text:
mutable struct State
    name
    control
    next
end

function merge_microcode(code_in)
    states = Dict()
    heads = []
    substates = []
    for (opcode, cycles) in Iterators.reverse(code_in)
        next = "INST_NOP"
        for (k, control) in Iterators.reverse(enumerate(cycles))
            control = get_bits(control)
            if k == length(cycles)
                control |= C_INST_DONE
            end

            key = (control, next)
            if haskey(states, key)
                curr = states[key]
            else
                name = k == 1 ? uppercase(string(opcode)) : uppercase(string(opcode,"_",k))
                curr = State(name, control, next)
                states[key] = curr
                if k==1; push!(heads, curr)
                else; push!(substates, curr)
                end
            end
            next = curr.name
        end
    end
    return reverse(vcat(substates, heads))
end

function emit_switch_module(io::IO, name::String, input_name::String, input_bits, output_name, output_bits, mapping)
    print(io, """
        module $(name)(
                input [$(input_bits)-1:0] $(input_name),
                output reg [$(output_bits)-1:0] $(output_name)
            );
            always @* begin
                case ($(input_name))
        """);
    for (key, value) in mapping
        println(io, "            $(key): $(output_name) = $(value);")
    end
    print(io, """
                endcase
            end
        endmodule
        """);
end

function emit_enum(io::IO, enum)
    for value in instances(enum)
        println(io, "`define $(value) $(Int(value))")
    end
end

function emit_constants(io::IO, mapping)
    for (name, value) in mapping
        println(io, "`define $(name) $(value)")
    end
end

function control2text(control::Int)
    strs = [result_sel_to_name[get_res(control)]]
    flags = get_flags(control)
    for i=0:max_control_bit_index
        if flags & (1<<i) != 0
            push!(strs, " | ", control_bit_to_name[i])
        end
    end
    return prod(strs)
end

function encode_instruction(inst::Instruction)
    opcode::Int = Int(inst.op)
    scale_code::Int = inst.context.n + 13
    addr::Int = inst.context.addr

    @assert 0 <= opcode <= 15
    @assert 0 <= scale_code <= 15
    @assert 0 <= addr <= 63

    return (opcode << 12) | (scale_code << 8) | addr;
end

result_sel_pairs = [(uppercase(string(key)), value) for (key, value) in result_selections]
control_flag_pairs = [(uppercase(string(key)), value) for (key, value) in control_flags]

result_sel_to_name = Dict((value, string("`", key)) for (key, value) in result_sel_pairs)
control_bit_to_name = Dict((round(Int, log2(value)), string("`", key)) for (key, value) in control_flags)
max_control_bit_index = maximum([key for (key, value) in control_bit_to_name])

# Sort so that the opcode's zero based index in the list corresponds to its value,
# so that the emitted C defines will be consistent with the emitted states
sorted_microcode = [(opcode => opcode_programs[opcode]) for opcode in instances(Opcode)]

inst_names = [uppercase(string(opcode)) for (opcode, cycles) in sorted_microcode]

states = merge_microcode(sorted_microcode)
state_names = [state.name for state in states]
next_states = [("`"*state.name, "`"*state.next) for state in states]
push!(next_states, ("default", "`INST_NOP"))
controls = [("`"*state.name, control2text(state.control)) for state in states]
push!(controls, ("default", controls[1][2]))

state_bits = ceil(Int,log2(length(state_names)))
control_bits = max_control_bit_index+1

function emit_defines(io::IO)
    # Key constants
    emit_constants(io, [("STATE_BITS", state_bits), ("CONTROL_BITS", control_bits), ("RESULT_SEL_BITS", result_sel_bits)])
    println(io, "\n// Result selector codes")
    emit_constants(io, [(key, string(const_sel_bits, "'d", value)) for (key, value) in result_sel_pairs]);
    println(io, "\n// Control flags")
    emit_constants(io, [(key, "(1 << $(round(Int, log2(value))))") for (key, value) in control_flags]);
    println(io, "\n// States")
    for (k, name) in enumerate(state_names)
        println(io, "`define $(name) $(k-1)")
    end
end

function emit_state_machine_code(io::IO)
    emit_switch_module(io, "fms_next_state_function", "state", "`STATE_BITS", "next_state", "`STATE_BITS", next_states)
    println(io)
    emit_switch_module(io, "fms_control_function", "state", "`STATE_BITS", "control", "`CONTROL_BITS", controls)
end

function emit_instruction_defines(io::IO, instructions)
    print(io, """
        // Generated by julia/codegen.jl -- rerun to regenerate
        #ifndef __instructions_generated_h
        #define __instructions_generated_h

        """)

    max_len = maximum(length(name) for (name, inst) in instructions)

    for (name, inst) in instructions
        inst_code = encode_instruction(inst)
        pad = " "^(max_len - length(name))
        name = "INST_" * uppercase(name)
        println(io, "#define $name$pad 0x$(string(inst_code, base=16, pad=4))")
    end

    print(io, """

        #endif // __instructions_generated_h
        """)
end

emit_defines(stdout)
emit_state_machine_code(stdout)

emit_instruction_defines(stdout, instructions)

if true
    open("../fmicrosynth_generated.vh", "w") do io
        print(io, """
            // Generated by julia/codegen.jl -- rerun to regenerate

            `ifndef __fmicrosynth_generated_vh
            `define __fmicrosynth_generated_vh

            """)
        emit_defines(io)
        println(io, "\n`endif // __fmicrosynth_generated_vh")
    end
    open("../fmicrosynth_generated.v", "w") do io
        print(io, """
            // Generated by julia/codegen.jl -- rerun to regenerate
            `include "fmicrosynth_generated.vh"

            """)
        emit_state_machine_code(io)
    end
end

open("../c/instructions_generated.h", "w") do io
    emit_instruction_defines(io, instructions)
end

max_stars_repo_name: toivoh/fmicrosynth

max_issues_repo_name: toivoh/fmicrosynth

max_forks_repo_name: toivoh/fmicrosynth


text:
@testset "ComesAfter" begin

	grammar = @cgrammar begin
		Real = 1 | 2 | 3
		Real = Real + Real
		Real = Real * Real
	end

	iter = ContextExpressionIterator(grammar, 2, :Real)
	elems = collect(iter)

	constraint = ComesAfter(3,5)
	addconstraint!(grammar, constraint)

	iter2 = ContextExpressionIterator(grammar, 2, :Real)
	elems2 = collect(iter2)

	@test length(elems) > length(elems2)
	@test length(elems) == 21
	@test length(elems2) == 15
end

@testset "Ordered rules" begin
	grammar = @cgrammar begin
		Real = 1 | 2 | 3
		Real = Real + Real
		Real = Real * Real
	end

	iter = ContextExpressionIterator(grammar, 2, :Real)
	elems = collect(iter)

	constraint = Ordered([1,2,3])
	addconstraint!(grammar, constraint)

	iter2 = ContextExpressionIterator(grammar, 2, :Real)
	elems2 = collect(iter2)

	@test length(elems) > length(elems2)
	@test length(elems) == 21
	@test length(elems2) == 5
end


text:
@testset "Conversion" begin
    @test convert(Node,:a) == Node(:a)
end

struct Start{S <: Number} <: EdgeModifier
    s::S
end

@testset "Promotion" begin
    @test eltype([Node(:a), :b]) <: Node
    @test eltype([Node(:a), :b, Start(1) * :c]) <: AbstractNode
    @test eltype([:b, Start(1) * :c]) <: AbstractNode
    @test eltype([Start(1) * :b, Start(1) * :c]) <: MetaNode
end


text:
module Skans

import Base: map!

using AbstractTrees:
    PreOrderDFS,
    StatelessBFS
using HTTP:
    Form,
    get,
    post
using Gumbo:
    HTMLDocument,
    HTMLElement,
    HTMLText,
    parsehtml,
    tag
using JSON:
    json,
    parse as jsonparse
using TOML:
    parse as tomlparse

const TRIPLEQUOTE = "\"\"\""

include("select.jl")
include("webpage.jl")
include("state.jl")
include("git.jl")
include("notify.jl")
include("skan.jl")

export body, clean
export WebPage
export urls
export GitHubRepo
export skan!

end # module


text:
#*******************************************************/
#* Copyright(c) 2018 by Artelys                        */
#* This source code is subject to the terms of the     */
#* MIT Expat License (see LICENSE.md)                  */
#*******************************************************/

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#  This example demonstrates how to use Knitro to solve the following
#  simple linear programming problem (LP).
#  (This example from "Numerical Optimization", J. Nocedal and S. Wright)
#
#     minimize     -4*x0 - 2*x1
#     subject to   x0 + x1 + x2        = 5
#                  2*x0 + 0.5*x1 + x3  = 8
#                 0 <= (x0, x1, x2, x3)
#  The optimal solution is:
#     obj=-17.333 x=[3.667,1.333,0,0]
#
#  The purpose is to illustrate how to invoke Knitro using the C
#  language API.
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

using KNITRO
using Test

function example_lp1(; verbose=true)
    # Create a new Knitro solver instance.
    kc = KNITRO.KN_new()

    # Illustrate how to override default options by reading from
    # the knitro.opt file.
    options = joinpath(dirname(@__FILE__), "..", "examples", "knitro.opt")
    KNITRO.KN_load_param_file(kc, options)

    # Initialize Knitro with the problem definition.

    # Add the variables and set their bounds.
    # Note: unset bounds assumed to be infinite.
    xIndices = KNITRO.KN_add_vars(kc, 4)
    for x in xIndices
        KNITRO.KN_set_var_lobnd(kc, x, 0.0)
    end

    # Add the constraints and set the rhs and coefficients.
    cons = KNITRO.KN_add_cons(kc, 2)
    KNITRO.KN_set_con_eqbnds_all(kc, [5.0, 8.0])
    # Add Jacobian structure and coefficients.
    # First constraint
    jacIndexCons = Int32[0, 0, 0]
    jacIndexVars = Int32[0, 1, 2]
    jacCoefs = [1.0, 1.0, 1.0]
    # Second constraint
    jacIndexCons = [jacIndexCons; Int32[1, 1, 1]]
    jacIndexVars = [jacIndexVars; Int32[0, 1, 3]]
    jacCoefs = [jacCoefs; [2.0, 0.5, 1.0]]
    KNITRO.KN_add_con_linear_struct(kc, 0, Int32[0, 1, 2], [1.0, 1.0, 1.0])
    KNITRO.KN_add_con_linear_struct(kc, 1, Int32[0, 1, 3], [2.0, 0.5, 1.0])

    # Set minimize or maximize (if not set, assumed minimize).
    KNITRO.KN_set_obj_goal(kc, KNITRO.KN_OBJGOAL_MINIMIZE)

    # Set the coefficients for the objective.
    objIndices = Int32[0, 1]
    objCoefs = [-4.0, -2.0]
    KNITRO.KN_add_obj_linear_struct(kc, 2, objIndices, objCoefs)

    kn_outlev = verbose ? KNITRO.KN_OUTLEV_ALL : KNITRO.KN_OUTLEV_NONE
    KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_OUTLEV, kn_outlev)

    # Solve the problem.
    #
    # Return status codes are defined in "kn_defines.jl" and described
    # in the Knitro manual.
    nStatus = KNITRO.KN_solve(kc)
    nStatus, objSol, x, lambda_ = KNITRO.KN_get_solution(kc)

    if verbose
        println("Knitro converged with final status = ", nStatus)
        println("  optimal objective value  = ", objSol)
        println("  optimal primal values x  = ", x)
        println("  feasibility violation    = ", KNITRO.KN_get_abs_feas_error(kc))
        println("  KKT optimality violation = ", KNITRO.KN_get_abs_opt_error(kc))
    end

    # Delete the Knitro solver instance.
    KNITRO.KN_free(kc)

    @testset "Example lp1" begin
        @test nStatus == 0
        @test objSol ≈ -17.333333 atol = 1e-5
        @test x ≈ [3.667, 1.333, 0, 0] atol = 1e-3
    end
end

example_lp1(; verbose=isdefined(Main, :KN_VERBOSE) ? KN_VERBOSE : true)


text:
function main()
	N = 100
	s = string( factorial(big(N)) )
	reduce(+, map(parseint,collect(s)))
end

main()

max_stars_repo_name: cmey/euler

max_issues_repo_name: cmey/euler

max_forks_repo_name: cmey/euler


text:
using Pkg
using SafeTestsets
using Test

@time @safetestset "Interface Tests" begin include("interface_tests.jl") end
@time @safetestset "Derivative Tests" begin include("derivative_tests.jl") end
@time @safetestset "Infinite Integral Tests" begin include("inf_integral_tests.jl") end


text:
# BigWig File Format
# ==================
#
# The BigWig/BigBed format is documented in
#   Kent, W. James, et al. "BigWig and BigBed: # enabling browsing of large
#   distributed datasets." Bioinformatics 26.17 (2010): 2204-2207.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

module BigWig

import BufferedStreams
import GenomicFeatures: GenomicFeatures, BBI, Interval
import Libz
using BioCore

include("header.jl")
include("reader.jl")
include("writer.jl")
include("record.jl")
include("overlap.jl")
include("stats.jl")

end


text:
module AHPQ
    using LinearAlgebra
    using SparseArrays               # For memory-efficient identity matrices
    using Random: shuffle!           # For random initialisation
    using CUDA: CuArray              # For GPU Processing
    using Zygote: gradient           # For fast gradient computation
    using Statistics: mean, norm     # For fast stats computation

    include("utils/configs.jl")
    include("utils/distances.jl")
    include("quantizer/quantizer.jl")
    include("quantizer/preclustering.jl")
    include("builder.jl")
    include("searcher.jl")
    include("utils/testmetrics.jl")

    export builder
    export MIPS

end


text:
using ModelingToolkit, DifferentialEquations
@variables t
D = Differential(t)
sts = @variables x(t) RHS(t)  # RHS is observed
ps = @parameters τ=0.3      # parameters
@named m = ODESystem([ RHS  ~ (1 - x)/τ, D(x) ~ RHS ], t, sts, ps)
ms = structural_simplify(m)
prob = ODEProblem(ms, [x => 1.1], (0.0,1.0), [])
sol = solve(prob, Tsit5())

sol(0.3, idxs=[RHS])

max_stars_repo_name: bgctw/MTKHelpers

max_issues_repo_name: bgctw/MTKHelpers

max_forks_repo_name: bgctw/MTKHelpers


text:
#### ExportFuncs

export NetCDFWriter
export export_state
export append_data
export init_data
export full_name

using NCDatasets
using OrderedCollections

#####
##### Diagnostic variable map
#####

struct DiagnosticVariable
    name::String
    attrib::OrderedDict
    DiagnosticVariable(name::String, attrib::OrderedDict = OrderedDict()) =
        new(name, attrib)
end

function var_attrib(
    units::String = "[units]",
    long_name::String = "[long_name_default]",
    standard_name::String = "[standard_name_default]",
    fill_value::Union{Nothing, Any} = nothing,
)
    attrib = OrderedDict(
        "units" => units,
        "long_name" => long_name,
        "standard_name" => standard_name,
    )
    if !isnothing(fill_value)
        attrib["_FillValue"] = fill_value
    end
    return attrib
end

#! format: off
const GLOBAL_VARIABLE_MAP = OrderedDict{String, DiagnosticVariable}(
    "w"      => DiagnosticVariable("w"     , var_attrib("m s^-1", "vertical wind", "upward_air_velocity")),
    "ρ"      => DiagnosticVariable("ρ"     , var_attrib("kg m^-3", "air density", "air_density")),
    "a_gm"   => DiagnosticVariable("ρ"     , var_attrib("kg m^-3", "air density", "air_density")),
    "ρ_0_gm" => DiagnosticVariable("ρ_0_gm", var_attrib("units", "ρ_0_gm", "ρ_0_gm")),
    "a_gm"   => DiagnosticVariable("a_gm"  , var_attrib("units", "a_gm", "a_gm")),
    "a_en"   => DiagnosticVariable("a_en"  , var_attrib("units", "a_en", "a_en")),
    "a_ud_1" => DiagnosticVariable("a_ud_1", var_attrib("units", "a_ud_1", "a_ud_1")),
    "w_gm"   => DiagnosticVariable("w_gm"  , var_attrib("units", "w_gm", "w_gm")),
    "w_en"   => DiagnosticVariable("w_en"  , var_attrib("units", "w_en", "w_en")),
    "w_ud_1" => DiagnosticVariable("w_ud_1", var_attrib("units", "w_ud_1", "w_ud_1")),
    "ϕ_gm"   => DiagnosticVariable("ϕ_gm"  , var_attrib("units", "ϕ_gm", "ϕ_gm")),
    "ϕ_en"   => DiagnosticVariable("ϕ_en"  , var_attrib("units", "ϕ_en", "ϕ_en")),
    "ϕ_ud_1" => DiagnosticVariable("ϕ_ud_1", var_attrib("units", "ϕ_ud_1", "ϕ_ud_1")),
    "ψ_gm"   => DiagnosticVariable("ψ_gm"  , var_attrib("units", "ψ_gm", "ψ_gm")),
    "ψ_en"   => DiagnosticVariable("ψ_en"  , var_attrib("units", "ψ_en", "ψ_en")),
    "ψ_ud_1" => DiagnosticVariable("ψ_ud_1", var_attrib("units", "ψ_ud_1", "ψ_ud_1")),
    "ρ_0_gm" => DiagnosticVariable("ρ_0_gm", var_attrib("units", "ρ_0_gm", "ρ_0_gm"))
)
#! format: on

#####
##### NetCDF writer
#####

abstract type AbstractWriter end

struct NetCDFWriter{M} <: AbstractWriter
    filename::String
    meta::M
end
NetCDFWriter(filename::String) = NetCDFWriter{Nothing}(filename, nothing)

full_name(writer::NetCDFWriter) = writer.filename * ".nc"

function get_attrib(var_name)
    if haskey(GLOBAL_VARIABLE_MAP, var_name)
        var = GLOBAL_VARIABLE_MAP[var_name]
        return var.attrib
    else
        return var_attrib()
    end
end

#####
##### Implementation
#####

function export_state(
    nc::NetCDFWriter,
    grid,
    sv::StateVec,
    simtime = 0;
    kwargs...,
)
    nc = init_data(nc, grid, sv; kwargs...)
    append_data(nc, sv, simtime)
end

function init_data(nc::NetCDFWriter, grid, sv::StateVec; kwargs...)
    # TODO: allow exporting onto cell faces
    vars = OrderedDict()
    include_ghost = true
    domain_range = include_ghost ? over_elems(grid) : over_elems_real(grid)
    z = grid.zc[domain_range]
    dims = OrderedDict("z" => (z, OrderedDict()))

    # For collecting var_names and adding to GLOBAL_VARIABLE_MAP:
    # @inbounds for name_id in var_names(sv)
    #     @inbounds for i in over_sub_domains(sv, name_id)
    #         var_name = var_string(sv, name_id, i)
    #         @show var_name
    #     end
    # end

    FT = eltype(grid)

    @inbounds for name_id in var_names(sv)
        @inbounds for i in over_sub_domains(sv, name_id)
            var_name = var_string(sv, name_id, i)
            attrib = get_attrib(var_name)
            vars[var_name] = (tuple(collect(keys(dims))...), FT, attrib)
        end
    end

    init_data(nc, dims, vars; kwargs...)
    return NetCDFWriter(nc.filename, domain_range)
end

function init_data(nc::NetCDFWriter, dims, vars; overwrite = true)
    nm = full_name(nc)
    if isfile(nm)
        if overwrite
            @warn "$(nm) exists and will be overwritten."
        else
            error("$(nm) exists and overwriting is forbidden.")
        end
    end
    Dataset(nm, "c") do ds
        # define spatial and time dimensions
        for (dn, (dv, da)) in dims
            defDim(ds, dn, length(dv))
        end
        defDim(ds, "time", Inf) # Inf sets UNLIMITED dimension

        # include dimensions as variables
        for (dn, (dv, da)) in dims
            defVar(ds, dn, dv, (dn,), attrib = da)
        end
        defVar(
            ds,
            "time",
            Float64,
            ("time",),
            attrib = OrderedDict(
                "units" => "seconds since 1900-01-01 00:00:00",
                "long_name" => "time",
            ),
        )

        # define variables
        for (vn, (vd, vt, va)) in vars
            defVar(ds, vn, vt, (vd..., "time"), attrib = va)
        end
    end
    return nothing
end

function append_data(nc::NetCDFWriter, sv::StateVec, simtime = 0)
    varvals = OrderedDict()
    domain_range = nc.meta
    @inbounds for name_id in var_names(sv)
        @inbounds for i in over_sub_domains(sv, name_id)
            var_name = var_string(sv, name_id, i)
            varvals[var_name] = [sv[name_id, k, i] for k in domain_range]
        end
    end
    append_data(nc, varvals, simtime)
end

function append_data(nc::NetCDFWriter, varvals, simtime)
    Dataset(full_name(nc), "a") do ds
        timevar = ds["time"]
        t = length(timevar) + 1
        timevar[t] = simtime
        for (vn, vv) in varvals
            dsvar = ds[vn]
            dsvar[ntuple(_ -> Colon(), ndims(vv))..., t] = vv
        end
    end
    return nothing
end


text:
"""
    bresenham(start::Tuple{<:T, <:T}, stop::Tuple{<:T, <:T}) where T<:Number

Implements the [Bresenham's line algorithm](https://en.wikipedia.org/wiki/Bresenham%27s_line_algorithm).

# Examples

```jldoctest
julia> bresenham()

```
"""
function bresenham(start::Tuple{T, T}, stop::Tuple{T, T}) where T<:Integer
    x0, y0 = start
    x1, y1 = stop
    if abs(y1 - y0) < abs(x1 - x0)
        if x0 > x1
            reverse!(low_line(stop, start))
        else
            low_line(start, stop)
        end
    else
        if y0 > y1
            reverse!(high_line(stop, start))
        else
            high_line(start, stop)
        end
    end
end

"""
    low_line(start::Tuple{T, T}, stop::Tuple{T, T}) where T<:Integer

# Examples

```jldoctest
julia> low_line()

```
"""
function low_line(start::Tuple{T, T}, stop::Tuple{T, T}) where T<:Integer
    x0, y0 = start
    x1, y1 = stop
    line = fill((zero(T), zero(T)), length(x0:x1))
    dx = x1 - x0
    dy = y1 - y0
    yi = 1
    if dy < 0
        yi = -1
        dy = -dy
    end
    D = 2*dy - dx
    y = y0

    for x in x0:x1
        line[x - x0 + 1] =  (x, y)
        if D > 0
            y = y + yi
            D = D - 2*dx
        end
        D = D + 2*dy
    end

    line
end

"""
    high_line(start::Tuple{T, T}, stop::Tuple{T, T}) where T<:Integer

# Examples

```jldoctest
julia> high_line()

```
"""
function high_line(start::Tuple{T, T}, stop::Tuple{T, T}) where T<:Integer
    x0, y0 = start
    x1, y1 = stop
    line = fill((zero(T), zero(T)), length(y0:y1))
    dx = x1 - x0
    dy = y1 - y0
    xi = 1
    if dx < 0
        xi = -1
        dx = -dx
    end
    D = 2*dx - dy
    x = x0

    for y in y0:y1
        line[y - y0 + 1] =  (x, y)
        if D > 0
            x = x + xi
            D = D - 2*dy
        end
        D = D + 2*dx
    end

    line
end


text:
# ==============================================================================
# Projete xTilde sur le polyedre X du SPA avec norme-L1
# version FP 2005

function Δ2SPA(A::Array{Int,2}, xTilde::Array{Int,1})

    nbctr = size(A,1)
    nbvar = size(A,2)
    idxTilde0, idxTilde1 = split01(xTilde)

    proj = Model(GLPK.Optimizer)
    @variable(proj, 0.0 <= x[1:length(xTilde)] <= 1.0 )
    @objective(proj, Min, sum(x[i] for i in idxTilde0) + sum((1-x[i]) for i in idxTilde1) )
    @constraint(proj, [i=1:nbctr],(sum((x[j]*A[i,j]) for j in 1:nbvar)) == 1)
    optimize!(proj)
    return objective_value(proj), value.(x)
end

# ==============================================================================
# Projete xTilde sur le polyedre X du SPA avec norme-L1
# version avec somme ponderee donnant la direction vers le generateur k

function Δ2SPAbis(A::Array{Int,2}, xTilde::Array{Int,1},
                  c1::Array{Int,1}, c2::Array{Int,1}, k::Int64, λ1::Vector{Float64}, λ2::Vector{Float64})

    nbctr = size(A,1)
    nbvar = size(A,2)
    idxTilde0, idxTilde1 = split01(xTilde)

    cλ = λ1[k].*c1 + λ2[k].*c2
    proj = Model(GLPK.Optimizer)
    @variable(proj, 0.0 <= x[1:length(xTilde)] <= 1.0 )
#    @objective(proj, Min, sum(λ1[k]*x[i] for i in idxTilde0) + sum(λ2[k]*(1-x[i]) for i in idxTilde1) )
    @objective(proj, Min, sum(cλ[i]*x[i] for i in idxTilde0) + sum(cλ[i]*(1-x[i]) for i in idxTilde1) )
    @constraint(proj, [i=1:nbctr],(sum((x[j]*A[i,j]) for j in 1:nbvar)) == 1)
    optimize!(proj)
    return objective_value(proj), value.(x)
end

# ==============================================================================
# projecte la solution entiere correspondant au generateur k et test d'admissibilite
function projectingSolution!(vg::Vector{tGenerateur}, k::Int64,
                             A::Array{Int,2}, c1::Array{Int,1}, c2::Array{Int,1},
                             λ1::Vector{Float64}, λ2::Vector{Float64},
                             d::tListDisplay)

    # --------------------------------------------------------------------------
    # Projete la solution entiere sur le polytope X

#    fPrj, vg[k].sPrj.x = Δ2SPA(A,vg[k].sInt.x)
    fPrj, vg[k].sPrj.x = Δ2SPAbis(A,vg[k].sInt.x,c1,c2,k,λ1,λ2)

    # Nettoyage de la valeur de vg[k].sPrj.x et calcul du point bi-objectif
    # reconditionne les valeurs 0 et 1 et arrondi les autres valeurs
    nettoyageSolution!(vg[k].sPrj.x)
#    verbose ? @printf("  %2dP : fPrj = %8.2f  ",k, round(fPrj, digits=2)) : nothing

    # recalcule la solution au regard des 2 objectifs
    vg[k].sPrj.y[1], vg[k].sPrj.y[2] = evaluerSolution(vg[k].sPrj.x, c1, c2)
    verbose ? @printf("  %2dP : [ %8.2f , %8.2f ] ",k, vg[k].sPrj.y[1], vg[k].sPrj.y[2]) : nothing

    # archive le point obtenu pour les besoins d'affichage
    if generateurVisualise == -1
        # archivage pour tous les generateurs
        push!(d.XProj, vg[k].sPrj.y[1])
        push!(d.YProj, vg[k].sPrj.y[2])
    elseif generateurVisualise == k
        # archivage seulement pour le generateur k
        push!(d.XProj, vg[k].sPrj.y[1])
        push!(d.YProj, vg[k].sPrj.y[2])
    end

    # ----------------------------------------------------------------
    # Teste si la projection est admissible

    if estAdmissible(vg[k].sPrj.x)

        # sauvegarde de la solution entiere admissible obtenue
        vg[k].sInt.x = deepcopy(vg[k].sPrj.x)
        vg[k].sInt.y[1] = vg[k].sPrj.y[1]
        vg[k].sInt.y[2] = vg[k].sPrj.y[2]
        vg[k].sFea = true
        @printf("→ Admissible "); print("                       ")

        # archive le point obtenu pour les besoins d'affichage
        if generateurVisualise == -1
            # archivage pour tous les generateurs
            push!(d.XFeas, vg[k].sPrj.y[1])
            push!(d.YFeas, vg[k].sPrj.y[2])
        elseif generateurVisualise == k
            # archivage seulement pour le generateur k
            push!(d.XFeas, vg[k].sPrj.y[1])
            push!(d.YFeas, vg[k].sPrj.y[2])
        end

    else

        vg[k].sFea = false
        @printf("→ x          "); print("                       ")
        # prepare pour l'iteration suivante
#        vg[k].xRlx = deepcopy(vg[k].sPrj.x) !!!!!!!!!!!!!

    end

end

max_stars_repo_name: JuanferM/GravityMachine

max_issues_repo_name: JuanferM/GravityMachine

max_forks_repo_name: JuanferM/GravityMachine


text:
module Dga
export Banjori, Corebot, Cryptolocker, Dircrypt, Kraken
export Lockyv2, Pykspa, Qakbot, Ramdo, Ramnit, Simda

include(string(@__DIR__) * "/Banjori.jl")
include(string(@__DIR__) * "/Corebot.jl")
include(string(@__DIR__) * "/Cryptolocker.jl")
include(string(@__DIR__) * "/Dircrypt.jl")
include(string(@__DIR__) * "/Kraken.jl")
include(string(@__DIR__) * "/Lockyv2.jl")
include(string(@__DIR__) * "/Pykspa.jl")
include(string(@__DIR__) * "/Qakbot.jl")
include(string(@__DIR__) * "/Ramdo.jl")
include(string(@__DIR__) * "/Ramnit.jl")
include(string(@__DIR__) * "/Simda.jl")
end # module


text:
include("EchoAlgorithm.jl")

struct ElectionIDMessage{T <: Integer} <: MessageContent
    id::T
end

struct Swallow <: StopContent end

mutable struct ElectionAlgorithmWinner{T <: Integer}
    echo_algorithm_instance::Union{EchoAlgorithm,Nothing}
    winner_process_id::T
    ElectionAlgorithmWinner{T}() where T <: Integer = new(nothing, 0)

end

struct ElectionAlgorithm{T <: Integer}
    graph::AbstractGraph{T}
    process_ids::AbstractArray{T} # unique identity p
    strongest_ids::AbstractArray{T} # local variable M_p per process
    strongest_id_callback::Function
    initiators::AbstractArray{T} # process id of initiators
    initiator_id_nodes_map::Dict{T,T}
    echo_algorithms::AbstractArray{EchoAlgorithm}
    winner::ElectionAlgorithmWinner{T}

    function ElectionAlgorithm{T}(graph::AbstractGraph{T}, process_ids::AbstractArray{T}, initiators::AbstractArray{T}) where T <: Integer
        1 ≤ length(initiators) ≤ length(process_ids) || throw(error("There must not be more initiators than processes and at least one initiator!"))
        length(process_ids) == nv(graph) || throw(error("Number of processes and vertices of the graph must be equal!"))
        allunique(process_ids) || throw(error("There must be no duplicate process ids!"))
        is_connected(graph) || throw(error("Graph must be connected!"))
        initiator_id_nodes_map = Dict{T,T}()
        for (i, id) ∈ enumerate(process_ids)
            if id ∈ initiators
                initiator_id_nodes_map[id] =  i
            end
        end
        length(keys(initiator_id_nodes_map)) == length(initiators) || throw(error("All initiators must have a process id in the array!"))
        echo_algorithms = map(i -> EchoAlgorithm{T}(graph, initiator_id_nodes_map[i]), initiators)
        strongest_ids = zeros(T, length(process_ids))
        strongest_id_callback = _elec_algorithm_echo_callback_gen(strongest_ids)
        new(graph, process_ids, strongest_ids, strongest_id_callback, initiators, initiator_id_nodes_map, echo_algorithms, ElectionAlgorithmWinner{T}())
    end
end
function _elec_algorithm_echo_callback_gen(strongest_ids::AbstractArray{T}) where {T <: Integer}
    swallow = Swallow()
    return function elec_content_callback(letter::Letter{T,ElectionIDMessage{T}}) where {T <: Integer}
        recipient_node, content = dst(letter.e), letter.content
        if strongest_ids[recipient_node] <= content.id
            strongest_ids[recipient_node] = content.id
            return content
        end
        return swallow
    end
end

function ElectionAlgorithm{T}(graph::AbstractGraph{T}) where {T <: Integer}
    processes = collect(T, vertices(graph))
    ElectionAlgorithm{T}(graph, processes, processes)
end

has_winner(elec_alg::ElectionAlgorithm) = elec_alg.winner.echo_algorithm_instance !== nothing
has_terminated(elec_alg::ElectionAlgorithm) = has_winner(elec_alg) && has_terminated(elec_alg.winner.echo_algorithm_instance) # TODO implement last round
is_dead(elec_alg::ElectionAlgorithm) = !isempty(elec_alg.echo_algorithms) && all(is_dead.(elec_alg.echo_algorithms))

function init_election_algorithm!(elec_alg::ElectionAlgorithm)
    for initiator ∈ elec_alg.initiators # set M_p to p
        elec_alg.strongest_ids[elec_alg.initiator_id_nodes_map[initiator]] = initiator
    end
    init_echo_algorithm!.(elec_alg.echo_algorithms, [ElectionIDMessage(id) for id in elec_alg.initiators])
    return elec_alg
end

function election_algorithm_step!(elec_alg::ElectionAlgorithm)
    if has_winner(elec_alg)
        return _winner_step!(elec_alg)
    else
        dead_echo_instances = Vector{Int64}()
        for (i, echo_alg) in enumerate(elec_alg.echo_algorithms)
            if has_terminated(echo_alg)
                return _establish_winner!(elec_alg, echo_alg)
            elseif is_dead(echo_alg)
                push!(dead_echo_instances, i)
            else
                echo_algorithm_step!(echo_alg, elec_alg.strongest_id_callback)
            end
        end
        deleteat!(elec_alg.echo_algorithms, dead_echo_instances)
        return elec_alg
    end
end

function _establish_winner!(elec_alg::ElectionAlgorithm, winner_echo_instance::EchoAlgorithm)
    elec_alg.winner.echo_algorithm_instance = winner_echo_instance
    elec_alg.winner.winner_process_id = elec_alg.process_ids[winner_echo_instance.initiator]
    empty!(elec_alg.echo_algorithms)
    return elec_alg
end

function _winner_step!(elec_alg::ElectionAlgorithm) # TODO implement last distribution of winner
    return elec_alg
end

max_stars_repo_name: tomtuamnuq/dist-algorithms

max_issues_repo_name: tomtuamnuq/dist-algorithms

max_forks_repo_name: tomtuamnuq/dist-algorithms


text:
import RationalApproximations: aaa_v1, aaa_v2, locate_extrema_idx
import RationalApproximations: RationalFunc, equioscillation!, minimax
using PyPlot

f = exp
S = range(-1, 1, length=201)
F = f.(S)
tol = 1e-6
max_m = 5
r, err = aaa_v1(F, S, max_m, tol)
m = length(r)

figure(1)
supp_pt = r[m].supp_pt
resid = F - r[m].(S)
plot(S, resid, supp_pt, f.(supp_pt)-r[m].(supp_pt), "o")
legend((L"$f-r$", "support points"))
title("Result using aaa_v1 (m = $m)")
grid(true)
xylims = axis()

figure(2)
r, err = aaa_v2(F, S, max_m, tol)
m = length(r)
supp_pt = r[m].supp_pt
resid = F - r[m].(S)
plot(S, resid, supp_pt, f.(supp_pt)-r[m].(supp_pt), "o")
legend((L"$f-r$", "support points"))
title("Result using aaa_v2 (m = $m)")
grid(true)
axis(xylims)

idx = locate_extrema_idx(resid, m)
x = S[idx]
figure(3)
plot(S, resid, x, f.(x)-r[m].(x), "o")
legend((L"$f-r$", "extrema"))
title("Locate extrema")
grid(true)
axis(xylims)

clamp = :no
α, β, λ = equioscillation!(f, x, supp_pt, clamp)
new_r = RationalFunc(α, β, supp_pt)
figure(4)
plot(S, f.(S) - new_r.(S),
     [-1.0, 1.0, NaN, -1.0, 1.0], [λ, λ, NaN, -λ, -λ], ":k",
     x, f.(x)-new_r.(x), "o")
grid(true)
title("Equioscillation procedure")
axis(xylims)

max_iterations = 4
x = S[idx]
r, zmin, zmax = minimax(f, x, supp_pt, (-1.0,1.0), max_iterations, clamp)
m = length(r)
figure(5)
plot(S, f.(S) - r[m].(S), x, f.(x) - r[m].(x), "o")
grid(true)
title("Result after $m Remez iterations")
axis(xylims)

max_stars_repo_name: billmclean/RationalApproximations

max_issues_repo_name: billmclean/RationalApproximations

max_forks_repo_name: billmclean/RationalApproximations


text:
using Oceananigans.Distributed: reconstruct_global_grid

function random_divergent_source_term(arch, grid)
    # Generate right hand side from a random (divergent) velocity field.
    Ru = XFaceField(arch, grid)
    Rv = YFaceField(arch, grid)
    Rw = ZFaceField(arch, grid)
    U = (u=Ru, v=Rv, w=Rw)

    Nx, Ny, Nz = size(grid)
    set!(Ru, (x, y, z) -> rand())
    set!(Rv, (x, y, z) -> rand())
    set!(Rw, (x, y, z) -> rand())

    fill_halo_regions!(Ru, arch)
    fill_halo_regions!(Rv, arch)
    fill_halo_regions!(Rw, arch)

    # Compute the right hand side R = ∇⋅U
    ArrayType = array_type(arch)
    R = zeros(Nx, Ny, Nz) |> ArrayType
    event = launch!(arch, grid, :xyz, divergence!, grid, U.u.data, U.v.data, U.w.data, R,
                    dependencies=Event(device(arch)))
    wait(device(arch), event)

    return R
end

function divergence_free_poisson_solution_triply_periodic(grid_points, ranks)
    topo = (Periodic, Periodic, Periodic)
    arch = MultiArch(CPU(), ranks=ranks, topology = topo)
    local_grid = RectilinearGrid(arch, topology=topo, size=grid_points, extent=(1, 2, 3))

    full_grid = reconstruct_global_grid(local_grid)
    solver = DistributedFFTBasedPoissonSolver(arch, full_grid, local_grid)

    R = random_divergent_source_term(child_architecture(arch), local_grid)
    first(solver.storage) .= R

    ϕc = first(solver.storage)
    solve!(ϕc, solver)

    p_bcs = FieldBoundaryConditions(local_grid, (Center, Center, Center))
    p_bcs = inject_halo_communication_boundary_conditions(p_bcs, arch.local_rank, arch.connectivity)

    ϕ   = CenterField(child_architecture(arch), local_grid, p_bcs)  # "pressure"
    ∇²ϕ = CenterField(child_architecture(arch), local_grid, p_bcs)

    interior(ϕ) .= real(first(solver.storage))
    compute_∇²!(∇²ϕ, ϕ, arch, local_grid)

    return R ≈ interior(∇²ϕ)
end

@testset "Distributed FFT-based Poisson solver" begin
    @info "  Testing distributed FFT-based Poisson solver..."
    @test divergence_free_poisson_solution_triply_periodic((16, 16, 1), (1, 4, 1))
    @test divergence_free_poisson_solution_triply_periodic((44, 44, 1), (1, 4, 1))
    @test divergence_free_poisson_solution_triply_periodic((44, 16, 1), (1, 4, 1))
    @test divergence_free_poisson_solution_triply_periodic((16, 44, 1), (1, 4, 1))
end


text:
module SignedDistanceFunctions

using Meshes

"""
Signed Distance Function. Returns the signed distance between an implicit surface defined as the root of this function, with a positive distance outside and a negative distance inside.
"""
abstract type SDF end

(p::SDF)(x) = p.f(x)

"""
Exact Signed Distance Function, valid outside and inside its implicit surface.
"""
struct TrueSDF <: SDF
    f::Function
    TrueSDF(f::Function) = new(f)
end

"""
Signed Distance Function validly defined outside its implicit surface. Inside the surface, only the sign of the returned value is correct.
"""
struct ExteriorSDF <: SDF
    f::Function
    ExteriorSDF(f::Function) = new(f)
end

struct BooleanSDF <: SDF
    f::Function
    BooleanSDF(f::Function) = new(f)
end

Base.union(sdfs::SDF...) = ExteriorSDF(x -> minimum(sdf.f(x) for sdf in sdfs))
Base.intersect(sdfs::SDF...) = BooleanSDF(x -> maximum(sdf.f(x) for sdf in sdfs))
Base.:-(sdf1::SDF, sdf2::SDF) = BooleanSDF(x -> max(sdf1.f(x), -sdf2.f(x)))

onion(thickness, sdf::SDF) = typeof(sdf)(x -> abs(sdf.f(x)) - thickness)
onion(thickness, obj) = onion(thickness, sdf(obj))
Base.round(radius, sdf::SDF) = typeof(sdf)(x -> sdf.f(x) - radius)

function smooth_union(x, y, factor)
    h = clamp(0.5 + 0.5 * (x - y)/factor, 0., 1.)
    coordinates(Segment(Point(x), Point(y))(h))[1] - factor * h * (1. - h)
end

smooth_union(sdf1::SDF, sdf2::SDF, factor) = BooleanSDF(x -> smooth_union(sdf1.f(x), sdf2.f(x), factor))

include("objects.jl")

export
    TrueSDF,
    ExteriorSDF,
    BooleanSDF,
    sdf,
    onion,
    smooth_union,

    Ellipse,
    HalfSpace

end


text:
using PyCall
using Makie

# Construct data
spt = pyimport_conda("scipy.spatial", "scipy")

# Generate random points
pnts = [Point(rand(2)...) for i in 1 : 100]

# Find the convex hull
hull = spt.ConvexHull(collect(hcat(collect.(pnts)...)') )
vtx =  [Point(val[1], val[2]) for val in eachrow(hull.points[hull.vertices .+ 1 , :])]

fig, ax, plt = scatter(pnts, color=:black, markersize=5)
lines!(vtx, linewidth=3)
lines!(vtx[[1, end]], linewidth=3)
fig


text:
using Test
using Dates
using ExtensibleScheduler
using TimeFrames

@testset "TimeFrameTrigger" begin
    tf = tf"5Min"
    trigger = Trigger(tf)

    n_executed = 0

    function sample_action()
        n_executed += 1
    end

    dt_start = DateTime(2010, 1, 3, 23, 58)
    dt = dt_start
    clock = SimClock(dt)
    sched = BlockingScheduler(clock=clock, delayfunc=NoSleep)
    action = Action(sample_action)

    add(sched, action, trigger)
    run_pending(sched)
    @test n_executed == 0

    for i in 1:13
        dt += Dates.Minute(1)
        #println("Clock set to $dt")
        set(clock, dt)
        run_pending(sched)
    end

    @test n_executed == 3

end


text:
# For this to work, you need to be able to run
#     using Images
#     include(joinpath(dirname(dirname(pathof(Images))), "test", "runtests.jl"))
# successfully. Even if you have the Images package, you may need to add packages.

using SnoopCompile

### Log the compiles
# This only needs to be run once (to generate "/tmp/images_compiles.csv")

SnoopCompile.@snoopc "/tmp/images_compiles.csv" begin
    using Images, Pkg
    include(joinpath(dirname(dirname(pathof(Images))), "test", "runtests.jl"))
end

### Parse the compiles and generate precompilation scripts
# This can be run repeatedly

data = SnoopCompile.read("/tmp/images_compiles.csv")

# Old versions of Images ran the tests inside a module ImagesTest, so all
# the precompiles would get credited to ImagesTest. Credit them to Images instead:
subst = ["ImagesTests"=>"Images", ]

# Blacklist can be used to help fix problems.
# For example, old versions of Julia output would MIME types with invalid syntax:
blacklist = ["MIME", ]

# Use these two lines if you want to create precompile functions for
# individual packages
pc = SnoopCompile.parcel(reverse!(data[2]), subst=subst, blacklist=blacklist)
SnoopCompile.write("/tmp/precompile", pc)

# Use these two lines if you want to add to your userimg.jl
pc = SnoopCompile.format_userimg(reverse!(data[2]), subst=subst, blacklist=blacklist)
SnoopCompile.write("/tmp/userimg_Images.jl", pc)


text:
#################################################################
#
#    Internal function testing
#
#    m = ReverseDiffSource module
#################################################################

@test m.isSymbol(:a)            == true
@test m.isSymbol(:(a[1]))       == false
@test m.isSymbol(:(a.b))        == false
@test m.isSymbol(:(exp(a)))     == false

@test m.isRef(:a)            == false
@test m.isRef(:(a[1]))       == true
@test m.isRef(:(a[x]))       == true
@test m.isRef(:(a.b))        == false
@test m.isRef(:(a.b[end]))   == false
@test m.isRef(:(a[end].b))   == false
@test m.isRef(:(exp(a)))     == false

@test m.isDot(:a)           == false
@test m.isDot(:(a[1]))      == false
@test m.isDot(:(a[x]))      == false
@test m.isDot(:(a.b))       == true
@test m.isDot(:(a.b[end]))  == false
@test m.isDot(:(a[end].b))  == false
@test m.isDot(:(exp(a)))    == false

@test m.dprefix("coucou")            == :dcoucou
@test m.dprefix(:tr)                 == :dtr

#### accumulator variable generation testing ####
    function zerocode(v)
        zex = m.zeronode( m.NConst(:abcd, [], [], v, false) )
        m.resetvar()
        m.tocode(zex)
    end

    @test zerocode(12)       == :(0.0;)
    @test zerocode(12.2)     == :(0.0;)
    @test zerocode(5:6)      == :(zeros(2); )
    @test zerocode(5:0.2:10) == :(zeros(2); )

    @test zerocode( [2., 2 , 3, 0] )    == :(zeros(size(tv)); )
    @test zerocode( [2.; [1,2 ]; 3.] )  == :(zeros(size(tv)); )

    @test zerocode( [false, true] )    == :(zeros(size(tv)); )

    @test zerocode(12.2-4im) == striplinenumbers( quote
                                    _tmp1 = cell(2)
                                    _tmp1[1] = 0.0
                                    _tmp1[2] = 0.0
                                    _tmp1
                                end )

    type Abcd
        a::Float64
        b
        c::Vector{Float64}
    end
    @test zerocode(Abcd(2., 3., [1,2 ]))  == striplinenumbers( quote
                                                    _tmp1 = cell(3)
                                                    _tmp1[1] = 0.0
                                                    _tmp1[2] = 0.0
                                                    _tmp1[3] = zeros(size(tv.c))
                                                    _tmp1
                                                end )

    @test zerocode([Abcd(2., 3., [1,2 ]), Abcd(2., 3., [1,2 ])])  ==
            striplinenumbers( quote
                                  _tmp1 = cell(size(tv))
                                  for i = 1:length(_tmp1)
                                      _tmp2 = cell(3)
                                      _tmp2[1] = 0.0
                                      _tmp2[2] = 0.0
                                      _tmp2[3] = zeros(size(tv[i].c))
                                      _tmp1[i] = _tmp2
                                  end
                                  _tmp1
                              end )

    @test zerocode( (2., 3., [1,2 ]) )    == striplinenumbers( quote
                                                    _tmp1 = cell(3)
                                                    _tmp1[1] = 0.0
                                                    _tmp1[2] = 0.0
                                                    _tmp1[3] = zeros(size(tv[3]))
                                                    _tmp1
                                                end )

    @test zerocode( Any[2., [1,2 ], 3.] )    == striplinenumbers( quote
                                                    _tmp1 = cell(3)
                                                    _tmp1[1] = 0.0
                                                    _tmp1[2] = zeros(size(tv[2]))
                                                    _tmp1[3] = 0.0
                                                    _tmp1
                                                end )

    @test zerocode( [2.+im, 3.-2im] )    == striplinenumbers( quote
                                                        _tmp1 = cell(size(tv))
                                                        for i = 1:length(_tmp1)
                                                            _tmp2 = cell(2)
                                                            _tmp2[1] = 0.0
                                                            _tmp2[2] = 0.0
                                                            _tmp1[i] = _tmp2
                                                        end
                                                        _tmp1
                                                    end )

### tmatch (naive multiple dispatch) testing  ###
    tts = Any[ Tuple{Real},
               Tuple{Real, Real},
               Tuple{Float64, Float64},
               Tuple{Float64, Int},
               Tuple{Float64, Int64},
               Tuple{Float64},
               Tuple{Int64},
               Tuple{AbstractString},
               Tuple{Any},
               Tuple{Any, Any} ]

    @test m.tmatch(Tuple{Float64}, tts)           == Tuple{Float64}
    @test m.tmatch(Tuple{Int64}, tts)             == Tuple{Int64}
    @test m.tmatch(Tuple{Float64,Int64}, tts)     == Tuple{Float64, Int64}
    @test m.tmatch(Tuple{Float64,Int32}, tts)     == Tuple{Real,Real}
    @test m.tmatch(Tuple{Any, Real}, tts)         == Tuple{Any,Any}
    @test m.tmatch(Tuple{Float64, Int, Int}, tts) == nothing
    @test m.tmatch(Tuple{AbstractString}, tts)    == Tuple{AbstractString}
    @test m.tmatch(Tuple{Float32}, tts)           == Tuple{Real}
    @test m.tmatch(Tuple{Float64,Real}, tts)      == Tuple{Real,Real}
    @test m.tmatch(Tuple{Float64,Float64}, tts)   == Tuple{Float64,Float64}
    @test m.tmatch(Tuple{Float64,Vector}, tts)    == Tuple{Any,Any}
    @test m.tmatch(Tuple{Real,Float64}, tts)      == Tuple{Real,Real}

### testing conversions for functions diff  ###
  ex = quote
    a = zeros(2)
    for i in 1:2
        a[i] = x
    end
    sum(a)
  end
  @test m.e2s(m.streamline(ex)) == "↑=→:a→↑call→:zeros→2↓↓↑for→↑=→:i→↑:→1→2↓↓" *
                               "→↑block→↑=→↑ref→:a→:i↓→:x↓↓↓↑call→:sum→:a↓"

  ex = quote
      a = zeros(2)
      c = 0
      for i in 1:2
          d = 3
          for j in 1:2
              a[i] = x
              a[i] += b[i+1]
              d = d * sin(2i)
          end
          c += d
      end
      z = c * sum(a)
      return z - sum(a)
  end
  @test m.e2s(m.streamline(ex)) == "↑=→:a→↑call→:zeros→2↓↓↑=→:c→0↓↑for→↑=→:i→↑:→1→2↓↓→↑block→" *
                               "↑=→:d→3↓→↑for→↑=→:j→↑:→1→2↓↓→↑block→↑=→↑ref→:a→:i↓→:x↓→" *
                               "↑+=→↑ref→:a→:i↓→↑ref→:b→↑call→:+→:i→1↓↓↓→↑=→:d→↑call→:*→:d→" *
                               "↑call→:sin→↑call→:*→2→:i↓↓↓↓↓↓→↑+=→:c→:d↓↓↓↑=→:z→↑call→:*→:c→" *
                               "↑call→:sum→:a↓↓↓↑return→↑call→:-→:z→↑call→:sum→:a↓↓↓"

  aex = Expr[ :(a = zeros(2)),
              striplinenumbers(:(for i = 1:2 ; a[i] = x ; end)),
              :(sum(a)) ]

  @test m.s2e("↑=→:a→↑call→:zeros→2↓↓↑for→↑=→:i→↑:→1→2↓↓" *
              "→↑block→↑=→↑ref→:a→:i↓→:x↓↓↓↑call→:sum→:a↓" ) == aex


text:
# Constants - from mathlink.h

@enum Packet begin
    PKT_ILLEGAL    =   0

    PKT_CALL       =   7
    PKT_EVALUATE   =  13
    PKT_RETURN     =   3

    PKT_INPUTNAME  =   8
    PKT_ENTERTEXT  =  14
    PKT_ENTEREXPR  =  15
    PKT_OUTPUTNAME =   9
    PKT_RETURNTEXT =   4
    PKT_RETURNEXPR =  16

    PKT_DISPLAY    =  11
    PKT_DISPLAYEND =  12

    PKT_MESSAGE    =   5
    PKT_TEXT       =   2

    PKT_INPUT      =   1
    PKT_INPUTSTR   =  21
    PKT_MENU       =   6
    PKT_SYNTAX     =  10

    PKT_SUSPEND    =  17
    PKT_RESUME     =  18

    PKT_BEGINDLG   =  19
    PKT_ENDDLG     =  20

    PKT_FIRSTUSER  = 128
    PKT_LASTUSER   = 255
end

@enum Token begin
    TK_ERROR    = 0           # /* bad token */

    TK_OLDINT   = Cint('I')   # /* 73 Ox49 01001001 */ # /* integer leaf node */
    TK_OLDREAL  = Cint('R')   # /* 82 Ox52 01010010 */ # /* real leaf node */
    TK_FUNC     = Cint('F')   # /* 70 Ox46 01000110 */ # /* non-leaf node */

    TK_STR      = Cint('"')   # /* 34 0x22 00100010 */
    TK_SYM      = Cint('#')   # /* 35 0x23 # 00100011 */

    TK_REAL     = Cint('*')   # /* 42 0x2A 00101010 */
    TK_INT      = Cint('+')   # /* 43 0x2B 00101011 */

    # /* The following defines are for internal use only */
    TK_PCTEND   = Cint(']')    # /* at end of top level expression */
    TK_END      = Cint('\n')
    TK_AEND     = Cint('\r')
    TK_SEND     = Cint(',')
    TK_CONT     = Cint('\\')
    TK_ELEN     = Cint(')')
    TK_NULL     = Cint('.')
    TK_OLDSYM   = Cint('Y')    # /* 89 0x59 01011001 */
    TK_OLDSTR   = Cint('S')    # /* 83 0x53 01010011 */
    TK_PACKED   = Cint('P')    # /* 80 0x50 01010000 */
    TK_ARRAY    = Cint('A')    # /* 65 0x41 01000001 */
    TK_DIM      = Cint('D')    # /* 68 0x44 01000100 */

    TK_INVALID  = 155
    TK_INT8     = 160
    TK_UINT8    = 161

    TK_INT16_BE    = 162
    TK_UINT16_BE   = 163
    TK_INT32_BE    = 164
    TK_UINT32_BE   = 165
    TK_INT64_BE    = 166
    TK_UINT64_BE   = 167

    TK_INT16_LE    = 226
    TK_UINT16_LE   = 227
    TK_INT32_LE    = 228
    TK_UINT32_LE   = 229
    TK_INT64_LE    = 230
    TK_UINT64_LE   = 231

    TK_FLOAT32_BE  = 180
    TK_FLOAT64_BE  = 182
    TK_FLOAT128_BE = 184

    TK_FLOAT32_LE  = 244
    TK_FLOAT64_LE  = 246
    TK_FLOAT128_LE = 248
end

@enum Error begin
    ERR_UNKNOWN         =   -1
    ERR_OK              =    0
    ERR_DEAD            =    1
    ERR_GBAD            =    2
    ERR_GSEQ            =    3
    ERR_PBTK            =    4
    ERR_PSEQ            =    5
    ERR_PBIG            =    6
    ERR_OVFL            =    7
    ERR_MEM             =    8
    ERR_ACCEPT          =    9
    ERR_CONNECT         =   10
    ERR_CLOSED          =   11
    ERR_DEPTH           =   12  # /* internal error */
    ERR_NODUPFCN        =   13  # /* stream cannot be duplicated */

    ERR_NOACK           =   15  # /* */
    ERR_NODATA          =   16  # /* */
    ERR_NOTDELIVERED    =   17  # /* */
    ERR_NOMSG           =   18  # /* */
    ERR_FAILED          =   19  # /* */

    ERR_GETENDEXPR      =   20
    ERR_PUTENDPACKET    =   21 # /* unexpected call of MLEndPacket */
    # /* currently atoms aren't
    # * counted on the way out so this error is raised only when
    # * MLEndPacket is called in the midst of an atom
    # */
    ERR_NEXTPACKET      =   22
    ERR_UNKNOWNPACKET   =   23
    ERR_GETENDPACKET    =   24
    ERR_ABORT           =   25
    ERR_MORE            =   26 # /* internal error */
    ERR_NEWLIB          =   27
    ERR_OLDLIB          =   28
    ERR_BADPARAM        =   29
    ERR_NOTIMPLEMENTED  =   30

    ERR_INIT            =   32
    ERR_ARGV            =   33
    ERR_PROTOCOL        =   34
    ERR_MODE            =   35
    ERR_LAUNCH          =   36
    ERR_LAUNCHAGAIN     =   37
    ERR_LAUNCHSPACE     =   38
    ERR_NOPARENT        =   39
    ERR_NAMETAKEN       =   40
    ERR_NOLISTEN        =   41
    ERR_BADNAME         =   42
    ERR_BADHOST         =   43
    ERR_RESOURCE        =   44  # /* a required resource was missing */
    ERR_LAUNCHFAILED    =   45
    ERR_LAUNCHNAME      =   46

    ERR_PDATABAD        =   47
    ERR_PSCONVERT       =   48
    ERR_GSCONVERT       =   49
    ERR_NOTEXE          =   50
    ERR_SYNCOBJECTMAKE  =   51
    ERR_BACKOUT         =   52

    ERR_TRACEON         =  996  # /* */
    ERR_TRACEOFF        =  997  # /* */
    ERR_DEBUG           =  998  # /* */
    ERR_ASSERT          =  999  # /* an internal assertion failed */
    ERR_USER            = 1000  # /* start of user defined errors */
end


text:
"""
    MetaDiGraph

A `MetaGraph` whose underlying graph is of type `Graphs.SimpleDiGraph`.
"""
const MetaDiGraph = MetaGraph{<:Any,<:Any,<:SimpleDiGraph}

Graphs.SimpleDiGraph(g::MetaDiGraph) = g.graph

Graphs.is_directed(::Type{<:MetaDiGraph}) = true
Graphs.is_directed(::MetaDiGraph) = true

function arrange(
    ::MetaDiGraph{<:Any,Label}, label_1::Label, label_2::Label, args...
) where {Label}
    return label_1, label_2
end


text:
module AnimationExport

export exportAnimation

using StaticArrays

import Modia3D
import Modia3D.Shapes
import Modia3D.Composition
import LinearAlgebra
import JSON
import UUIDs

include("exportAnimation.jl")

end


text:
function search_and_score!(r::Robot)
    score = 0
    count = 0
    side = Ost
    while true
        while !isborder(r, side)
            if ismarker(r)
                score+=temperature(r)
                count+=1
            end
            move!(r, side)
        end
        if ismarker(r)
            score+=temperature(r)
            count+=1
        end
        if isborder(r, Sud)
            break
        end
        move!(r, Sud)
        side = reversSide(side)
    end
    println("Average score: ", score/count, "\nQuantity: ", count)
end


text:
"""
    x, randx, logev = noisyintersectionvi(; posterior = posterior, g = g, [σ2 = 1e-8*ones(length(g))])

`posterior` is the posterior distribution returned by bpca, of type ```Distributions.MvNormal```.

`g` ∈ ℜⁿ₊ is the candidate direction vector that is being tested, of type ```Array{Float64, 1}```.

The function returns:

(1) the average intersection point `x`.

(2) a function ```randx()``` that draws a sample intersection point.

(3) the ```logevidence```, a measure of how good an intersection the candidate direction vector is.
"""
function noisyintersectionvi(; posterior = posterior, g = g, σ2 = 1e-6*ones(length(g)))

    # check dimensions of arguments

    @assert(length(g) == round(Int, length(posterior)/2))

    # check that components of vector are positive

    @assert(all(g .> 0.0))

    # work with normalised vector

    gunit = g / norm(g)

    linemeansigma = marginaline(posterior, σ2)

    function logp(p)

        local μ, Σ = linemeansigma(p[2])

        logpdf(MvNormal(μ, Σ), p[1]*gunit)

    end

    # Initialise with ML estimate

    bestsolution = noisyintersectionML(; posterior = posterior, g = gunit, σ2 = σ2)

    # Call variational inference

    postp, logevidence = VI(x->logp(x), [bestsolution.+(i/10)*(randn(length(bestsolution))) for i=1:30], optimiser = NelderMead(), S = 750, inititerations = 50, iterations = 1000000, show_every = 25)

    pdfscale = Normal(mean(postp)[1], sqrt(cov(postp)[1,1]))

    randx = ()->rand(pdfscale)*gunit

    return mean(pdfscale)*gunit, randx, logevidence

end


text:
### road_constraints.jl

flow_var(G, i, j) = get_prop(G, i, j, :flow_var)

"""
Standard flow conservation constraints, not applied at dead-ends
"""
function cs_conserve_flow!(m, M::RoadModel)
    G = M.G
    for v in vertices(G)
		if length(Set(all_neighbors(G, v))) <= 1
			continue
		end

		expr = 0
		for i in inneighbors(G, v)
			expr += flow_var(G, i, v)
		end

		for o in outneighbors(G, v)
			expr -= flow_var(G, v, o)
		end

		push!(M.constraints, @constraint(m, expr == 0))
	end
end

function cs_inflow!(m, M::RoadModel)
    G = M.G
	for v in vertices(G)
		for iv in inneighbors(G, v)
			expr = -flow_var(G, iv, v)

			ov_count = 0
			for ov in outneighbors(G, v)
				if ov == iv
					continue
				end
				expr += flow_var(G, v, ov)
				ov_count += 1
			end

			if ov_count == 0
				continue # dead-end
			end

			push!(M.constraints, @constraint(m, expr >= 0))
		end
	end
end

"""
Introduces additional variables for which direction flow takes at each junction. Can be used to model allowed turning directions
"""
function cs_junction!(m, M::RoadModel)
	G = M.G
	for v in vertices(G)
		if length(outneighbors(G, v)) == 0 || length(inneighbors(G, v)) == 0
			continue
		end

		out_exprs = Dict(ov => -flow_var(G, v, ov) for ov in outneighbors(G, v))
		for iv in inneighbors(G, v)
			in_expr = -flow_var(G, iv, v)
			for ov in outneighbors(G, v)
				if ov == iv
					continue
				end

				# declare anonymous junction variable
				iv_ov_var = @variable(m)
				set_lower_bound(iv_ov_var, 0)
				# set_upper_bound(iv_ov_var, 500) # maybe upper bound speeds up solver? (nope)

				out_exprs[ov] += iv_ov_var
				in_expr += iv_ov_var
			end
			push!(M.constraints, @constraint(m, in_expr == 0))
		end

		for oc in values(out_exprs)
			push!(M.constraints, @constraint(m, oc == 0))
		end
	end
end

"""
Impose "max flow creation" constraints on larger basic cycles. Not feasible to do this for *all* cycles (unless using lazy constraints?)
"""
function cs_anticycle_local!(m, M::RoadModel)
    # TODO: Waiting on cycle area implementation in `road_utils`
end

"""
Form connected components of "bad/small" cycles and bound all contained edges by the net flow entering the component
"""
function cs_anticycle_nonlocal!(m, M::RoadModel)
	# Not implemented
end

function cs_sensor!(m, M::RoadModel)
    for (u, v, prop) in edge_prop_triples(M.G)
        push!(M.constraints, @constraint(m, prop[:flow_var] <= prop[:capacity] * (1 - M.MP.sensor_cap_ratio * prop[:sensor_var])))
    end

    if M.MP.force_sensor_bidirectional
		for (u, v, prop) in edge_prop_triples(M.G)
			# sort edge components so that constraint is only added once
			if has_edge(M.G, v, u) && u < v
				push!(M.constraints, @constraint(m, prop[:sensor_var] == get_prop(M.G, v, u, :sensor_var)))
			end
		end
    end
end

export flow_var, cs_conserve_flow!, cs_inflow!, cs_junction!, cs_anticycle_local!, cs_anticycle_nonlocal!, cs_sensor!


text:
import Pkg

root_directory = dirname(@__DIR__)

nnlib = Pkg.PackageSpec(path = root_directory)
Pkg.develop(nnlib)
Pkg.precompile()


text:
# GPUArrays.jl interface

#
# Device functionality
#

## device properties

GPUArrays.threads(dev::CUDAdrv.CuDevice) =
    CUDAdrv.attribute(dev, CUDAdrv.DEVICE_ATTRIBUTE_MAX_THREADS_PER_BLOCK)

## execution

struct CuArrayBackend <: AbstractGPUBackend end

struct CuKernelContext <: AbstractKernelContext end

function GPUArrays.gpu_call(::CuArrayBackend, f, args...; threads::Int, blocks::Int)
    @cuda threads=threads blocks=blocks f(CuKernelContext(), args...)
end

GPUArrays.synchronize(A::CuArray) = CUDAdrv.synchronize()

## on-device

# indexing

GPUArrays.blockidx(ctx::CuKernelContext) = CUDAnative.blockIdx().x
GPUArrays.blockdim(ctx::CuKernelContext) = CUDAnative.blockDim().x
GPUArrays.threadidx(ctx::CuKernelContext) = CUDAnative.threadIdx().x
GPUArrays.griddim(ctx::CuKernelContext) = CUDAnative.gridDim().x

# memory

@inline function GPUArrays.LocalMemory(::CuKernelContext, ::Type{T}, ::Val{dims}, ::Val{id}
                                      ) where {T, dims, id}
    ptr = CUDAnative._shmem(Val(id), T, Val(prod(dims)))
    CuDeviceArray(dims, DevicePtr{T, CUDAnative.AS.Shared}(ptr))
end

# synchronization

@inline GPUArrays.synchronize_threads(::CuKernelContext) = CUDAnative.sync_threads()

#
# Host abstractions
#

GPUArrays.device(A::CuArray) = CUDAdrv.device(CUDAdrv.CuCurrentContext())

GPUArrays.backend(::Type{<:CuArray}) = CuArrayBackend()

GPUArrays.unsafe_reinterpret(::Type{T}, A::CuArray, size::NTuple{N, Integer}) where {T, N} =
  CuArray{T,N}(convert(CuPtr{T}, A.ptr), size, A)


text:
basedir = realpath(joinpath(@__DIR__, ".."))
basedir = realpath("/home/raf/.julia/dev/SpottedWingPaper")
@show basedir
floattype = Float64
# Set the input data used in human dispersal models
humanactivity = :gdp
# humanactivity = :pop

# Opimisation and validation settings
cpu_cores = :multi
# cpu_cores = :single
ngroups = 5 # should be your physical cpu cores minus 1
# groupsize = 50 # Used in the paper - for 250 replicates when ngroups == 5
groupsize = 1 # Fast - ngroups total replicates, to check the scripts

### Load data and packages ###

include(joinpath(basedir, "src", "spreadsetup.jl"))
using Optim, LossFunctions, JLD2, FileIO, DataStructures, LabelledArrays, DynamicGridsGtk

### Load common sctipt functions ###

# Main parametrisation method
function parametrise_ruleset!(optimresults, parametriser, key)
    println("\n\n==== Model: ", key, now(), " ===================================================")
    ruleset = parametriser.ruleset
    initpars = flatten(rules(ruleset))
    println(initpars)
    # Make a labelled array so we can ignore the order
    parnames = fieldnameflatten(rules(ruleset), Number)
    # Get the lower and upper bounds for params with flatten
    bnds = metaflatten(rules(ruleset), FieldMetadata.bounds, Number)
    lower = first.(bnds)
    upper = last.(bnds)
    names = fieldnameflatten(rules(ruleset), Number)
    namedpars = @LArray [initpars...] parnames
    display(DataFrame(names=[names...], pars=[initpars...], lower=[lower...], upper=[upper...]))
    # Zero params to make sure parametrizer is working
    output.running = false
    ruleset.rules = reconstruct(ruleset.rules, initpars .* 0)
    println("Running optimiser...")
    @time res = Optim.optimize(
        parametriser, [lower...], [upper...], namedpars, SAMIN(),
        Optim.Options(iterations=iterations, show_trace=true, store_trace=true)
    )
    optimresults[key] = res
    println("Saving JLD2 file")
    save("output/optimresults_$key.jld2", "optimresults", optimresults)
    println("Calculating losses...")
    optimresults
end

# Method to add results to a dataframe table
function ruleset_results!(resultdf, rulesets, parametriser, params, location, scenario=:regular)
    for key in keys(params)
        println("\n==== Model: ", key, now(), " ===================================================")
        ruleparams = params[key]
        # Set the parameters to zero, to make sure they are updated in the optimizer
        ruleset = deepcopy(rulesets[key])
        ruleset.rules = reconstruct(ruleset.rules, ruleparams .* 0)
        parametriser = @set parametriser.ruleset = ruleset
        loss = parametriser(ruleparams)
        # Standard deviation of accuracy
        std_ = std(accuracy.(Ref(Dispersal.targets(parametriser.objective)),
                             vcat(parametriser.results...)))
        row = (
            location=location,
            model=key,
            scenario=scenario,
            loss=loss,
            std=std_,
            accuracy=accuracy(Dispersal.targets(parametriser.objective), loss)
        )
        push!(resultdf, row)
    end
end
accuracy(target, loss) = one(loss) - loss/length(target)

# Plot accuracy heatmap
function lossplot(output::RegionOutput, ruleset)
    sim!(output, ruleset)
    diff = (transformfunc.(output[1]) - transformfunc.(output.objective.occurrence))
    sleep(0.5)
    println("Accuracy: $(1 - (sum(abs.(diff)) * 0.5 / length(diff)))")
    heatmap(diff; clims=(-1,1))
end

# Use ColorRegionFit in GtkOutput
function regionsim(ruleset, objective::RegionObjective,
        mask, popgrid, growthrates, tspan, carrycap; fps=20, outputtype=GtkOutput, filename="output/regionsim.gif")
    truescheme = Greyscale()
    falsescheme = reverse(ColorSchemes.Reds_3)
    truezerocolor =  (0.5, 0.5, 0.5)
    falsezerocolor = falsescheme[1]
    maskcolor =      (0.0, 1.0, 1.0)
    region_processor = ColorRegionFit(
        objective, truescheme, falsescheme,
        truezerocolor, falsezerocolor, maskcolor
    )
    gtkoutput = outputtype((population=popgrid,); filename=filename, mask=BitArray(collect(mask)), aux=(growthrates=growthrates,),
        store=true,
        tspan=tspan,
        processor=region_processor,
        minval=zero(carrycap),
        maxval=carrycap,
        fps=fps,
    )
    sim!(gtkoutput, ruleset)
    gtkoutput
end

# Sensitivity routines
function growthrate_sensitivity!(resultdf, rulesets, parametriser, params, stressors, change, selectors, location)
    output = parametriser.output
    crs_ = if output isa Output
        DynamicGrids.aux(output)[:growthrates]
    elseif output isa Vector{<:Output}
        DynamicGrids.aux(first(output))[:growthrates]
    end |> crs
    for s in stressors, c in change
        scenario = Symbol(s, :_, c)
        println(scenario)
        # Load growthrates file matching this scenario
        growthratesfilepath = joinpath(basedir, "output/sensitivity/growthrates_$(scenario).ncd")
        isfile(growthratesfilepath) || error("$growthratesfilepath not found")
        growthrates = NCDarray(growthratesfilepath; crs=crs_, dimcrs=EPSG(4326))[selectors...]  |>
            x -> setdims(x, (@set dims(x, Ti).mode.span = Regular(Month(1)))) |>
            x -> convertmode(Projected, x) |>
            x -> permutedims(x, (Lat, Lon, Ti)) |>
            x -> reverse(x; dims=Lat) |>
            x -> replace_missing(x, 0)
        # set the output aux to the new growthrates
        crs_ = if output isa Output
            @set! parametriser.output.extent.aux = (growthrates=growthrates,)
        elseif output isa Vector{<:Output}
            newoutputs = map(output) do o
                @set! o.extent.aux = (growthrates=growthrates,)
            end
            @set! parametriser.output = newoutputs
        end
        println("\n\n\n######### scenario: $scenario ########################################")
        ruleset_results!(resultdf, rulesets, parametriser, params, location, scenario)
    end
end

function param_sensitivity!(resultdf, rulesets, parametrizer, params, changekeys, change, location)
    # Rebuild all rulesets with adjusted params
    for paramkey in changekeys, c in change
        scenario = Symbol(paramkey, :_, c)
        modified_rulesets = Dict()
        for rulesetkey in keys(rulesets)
            rs = rulesets[rulesetkey]
            names = fieldnameflatten(rules(rs), Number)
            pars = NamedTuple{names}(flatten(rules(rs), Number))
            if paramkey in names
                ruleset = deepcopy(rs)
                # multiply the parameter by the amount of change
                pars = @set pars[paramkey] = pars[paramkey] * c
                # rebuild
                ruleset.rules = reconstruct(rules(ruleset), pars, Number)
                modified_rulesets[rulesetkey] = ruleset
            end
        end
        println("\n\n\n######### scenario: $scenario% ########################################")
        ruleset_results!(resultdf, modified_rulesets, parametrizer, params, location, scenario)
    end
end

function years_detected(df, years, nstates)
    occurrence = zeros(Bool, nstates, length(years))
    for y in eachindex(years), state in 1:nstates
        year = years[y]
        occurrence[state, y] = year >= df.year_detected[state]
    end
    occurrence
end

# Load US data to optimise the model:

timestep = Month(1)
us_df = CSV.read(joinpath(basedir, "data/region_year_detected_US.csv"))
us_years = Dates.year(first(us_tspan)):Dates.year(last(us_tspan))
us_nstates = round(Int, maximum(us_states))
us_occurrence = years_detected(us_df, us_years, us_nstates)

heatmap(us_occurrence; title="US occurrence", xlab="Year", ylab="State ID")

# Define model combinations for comparison

kwargs = (
    timestep=timestep,
    opt=SparseOpt(),
)
us_stochastic_rulesets = (
    full=Ruleset(us_humandisp, Chain(localdisp, allee, growth); kwargs...),
    nolocal=Ruleset(us_humandisp, Chain(allee, growth); kwargs...),
    noallee=Ruleset(us_humandisp, Chain(localdisp, growth); kwargs...),
    noclimate = Ruleset(us_humandisp, Chain(localdisp, allee, constant_growth); kwargs...),
)
null_rulesets = (
    null=Ruleset(; kwargs...),
)
null_params = (null=(),)

## Set up the optimiser

framesperstep = 12
startmonth = 5
transformfunc = x -> 2x - 1 # Transform to -1, +1 instead of 0, 1
lossfunc = ZeroOneLoss()
iterations = 1000
detectionthreshold = floattype(1e7)
threading = if cpu_cores == :multi
    ThreadedReplicates()
else
    SingleCoreReplicates()
end

us_objective = RegionObjective(
    detectionthreshold,
    parent(us_states),
    us_occurrence,
    framesperstep,
    startmonth,
)
us_aux = (growthrates=us_growthrates_zeromissing,)
us_ruleset = us_stochastic_rulesets[:full]
#us_ruleset = us_deterministic_rulesets[:nomodel]
us_output = RegionOutput((population=parent(us_populationgrid),);
    aux=us_aux,
    mask=BitArray(us_boolmask),
    tspan=us_tspan,
    objective=us_objective,
)

# rs = us_full_parametrized # If you have loaded the parametrised version
# rs = us_stochastic_rulesets[:full] # The initial version
# using DynamicGridsInteract
# InteractOutput((population=parent(us_populationgrid),);
    # ruleset=us_full_parametrized,
    # aux=us_aux,
    # mask=BitArray(us_boolmask),
    # tspan=us_tspan,
    # maxval=carrycap,
    # processor=inferno,
# ) |> display

#' ## Profile model performance
#'
#' It's good to check the everythign is running well
#' before starting a 24 hour simulation:
#'
#+

# simdata = DynamicGrids.SimData(us_output.extent, us_ruleset)
# using BenchmarkTools
# @btime sim!($us_output, $us_ruleset)
# @profiler for i = 1:10 sim!(us_output, us_ruleset; simdata=simdata) end
# isinferred(us_output, us_ruleset)

## Run the optimizer

# **WARNING** this might take a few days if there is no params file already saved

stochparamsfile = "output/stochparams_$humanactivity.jld2"
optimresultsfile = "output/optimresults_$humanactivity.jld2"
optimresults = OrderedDict()
if isfile(stochparamsfile)
    stochparams = load(stochparamsfile, "stochparams_$humanactivity")
    stochparams = load(stochparamsfile, "stochparams_$humanactivity")
    stochrulesetkeys = keys(stochparams)
    stochparamvals = values(stochparams)
else
    for key in keys(us_stochastic_rulesets)
        us_parametriser = Parametriser(
            us_stochastic_rulesets[key], us_output, us_objective, transformfunc, lossfunc, ngroups, groupsize, threading
        )
        parametrise_ruleset!(optimresults, us_parametriser, key)
        # Save as we go
        stochrulesetkeys = Tuple(keys(optimresults))
        stochparamvals = ((Optim.minimizer(optimresults[k]) for k in stochrulesetkeys)...,)
        stochparams = NamedTuple{stochrulesetkeys}(stochparamvals)
        save(optimresultsfile, "optimresults", optimresults)
        save("output/stochparams_$humanactivity.jld2", "stochparams_$humanactivity", stochparams)
    end
    optimresults
end

# Test the output:

us_full_paremetrized = reconstruct(us_stochastic_rulesets[:full], stochparams[:full])
lossplot(us_output, us_full_paremetrized)
# regionoutput = regionsim(us_full_paremetrized, us_objective, us_boolmask, us_populationgrid,
    # us_growthrates_zeromissing, us_tspan, carrycap)

# Add US results to dataframes

# Add parameters to DataFrame
allparamkeys = reduce(union, map(keys, stochparams); init=Symbol[])
colnames = [:lower, :upper, stochrulesetkeys...]
missingparams = [Any[missing for pk in allparamkeys] for rk in colnames]
paramdf = DataFrame([allparamkeys, missingparams...], [:param, colnames...])
fullbounds = metaflatten(us_stochastic_rulesets[:full].rules, FieldMetadata.bounds, Union{Real,Quantity})
fullkeys = fieldnameflatten(us_stochastic_rulesets[:full].rules, Union{Real,Quantity})
for key in stochrulesetkeys
    stochrulesetresults = stochparams[key]
    for (row, rowkey) in enumerate(allparamkeys)
        if rowkey in keys(stochparams[key])
            paramdf[key][row] = stochrulesetresults[rowkey]
        end
    end
end
# Add param bounds, using :full model
for (row, rowkey) in enumerate(allparamkeys)
    if rowkey in fullkeys
        paramdf[:lower][row], paramdf[:upper][row] = fullbounds[row]
    end
end
# Set intrinsicrate bounds from ExactLogisticGrowth model
select = paramdf.param .== :intrinsicrate
paramdf[select, :lower], paramdf[select, :upper] =
    FieldMetadata.bounds(ExactLogisticGrowth, :intrinsicrate)

# Add parameter names
stochparamnames = keys(stochparams[:full])

### Get the results in a dataframe: ###

# Loss/accuracy results dataframe
resultdf = DataFrame(location=[], model=[], scenario=[], loss=Float64[], std=Float64[], accuracy=Float64[])
us_parametriser = Parametriser(
    first(values(us_stochastic_rulesets)), us_output, us_objective, transformfunc,
    lossfunc, ngroups, groupsize, threading
)
us_null_parametriser = Parametriser(
    first(values(null_rulesets)), us_output, us_objective, transformfunc,
    lossfunc, ngroups, groupsize, threading
)
us_location = :US

# Get the results for all stochastic rulsets, and the null rulesets:

ruleset_results!(resultdf, us_stochastic_rulesets, us_parametriser, stochparams, us_location)
ruleset_results!(resultdf, null_rulesets, us_null_parametriser, null_params, us_location)

# Show the output in the plot pane:

display(paramdf)
display(resultdf)
showtable(resultdf)

stressors = (:heatstress, :coldstress, :wiltstress, :wetstress)
paramnames = (:human_exponent, :dist_exponent, :dispersalperpop, :max_dispersers)
change = 80, 120
selectors = us
growthrate_sensitivity!(resultdf, us_stochastic_rulesets, us_parametriser, stochparams, stressors, change, selectors, us_location)
param_sensitivity!(resultdf, us_stochastic_rulesets, us_parametriser, stochparams, paramnames, change, us_location)
showtable(resultdf)

# Load data for validation against the EU

eu_stochastic_rulesets = (
    full=Ruleset(eu_humandisp, Chain(localdisp, allee, growth); kwargs...),
    nolocal=Ruleset(eu_humandisp, Chain(allee, growth); kwargs...),
    noallee=Ruleset(eu_humandisp, Chain(localdisp, growth); kwargs...),
    noclimate = Ruleset(eu_humandisp, Chain(localdisp, allee, constant_growth); kwargs...)
)
eu_rulesets = map((r, p) -> reconstruct(r, p, Number), eu_stochastic_rulesets, stochparams)
eu_df = CSV.read(joinpath(basedir, "data/region_year_detected_EU.csv"))
eu_years = Dates.year(first(eu_tspan)):Dates.year(last(eu_tspan))
eu_occurrence = years_detected(eu_df, eu_years, maximum(eu_states))
eu_objective = RegionObjective(detectionthreshold, eu_states, eu_occurrence, framesperstep, startmonth)
eu_aux = (growthrates=eu_growthrates_zeromissing,)
eu_init = (population=eu_populationgrid,)
eu_output = RegionOutput(eu_init;
    aux=eu_aux,
    mask=eu_boolmask,
    tspan=eu_tspan,
    objective=eu_objective,
)
eu_parametriser = Parametriser(
    first(values(eu_rulesets)), eu_output, eu_objective,
    transformfunc, lossfunc, ngroups, groupsize, threading
)
eu_null_parametriser = Parametriser(
    first(values(null_rulesets)), eu_output, eu_objective, transformfunc,
    lossfunc, ngroups, groupsize, threading
)
eu_location = :EU

# Test output

heatmap(eu_occurrence; title="EU occurrence", xlab="Year", ylab="State ID")
eu_output.extent.init.population |> plot
lossplot(eu_output, null_rulesets[:null])
savefig("output/eu_null_loss.png")

eu_full = reconstruct(eu_stochastic_rulesets[:full], stochparams[:full])
lossplot(eu_output, eu_full)
savefig("output/eu_full_loss.png")
# regionoutput = regionsim(eu_full, eu_objective, eu_boolmask, eu_populationgrid,
          # eu_growthrates_zeromissing, eu_tspan, carrycap; fps=10)#, outputtype=GifOutput, filename="eu_gdp.gif")

### Validate US/EU ###

selectors = eu
ruleset_results!(resultdf, null_rulesets, eu_null_parametriser, null_params, eu_location)
ruleset_results!(resultdf, eu_stochastic_rulesets, eu_parametriser, stochparams, eu_location)
growthrate_sensitivity!(resultdf, eu_stochastic_rulesets, eu_parametriser, stochparams, stressors, change, selectors, eu_location)
param_sensitivity!(resultdf, eu_stochastic_rulesets, eu_parametriser, stochparams, paramnames, change, eu_location)
display(resultdf)
@save "output/resultdf_$(humanactivity).jld2" resultdf
CSV.write("output/sensitivity_$(humanactivity).csv", resultdf)
CSV.write("output/params_$(humanactivity).csv", paramdf)

### NoHuman model ###

# This needs to be separated out so that we can
# optimise the no-human model with a free λ parameter.

println("==== Running nohuman model ==============================================")

λ = floattype(0.0125)
radius = 6

@time hood = DispersalKernel{radius}(
    formulation=ExponentialKernel(λ),
    distancemethod=AreaToArea(30),
    cellsize=floattype(1.0)
)

free_localdisp = InwardsPopulationDispersal{:population,:population}(hood)
log.(hood.kernel) |> heatmap

nohuman = Ruleset(Chain(free_localdisp, allee, growth); kwargs...)
nohuman_rulesets = (nohuman=nohuman,)
# Parametrise λ - it's now a general spread model, not local movement
@flattenable @bounds ExponentialKernel begin
    λ | true | (0.0, 100.0)
end

# This model is deteministic so we don't need replicates
ngroups = 1
groupsize = 1
threading = SingleCoreReplicates()
us_nohuman_parametriser = Parametriser(
    nohuman, us_output, us_objective, transformfunc, lossfunc, ngroups, groupsize, threading
)
eu_nohuman_parametriser = Parametriser(
    nohuman, eu_output, eu_objective, transformfunc, lossfunc, ngroups, groupsize, threading
)

# Parametrise unless an output is already saved
nohumanfile = joinpath(basedir, "output/params_$(humanactivity)_nohuman.jld2")
if isfile(nohumanfile)
    nohuman_params = load(nohumanfile, "nohuman_params")
    nohuman_paramvals = values(nohuman_params)
else
    parametrise_ruleset!(optimresults, us_nohuman_parametriser, :nohuman)
    nohuman_rulesetkeys = (:nohuman,)
    nohuman_paramvals = Optim.minimizer(optimresults[:nohuman])
    nohuman_params = NamedTuple{(:nohuman,)}((nohuman_paramvals,))
end

parametrized_nohuman = reconstruct(nohuman, nohuman_paramvals, Number)
lossplot(us_output, parametrized_nohuman)
# regionoutput = regionsim(
    # parametrized_nohuman, us_objective, us_boolmask, us_populationgrid,
    # us_growthrates_zeromissing, us_tspan, carrycap;
    #outputtype=GifOutput, filename="output/us_nohuman_regionsim.gif", fps=10,
# )
lossplot(eu_output, nohuman)
# regionoutput = regionsim(
    # parametrized_nohuman, eu_objective, eu_boolmask, eu_populationgrid,
    # eu_growthrates_zeromissing, eu_tspan, carrycap;
    #outputtype=GifOutput, filename="output/eu_nohuman_regionsim.gif", fps=10,
# )

# Regular results
ruleset_results!(resultdf, nohuman_rulesets, us_nohuman_parametriser, nohuman_params, us_location)
ruleset_results!(resultdf, nohuman_rulesets, eu_nohuman_parametriser, nohuman_params, eu_location)
# Growth rate sensitivity
growthrate_sensitivity!(resultdf, nohuman_rulesets, us_nohuman_parametriser, nohuman_params, stressors, change, selectors, us_location)
growthrate_sensitivity!(resultdf, nohuman_rulesets, eu_nohuman_parametriser, nohuman_params, stressors, change, selectors, eu_location)
# Dispersal parameter sensitivity
nohuman_paramnames = keys(nohuman_params[:nohuman])
param_sensitivity!(resultdf, nohuman_rulesets, us_nohuman_parametriser, nohuman_params, nohuman_paramnames, change, us_location)
param_sensitivity!(resultdf, nohuman_rulesets, eu_nohuman_parametriser, nohuman_params, nohuman_paramnames, change, eu_location)

nohumancol = deepcopy(paramdf.full) .= missing
insert!(paramdf, 8, nohumancol, :nohuman)
select = paramdf.param .== :minfounders
paramdf[select, :nohuman] = nohuman_params[1].minfounders
newrow = [:λ, 0.0, 100.0, missing, missing, missing, missing, nohuman_params[1].λ]
push!(paramdf, newrow)
showtable(paramdf)
showtable(resultdf)
display(paramdf)
display(resultdf)

# Save results as a csv file

@save "output/params_$(humanactivity)_nohuman.jld2" nohuman_params
CSV.write(joinpath(basedir, "output/sensitivity_$(humanactivity)_nohuman.csv"), resultdf)
CSV.write(joinpath(basedir, "output/params_$(humanactivity)_nohuman.csv"), paramdf)

max_stars_repo_name: rafaqz/DSuzukiiInvasionPaper

max_issues_repo_name: rafaqz/DSuzukiiInvasionPaper


text:
using LinearAlgebra, Plots, ECOS, Convex, DelimitedFiles, Test, PandaMPC

import PandaMPC

#Initial conditions
x0 = -1;
y0 = 0;
z0 = 0;
v0x = 0.9;
v0y = 0.2;
v0z = 4;
cond0 = [x0; y0; z0; v0x; v0y; v0z];
simtimeball = 1; # Simulation time

# Here we get the trajectory of the ball
t, posdrag, posnodrag = ball(cond0, simtimeball)
xplane = 0.4; #distance of grabbing to the ball from the base
xcapture, ycapture, zcapture, tcapture = intercplane(t, posdrag, xplane)

pesoQ = 1 #state penalization weight
pesoR = 0.001 #input penalization weight
m = 1 #number of joints

Q = [0 0 0; 0 1 0; 0 0 1]
R = pesoR
@test (Q,R) == weights(pesoQ, pesoR, m)
@test (xcapture, ycapture, zcapture, tcapture) == intercplane(t, posdrag, xplane)

span = 1
resample = 0.02
final_z =  [1; 0; 0; 1; 0; 0; 1; 0; 0; 1; 0; 0; 1; 0; 0; 1; 0; 0 ; 1; 0; 0]
initial_x = [0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0]
initial_u = [0; 0; 0; 0; 0; 0; 0]
initial_z = [0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0]
pesoQ = 1
pesoR = 0.001
m = 7;  n = 21; l = 21
T = 15+1 # Number of timesteps/ Prediction (and Control) Horizon

tfin, zfin, ufin = go7(span, final_z, initial_x, initial_u, initial_z, resample, T, n, l, m, pesoQ, pesoR)

a = (abs.(zfin[:,end]).-final_z).<0.01

@test a[1]
@test a[2]
@test a[3]
@test a[4]
@test a[5]
@test a[6]
@test a[7]
@test a[8]
@test a[9]
@test a[10]
@test a[11]
@test a[12]
@test a[13]
@test a[14]
@test a[15]
@test a[16]
@test a[17]
@test a[18]
@test a[19]
@test a[20]
@test a[21]


text:
using Flux
using Flux: @treelike

struct ResidualBlock{L,S}
  layers::L
  shortcut::S
end

struct ConvNorm{C, N}
  conv::C
  norm::N
end
(cn::ConvNorm)(value) = cn.norm(cn.conv(value))

@treelike ResidualBlock
@treelike ConvNorm

function ResidualBlock(filters, kernels::Array{Tuple{Int,Int}}, pads::Array{Tuple{Int,Int}}, strides::Array{Tuple{Int,Int}}, shortcut = identity)
  local layers = []
  for i in 2:length(filters)
    push!(layers, ConvNorm(
        Conv(kernels[i-1], filters[i-1]=>filters[i], pad = pads[i-1], stride = strides[i-1]),
        BatchNorm(filters[i])))
  end
  ResidualBlock(Tuple(layers),shortcut)
end

function ResidualBlock(filters, kernels::Array{Int}, pads::Array{Int}, strides::Array{Int}, shortcut = identity)
  @assert length(filters) == 4
  ResidualBlock(filters, [(i,i) for i in kernels], [(i,i) for i in pads], [(i,i) for i in strides], shortcut)
end

# This is a hack, to avoid Zygote's trouble with destructuring
rb_apply(t::Tuple{Any}, value, shortcut) = relu.(first(t)(value) + shortcut)
rb_apply(t::Tuple, value, shortcut) = rb_apply(Base.tail(t), relu.(first(t)(value)), shortcut)
function (block::ResidualBlock)(input)
  rb_apply(block.layers, copy(input), block.shortcut(input))
end

function Bottleneck(filters::Int, downsample::Bool = false, res_top::Bool = false)
  if(!downsample && !res_top)
    return ResidualBlock([4 * filters, filters, filters, 4 * filters], [1,3,1], [0,1,0], [1,1,1])
  elseif(downsample && res_top)
    return ResidualBlock([filters, filters, filters, 4 * filters], [1,3,1], [0,1,0], [1,1,1], Chain(Conv((1,1), filters=>4 * filters, pad = (0,0), stride = (1,1)), BatchNorm(4 * filters)))
  else
    shortcut = Chain(Conv((1,1), 2 * filters=>4 * filters, pad = (0,0), stride = (2,2)), BatchNorm(4 * filters))
    return ResidualBlock([2 * filters, filters, filters, 4 * filters], [1,3,1], [0,1,0], [1,1,2], shortcut)
  end
end

function resnet50()
  local layers = [3, 4, 6, 3]
  local layer_arr = []

  push!(layer_arr, Conv((7,7), 3=>64, pad = (3,3), stride = (2,2)))
  push!(layer_arr, x -> maxpool(x, (3,3), pad = (1,1), stride = (2,2)))

  initial_filters = 64
  for i in 1:length(layers)
    push!(layer_arr, Bottleneck(initial_filters, true, i==1))
    for j in 2:layers[i]
      push!(layer_arr, Bottleneck(initial_filters))
    end
    initial_filters *= 2
  end

  push!(layer_arr, x -> meanpool(x, (7,7)))
  push!(layer_arr, x -> reshape(x, :, size(x,4)))
  push!(layer_arr, (Dense(2048, 1000)))

  # Disable softmax, use logitcrossentropy instead
  #push!(layer_arr, softmax)

  Chain(layer_arr...)
end