text:
# using LibGit2
import Pkg
using Test
using Pluto.Configuration: CompilerOptions
using Pluto.WorkspaceManager: _merge_notebook_compiler_options
import Pluto: update_save_run!, update_run!, WorkspaceManager, ClientSession, ServerSession, Notebook, Cell, project_relative_path, SessionActions, load_notebook
import Pluto.PkgUtils
import Pluto.PkgCompat
import Distributed

const pluto_test_registry_spec = Pkg.RegistrySpec(;
    url="https://github.com/JuliaPluto/PlutoPkgTestRegistry",
    uuid=Base.UUID("96d04d5f-8721-475f-89c4-5ee455d3eda0"),
    name="PlutoPkgTestRegistry",
)

@testset "Built-in Pkg" begin

    # Pkg.Registry.rm("General")
    Pkg.Registry.add(pluto_test_registry_spec)

    @testset "Basic" begin
        fakeclient = ClientSession(:fake, nothing)
        🍭 = ServerSession()
        🍭.connected_clients[fakeclient.id] = fakeclient

        # See https://github.com/JuliaPluto/PlutoPkgTestRegistry

        notebook = Notebook([
            Cell("import PlutoPkgTestA"), # cell 1
            Cell("PlutoPkgTestA.MY_VERSION |> Text"),
            Cell("import PlutoPkgTestB"), # cell 3
            Cell("PlutoPkgTestB.MY_VERSION |> Text"),
            Cell("import PlutoPkgTestC"), # cell 5
            Cell("PlutoPkgTestC.MY_VERSION |> Text"),
            Cell("import PlutoPkgTestD"), # cell 7
            Cell("PlutoPkgTestD.MY_VERSION |> Text"),
            Cell("import Dates"),
            # eval to hide the import from Pluto's analysis
            Cell("eval(:(import DataFrames))")
        ])
        fakeclient.connected_notebook = notebook

        @test !notebook.nbpkg_ctx_instantiated

        update_save_run!(🍭, notebook, notebook.cells[[1, 2, 7, 8]]) # import A and D
        @test notebook.cells[1].errored == false
        @test notebook.cells[2].errored == false
        @test notebook.cells[7].errored == false
        @test notebook.cells[8].errored == false

        @test notebook.nbpkg_ctx !== nothing
        @test notebook.nbpkg_restart_recommended_msg === nothing
        @test notebook.nbpkg_restart_required_msg === nothing
        @test notebook.nbpkg_ctx_instantiated
        @test notebook.nbpkg_busy_packages |> isempty

        terminals = notebook.nbpkg_terminal_outputs

        @test haskey(terminals, "PlutoPkgTestA")
        @test haskey(terminals, "PlutoPkgTestD")
        # they were installed in one batch, so their terminal outputs should be the same
        @test terminals["PlutoPkgTestA"] == terminals["PlutoPkgTestD"]

        @test notebook.cells[2].output.body == "0.3.1" # A
        @test notebook.cells[8].output.body == "0.1.0" # D

        @test PkgCompat.get_manifest_version(notebook.nbpkg_ctx, "PlutoPkgTestA") == v"0.3.1"
        @test PkgCompat.get_manifest_version(notebook.nbpkg_ctx, "PlutoPkgTestD") == v"0.1.0"

        old_A_terminal = deepcopy(terminals["PlutoPkgTestA"])
        @show old_A_terminal

        update_save_run!(🍭, notebook, notebook.cells[[3, 4]]) # import B

        @test notebook.cells[3].errored == false
        @test notebook.cells[4].errored == false

        @test notebook.nbpkg_ctx !== nothing
        @test notebook.nbpkg_restart_recommended_msg === nothing
        @test notebook.nbpkg_restart_required_msg === nothing

        @test haskey(terminals, "PlutoPkgTestB")
        @test terminals["PlutoPkgTestA"] == terminals["PlutoPkgTestD"] == old_A_terminal

        @test terminals["PlutoPkgTestA"] != terminals["PlutoPkgTestB"]

        @test notebook.cells[4].output.body == "1.0.0" # B

        # running the 5th cell will import PlutoPkgTestC, putting a 0.2 compatibility bound on PlutoPkgTestA. This means that a notebook restart is required, since PlutoPkgTestA was already loaded at version 0.3.1.
        update_save_run!(🍭, notebook, notebook.cells[[5, 6]])

        @test notebook.cells[5].errored == false
        @test notebook.cells[6].errored == false

        @test notebook.nbpkg_ctx !== nothing
        @test (
            notebook.nbpkg_restart_recommended_msg !==  nothing || notebook.nbpkg_restart_required_msg !== nothing
        )
        @test notebook.nbpkg_restart_required_msg !== nothing

        # running cells again should persist the restart message

        update_save_run!(🍭, notebook, notebook.cells[1:8])
        @test notebook.nbpkg_restart_required_msg !== nothing

        Pluto.response_restrart_process(Pluto.ClientRequest(
            session=🍭,
            notebook=notebook,
        ); run_async=false)

        # @test_nowarn SessionActions.shutdown(🍭, notebook; keep_in_session=true, async=true)
        # @test_nowarn update_save_run!(🍭, notebook, notebook.cells[1:8]; , save=true)

        @test notebook.cells[1].errored == false
        @test notebook.cells[2].errored == false
        @test notebook.cells[3].errored == false
        @test notebook.cells[4].errored == false
        @test notebook.cells[5].errored == false
        @test notebook.cells[6].errored == false
        @test notebook.cells[7].errored == false
        @test notebook.cells[8].errored == false

        @test notebook.nbpkg_ctx !== nothing
        @test notebook.nbpkg_restart_recommended_msg === nothing
        @test notebook.nbpkg_restart_required_msg === nothing

        @test notebook.cells[2].output.body == "0.2.2"
        @test notebook.cells[4].output.body == "1.0.0"
        @test notebook.cells[6].output.body == "1.0.0"
        @test notebook.cells[8].output.body == "0.1.0"

        update_save_run!(🍭, notebook, notebook.cells[9])

        @test notebook.cells[9].errored == false
        @test notebook.nbpkg_ctx !== nothing
        @test notebook.nbpkg_restart_recommended_msg === nothing
        @test notebook.nbpkg_restart_required_msg === nothing

        # we should have an isolated environment, so importing DataFrames should not work, even though it is available in the parent process.
        update_save_run!(🍭, notebook, notebook.cells[10])
        @test notebook.cells[10].errored == true

        ptoml_file() = PkgCompat.project_file(notebook)
        mtoml_file() = PkgCompat.manifest_file(notebook)

        ptoml_contents() = read(ptoml_file(), String)
        mtoml_contents() = read(mtoml_file(), String)

        nb_contents() = read(notebook.path, String)

        @testset "Project & Manifest stored in notebook" begin

            @test occursin(ptoml_contents(), nb_contents())
            @test occursin(mtoml_contents(), nb_contents())

            @test occursin("PlutoPkgTestA", mtoml_contents())
            @test occursin("PlutoPkgTestB", mtoml_contents())
            @test occursin("PlutoPkgTestC", mtoml_contents())
            @test occursin("PlutoPkgTestD", mtoml_contents())
            @test occursin("Dates", mtoml_contents())
            @test count("PlutoPkgTestA", ptoml_contents()) == 2 # once in [deps], once in [compat]
            @test count("PlutoPkgTestB", ptoml_contents()) == 2
            @test count("PlutoPkgTestC", ptoml_contents()) == 2
            @test count("PlutoPkgTestD", ptoml_contents()) == 2
            @test count("Dates", ptoml_contents()) == 1 # once in [deps], but not in [compat] because it is a stdlib

            ptoml = Pkg.TOML.parse(ptoml_contents())

            @test haskey(ptoml["compat"], "PlutoPkgTestA")
            @test haskey(ptoml["compat"], "PlutoPkgTestB")
            @test haskey(ptoml["compat"], "PlutoPkgTestC")
            @test haskey(ptoml["compat"], "PlutoPkgTestD")
            @test !haskey(ptoml["compat"], "Dates")
        end

        ## remove `import Dates`
        setcode(notebook.cells[9], "")
        update_save_run!(🍭, notebook, notebook.cells[9])

        # removing a stdlib does not require a restart
        @test notebook.cells[9].errored == false
        @test notebook.nbpkg_ctx !== nothing
        @test notebook.nbpkg_restart_recommended_msg === nothing
        @test notebook.nbpkg_restart_required_msg === nothing

        @test count("Dates", ptoml_contents()) == 0

        ## remove `import PlutoPkgTestD`
        setcode(notebook.cells[7], "")
        update_save_run!(🍭, notebook, notebook.cells[7])

        @test notebook.cells[7].errored == false
        @test notebook.nbpkg_ctx !== nothing
        @test notebook.nbpkg_restart_recommended_msg !== nothing # recommend restart
        @test notebook.nbpkg_restart_required_msg === nothing

        @test count("PlutoPkgTestD", ptoml_contents()) == 0

        WorkspaceManager.unmake_workspace((🍭, notebook))
    end

    simple_import_path = joinpath(@__DIR__, "simple_import.jl")
    simple_import_notebook = read(simple_import_path, String)

    @testset "Manifest loading" begin
        fakeclient = ClientSession(:fake, nothing)
        🍭 = ServerSession()
        🍭.connected_clients[fakeclient.id] = fakeclient

        dir = mktempdir()
        path = joinpath(dir, "hello.jl")
        write(path, simple_import_notebook)

        notebook = SessionActions.open(🍭, path; run_async=false)
        fakeclient.connected_notebook = notebook

        @test num_backups_in(dir) == 0

        @test notebook.nbpkg_ctx !== nothing
        @test notebook.nbpkg_restart_recommended_msg === nothing
        @test notebook.nbpkg_restart_required_msg === nothing

        @test notebook.cells[1].errored == false
        @test notebook.cells[2].errored == false

        @test notebook.cells[2].output.body == "0.2.2"

        WorkspaceManager.unmake_workspace((🍭, notebook))
    end

    @testset "Pkg cell" begin
        fakeclient = ClientSession(:fake, nothing)
        🍭 = ServerSession()
        🍭.connected_clients[fakeclient.id] = fakeclient

        notebook = Notebook([
            Cell("1"),
            Cell("2"),
            Cell("3"),
            Cell("4"),
            Cell("5"),
            Cell("6"),
        ])
        fakeclient.connected_notebook = notebook

        update_save_run!(🍭, notebook, notebook.cells)

        # not necessary since there are no packages:
        # @test has_embedded_pkgfiles(notebook)

        setcode(notebook.cells[1], "import Pkg")
        update_save_run!(🍭, notebook, notebook.cells[1])
        setcode(notebook.cells[2], "Pkg.activate(mktempdir())")
        update_save_run!(🍭, notebook, notebook.cells[2])

        @test notebook.cells[1].errored == false
        @test notebook.cells[2].errored == false
        @test notebook.nbpkg_ctx === nothing
        @test notebook.nbpkg_restart_recommended_msg === nothing
        @test notebook.nbpkg_restart_required_msg === nothing
        @test !has_embedded_pkgfiles(notebook)

        setcode(notebook.cells[3], "Pkg.add(\"JSON\")")
        update_save_run!(🍭, notebook, notebook.cells[3])
        setcode(notebook.cells[4], "using JSON")
        update_save_run!(🍭, notebook, notebook.cells[4])
        setcode(notebook.cells[5], "using Dates")
        update_save_run!(🍭, notebook, notebook.cells[5])

        @test notebook.cells[3].errored == false
        @test notebook.cells[4].errored == false
        @test notebook.cells[5].errored == false

        @test !has_embedded_pkgfiles(notebook)

        setcode(notebook.cells[2], "2")
        setcode(notebook.cells[3], "3")
        update_save_run!(🍭, notebook, notebook.cells[2:3])

        @test notebook.nbpkg_ctx !== nothing
        @test notebook.nbpkg_restart_required_msg !== nothing
        @test has_embedded_pkgfiles(notebook)

        WorkspaceManager.unmake_workspace((🍭, notebook))
    end

    pre_pkg_notebook = read(joinpath(@__DIR__, "old_import.jl"), String)

    local post_pkg_notebook = nothing

    @testset "Backwards compat" begin
        fakeclient = ClientSession(:fake, nothing)
        🍭 = ServerSession()
        🍭.connected_clients[fakeclient.id] = fakeclient

        dir = mktempdir()
        path = joinpath(dir, "hello.jl")
        write(path, pre_pkg_notebook)

        @test num_backups_in(dir) == 0

        notebook = SessionActions.open(🍭, path; run_async=false)
        fakeclient.connected_notebook = notebook
        nb_contents() = read(notebook.path, String)

        @test num_backups_in(dir) == 0
        # @test num_backups_in(dir) == 1

        post_pkg_notebook = read(path, String)

        # test that pkg cells got added
        @test length(post_pkg_notebook) > length(pre_pkg_notebook) + 50
        @test has_embedded_pkgfiles(notebook)

        @test notebook.nbpkg_ctx !== nothing
        @test notebook.nbpkg_restart_recommended_msg === nothing
        @test notebook.nbpkg_restart_required_msg === nothing

        WorkspaceManager.unmake_workspace((🍭, notebook))
    end

    @testset "Forwards compat" begin
        # Using Distributed, we will create a new Julia process in which we install Pluto 0.14.7 (before PlutoPkg). We run the new notebook file on the old Pluto.
        p = Distributed.addprocs(1) |> first

        @test post_pkg_notebook isa String

        Distributed.remotecall_eval(Main, p, quote
            path = tempname()
            write(path, $(post_pkg_notebook))
            import Pkg
            # optimization:
            if isdefined(Pkg, :UPDATED_REGISTRY_THIS_SESSION)
                Pkg.UPDATED_REGISTRY_THIS_SESSION[] = true
            end

            Pkg.activate(mktempdir())
            Pkg.add(Pkg.PackageSpec(;name="Pluto",version=v"0.14.7"))
            import Pluto
            @assert Pluto.PLUTO_VERSION == v"0.14.7"

            s = Pluto.ServerSession()
            s.options.evaluation.workspace_use_distributed = false

            nb = Pluto.SessionActions.open(s, path; run_async=false)

            nothing
        end)

        # Cells that use Example will error because the package is not installed.

        # @test Distributed.remotecall_eval(Main, p, quote
        #     nb.cells[1].errored == false
        # end)
        @test Distributed.remotecall_eval(Main, p, quote
            nb.cells[2].errored == false
        end)
        # @test Distributed.remotecall_eval(Main, p, quote
        #     nb.cells[3].errored == false
        # end)
        # @test Distributed.remotecall_eval(Main, p, quote
        #     nb.cells[3].output.body == "25"
        # end)

        Distributed.rmprocs([p])
    end

    @testset "PkgUtils -- reset" begin
        dir = mktempdir()
        f = joinpath(dir, "hello.jl")

        write(f, simple_import_notebook)

        @test num_backups_in(dir) == 0
        Pluto.reset_notebook_environment(f)

        @test num_backups_in(dir) == 1
        @test !has_embedded_pkgfiles(read(f, String))
    end

    @testset "PkgUtils -- update" begin
        dir = mktempdir()
        f = joinpath(dir, "hello.jl")

        write(f, simple_import_notebook)
        @test !occursin("0.3.1", read(f, String))

        @test num_backups_in(dir) == 0
        Pluto.update_notebook_environment(f)

        @test num_backups_in(dir) == 1
        @test has_embedded_pkgfiles(read(f, String))
        @test !Pluto.only_versions_differ(f, simple_import_path)
        @test occursin("0.3.1", read(f, String))

        Pluto.update_notebook_environment(f)
        @test_skip num_backups_in(dir) == 1
    end

    @testset "Bad files" begin
        @testset "$(name)" for name in ["corrupted_manifest", "unregistered_import"]

            original_path = joinpath(@__DIR__, "$(name).jl")
            original_contents = read(original_path, String)

            fakeclient = ClientSession(:fake, nothing)
            🍭 = ServerSession()
            🍭.connected_clients[fakeclient.id] = fakeclient

            dir = mktempdir()
            path = joinpath(dir, "hello.jl")
            write(path, original_contents)

            @test num_backups_in(dir) == 0

            notebook = SessionActions.open(🍭, path; run_async=false)
            fakeclient.connected_notebook = notebook
            nb_contents() = read(notebook.path, String)

            should_restart = (
                notebook.nbpkg_restart_recommended_msg !==  nothing || notebook.nbpkg_restart_required_msg !== nothing
            )

            # if name == "corrupted_manifest"
            #     @test !should_restart
            # end

            # this breaks julia for somee reason:
            #     # we don't want to recommend restart right after launch, but it's easier for us
            #     @test_broken !should_restart
            # end

            if should_restart
                Pluto.response_restrart_process(Pluto.ClientRequest(
                    session=🍭,
                    notebook=notebook,
                ); run_async=false)
            end

            if name != "unregistered_import"
                @test notebook.cells[1].errored == false
                @test notebook.cells[2].errored == false
                @test notebook.cells[2].output.body == "0.2.2" # the Project.toml remained, so we did not lose our compat bound.
                @test has_embedded_pkgfiles(notebook)
            end

            @test !Pluto.only_versions_differ(notebook.path, original_path)

            @test notebook.nbpkg_ctx !== nothing
            @test notebook.nbpkg_restart_recommended_msg === nothing
            @test notebook.nbpkg_restart_required_msg === nothing

            setcode(notebook.cells[2], "1 + 1")
            update_save_run!(🍭, notebook, notebook.cells[2])
            @test notebook.cells[2].output.body == "2"

            setcode(notebook.cells[2], """
            begin
                import PlutoPkgTestD
                PlutoPkgTestD.MY_VERSION |> Text
            end
            """)
            update_save_run!(🍭, notebook, notebook.cells[2])
            @test notebook.cells[2].output.body == "0.1.0"

            @test has_embedded_pkgfiles(notebook)

            WorkspaceManager.unmake_workspace((🍭, notebook))
        end

    end

    # @test false

    # @testset "File format" begin
    #     notebook = Notebook([
    #         Cell("import PlutoPkgTestA"), # cell 1
    #         Cell("PlutoPkgTestA.MY_VERSION |> Text"),
    #         Cell("import PlutoPkgTestB"), # cell 3
    #         Cell("PlutoPkgTestB.MY_VERSION |> Text"),
    #         Cell("import PlutoPkgTestC"), # cell 5
    #         Cell("PlutoPkgTestC.MY_VERSION |> Text"),
    #         Cell("import PlutoPkgTestD"), # cell 7
    #         Cell("PlutoPkgTestD.MY_VERSION |> Text"),
    #         Cell("import PlutoPkgTestE"), # cell 9
    #         Cell("PlutoPkgTestE.MY_VERSION |> Text"),
    #         Cell("eval(:(import DataFrames))")
    #     ])

    #     file1 = tempname()
    #     notebook.path = file1

    #     save_notebook()

    #     save_notebook
    # end

    Pkg.Registry.rm(pluto_test_registry_spec)
    # Pkg.Registry.add("General")
end

# reg_path = mktempdir()
# repo = LibGit2.clone("https://github.com/JuliaRegistries/General.git", reg_path)

# LibGit2.checkout!(repo, "aef26d37e1d0e8f8387c011ccb7c4a38398a18f6")


text:
##
# Based of https://github.com/gpudirect/libmp/blob/master/examples/mp_sendrecv_kernel.cu
##
import MPI

using Hydrozoa
using Hydrozoa.Device
using CUDAnative

const MAX_SIZE = 64*1024
const ITER_COUNT_SMALL = 20
const ITER_COUNT_LARGE = 10
const WINDOW_SIZE = 64

struct Comms{N}
    tx::NTuple{N, SendDesc}
    tx_wait::NTuple{N, WaitDesc}
    rx_wait::NTuple{N, WaitDesc}
end

function exchange_kernel(rank, descs::Comms, N)
    @cuassert(gridDim().x == 1)
    for i in 1:N
        if rank == 0
            if threadIdx().x == 1
                @cuprintf("i=%ld, send+recv", i)
                threadfence()
                send(descs.tx[i])
                wait(descs.tx_wait[i])
                signal(descs.tx_wait[i])
                wait(descs.rx_wait[i])
                signal(descs.rx_wait[i])
            end
            sync_threads()
        else
            if threadIdx().x == 1
                @cuprintf("i=%ld, recv+send", i)
                threadfence()
                wait(descs.rx_wait[i])
                signal(descs.rx_wait[i])
                send(descs.tx[i])
                wait(descs.tx_wait[i])
            end
            sync_threads()
        end
    end
end

function exchange(comm, size, iter_count, validate)
    @info "implementation not yet done"

    # @cuda blocks=1, threads=16, stream=stream exchange_kernel(rank, descs, N)
end

function main(validate = true)
    MPI.Init()
    comm = MPI.COMM_WORLD

    if MPI.Comm_size(comm) != 2
        @error "This test requires exactly two processes"
        MPI.Abort(comm, -1)
        return
    end

    rank = MPI.Comm_rank(comm)

    # init gpus
    # init MP
    iter_count = ITER_COUNT_SMALL

    size = 1
    while size <= MAX_SIZE
        if size > 1024
            iter_count = ITER_COUNT_LARGE
        end

        exchange(comm, size, iter_count, validate)

        if rank == 0
            @info "SendRecv test passed validation with message size" size
        end

        size *= 2
    end

    # mp_finalize
    MPI.Barrier(comm)
    MPI.finalize()
end

ext: jl


text:
struct RunDump <: Output
  path_prefix::String
  RunDump(path_prefix::AbstractString) = new(path_prefix)
end
RunDump(path_prefix::AbstractString, ::Integer) = RunDump(path_prefix)

function pushtrajectory!(d::RunDump, trajectory_id::Integer, writelock::Base.AbstractLock, state::MocosSim.SimState, ::MocosSim.SimParams, callback::DetectionCallback)
  try lock(writelock)
    f = jldopen(d.path_prefix*"_$trajectory_id.jld2", "w", compress=true)
    try
      MocosSim.saveparams(f, state)
      saveparams(f, callback)
    finally
      close(f)
    end
  finally
    unlock(writelock)
  end
  nothing
end

max_stars_repo_name: Novirmor/MocosSimLauncher

max_issues_repo_name: Novirmor/MocosSimLauncher

max_forks_repo_name: Novirmor/MocosSimLauncher


text:
include("../plot/plot_plots.jl")
# include("../plot/plot_makie.jl")

using LinearAlgebra

"""
matrix = cross_matrix(vector)

2つのベクトルの外積cross(a, v)のvを行列化し、v*aで外積を計算できるようにする

# Argments
 - `vector`：外積するベクトル

# Return
 - `matrix`：行列化したベクトル
"""
function cross_matrix(vector)
	matrix = zeros(3,3)
	matrix[1,1] = 0
	matrix[1,2] = vector[3]
	matrix[1,3] = -vector[2]
	matrix[2,1] = -vector[3]
	matrix[2,2] = 0
	matrix[2,3] = vector[1]
	matrix[3,1] = vector[2]
	matrix[3,2] = -vector[1]
	matrix[3,3] = 0
	return matrix
end

"""
t_req  = lyapunov_torque(sat_tar, sat_att, kp, kr, ω)

リアプノフ関数から目標姿勢までに必要なトルクを求める

# Argments
 - `sat_tar`：目標姿勢クォータニオン
 - `sat_att`：現在の衛星姿勢クォータニオン
 - `kp`：ポイントゲイン
 - `kr`：レートゲイン
 - `ω`：位置ベクトルの角速度（ベクトル）

# Return
 - `t_req`：必要トルク
"""
function lyapunov_torque(sat_tar, sat_att, kp, kr, ω)
	q_e = zeros(3)
	q = zeros(3)
	q_t = zeros(3)
	q0 = sat_att.q0
	q[1] = sat_att.q1
	q[2] = sat_att.q2
	q[3] = sat_att.q3
	q0_t = sat_tar.q0
	q_t[1] = sat_tar.q1
	q_t[2] = sat_tar.q2
	q_t[3] = sat_tar.q3
	#println("q", q)
	#println("qt", q_t)
	q_e = - q0*q + q0*q_t - cross(q,q_t)
	# q0_e = q0 * q0_t + dot(q, q_t)
	#println("qe", q_e)
	return kp*q_e - kr*ω
end

"""
m = B_dot()

B-dot法により
姿勢制御する際に必要な磁気モーメントを求める

# Argments
 - `B`：地磁気ベクトル@ECEF
 - `ω`：位置ベクトルの角速度（ベクトル)

# Return
- `m`：必要磁気モーメント
"""
function B_dot(B, ω, ω_b)
	k = zeros(3)

	k = [40000, 40000, 40000]

	m = -1 * k .* cross(B, ω)

	return m
end

"""
m = cross_product(sat_tar, sat_att, kp, kr, ω, B)

cross_product法により
姿勢制御する際に必要な磁気モーメントを求める

# Argments
 - `sat_tar`：目標姿勢クォータニオン
 - `sat_att`：現在の衛星姿勢クォータニオン
 - `kp`：ポイントゲイン
 - `kr`：レートゲイン
 - `ω`：位置ベクトルの角速度（ベクトル）
 - `B`：地磁気ベクトル@ECEF

# Return
 - `m`：必要磁気モーメント
"""
function cross_product(sat_tar, sat_att, kp, kr, ω, B)
	m = zeros(3)
	#必要トルクを求める
	t_req = lyapunov_torque(sat_tar, sat_att, kp, kr, ω)
	#println("t_req", t_req)
	#必要磁気モーメントを求める
	m = -(cross(t_req, B/norm(B)))/(norm(B))
	return t_req, m
end

"""
m = pseudo_inverse(sat_tar, sat_att, kp, kr, ω, B)

擬似逆行列により
姿勢制御する際に必要な磁気モーメントを求める

# Argments
 - `sat_tar`：目標姿勢クォータニオン
 - `sat_att`：現在の衛星姿勢クォータニオン
 - `kp`：ポイントゲイン
 - `kr`：レートゲイン
 - `ω`：位置ベクトルの角速度（ベクトル）
 - `B`：地磁気ベクトル@ECEF
 - `ω`：位置ベクトルの角速度

# Return
- `m`：必要磁気モーメント
"""
function pseudo_inverse(sat_tar, sat_att, kp, kr, ω, B)
	m = zeros(3,3)
	B_mat = zeros(3,3)
	B_mat = cross_matrix(B)
	#B_matの擬似逆行列
	B_pse = pinv(B_mat)
	#必要トルクを求める
	t_req = lyapunov_torque(sat_tar, sat_att, kp, kr, ω)
	m = B_pse*t_req
	return t_req, m
end

"""
out_cur = attitude_control()

目標姿勢まで姿勢制御する

# Argments
 - `sat_att`：現在の衛星姿勢クォータニオン
 - `sat_tar`：目標姿勢クォータニオン
 - `ω`：位置ベクトルの角速度（ベクトル）
 - `I`：慣性テンソル

# Return
- `out_cur`：加える電流値
"""
function attitude_control(sat_att, sat_tar, B, ω, I)

end

"""
i = mm2current_theory(M)

出力したい磁気モーメントから磁気トルカに流す電流値（理論値）を求める

#Argments
M : 必要とされる磁気モーメント [A/m2]

#Return
i : 磁気トルカに流す電流 [A]
"""
function mm2current_theory(M)
	n = 400  # 巻き数　[巻]
	L = 80  # コア長さ [mm]
	D = 10  # コア直径 [mm]
	μ = 5000  #コア初期比透磁率

	S = π*(D/1000)^2/4  # 断面積 [m]
	p = L/D
	μeff = 1 / (1/μ + (log(p)-1)/(p^2))  # コア実効比透磁率

	i = M / (μeff*n*S)

	return i
end

"""
i = mm2current_measure(M)

出力したい磁気モーメントから磁気トルカに流す電流値（実験値）を求める

#Argments
M : 必要とされる磁気モーメント [A/m2]

#Return
i : 磁気トルカに流す電流 [A]
"""
function mm2current_measure(M)
	slope = [0.92, 0.92, 0.92]
	intercept = [0.001, 0.001, 0.001]  # ix = slope[1] * Mx + intercept[1]となる様な切片、傾き

	i = slope .* M + intercept

	return i
end

"""
crossproduct_adj = target_adjustment(targetqua)

"""

function crossproduct_adj(sat_tar, sat_att, kp, kr, ω, B, B2)
	m = zeros(3)
	n = 0
	#必要トルクを求める
	t_req = lyapunov_torque(sat_tar, sat_att, kp, kr, ω)

	t_req1 = [t_req[1], t_req[2], 0.]
	t_req2 = [0., t_req[2], t_req[3]]

	if norm(B2[1]) < norm(B2[3])
		Treq = t_req1
	else
		Treq = t_req2
		n = 1
	end

	M = -(cross(Treq, B/norm(B)))/(norm(B))

	return Treq, M, n
end

"""
t_req  = lyapunov_torque_addIntg(sat_tar, sat_att, kp, kr, ki, ω, errIntg)

リアプノフ関数から目標姿勢までに必要なトルクを求める

# Argments
 - `sat_tar`：目標姿勢クォータニオン
 - `sat_att`：現在の衛星姿勢クォータニオン
 - `kp`：ポイントゲイン
 - `kr`：レートゲイン
 - `ki`：積分ゲイン
 - `ω`：位置ベクトルの角速度（ベクトル）
 - `errIntg`：蓄積誤差量

# Return
 - `t_req`：必要トルク
"""
function lyapunov_torque_addIntg(sat_tar, sat_att, kp, kr, ki, ω, errIntg)
	q_e = zeros(3)
	q = zeros(3)
	q_t = zeros(3)
	q0 = sat_att.q0
	q[1] = sat_att.q1
	q[2] = sat_att.q2
	q[3] = sat_att.q3
	q0_t = sat_tar.q0
	q_t[1] = sat_tar.q1
	q_t[2] = sat_tar.q2
	q_t[3] = sat_tar.q3
	#println("q", q)
	#println("qt", q_t)
	q_e = - q0*q + q0*q_t - cross(q,q_t)
	q0_e = q0 * q0_t + dot(q, q_t)
	errIntg = errIntg + q0_e
	#println("qe", q_e)
	return kp*q_e - kr*ω .+ ki * errIntg, errIntg
end

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
using NearestNeighbors
using ReverseGeocode

"""
Prepares smaller dataset for testing.
Run this script after changing the testcases and
after downloading and processing the dataset.
"""
function prepare_testing_dataset(points;
        full_data="$(ReverseGeocode.DATA_DIR)/$(ReverseGeocode.GEO_FILE).csv")
    gc = Geocoder()

    # For each tested point, use 10 closest points
    idxs, dist = knn(gc.tree, points, 4)

    # Find the corresponding lines and save them
    lns = unique([1, vcat(vec(idxs)...).+1...]) # 1 for the header

    full_data_lines = readlines(full_data)
    open("./test/data/test_cities.csv", "w") do io
        for ln in lns
            line = full_data_lines[ln]
            write(io, "$line\n")
        end
    end
end

println("Preparing test data...")
points = [l.loc for l in test_locs] # test_locs is variable used in tests.
prepare_testing_dataset(points)
println("Done.")


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************

function plot_electric_dispatch(d::Dict; title="Electric Systems Dispatch", save_html=false)

    traces = PlotlyJS.GenericTrace[]

    layout = PlotlyJS.Layout(
        title_text = title,
        xaxis_title_text = "time step",
        yaxis_title_text = "kW"
    )

    eload = d["ElectricLoad"]["load_series_kw"]
    T = length(eload)

    push!(traces, PlotlyJS.scatter(
        name = "total load",
        x = 1:T,
        y = d["ElectricLoad"]["load_series_kw"],
        fill = "none",
        line = PlotlyJS.attr(
            width = 3
        ),
    ))

    push!(traces, PlotlyJS.scatter(
        name = "grid supply",
        x = 1:T,
        y = d["ElectricUtility"]["year_one_to_load_series_kw"],
        fill = "tozeroy",
        marker = PlotlyJS.attr(
            color="rgb(12,12,12)",
        ),
        line = PlotlyJS.attr(
            width = 0
        ),
    ))
    # invisible line for stacking
    push!(traces, PlotlyJS.scatter(
        name = "invisible",
        x = 1:T,
        y = d["ElectricUtility"]["year_one_to_load_series_kw"],
        fill = Nothing,
        line = PlotlyJS.attr(
            width = 0
        ),
        showlegend = false,
        hoverinfo = "skip",
    ))

    pv_to_load = zeros(T)
    if "PV" in keys(d)  # TODO multiple PVs
        pv_to_load = d["PV"]["year_one_to_load_series_kw"]
        push!(traces, PlotlyJS.scatter(
            name = "PV+grid supply",
            x = 1:T,
            y = pv_to_load .+ d["ElectricUtility"]["year_one_to_load_series_kw"],
            fill = "tonexty",
            marker = PlotlyJS.attr(
                color="rgb(255, 127, 14)",
            ),
            line = PlotlyJS.attr(
                width = 0
            ),
        ))
    end

    layout = PlotlyJS.Layout()

    if "ElectricStorage" in keys(d)
        # invisible line for stacking
        push!(traces, PlotlyJS.scatter(
            name = "invisible",
            x = 1:T,
            y = d["ElectricUtility"]["year_one_to_load_series_kw"] .+ pv_to_load,
            fill = Nothing,
            line = PlotlyJS.attr(
                width = 0
            ),
            showlegend = false,
            hoverinfo = "skip",
        ))
        push!(traces, PlotlyJS.scatter(
            name = "battery+PV+grid supply",
            x = 1:T,
            y = d["ElectricUtility"]["year_one_to_load_series_kw"] .+ pv_to_load .+ d["ElectricStorage"]["year_one_to_load_series_kw"],
            fill = "tonexty",
            marker = PlotlyJS.attr(
                color="rgb(44, 160, 44)",
            ),
            line = PlotlyJS.attr(
                width = 0
            ),
        ))
        push!(traces, PlotlyJS.scatter(
            name = "ElectricStorage SOC",
            x = 1:T,
            y = d["ElectricStorage"]["year_one_soc_series_pct"],
            yaxis="y2",
            marker = PlotlyJS.attr(
                color="rgb(100,100,100)",
            ),
        ))

        layout = PlotlyJS.Layout(
            title_text = title,
            xaxis_title_text = "time step",
            yaxis_title_text = "kW",
            yaxis2 = PlotlyJS.attr(
                title = "SOC",
                overlaying = "y",
                side = "right"
            )
        )
    end
    p = PlotlyJS.plot(traces, layout)

    if save_html
        PlotlyJS.savefig(p, replace(title, " " => "_") * ".html")
    end

    PlotlyJS.plot(traces, layout)  # will not produce plot in a loop

end


text:
using Flux
using FluxGAN
using FileIO
using MLDatasets

const n = parse(Int, ARGS[1])    # iterations
const m = parse(Int, ARGS[2])    # minibatch size
const skip = parse(Int, ARGS[3]) # verbose frequency

# dataset
train_tensor, train_labels = CIFAR10.traindata(Float32)

animal_classnames = [
    "bird",
    "cat",
    "deer",
    "dog",
    "frog",
    "horse"
]

classnames = CIFAR10.classnames()

# removing any images that aren't animals
animal_indices = findall(lab -> classnames[lab + 1] in animal_classnames, train_labels)
animals_tensor = train_tensor[:,:,:,animal_indices]

const img_size = size(animals_tensor)[1:end-1]
const img_num  = size(animals_tensor)[end]
const img_dim  = *(img_size...)

# centering pixels and flattening images
animals_tensor = reshape((@. 2f0 * animals_tensor - 1f0), (img_dim, img_num))

const d = 100

generator = Chain(
    Dense(d, 8000, relu),
    Dense(8000, 8000, tanh_fast),
    Dense(8000, img_dim, tanh_fast)
)

discriminator = Chain(
    Maxout(() -> Dense(img_dim, 1600), 5),
    Maxout(() -> Dense(1600, 1600), 5),
    Dense(1600, 1)
)

# learning rate
const η = 0.0001f0

# minibatch size

model = GAN(generator, discriminator;
    η_gen=η,
    η_dscr=η,
    latent_dim=d,
    minibatch=m,
    img_size=img_size
)

train!(model, animals_tensor, skip=skip, iterations=n)

# uncomment below to save model after training

println("\nsaving model...")
save("models/cifar10_goodfellow_MLP_n_$n.jld2", Dict("model" => model))
println("saved!\n")

# println("plotting generated images...")

# layout = (5, 5)

# output_dir = "images/CIFAR10"

# info = ["animals", "MLP", "n", n, "m", m]

# save_image_grid(model, output_dir;
#     layout=layout,
#     file_info=info,
#     img_res=150
# )

println("finished!\n")

max_stars_repo_name: aarontrowbridge/GAN

max_issues_repo_name: aarontrowbridge/GAN

max_forks_repo_name: aarontrowbridge/GAN


text:
# Drive train example
# ============================

using HybridSystems, MathematicalSystems, LazySets, Plots
import SX:is_halfspace, is_hyperplane, convert, parse_sxmath
import LazySets.HalfSpace

expr_p = parse_sxmath("0.2<=x & x<=0.3 & -0.1<=y & y<=0.1 & z==0 & x1==0 & x2==0 & x3==0")
#expr = [x,y,x1,x2,x3,z]
x0sets = []
vars = Basic[:x, :y, :x1, :x2, :x3, :z]
for expr_i in expr_p
    if is_halfspace(expr_i)
        push!(x0sets, convert(HalfSpace, expr_i, vars=vars))
    elseif is_hyperplane(expr_i)
        push!(x0sets, convert(Hyperplane, expr_i, vars=vars))
    end
end

AFFINE_SYSTEM = ConstrainedLinearControlContinuousSystem
HS = readsxmodel("models/SpaceEx/FilteredOscillator/filtered_oscillator_flattened.xml", ST=AFFINE_SYSTEM)

# calculate reachable states up to time T
prob = InitialValueProblem(HS, x0sets)
input_options = Options(:mode=>"reach");

problem_options = Options(:vars=>1:system_dimension, :T=>5.0, :δ=>0.05,
                          :plot_vars=>[1, 2], :max_jumps=>1, :verbosity=>1, :init_loc=>3);
options_input = merge(problem_options, input_options)
using Polyhedra
sol = solve(prob, options_input);

plot(sol, indices=1:2:length(sol.Xk))

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
mutable struct gps_satellite_info_t <: LCMType
    used_for_nav::Bool
    differential_correction_avail::Bool
    orbit_info_avail::Bool
    ephemeris::Bool
    healthy::Bool

    signal_quality::Int16  # 0-idle,(1,2)-searching,3-detected but unstable,4-lock on signal,(5,6)-code and carrier locked, 7-6+receiving @ 50bps
    carrier_to_noise::Int16
    azimuth::Float64
    elevation::Float64
end

@lcmtypesetup(gps_satellite_info_t)


text:
abstract type AbstractDataGenerator end
function launch(launchable::AbstractDataGenerator) error("unimplemented") end

abstract type  AbstractContext end
abstract type  AbstractSchedule end

"""
    AbstractRunnable

an abstract type interface for all decisional and utilitary processes (TSO, Market and assessment steps).
Such structure might need to overload the following functions :
- `compute_firmness(runnable::AbstractRunnable, ech, next_ech, TS, context::AbstractContext)`
- `run(runnable::AbstractRunnable, ech, firmness, TS, context::AbstractContext)`
- `affects_market_schedule(runnable::AbstractRunnable)`
- `update_market_schedule!(context::AbstractContext, ech, result, firmness, runnable::AbstractRunnable)`
- `affects_tso_schedule(runnable::AbstractRunnable)`
- `update_tso_schedule!(context::AbstractContext, ech, result, firmness, runnable::AbstractRunnable)`
- `affects_tso_actions(runnable::AbstractRunnable)`
- `update_tso_actions!(context::AbstractContext, ech, result, firmness, runnable::AbstractRunnable)`
"""
abstract type  AbstractRunnable end
"""
Determines the firmnes of commitment and powerlevel decisions for each timestep
"""
function compute_firmness(runnable::AbstractRunnable, ech, next_ech, TS, context::AbstractContext) error("unimplemented") end
"""
Launches the runnable :
    get relevant data from context
    build the optimization model
    solve the model
    return a result object to be exploited by the updating methods
"""
function run(runnable::AbstractRunnable, ech, firmness, TS, context::AbstractContext) error("unimplemented") end
"""
If returns True, will launch the market_schedule related updates during sequence execution
"""
function affects_market_schedule(runnable::AbstractRunnable) return false end
"""
Updates the market schedule using the result of run()
"""
function update_market_schedule!(context::AbstractContext, ech, result, firmness, runnable::AbstractRunnable) error("unimplemented") end
"""
If returns True, will launch the tso_schedule related updates during sequence execution
"""
function affects_tso_schedule(runnable::AbstractRunnable) return false end
"""
Updates the TSO schedule using the result of run()
"""
function update_tso_schedule!(context::AbstractContext, ech, result, firmness, runnable::AbstractRunnable) error("unimplemented") end
"""
If returns True, will launch the tso actions related updates during sequence execution
"""
function affects_tso_actions(runnable::AbstractRunnable) return false end
"""
Updates the TSO actions using the result of run()
"""
function update_tso_actions!(context::AbstractContext, ech, result, firmness, runnable::AbstractRunnable) end
#FIXME are these compatible with the Assessment step ?

abstract type  AbstractTSO <: AbstractRunnable  end
function affects_tso_schedule(runnable::AbstractTSO) return true end
function affects_tso_actions(runnable::AbstractTSO) return true end

abstract type  AbstractMarket <: AbstractRunnable  end
function affects_market_schedule(runnable::AbstractMarket) return true end

abstract type DeciderType end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
module base_things
    mutable struct 변수
        input_num::Float64
        Weight::Float64
        real_num::Float64
        abs_error::Float64
        minus_weight ::Float64
    end

    function abs_error_cal(var)::Float64
        var.abs_error = abs(var.input_num - var.real_num)
        return var.abs_error
    end

    function forward(var,sum_of_W)
        percent_of_W = var.Weight/sum_of_W
        out = var.input_num*percent_of_W
        out
    end

    function graph(add,max,avg)
        if  add > avg
            return -(1.5+(add*avg))
        else
            y = cos((10.0*add)/(π*max)) - avg*add
        return y
        end
    end

    function backward(var,new_w)
        new_v = var.Weight + var.minus_weight + new_w
        if new_v < 0
            var.minus_weight = new_v
            var.Weight = 0
            return var.Weight , var.minus_weight
        else
            return new_v , 0
        end
    end

end


text:
# ---
# title: 1523. Count Odd Numbers in an Interval Range
# id: problem1523
# author: zhwang
# date: 2022-03-18
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/count-odd-numbers-in-an-interval-range/description/>
# hidden: true
# ---
#
# Given two non-negative integers `low` and `high`. Return the _count of odd
# numbers between_`low` _and_` high` _  (inclusive)_.
#
#
#
# **Example 1:**
#
#
#
#     Input: low = 3, high = 7
#     Output: 3
#     Explanation: The odd numbers between 3 and 7 are [3,5,7].
#
# **Example 2:**
#
#
#
#     Input: low = 8, high = 10
#     Output: 1
#     Explanation: The odd numbers between 8 and 10 are [9].
#
#
#
# **Constraints:**
#
#   * `0 <= low <= high <= 10^9`
#
#
## @lc code=start
using LeetCode

count_odd_numbers(low::Int, high::Int) = (high - low + 1 + (low & 1)) >> 1

## @lc code=end


text:
using WindMixing

FILE_NAME = "NDE_training_mpp_10sim_windcooling_windheating_diffusivity_1e-1_Ri_1e-1_weights_divide1f5_gradient_smallNN_scale_5e-3_rate_1e-4"

test_files = [
    "wind_-1e-3_heating_-4e-8",
    "wind_-1e-3_heating_-1e-8",
    "wind_-5e-4_heating_-3e-8",
    "wind_-2e-4_heating_-5e-8",
    "wind_-2e-4_heating_-1e-8",
    "wind_-1e-3_cooling_4e-8",
    "wind_-2e-4_cooling_1e-8",
    "wind_-1e-3_cooling_2e-8",
    "wind_-2e-4_cooling_5e-8",
    "wind_-5e-4_cooling_3e-8",
]

OCEANANIGANS_OUTPUT_DIR = joinpath("D:\\Output_o", FILE_NAME)
animate_training_results_oceananigans(test_files, FILE_NAME, OCEANANIGANS_OUTPUT_DIR, timestep=60)


text:
using Markdown
#data = read("learn_markdown.md", String)
d = Markdown.parse_file("learn_markdown.md")
h = html(d)
open("learn_markdown.html", "w") do io
    println(io, h)
end


text:
struct SSHBinner
    thetas::Vector{Float64}
    phi_divisions::Vector{Int64}

    # bin_many helpers
    zs::Vector{Float64}
    phi_N_over_2pi::Vector{Float64}
    cumsum::Vector{Int64}
    N_thetas::Float64

    bins::Vector{Int64}
end

################################################################################
### Constructor
################################################################################

function SSHBinner(thetas::Vector{Float64}, phi_divisions::Vector{Int64})
    # Precalculate N/2pi for bin_many!
    # use slightly more than 2pi to guarantee phi * phi_N_over_pi ∈ [0, 1)
    PI_PLUS = 2.0pi + 1e-14
    phi_N_over_2pi = phi_divisions ./ PI_PLUS

    # Precalculate z(θ) for bin_many!
    # move slightly to negative values to guarantee that for an input value
    # value[3] >= minimum(zs)
    zs = cos.(thetas) #.- 1e-14 # 1...-1
    zs[end] -= 1e-14

    # Also precalculate cumulative number of bins for each theta (from 0..pi)
    _cumsum = cumsum(phi_divisions)

    # to avoid conversions to Float64
    N_thetas = Float64(length(thetas)-1)

    N_bins = sum(phi_divisions)
    bins = zeros(Int64, N_bins)
    SSHBinner(thetas, phi_divisions, zs, phi_N_over_2pi, _cumsum, N_thetas, bins)
end

"""
    SSHBinner(N_bins)

Generate a Sphere Surface Histogram Binner with approximately N_bins.
"""
SSHBinner(N_bins::Real; kwargs...) = SSHBinner(round(Int64, N_bins); kwargs...)
function SSHBinner(N_bins::Int64; method=partition_sphere2)
    thetas, phi_divisions = partition_sphere_optim(4pi/N_bins, method=method)
    SSHBinner(thetas, phi_divisions)
end

"""
    SSHBinner(bins)

Generate a Sphere Surface Histogram Binner using the `bins` of another
SSHBinner. This effectively creates a copy of the other source `SSHBinner`.
"""
function SSHBinner(bins::Vector; method=partition_sphere2)
    N_bins = length(bins)
    thetas, phi_divisions = method(4pi/N_bins)

    PI_PLUS = 2.0pi + 1e-14
    phi_N_over_2pi = phi_divisions ./ PI_PLUS
    zs = cos.(thetas)
    zs[end] -= 1e-14
    _cumsum = cumsum(phi_divisions)
    N_thetas = Float64(length(thetas)-1)

    SSHBinner(thetas, phi_divisions, zs, phi_N_over_2pi, _cumsum, N_thetas, bins)
end

"""
    partition_sphere2(dA)

Partitions a sphere into bins of size `dA` where number of divisions in ϕ
(0..2π) is restricted to a power of 2.

Note that for a random `dA` this will not produce equally sized bins.
(Specifically bins at z ≈ 0 will be wrong.) The algorithm relies on
`partition_sphere_optim` to optimize `dA` to a value fitting on the sphere
surface which happens during the creation of an `SSHBinner`.
"""
function partition_sphere2(dA)
    # Use symmetry - only 0..pi/2 needed (± some delta if center of dA @ pi/2)

    # Starting at some point theta1 with width dphi,
    # computes theta2 for an area dA
    next_theta(dA, theta1, dphi) = acos(max(min(cos(theta1) - dA / dphi, 1.0), -1.0))

    # First area
    # This will now be our goal theta-width
    dtheta_goal = next_theta(dA, 0.0, 2pi)

    # Devide 2pi by phi_div, such that phi = 2pi/n with n integer
    phi_divs = Int64[1]
    thetas = Float64[0.0, dtheta_goal]

    # Currently radius of circle, should be sidelength of square
    # pi r^2 = x^2 => sqrt(pi)r = x
    dtheta_goal *= sqrt(pi)

    while true
        phi_div = last(phi_divs)
        theta1 = last(thetas)

        # Find best dtheta
        dtheta = next_theta(dA, theta1, 2pi / phi_div) - theta1
        best_dtheta = dtheta
        while true
            phi_div *= 2
            dtheta = next_theta(dA, theta1, 2pi / phi_div) - theta1
            if abs(dtheta - dtheta_goal) < abs(best_dtheta - dtheta_goal)
                # New dtheta better? Keep it and keep trying
                best_dtheta = dtheta
            else
                # new dtheta worse? undo last phi_div increment and stop
                phi_div = div(phi_div, 2)
                break
            end
        end

        if (theta1 + best_dtheta / 2.0) > pi/2
            break
        else
            push!(thetas, theta1 + best_dtheta)
            push!(phi_divs, phi_div)
        end
    end

    dtheta = thetas[end] - thetas[end-1]
    theta_fourth = thetas[end-1] + 3dtheta/4
    if theta_fourth > pi/2
        # center should be flush
        thetas[end] = pi - thetas[end-1]
        for i in length(thetas)-2:-1:1
            push!(thetas, pi - thetas[i])
            push!(phi_divs, phi_divs[i])
        end
    else
        # last(thetas) should be flush
        thetas[end] = pi/2
        for i in length(thetas)-1:-1:1
            push!(thetas, pi - thetas[i])
            push!(phi_divs, phi_divs[i])
        end
    end

    return thetas, phi_divs
end

"""
    partition_sphere1(dA)

Partitions a sphere into bins of size `dA` where number of divisions in ϕ
(0..2π) is not restricted, i.e. it can be any integer.

Note that for a random `dA` this will not produce equally sized bins.
(Specifically bins at z ≈ 0 will be wrong.) The algorithm relies on
`partition_sphere_optim` to optimize `dA` to a value fitting on the sphere
surface which happens during the creation of an `SSHBinner`.
"""
function partition_sphere1(dA,)
    next_theta(dA, theta1, dphi) = acos(clamp(cos(theta1) - dA / dphi, -1.0, 1.0))

    # First area
    # This will now be our goal theta-width
    dtheta_goal = next_theta(dA, 0.0, 2pi)

    # Devide 2pi by phi_div, such that phi = 2pi/n with n integer
    phi_divs = Int64[1]
    thetas = Float64[0.0, dtheta_goal]

    # Currently radius of circle, should be sidelength of square
    # pi r^2 = x^2 => sqrt(pi)r = x
    dtheta_goal *= sqrt(pi)

    while true
        phi_div = last(phi_divs)
        theta1 = last(thetas)

        # Find best dtheta
        dtheta = next_theta(dA, theta1, 2pi / phi_div) - theta1
        best_dtheta = dtheta
        while true
            phi_div += 1
            dtheta = next_theta(dA, theta1, 2pi / phi_div) - theta1
            if abs(dtheta - dtheta_goal) < abs(best_dtheta - dtheta_goal)
                # New dtheta better? Keep it and keep trying
                best_dtheta = dtheta
            else
                # new dtheta worse? undo last phi_div increment and stop
                phi_div -= 1
                break
            end
        end

        if (theta1 + best_dtheta / 2.0) > pi/2
            break
        else
            push!(thetas, theta1 + best_dtheta)
            push!(phi_divs, phi_div)
        end
    end

    dtheta = thetas[end] - thetas[end-1]
    theta_fourth = thetas[end-1] + 3dtheta/4
    if theta_fourth > pi/2
        # center should be flush
        thetas[end] = pi - thetas[end-1]
        for i in length(thetas)-2:-1:1
            push!(thetas, pi - thetas[i])
            push!(phi_divs, phi_divs[i])
        end
    else
        # last(thetas) should be flush
        thetas[end] = pi/2
        for i in length(thetas)-1:-1:1
            push!(thetas, pi - thetas[i])
            push!(phi_divs, phi_divs[i])
        end
    end

    return thetas, phi_divs
end

"""
    partition_sphere_optim(dA[; max_iter = 10, method = partition_sphere2])

Optimizes the sphere partitioning by modifying dA (~ the number of bins). After
a few iterations, dA becomes consistent for all bins.
"""
function partition_sphere_optim(dA; max_iter = 10, method=partition_sphere2)
    # NOTE this converges really fast

    new_dA = dA
    N = round(Int64, 4pi/dA)
    _iterations = 0
    for i in 1:max_iter
        _iterations = i
        thetas, phi_divs = method(new_dA)
        next_N = sum(phi_divs)
        if N == next_N
            break
        else
            N = next_N
            new_dA = 4pi / N
        end
    end

    method(new_dA)
end

################################################################################
### binning (push, append)
################################################################################

@fastmath @inline function fast_theta_index_approximation(x, l::Float64)
    # trunc faster than floor
    # Brackets cause llvm code to be `fmul` rather than `afoldl`, much
    # faster (~100x according to BenchmarkTools)
    # equivalent to
    # trunc(Int64, l*(0.5 - (7x + x^5 + x^9 + x^17 + x^31) / 22)
    # which is ≈ trunc(Int64, l*acos(x))
    trunc(Int64,
        l * (0.5 - x*(
            0.3181818181818182 + ((((x*x)*x)*x)*0.045454545454545456)*(
                1.0 + (((x*x)*x)*x) * (
                    1.0 + (((((((x*x)*x)*x)*x)*x)*x)*x)*(
                        1.0 + (((((((((((((x*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)
                    )
                )
            )
        ))
    ) + 1
end

# Tested for 10, 1000, 10_000_000 bins, each with 1_000_000 random spins pushed
@inline function fast_theta_index_search(B::SSHBinner, value)
    index = fast_theta_index_approximation(value, B.N_thetas)
    if value > B.zs[index]
        while true
            index -= 1
            value <= B.zs[index] && (return index)
        end
    elseif B.zs[index] >= value
        while true
            index += 1
            B.zs[index] < value && (return index-1)
        end
    end
    return index
end

"""
    push!(B::SSHBinner, value)

Bins a single value (three dimensional unit vector).
"""
function Base.push!(B::SSHBinner, value)
    theta_index = fast_theta_index_search(B, value[3])
    phi_index = if B.phi_divisions[theta_index] == 1
        1
    else
        # Normalize xy-component
        # y >= 0 for 0..pi, else y < 0
        # acos(x) well defined for 0..pi
        R = sqrt(value[1]^2 + value[2]^2)
        phi = if value[2] >= 0.0
            acos(value[1]/R)
        else
            2.0pi - acos(value[1]/R)
        end

        trunc(Int64, phi * B.phi_N_over_2pi[theta_index]) + 1
    end

    # Calculate position in bins
    l = theta_index > 1 ? B.cumsum[theta_index-1] : 0
    B.bins[l+phi_index] += 1
    nothing
end

"""
    append!(B::SSHBinner, values)

Bins an array of values (three dimensional unit vectors).
"""
function Base.append!(B::SSHBinner, values)
    for v in values
        push!(B, v)
    end
    nothing
end

"""
    unsafe_append!(B::SSHBinner, values)

Bins an arra of values (three dimensional unit vector) without boundschecks.
"""
function unsafe_append!(B::SSHBinner, values)
    for v in values
        @inbounds push!(B, v)
    end
    nothing
end


text:
#####
##### More Advanced Policy
#####

# TODO: Mutable Binary Heaps from DataStructurtes.jl aren't really that great.
# Consider using a different data structure that doesn't suffer the memory blowup of
# the DataStructures.jl implementation.

# Type parameter: number of bins
mutable struct OptaneTracker{N}
    count::Int

    # The minimum size of each bin.
    bins::NTuple{N,Int}

    # Objects that live in the local heap that haven't been marked as easily evictable.
    regular_objects::NTuple{N,LRU{Block}}
    evictable_objects::NTuple{N,LRU{Block}}

    # For critical regions with concurrent allocations
    movement_enabled::Bool
end

function OptaneTracker(bins::NTuple{N,Int}) where {N}
    count = 0
    regular_objects = ntuple(_ -> LRU{Block}(), Val(N))
    evictable_objects = ntuple(_ -> LRU{Block}(), Val(N))

    return OptaneTracker{N}(
        count,
        bins,
        regular_objects,
        evictable_objects,
        true,
    )
end

# Almost exactly like `Base.findlast`, but instead of returning `N+1` if the predicate
# is never true, return `N`.
@inline function findbin(bins::NTuple{N,Int}, x; inbounds = true) where {N}
    for i in Base.OneTo(N)
        x <= @inbounds(bins[i]) && return i
    end
    return inbounds ? N : (N + 1)
end

function increment!(policy::OptaneTracker)
    x = policy.count
    policy.count += 1
    return x
end

#####
##### Policy API
#####

function Base.push!(
    policy::OptaneTracker,
    block::Block,
    pool = getpool(block);
    len = length(block),
)
    # Don't explicitly track remote blocks.
    pool == Remote && return block
    bin = findbin(policy.bins, len; inbounds = true)
    lru = policy.regular_objects[bin]
    push!(lru, block, increment!(policy))
    return block
end

function Base.delete!(policy::OptaneTracker, block::Block; len = length(block))
    getpool(block) == Remote && return true

    # TODO: Metadata in block so we don't have to do two searches.
    bin = findbin(policy.bins, len; inbounds = true)
    lru = policy.regular_objects[bin]
    if in(block, lru)
        delete!(lru, block)
        return true
    end

    lru = policy.evictable_objects[bin]
    if in(block, lru)
        delete!(lru, block)
        return true
    end
    return false
end

# TODO
update!(policy::OptaneTracker, _) = nothing

function softevict!(policy::OptaneTracker, manager, block)
    # Is this block in any of the live heaps.
    bin = findbin(policy.bins, length(block); inbounds = true)
    lru = policy.regular_objects[bin]
    if in(block, lru)
        delete!(lru, block)
        push!(policy.evictable_objects[bin], block, increment!(policy))
    end
    return nothing
end

#####
##### Policy <---> CacheManager bridge
#####

### policy_new_alloc
function policy_new_alloc(
    policy::OptaneTracker{N},
    manager,
    bytes,
    id,
    priority::AllocationPriority,
) where {N}
    # Can we try to allocate locally?
    if priority != ForceRemote
        @return_if_exists ptr = _try_alloc_local(policy, manager, bytes, id, priority)
    end
    @return_if_exists ptr = unsafe_alloc_direct(RemotePool(), manager, bytes, id)

    # Failed so far - try again if we have enough memory locally to defrag.
    allocated, total = getstate(getheap(manager, LocalPool()))
    if (total - allocated) > bytes
        defrag!(manager, policy)
    end
    @return_if_exists ptr = _try_alloc_local(policy, manager, bytes, id, priority)

    @show bytes, manager
    error("Ran out of memory!")
end

function defrag!(manager, policy::OptaneTracker)
    cb = () -> unsafe_cleanup!(manager)
    defrag!(getheap(manager, LocalPool()); queued_callback = cb) do _, newblock, oldblock
        # First, we need to check if this is even the primary block for the corresponding
        # object.
        #
        # N.B.: No need to handle the other case since the "defrag!" routine for the heap
        # will preserve siblings.
        #
        # This may change ...
        primary = getprimary(manager, oldblock)
        if primary === oldblock
            _ = unsafe_setprimary!(manager, oldblock, newblock; unsafe = true)
            @check delete!(policy, oldblock; len = length(oldblock))
            push!(policy, newblock, pool; len = length(oldblock))
        end
    end
end

function prefetch!(block::Block, policy::OptaneTracker, manager; readonly = false)
    policy.movement_enabled || return nothing
    block.size > sizeof(getheap(manager, LocalPool())) && return nothing

    if getpool(block) == Local
        return nothing
    end

    @check getsibling(block) === nothing
    bytes = length(block)
    if !canalloc(getheap(manager, LocalPool()), bytes)
        _eviction!(policy, manager, bytes)
    end

    if isqueued(block)
        return nothing
    end
    #safeprint("Prefetching block $(getid(block)) with size $(block.size)"; force = true)

    # Allocate and move.
    # Don't free the old block since we're functioning as a cache.
    ptr = unsafe_alloc_direct(LocalPool(), manager, length(block), getid(block))
    ptr === nothing && return nothing
    newblock = unsafe_block(ptr)
    copyto!(newblock, block, manager)
    result = unsafe_setprimary!(manager, block, newblock)
    if result === nothing
        unsafe_free_direct(manager, newblock)
    else
        # Need to update the local policy tracking to know that this is now local.
        push!(policy, newblock, Local)
        link!(newblock, block)

        # If this is not a read-only prefetch, then we need to be conservative and mark the
        # prefetched block as dirty.
        setdirty!(newblock, !readonly)
    end

    return nothing
end

function evict!(A, policy::OptaneTracker, manager)
    policy.movement_enabled || return nothing
    if getpool(metadata(A)) != Remote
        eviction_callback(manager, policy, metadata(A))
    end
    return nothing
end

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

# TODO: Time since last GC?
function _try_alloc_local(policy::OptaneTracker, manager, bytes, id, priority)
    # allocated, total = getstate(getheap(manager, LocalPool()))
    # if allocated / total >= 0.90
    #     # Trigger full GC and try to get pending finalizers to run.
    #     GC.gc(false)
    # end

    # If allocation is successful, good!
    @return_if_exists ptr = unsafe_alloc_direct(LocalPool(), manager, bytes, id)

    # If we're here, and we really need local, then we have to start evicting.
    if policy.movement_enabled && priority == ForceLocal
        _eviction!(policy, manager, bytes)
        ptr = unsafe_alloc_direct(LocalPool(), manager, bytes, id)
    end
    return ptr
end

function _cleanup!(manager, id)
    cleaned = unsafe_cleanup!(manager, id)
    @check cleaned
    return nothing
end

function eviction_callback(manager, policy, block::Block)
    # Copy back to sibling if one already exists.
    id = getid(block)
    if isqueued(block)
        _cleanup!(manager, id)
        return false
    end

    sibling = getsibling(block)
    if sibling !== nothing
        isdirty(block) && copyto!(sibling, block, manager)
        result = unsafe_setprimary!(manager, block, sibling)
        if result === nothing
            _cleanup!(manager, id)
        else
            unlink!(block, sibling)
            @check delete!(policy, block)
            unsafe_free_direct(manager, block)
        end
        return false
    end

    # Rare case where a pointer may be null.
    # In this case, we need to trigger the abort mechanism in the eviciton process so we
    # don't leave the heap in an undefined state.
    ptr = unsafe_alloc_direct(RemotePool(), manager, length(block), getid(block))
    ptr === nothing && return true
    #safeprint("Evicting block $(getid(block)) with size $(block.size)"; force = true)

    newblock = unsafe_block(ptr)
    copyto!(newblock, block, manager)
    result = unsafe_setprimary!(manager, block, newblock)

    # Setting primary failed because the current block is queued for freeing.
    # Roll back our change.
    #
    # Otherwise, setting the primary was successful, so we free the old block.
    if result === nothing
        unsafe_free_direct(manager, newblock)
        _cleanup!(manager, id)
    else
        # For now, make sure that something was actually deleted because it should ALWAYS
        # be tracked by the policy.
        @check delete!(policy, block)
        unsafe_free_direct(manager, block)
    end
    return false
end

function _eviction!(policy::OptaneTracker{N}, manager, bytes) where {N}
    # Get a bin number for this allocation.
    # If we have an easily evictable object of the correct size, then use that.
    # Pass the "inbounds = true" flag because if we don't have a block big
    # enough, then we want to evict from the largest bin.
    bin = findbin(policy.bins, bytes; inbounds = true)

    # TODO: Cleanup
    # canabort = manager.abort_callback
    # canabort.bytes = bytes
    cb = x -> eviction_callback(manager, policy, x)
    localheap = getheap(manager, LocalPool())

    # Check this bin and all higher bins.
    for i in bin:N
        lru = policy.evictable_objects[i]
        if !isempty(lru)
            block = first(lru)
            @check getpool(block) == Local
            evictfrom!(localheap, block, bytes; cb)
            return nothing
        end
    end

    # Try evicting not from the easily evictable trackers.
    for i in bin:N
        lru = policy.regular_objects[i]
        if !isempty(lru)
            block = first(lru)
            @check getpool(block) == Local
            evictfrom!(localheap, block, bytes; cb)
            return nothing
        end
    end
    return nothing
end


text:
using ModelingToolkit, Test

@testset "Parsing Test" begin include("variable_parsing.jl") end
@testset "Differentiation Test" begin include("derivatives.jl") end
@testset "Simplify Test" begin include("simplify.jl") end
@testset "Direct Usage Test" begin include("direct.jl") end
@testset "System Construction Test" begin include("system_construction.jl") end
@testset "Domain Test" begin include("domains.jl") end
@testset "Constraints Test" begin include("constraints.jl") end
@testset "PDE Construction Test" begin include("pde.jl") end
@testset "Distributed Test" begin include("distributed.jl") end
#@testset "Latexify recipes Test" begin include("latexify.jl") end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, UUIDs, Random

u1 = uuid1()
u4 = uuid4()
u5 = uuid5(u1, "julia")
@test uuid_version(u1) == 1
@test uuid_version(u4) == 4
@test uuid_version(u5) == 5
@test u1 == UUID(string(u1)) == UUID(GenericString(string(u1)))
@test u4 == UUID(string(u4)) == UUID(GenericString(string(u4)))
@test u5 == UUID(string(u5)) == UUID(GenericString(string(u5)))
@test u1 == UUID(UInt128(u1))
@test u4 == UUID(UInt128(u4))
@test u5 == UUID(UInt128(u5))
@test uuid4(MersenneTwister(0)) == uuid4(MersenneTwister(0))
@test_throws ArgumentError UUID("550e8400e29b-41d4-a716-446655440000")
@test_throws ArgumentError UUID("550e8400e29b-41d4-a716-44665544000098")
@test_throws ArgumentError UUID("z50e8400-e29b-41d4-a716-446655440000")

# results similar to Python builtin uuid
# To reproduce the sequence
#=
import uuid
uuids = [uuid.UUID("22b4a8a1-e548-4eeb-9270-60426d66a48e")]
for _ in range(5):
    uuids.append(uuid.uuid5(uuids[-1], "julia"))
=#

const following_uuids = [
    UUID("22b4a8a1-e548-4eeb-9270-60426d66a48e"),
    UUID("30ea6cfd-c270-569f-b4cb-795dead63686"),
    UUID("31099374-e3a0-5fde-9482-791c639bf29b"),
    UUID("6b34b357-a348-53aa-8c71-fb9b06c3a51e"),
    UUID("fdbd7d4d-c462-59cc-ae6a-0c3b010240e2"),
    UUID("d8cc6298-75d5-57e0-996c-279259ab365c"),
]

for (idx, init_uuid) in enumerate(following_uuids[1:end-1])
    next_id = uuid5(init_uuid, "julia")
    @test next_id == following_uuids[idx+1]
end

# Python-generated UUID following each of the standard namespaces
const standard_namespace_uuids = [
    (UUIDs.namespace_dns,  UUID("00ca23ad-40ef-500c-a910-157de3950d07")),
    (UUIDs.namespace_oid,  UUID("b7bf72b0-fb4e-538b-952a-3be296f07f6d")),
    (UUIDs.namespace_url,  UUID("997cd5be-4705-5439-9fe6-d77b18d612e5")),
    (UUIDs.namespace_x500, UUID("993c6684-82e7-5cdb-bd46-9bff0362e6a9")),
]

for (init_uuid, next_uuid) in standard_namespace_uuids
    result = uuid5(init_uuid, "julia")
    @test next_uuid == result
end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
module TheoryOfCISS

using InlineExports

# using AtomicSimulationEnvironment
# @export const ase = AtomicSimulationEnvironment

using BSON: @save, @load
using Dates: day, month, now, today, year
using LinearAlgebra
using Molecules
using Molecules: Molecule, Bond,
    Hydrogen, Carbon, Nitrogen, Oxygen,
    angularmomentum, hamiltonian, makeaxes, mass, nB, overlap
using ProgressMeter
using Random
using SlaterKoster: SlaterKoster, loaddir
using SparseArrays
using StaticArrays
using Statistics

@export const σ0 = [1.0 0.0; 0.0 1.0]
@export const σ = (
    x = [0.0 1.0; 1.0 0.0],
    y = [0.0 -1.0im; 1.0im 0.0],
    z = [1.0 0.0; 0.0 -1.0]
)
@export const ⊗ = kron

vec(x, y, z) = SVector(x, y, z)

abstract type Simulation end

include("globals.jl")

include("mol.jl")
include("continuous.jl")

include("lite.jl")
include("lite_datagen.jl")
include("plots_recipes.jl")

include("../extra/scratch.jl")

# include("basic.jl")

end # module


text:
#!/usr/bin/env julia
# Copyright 2017 Christian Henning
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

workspace() # In case, this code block is copied to a notebook.

using PyPlot

include("network.jl")

# Specify dataset.
function linear_labels(p_x::Array{Float64, 2}, out_size::Int)
    w = ones(size(p_x, 1))
    w /= norm(w)

    #b = norm(ones(size(p_x, 1)) * 0.5)
    b = 0

    y = zeros(Bool, out_size, size(p_x, 2))
    for i = 1:size(p_x, 2)
        if (dot(w, p_x[:,i]) - b) > 0
            y[:,i] = 1
        end
    end

    y
end

N_train = 100 # Size of training set.
N_test = 100 # Size of test set.
n_in = 10 # Number of input units.
n_out = 20 # Number of output units.

train_in = rand(n_in, N_train) .- 0.5
test_in = rand(n_in, N_test) .- 0.5

# Uncomment, to work with static sample.
#sample = rand(n_in)
#for i = 1:n_in
#    train_in[i,:] = sample[i]
#    test_in[i,:] = sample[i]
#end

train_out = linear_labels(train_in, n_out)
test_out = linear_labels(test_in, n_out)

in_size = size(train_in)[1]
out_size = size(train_out)[1]

# Define hyperparameters.
hparams = Hyperparameters()
hparams.λ = 1
hparams.γ = 0.99
hparams.αₛ = 0.5
hparams.αₚ = 0.5
hparams.rᵦ = 0.0
hparams.η = 0.9
hparams.num_epochs = 1000
hparams.batch_size = 1
hparams.sample_size = 100
hparams.clip_gradients = false
hparams.clipping_value = 0.1

# Weight initializer.
initializer(dims...) = 1e-1 * (2*rand(dims) - 1)

# Build network.
network = Network(hparams = hparams)

loss_params = LossParameters()
loss_params.sparsity = 0.5

layer0 = input_layer("in", in_size)
add_layer(network, layer0)

layer1 = fully_connected_layer("fc1", layer0, 4, loss_params=loss_params,
                               initializer=initializer)
add_regularizer(layer1, sparsity_regularizer)
add_regularizer(layer1, participation_regularizer)
add_layer(network, layer1)

layer2 = fully_connected_layer("out", layer1, out_size, loss_params=loss_params,
                               initializer=initializer)
add_layer(network, layer2)

# Train Network.
T, E = train_network(network, train_in, train_out = train_out,
                    test_in = test_in,
                    test_out = test_out,
                    eval_interval = 10,
                    eval_accuracy = false)

# Plot training development.
title("Development of training loss")
xlabel("Epoch")
ylabel("Train Loss")
plot(1:length(T), T, label = "Train Loss")
legend()
show()

title("Development of global test loss (without regularizers) during training")
xlabel("Epoch")
ylabel("Test Eval")
plot(10:10:length(T), E, label="Test Evaluation")
axhline(y=0.0, color="k", linestyle="--")
legend()
ax = gca()
ax[:set_xlim]([0, length(T)])
show()


text:
# Top level of the Push interpreter.
run(p::String) = run(p, default_configuration)
run(p::String, cfg::String) = run(p, load_configuration(cfg))
run(p::String, cfg::Configuration) = run(Parser.to_push(p), configure(cfg))
function run(p::Any, s::State)
  s.parameters.top_level_push_code && push!(s.code, p)
  push!(s.exec, p)
  execute(s)
  s.parameters.top_level_pop_code && pop!(s.code)
  return s
end

# No termination limit for now.
execute(s::State) = while !isempty(s.exec)
  if s.flag_quote_code
    s.flag_quote_code = false
    push!(s.code, pop!(s.exec))
  else
    execute(s, pop!(s.exec))
  end
end

# Could be made faster...
execute(s::State, v::Vector{Any}) = for elem in reverse(v)
  push!(s.exec, elem)
end

execute(s::State, v::Bool) =
  push!(s.boolean, v)

execute(s::State, v::Int32) =
  push!(s.integer, v)

execute(s::State, f::Float32) =
  push!(s.float, f)

function execute(s::State, v::Symbol)
  # If name quoting is enabled, push this name onto the NAME stack and disable quoting.
  if s.flag_quote_name
    s.flag_quote_name = false
    push!(s.name, v)

  # Check if this name refers to a "built-in" instruction.
  elseif haskey(s.instructions, v)
    s.instructions[v](s)

  # Check if the name refers to a stored definition.
  elseif haskey(s.definitions, v)
    push!(s.exec, s.definitions[v])

  # If neither of the above, treat the symbol as a literal name
  # and add it to the NAME stack.
  else
    push!(s.name, v)
  end
end


text:
export Domain,IntervalDomain,PeriodicDomain,tocanonical,fromcanonical,fromcanonicalD,∂
export chebyshevpoints,fourierpoints,isambiguous,arclength

# T is the numeric type used to represent the domain
# d is the dimension
abstract Domain{T,d}
typealias UnivariateDomain{T} Domain{T,1}
typealias BivariateDomain{T} Domain{T,2}

Base.eltype{T}(::Domain{T}) = T
Base.eltype{T,d}(::Type{Domain{T,d}}) = T
Base.isreal{T<:Real}(::Domain{T}) = true
Base.isreal{T}(::Domain{T}) = false
dimension{T,d}(::Domain{T,d}) = d
dimension{T,d}(::Type{Domain{T,d}}) = d
dimension{DT<:Domain}(::Type{DT}) = dimension(supertype(DT))

Base.length(d::Domain) = 1

# prectype gives the precision, including for Vec
prectype(d::Domain) = eltype(eltype(d))

#TODO: bivariate AnyDomain
immutable AnyDomain <: Domain{UnsetNumber} end
immutable EmptyDomain <: Domain{UnsetNumber} end

isambiguous(::AnyDomain) = true
dimension(::AnyDomain) = 1

complexlength(::AnyDomain) = NaN
arclength(::AnyDomain) = NaN

Base.reverse(a::Union{AnyDomain,EmptyDomain}) = a

canonicaldomain(a::Union{AnyDomain,EmptyDomain}) = a

Base.in(x::Domain,::EmptyDomain) = false

##General routines

Base.isempty(::EmptyDomain) = true
Base.isempty(::Domain) = false
Base.intersect(a::Domain,b::Domain) = a==b?a:EmptyDomain()
\(a::Domain,b::Domain) = setdiff(a,b)

## Interval Domains

abstract IntervalDomain{T} <: UnivariateDomain{T}

canonicaldomain(d::IntervalDomain) = Interval{real(eltype(eltype(d)))}()

Base.isapprox(a::Domain,b::Domain) = a==b
domainscompatible(a,b) = domainscompatible(domain(a),domain(b))
domainscompatible(a::Domain,b::Domain) = isambiguous(a) || isambiguous(b) ||
                    isapprox(a,b) || isapprox(a,reverse(b))

function chebyshevpoints{T<:Number}(::Type{T},n::Integer;kind::Integer=1)
    if kind == 1
        T[sinpi((n-2k-one(T))/2n) for k=n-1:-1:0]
    elseif kind == 2
        if n==1
            zeros(T,1)
        else
            T[cospi(k/(n-one(T))) for k=n-1:-1:0]
        end
    end
end
chebyshevpoints(n::Integer;kind::Integer=1) = chebyshevpoints(Float64,n;kind=kind)

##TODO: Should fromcanonical be fromcanonical!?

points{T}(d::IntervalDomain{T},n::Integer) = fromcanonical(d,chebyshevpoints(real(eltype(eltype(T))),n))  # eltype to handle point

points(d::AbstractVector,n::Integer)=points(convert(Domain,d),n)
bary(v::AbstractVector{Float64},d::IntervalDomain,x::Float64)=bary(v,tocanonical(d,x))

#TODO consider moving these
Base.first{T}(d::IntervalDomain{T})=fromcanonical(d,-one(T))
Base.last{T}(d::IntervalDomain{T})=fromcanonical(d,one(T))

Base.in(x,::AnyDomain)=true
function Base.in(x,d::IntervalDomain)
    T=real(eltype(eltype(eltype(d))))
    y=tocanonical(d,x)
    ry=real(y)
    iy=imag(y)
    sc=norm(fromcanonicalD(d,ry<-1?-1:(ry>1?1:ry)))  # scale based on stretch of map on projection to interal
    isapprox(fromcanonical(d,y),x) &&
        -one(T)-100eps(T)/sc≤ry≤one(T)+100eps(T)/sc &&
        -100eps(T)/sc≤iy≤100eps(T)/sc
end

pieces(d::Domain)=[d]
issubcomponent(a::Domain,b::Domain)=a in pieces(b)

###### Periodic domains

abstract PeriodicDomain{T} <: UnivariateDomain{T}

canonicaldomain(::PeriodicDomain)=PeriodicInterval()

points{T}(d::PeriodicDomain{T},n::Integer) = fromcanonical(d, fourierpoints(real(eltype(eltype(T))),n))

fourierpoints(n::Integer) = fourierpoints(Float64,n)
fourierpoints{T<:Number}(::Type{T},n::Integer)= convert(T,π)*collect(-n:2:n-2)/n

function Base.in{T}(x,d::PeriodicDomain{T})
    y=tocanonical(d,x)
    if !isapprox(fromcanonical(d,y),x)
        return false
    end

    l=arclength(d)
    if isinf(l)
        abs(imag(y))<20eps(T) && -π-2eps(T)<real(y)<π+2eps(T)
    else
        abs(imag(y))/l<20eps(T) && -π-2l*eps(T)<real(y)<π+2l*eps(T)
    end
end

Base.issubset(a::Domain,b::Domain)=a==b

Base.first(d::PeriodicDomain)=fromcanonical(d,-π)
Base.last(d::PeriodicDomain)=fromcanonical(d,π)

immutable AnyPeriodicDomain <: PeriodicDomain{UnsetNumber} end
isambiguous(::AnyPeriodicDomain)=true

Base.convert{D<:PeriodicDomain}(::Type{D},::AnyDomain)=AnyPeriodicDomain()

## conveninece routines

Base.ones(d::Domain)=ones(prectype(d),Space(d))
Base.zeros(d::Domain)=zeros(prectype(d),Space(d))

function commondomain(P::AbstractVector)
    ret = AnyDomain()

    for op in P
        d = domain(op)
        @assert ret == AnyDomain() || d == AnyDomain() || ret == d

        if d != AnyDomain()
            ret = d
        end
    end

    ret
end

commondomain{T<:Number}(P::AbstractVector,g::Array{T})=commondomain(P)
commondomain(P::AbstractVector,g)=commondomain([P;g])

domain(::Number)=AnyDomain()

## rand

Base.rand(d::IntervalDomain,k...)=fromcanonical(d,2rand(k...)-1)
Base.rand(d::PeriodicDomain,k...)=fromcanonical(d,2π*rand(k...)-π)

checkpoints(d::IntervalDomain)=fromcanonical(d,[-0.823972,0.3273484])
checkpoints(d::PeriodicDomain)=fromcanonical(d,[1.223972,-2.83273484])

## boundary

∂(d::Domain)=EmptyDomain()   # This is meant to be overriden
∂(d::IntervalDomain)=[first(d),last(d)] #TODO: Points domain
∂(d::PeriodicDomain)=EmptyDomain()

## map domains
# we auto vectorize arguments
tocanonical(d::Domain,x,y,z...)=tocanonical(d,Vec(x,y,z...))
fromcanonical(d::Domain,v::AbstractArray) =
    eltype(d)[fromcanonical(d,vk) for vk in v]

mappoint(d1::Domain,d2::Domain,x...)=fromcanonical(d2,tocanonical(d1,x...))
invfromcanonicalD(d::Domain,x...)=1./fromcanonicalD(d,x...)

## domains in higher dimensions

fromcanonical{V<:Vec}(d::Domain{V},p::AbstractArray) =
    V[fromcanonical(d,x) for x in p]

# points{T<:Array}(d::IntervalDomain{T},n::Integer) = T[fromcanonical(d,x) for x in chebyshevpoints(real(eltype(T)),n)]
# points{T<:Array}(d::PeriodicDomain{T},n::Integer) = T[fromcanonical(d,x) for x in fourierpoints(real(eltype(T)),n)]

## sorting
# we sort spaces lexigraphically by default

for OP in (:<,:(<=),:>,:(>=),:(Base.isless))
    @eval $OP(a::Domain,b::Domain)=$OP(string(a),string(b))
end

## Other special domains

immutable PositiveIntegers <: Domain{Int,0} end
immutable Integers <: Domain{Int,0} end

const ℕ = PositiveIntegers()
const ℤ = Integers()


text:
# general helper functions for integrators.

"""
Create nonlinear solver object for a system of `N` equations with data type `DT`.
The function ``f(x)=0`` to be solved for is determined by a julia function
`function_stages!(x, b, params)`, where `x` is the current solution and `b` is
the output vector, s.th. ``b = f(x)``. `params` are a set of parameters depending
on the equation and integrator that is used.
The solver type is obtained from the config dictionary (`:nls_solver`).
"""
function create_nonlinear_solver(DT, N, params)
    # create solution vector for nonlinear solver
    x = zeros(DT, N)

    # create wrapper function f(x,b) that calls `function_stages!(x, b, params)`
    # with the appropriate `params`
    f = (x,b) -> function_stages!(x, b, params)

    # create nonlinear solver with solver type obtained from config dictionary
    s = get_config(:nls_solver)(x, f)
end

"""
Create a solution vector of type `TwicePrecision{DT}` for a problem with `D` dimensions
and `M` independent initial conditions.
"""
function create_solution_vector(DT, D, M)
    x = Array{Vector{TwicePrecision{DT}}}(undef, M)

    for i in 1:M
        x[i] = zeros(TwicePrecision{DT}, D)
    end

    return x
end

"""
Create a solution vector of type `TwicePrecision{DT}` for a problem with `D` dimensions,
`NS' sample paths, and `NI` independent initial conditions.
"""
function create_solution_vector(DT, D, NS, NI)
    x = Array{Vector{TwicePrecision{DT}}}(undef, NS, NI)

    for i in 1:NS
        for j in 1:NI
            x[i,j] = zeros(TwicePrecision{DT}, D)
        end
    end

    return x
end

"""
Create a vector of S solution vectors of type DT to store the solution of S
internal stages for a problem with `D` dimensions.
"""
function create_internal_stage_vector(DT, D, S)
    a = Array{Vector{DT}}(undef, S)

    for i in 1:S
        a[i] = zeros(DT,D)
    end

    return a
end

"""
Create a vector of S solution matrices of type DT to store the solution of S
internal stages for a problem with `DxM` dimensions.
"""
function create_internal_stage_vector(DT, D, M, S)
    a = Array{Matrix{DT}}(undef, S)

    for i in 1:S
        a[i] = zeros(DT,D,M)
    end

    return a
end

"""
Create a vector of S+1 solution vectors of type DT to store the solution of S
internal stages and the solution of the previous timestep for a problem with `D`
    dimensions.
"""
function create_internal_stage_vector_with_zero(DT, D, S)
    a = OffsetArray{Vector{DT}}(undef, 0:S)

    for i in 0:S
        a[i] = zeros(DT,D)
    end

    return a
end

function check_solution_dimension_asserts(sol::Solution, m::Int, n::Int)
    @assert m ≥ 1
    @assert m ≤ sol.ni

    @assert n ≥ 1
    @assert n ≤ sol.ntime
end

function update_solution!(x::Vector{T}, xₑᵣᵣ::Vector{T}, ẋ::Matrix{T}, b::Vector{T}, Δt::T) where {T}
    @assert length(x) == length(xₑᵣᵣ)
    @assert length(x) == size(ẋ, 1)
    @assert length(b) == size(ẋ, 2)

    for k in axes(ẋ, 1)
        for i in axes(ẋ, 2)
            x[k], xₑᵣᵣ[k] = compensated_summation(Δt * b[i] * ẋ[k,i], x[k], xₑᵣᵣ[k])
        end
    end
end

function update_solution!(x::Vector{T}, xₑᵣᵣ::Vector{T}, ẋ::Vector{Vector{T}}, b::Vector{T}, Δt::T) where {T}
    @assert length(b) == length(ẋ)
    @assert length(x) == length(ẋ[1])
    @assert length(x) == length(xₑᵣᵣ)

    for i in eachindex(ẋ)
        for k in eachindex(ẋ[i])
            x[k], xₑᵣᵣ[k] = compensated_summation(Δt * b[i] * ẋ[i][k], x[k], xₑᵣᵣ[k])
        end
    end
end

function update_solution!(x::Union{Vector{T}, Vector{TwicePrecision{T}}}, ẋ::Matrix{T}, b::Vector{T}, Δt::T) where {T}
    @assert length(x) == size(ẋ, 1)
    @assert length(b) == size(ẋ, 2)

    local Δx::eltype(x)

    for k in axes(ẋ, 1)
        Δx = 0
        for i in axes(ẋ, 2)
            Δx += b[i] * ẋ[k,i]
        end
        x[k] += Δt * Δx
    end
end

# For stochastic Runge-Kutta methods (SIRK and WIRK)
# x - the solution vector to be updated
# Vx - the matrix containing the drift vector evaluated at the internal stages v(Q_i) (SIRK) or v(Q0_i) (WIRK)
# Bx - the array containing the diffusion matrix evaluated at the internal stages B(Q_i) (SIRK) or B(Q1^(l)_i) (WIRK)
# bdrift - the Runge-Kutta coefficients for the drift part
# bdiff - the Runge-Kutta coefficients for the diffusion part
# Δt - the time step
# ΔW - the increments of the Brownian motion (SFIRK) or the increments represented by the random variables \hat I^(k) (WFIRK)
function update_solution!(x::Union{Vector{T}, Vector{TwicePrecision{T}}}, Vx::Vector{Vector{T}}, Bx::Vector{Matrix{T}}, bdrift::Vector{T}, bdiff::Vector{T}, Δt::T, ΔW::Vector{T}) where {T}

    @assert length(bdrift) == length(bdiff) == length(Vx) == length(Bx)

    for i in eachindex(Vx, Bx)
        @assert length(x) == length(Vx[i]) == size(Bx[i], 1)
        @assert length(ΔW)== size(Bx[i], 2)
    end

    local Δx::eltype(x)

    # Contribution from the drift part
    for k in eachindex(x)
        Δx = 0.
        for i in eachindex(Vx)
            Δx += bdrift[i] * Vx[i][k]
        end
        x[k] += Δt * Δx
    end

    local Δy = zero(ΔW)

    # Contribution from the diffusion part
    for k in eachindex(x)
        Δy .= 0
        for i in eachindex(Bx)
            Δy .+= bdiff[i] .* Bx[i][k,:]
        end
        x[k] += dot(Δy,ΔW)
    end
end

# For stochastic Runge-Kutta methods
# x - the solution vector to be updated
# Vx - the matrix containing the drift vector evaluated at the internal stages v(Q_i)
# Bx - the array containing the diffusion matrix evaluated at the internal stages B(Q_i)
# bdrift - the Runge-Kutta coefficients for the drift part
# bdiff - the Runge-Kutta coefficients for the ΔW terms of the diffusion part
# bdiff2- the Runge-Kutta coefficients for the ΔZ terms of the diffusion part
# Δt - the time step
# ΔW - the increments of the Brownian motion
# ΔZ - the integrals of the increments of the Brownian motion
function update_solution!(x::Union{Vector{T}, Vector{TwicePrecision{T}}}, Vx::Vector{Vector{T}}, Bx::Vector{Matrix{T}}, bdrift::Vector{T}, bdiff::Vector{T}, bdiff2::Vector{T}, Δt::T, ΔW::Vector{T}, ΔZ::Vector{T}) where {T}

    @assert length(bdrift) == length(bdiff) == length(bdiff2) == length(Vx) == length(Bx)

    for i in eachindex(Vx, Bx)
        @assert length(x) == length(Vx[i]) == size(Bx[i], 1)
        @assert length(ΔW) == length(ΔZ) == size(Bx[i], 2)
    end

    local Δx::eltype(x)

    # Contribution from the drift part
    for k in eachindex(x)
        Δx = 0.
        for i in eachindex(Vx)
            Δx += bdrift[i] * Vx[i][k]
        end
        x[k] += Δt * Δx
    end

    local Δy = zero(ΔW)

    # Contribution from the diffusion part (ΔW terms)
    for k in eachindex(x)
        Δy .= 0
        for i in eachindex(Bx)
            Δy .+= bdiff[i] .* Bx[i][k,:]
        end
        x[k] += dot(Δy,ΔW)
    end

    # Contribution from the diffusion part (ΔZ terms)
    for k in eachindex(x)
        Δy .= 0
        for i in eachindex(Bx)
            Δy .+= bdiff2[i] .* Bx[i][k,:]
        end
        x[k] += dot(Δy,ΔZ)/Δt
    end
end

# For stochastic partitioned Runge-Kutta methods
# q, p - the solution vector to be updated
# Vqp, Fqp - the matrix containing the drift vectors evaluated at the internal stages v(Q_i), f(Q_i)
# Bqp, Gqp - the array containing the diffusion matrices evaluated at the internal stages B(Q_i), G(Q_i)
# bqdrift, bpdrift - the Runge-Kutta coefficients for the drift parts of the q and p equations
# bqdiff, bpdiff - the Runge-Kutta coefficients for the diffusion parts of the q and p equations
# Δt - the time step
# ΔW - the increments of the Brownian motion
function update_solution!(q::Vector{VT}, p::Vector{VT},
                          Vqp::Vector{Vector{T}}, Fqp::Vector{Vector{T}},
                          Bqp::Vector{Matrix{T}}, Gqp::Vector{Matrix{T}},
                          bqdrift::Vector{T}, bqdiff::Vector{T},
                          bpdrift::Vector{T}, bpdiff::Vector{T},
                          Δt::T, ΔW::Vector{T}) where {T, VT <: Union{T, TwicePrecision{T}}}

    @assert length(bqdrift) == length(bqdiff) == length(bpdrift) == length(bpdiff) == length(Vqp) == length(Fqp) == length(Bqp) == length(Gqp)

    for i in eachindex(Vqp, Fqp, Bqp, Gqp)
        @assert length(q) == length(p) == length(Vqp[i]) == length(Fqp[i]) == size(Bqp[i], 1) == size(Gqp[i], 1)
        @assert length(ΔW) == size(Bqp[i], 2) == size(Gqp[i], 2)
    end

    local Δq::eltype(q)
    local Δp::eltype(p)

    # Contribution from the drift part
    for k in eachindex(q, p)
        Δq = 0.
        Δp = 0.
        for i in eachindex(Vqp, Fqp)
            Δq += bqdrift[i] * Vqp[i][k]
            Δp += bpdrift[i] * Fqp[i][k]
        end
        q[k] += Δt * Δq
        p[k] += Δt * Δp
    end

    local Δy = zero(ΔW)
    local Δz = zero(ΔW)

    # Contribution from the diffusion part
    for k in eachindex(q, p)
        Δy .= 0
        Δz .= 0

        for i in eachindex(Bqp, Gqp)
            Δy .+= bqdiff[i] .* Bqp[i][k,:]
            Δz .+= bpdiff[i] .* Gqp[i][k,:]
        end

        q[k] += dot(Δy,ΔW)
        p[k] += dot(Δz,ΔW)
    end
end

# For stochastic split partitioned Runge-Kutta methods
# q, p - the solution vector to be updated
# Vqp, Fqp1, Fqp2 - the matrix containing the drift vectors evaluated at the internal stages v(Q_i), fi(Q_i)
# Bqp, Gqp1, Gqp2 - the array containing the diffusion matrices evaluated at the internal stages B(Q_i), Gi(Q_i)
# bqdrift, bpdrift1, bpdrift2 - the Runge-Kutta coefficients for the drift parts of the q and p equations
# bqdiff, bpdiff1, bpdiff2 - the Runge-Kutta coefficients for the diffusion parts of the q and p equations
# Δt - the time step
# ΔW - the increments of the Brownian motion
function update_solution!(q::Union{Vector{T}, Vector{TwicePrecision{T}}}, p::Union{Vector{T}, Vector{TwicePrecision{T}}},
                            Vqp::Vector{Vector{T}}, Fqp1::Vector{Vector{T}}, Fqp2::Vector{Vector{T}},
                            Bqp::Vector{Matrix{T}}, Gqp1::Vector{Matrix{T}}, Gqp2::Vector{Matrix{T}},
                            bqdrift::Vector{T}, bqdiff::Vector{T},
                            bpdrift1::Vector{T}, bpdrift2::Vector{T},
                            bpdiff1::Vector{T}, bpdiff2::Vector{T},
                            Δt::T, ΔW::Vector{T}) where {T}

   @assert length(bqdrift) == length(bqdiff) == length(bpdrift1) == length(bpdrift2) == length(bpdiff1) == length(bpdiff2) == length(Vqp) == length(Fqp1) == length(Fqp2) == length(Bqp) == length(Gqp1) == length(Gqp2)

   for i in eachindex(Vqp, Fqp1, Fqp2, Bqp, Gqp1, Gqp2)
       @assert length(q) == length(p) == length(Vqp[i]) == length(Fqp1[i]) == length(Fqp2[i]) == size(Bqp[i], 1) == size(Gqp1[i], 1) == size(Gqp2[i], 1)
       @assert length(ΔW) == size(Bqp[i], 2) == size(Gqp1[i], 2) == size(Gqp2[i], 2)
   end

   local Δq::eltype(q)
   local Δp::eltype(p)

   # Contribution from the drift part
   for k in eachindex(q, p)
       Δq = 0.
       Δp = 0.
       for i in eachindex(Vqp, Fqp1, Fqp2)
           Δq += bqdrift[i] * Vqp[i][k]
           Δp += bpdrift1[i] * Fqp1[i][k] + bpdrift2[i] * Fqp2[i][k]
       end
       q[k] += Δt * Δq
       p[k] += Δt * Δp
   end

   local Δy = zero(ΔW)
   local Δz = zero(ΔW)

   # Contribution from the diffusion part
   for k in eachindex(q, p)
       Δy .= 0
       Δz .= 0

       for i in eachindex(Bqp, Gqp1, Gqp2)
           Δy .+= bqdiff[i]  .* Bqp[i][k,:]
           Δz .+= bpdiff1[i] .* Gqp1[i][k,:] .+ bpdiff2[i] .* Gqp2[i][k,:]
       end

       q[k] += dot(Δy,ΔW)
       p[k] += dot(Δz,ΔW)
   end
end

# For weak Runge-Kutta methods WERK
# x - the solution vector to be updated
# Vx - the matrix containing the drift vector evaluated at the internal stages v(Q_i)
# Bx1 - the array containing the diffusion matrix evaluated at the internal stages H^(l)_i, such that Bx1[:,l,i] is evaluated at H^(l)_i
# Bx2 - the array containing the diffusion matrix evaluated at the internal stages \hat H^(l)_i, such that Bx2[:,l,i] is evaluated at \hat H^(l)_i
# alpha - the Runge-Kutta coefficients for the drift part
# beta1 - the Runge-Kutta coefficients for the diffusion term with the random increments
# beta2 - the Runge-Kutta coefficients for the second diffusion term
# Δt - the time step
# ΔW - the increments of the Brownian motion represented by the random variables \hat I^(k)
function update_solution!(x::Union{Vector{T}, Vector{TwicePrecision{T}}}, Vx::Vector{Vector{T}}, Bx1::Vector{Matrix{T}}, Bx2::Vector{Matrix{T}}, alpha::Vector{T}, beta1::Vector{T}, beta2::Vector{T}, Δt::T, ΔW::Vector{T}) where {T}
    @assert length(alpha) == length(beta1) == length(beta2) == length(Vx) == length(Bx1) == length(Bx2)

    for i in eachindex(Vx, Bx1, Bx2)
        @assert length(x) == length(Vx[i]) == size(Bx1[i], 1) == size(Bx2[i], 1)
        @assert length(ΔW)== size(Bx1[i], 2) == size(Bx2[i], 2)
    end

    local Δx::eltype(x)

    # Contribution from the drift part
    for k in eachindex(x)
        Δx = 0.
        for i in eachindex(Vx)
            Δx += alpha[i] * Vx[i][k]
        end
        x[k] += Δt * Δx
    end

    local Δy = zero(ΔW)

    # Contribution from the diffusion term with the random variables I^(k)_i
    for k in eachindex(x)
        Δy .= 0
        for i in eachindex(Bx1)
            Δy .+= beta1[i] .* Bx1[i][k,:]
        end
        x[k] += dot(Δy,ΔW)
    end

    # Contribution from the second diffusion term
    for k in eachindex(x)
        Δx= 0.
        for i in eachindex(Bx2)
            for l in axes(Bx2[i], 2)
                Δx += beta2[i] * Bx2[i][k,l]
            end
        end
        x[k] += sqrt(Δt)*Δx
    end
end

function update_solution!(x::Union{Vector{T}, Vector{TwicePrecision{T}}}, ẋ::Vector{Vector{T}}, b::Vector{T}, Δt::T) where {T}
    @assert length(b) == length(ẋ)
    @assert length(x) == length(ẋ[1])

    for i in eachindex(ẋ)
        for k in eachindex(ẋ[i])
            x[k] += Δt * b[i] * ẋ[i][k]
        end
    end
end

function update_solution!(x::Vector{T}, xₑᵣᵣ::Vector{T}, ẋ::Union{Matrix{T},Vector{Vector{T}}}, b::Vector{T}, b̂::Vector, Δt::T) where {T}
    update_solution!(x, xₑᵣᵣ, ẋ, b, Δt)
    update_solution!(x, xₑᵣᵣ, ẋ, b̂, Δt)
end

function update_solution!(x::Union{Vector{T}, Vector{TwicePrecision{T}}}, ẋ::Union{Matrix{T},Vector{Vector{T}}}, b::Vector{T}, b̂::Vector, Δt::T) where {T}
    update_solution!(x, ẋ, b, Δt)
    update_solution!(x, ẋ, b̂, Δt)
end

# For stochastic Runge-Kutta methods
function update_solution!(x::Union{Vector{T}, Vector{TwicePrecision{T}}}, Vx::Vector{Vector{T}}, Bx::Vector{Matrix{T}},
                            bdrift::Vector{T}, b̂drift::Vector, bdiff::Vector{T}, b̂diff::Vector, Δt::T, ΔW::Vector{T}) where {T}
    update_solution!(x, Vx, Bx, bdrift, bdiff, Δt, ΔW)
    update_solution!(x, Vx, Bx, b̂drift, b̂diff, Δt, ΔW)
end

# For stochastic partitioned Runge-Kutta methods
function update_solution!(q::Union{Vector{T}, Vector{TwicePrecision{T}}}, p::Union{Vector{T}, Vector{TwicePrecision{T}}},
                            Vqp::Vector{Vector{T}}, Fqp::Vector{Vector{T}}, Bqp::Vector{Matrix{T}}, Gqp::Vector{Matrix{T}},
                            bqdrift::Vector{T}, b̂qdrift::Vector, bqdiff::Vector{T}, b̂qdiff::Vector,
                            bpdrift::Vector{T}, b̂pdrift::Vector, bpdiff::Vector{T}, b̂pdiff::Vector, Δt::T, ΔW::Vector{T}) where {T}
    update_solution!(q, p, Vqp, Fqp, Bqp, Gqp, bqdrift, bqdiff, bpdrift, bpdiff, Δt, ΔW)
    update_solution!(q, p, Vqp, Fqp, Bqp, Gqp, b̂qdrift, b̂qdiff, b̂pdrift, b̂pdiff, Δt, ΔW)
end

# For stochastic split partitioned Runge-Kutta methods
function update_solution!(q::Union{Vector{T}, Vector{TwicePrecision{T}}}, p::Union{Vector{T}, Vector{TwicePrecision{T}}},
                            Vqp::Vector{Vector{T}}, Fqp1::Vector{Vector{T}}, Fqp2::Vector{Vector{T}},
                            Bqp::Vector{Matrix{T}}, Gqp1::Vector{Matrix{T}}, Gqp2::Vector{Matrix{T}},
                            bqdrift::Vector{T}, b̂qdrift::Vector, bqdiff::Vector{T}, b̂qdiff::Vector,
                            bpdrift1::Vector{T}, b̂pdrift1::Vector, bpdrift2::Vector{T}, b̂pdrift2::Vector,
                            bpdiff1::Vector{T}, b̂pdiff1::Vector, bpdiff2::Vector{T}, b̂pdiff2::Vector,
                            Δt::T, ΔW::Vector{T}) where {T}
    update_solution!(q, p, Vqp, Fqp1, Fqp2, Bqp, Gqp1, Gqp2, bqdrift, bqdiff, bpdrift1, bpdrift2, bpdiff1, bpdiff2, Δt, ΔW)
    update_solution!(q, p, Vqp, Fqp1, Fqp2, Bqp, Gqp1, Gqp2, b̂qdrift, b̂qdiff, b̂pdrift1, b̂pdrift2, b̂pdiff1, b̂pdiff2, Δt, ΔW)
end

 # For weak Runge-Kutta methods WERK
 function update_solution!(x::Union{Vector{T}, Vector{TwicePrecision{T}}}, Vx::Vector{Vector{T}}, Bx1::Vector{Matrix{T}}, Bx2::Vector{Matrix{T}},
                             alpha::Vector{T}, âlpha::Vector, beta1::Vector{T}, b̂eta1::Vector, beta2::Vector{T}, b̂eta2::Vector, Δt::T, ΔW::Vector{T}) where {T}
     update_solution!(x, Vx, Bx1, Bx2, alpha, beta1, beta2, Δt, ΔW)
     update_solution!(x, Vx, Bx1, Bx2, âlpha, b̂eta1, b̂eta2, Δt, ΔW)
 end

function update_multiplier!(λ::SolutionVector{T}, Λ::Vector{Vector{T}}, b::Vector{T}) where {T}
    @assert length(λ) == size(Λ, 1)
    local t::T
    @inbounds for i in eachindex(λ)
        t = zero(T)
        for j=eachindex(b)
            t += b[j] * Λ[j][i]
        end
        λ[i] = t
    end
    nothing
end

function cut_periodic_solution!(x::Vector{T}, periodicity::Vector{T}) where {T}
    @assert length(x) == length(periodicity)

    for k in eachindex(x, periodicity)
        if periodicity[k] ≠ 0
            while x[k] < 0
                x[k] += periodicity[k]
            end
            while x[k] ≥ periodicity[k]
                x[k] -= periodicity[k]
            end
        end
    end
end

function cut_periodic_solution!(x::Vector{TwicePrecision{T}}, periodicity::Vector{T}) where {T}
    @assert length(x) == length(periodicity)

    for k in eachindex(x, periodicity)
        if periodicity[k] ≠ 0
            while x[k].hi < 0
                x[k] += periodicity[k]
            end
            while x[k].hi ≥ periodicity[k]
                x[k] -= periodicity[k]
            end
        end
    end
end

function cut_periodic_solution!(x::Vector{T}, periodicity::Vector{T}, shift::Vector{T}) where {T}
    @assert length(x) == length(periodicity)
    shift .= 0
    for k in eachindex(x, periodicity, shift)
        if periodicity[k] ≠ 0
            while x[k] + shift[k] < 0
                shift[k] += periodicity[k]
            end
            while x[k] + shift[k] ≥ periodicity[k]
                shift[k] -= periodicity[k]
            end
        end
    end
end

function cut_periodic_solution!(x::Vector{TwicePrecision{T}}, periodicity::Vector{T}, shift::Vector{T}) where {T}
    @assert length(x) == length(periodicity)
    shift .= 0
    for k in eachindex(x, periodicity, shift)
        if periodicity[k] ≠ 0
            while x[k].hi + shift[k] < 0
                shift[k] += periodicity[k]
            end
            while x[k].hi + shift[k] ≥ periodicity[k]
                shift[k] -= periodicity[k]
            end
        end
    end
end


text:
immutable ODEModel{P<:AbstractParameter,D<:AbstractData,N<:AbstractNoiseModel} <: AbstractModel
    #Generic model specs
    name::AbstractString
    parameters::Vector

    #Data + noise
    measurements::D
    noisemodel::N

    #ODE specific
    ode::Function
    initial::Vector
    numstates::Int
    observed::Vector{Int}
    unobserved::Vector{Int}
    abstol::Real
    reltol::Real

    ODEModel(name::AbstractString,parameters::Vector{P},measurements::D,noisemodel::N,
             ode::Function,initial::Vector,numstates::Int,observed::Vector{Int},unobserved::Vector{Int},abstol::Real,reltol::Real) =
        new(name,parameters,measurements,noisemodel,ode,initial,numstates,observed,unobserved,abstol,reltol)
end

###Factory function with default values for tolerance
function _model(::Type{Val{:ode}},parameters::Vector,data::AbstractData,noise::AbstractNoiseModel,
       ode::Function,initial::Vector,numstates::Int,observed::Vector{Int};name ="ODE",abstol =1e-6,reltol =1e-6)
    P = eltype(parameters)
    D = typeof(data)
    N = typeof(noise)
    ODEModel{P,D,N}(name,parameters,data,noise,ode,initial,numstates,observed,setdiff(1:numstates,observed),abstol,reltol)
end

###Evaluate the model for the given parameter values
function evaluate!(m::ODEModel,vals::AbstractVector)
    o(t,y,ydot) = m.ode(t,y,ydot,vals)
    sub(Sundials.cvode(o,m.initial,dataindex(m.measurements);reltol=m.reltol,abstol=m.abstol),:,m.observed)
end

###Utility functions used in generic implementations in AbstractModel
@inline dataindex(m::ODEModel) = dataindex(m.measurements)
@inline measurements(m::ODEModel) = datavalues(m.measurements)
@inline noisemodel(m::ODEModel) = m.noisemodel

###Base functionality
function show(io::IO,m::ODEModel)
    println(io,"ODEModel ",m.name)
    print(io,"parameters: ") ; show(io,m.parameters)
    print(io,"measurements: ") ; show(io,m.measurements)
    print(io,"noisemodel: ") ; show(io,m.noisemodel)
    println(io,"ode: ",m.ode)
    println(io,"initial: ",m.initial)
    println(io,"numstates: ",m.numstates)
    println(io,"observed: ",m.observed)
    if ~isempty(m.unobserved)
        println(io,"unobserved: ",m.unobserved)
    end
    println(io,"abstol: ",m.abstol)
    println(io,"reltol: ",m.reltol)
    nothing
end