text:
using Revise
using Base.Test
using DataStructures: OrderedSet
using Compat

to_remove = String[]

@testset "Revise" begin

    function collectexprs(ex::Revise.RelocatableExpr)
        exs = Revise.RelocatableExpr[]
        for item in Revise.LineSkippingIterator(ex.args)
            push!(exs, item)
        end
        exs
    end

    @static if is_apple()
        yry() = (sleep(1.1); revise(); sleep(1.1))
    else
        yry() = (sleep(0.1); revise(); sleep(0.1))
    end

    @testset "LineSkipping" begin
        ex = Revise.relocatable!(quote
                                 f(x) = x^2
                                 g(x) = sin(x)
                                 end)
        @test length(ex.args) == 4  # including the line number expressions
        exs = collectexprs(ex)
        @test length(exs) == 2
        @test isequal(exs[1], Revise.relocatable!(:(f(x) = x^2)))
        @test !isequal(exs[2], Revise.relocatable!(:(f(x) = x^2)))
        @test isequal(exs[2], Revise.relocatable!(:(g(x) = sin(x))))
        @test !isequal(exs[1], Revise.relocatable!(:(g(x) = sin(x))))
    end

    @testset "Parse errors" begin
        warnfile = joinpath(tempdir(), randstring(10))
        open(warnfile, "w") do io
            redirect_stderr(io) do
                md = Revise.ModDict(Main=>OrderedSet{Revise.RelocatableExpr}())
                @test !Revise.parse_source!(md, """
f(x) = 1
g(x) = 2
h{x) = 3  # error
k(x) = 4
""",
                                            :test, 1, Main, tempdir())
                @test convert(Revise.RelocatableExpr, :(g(x) = 2)) ∈ md[Main]
            end
        end
        @test contains(read(warnfile, String), "parsing error near line 3")
        rm(warnfile)
    end

    @testset "Comparison" begin
        fl1 = joinpath(@__DIR__, "revisetest.jl")
        fl2 = joinpath(@__DIR__, "revisetest_revised.jl")
        fl3 = joinpath(@__DIR__, "revisetest_errors.jl")
        include(fl1)  # So the modules are defined
        # test the "mistakes"
        @test ReviseTest.cube(2) == 16
        @test ReviseTest.Internal.mult3(2) == 8
        oldmd = Revise.parse_source(fl1, Main, dirname(fl1))
        newmd = Revise.parse_source(fl2, Main, nothing)
        revmd = Revise.revised_statements(newmd, oldmd)
        @test length(revmd) == 2
        @test haskey(revmd, ReviseTest) && haskey(revmd, ReviseTest.Internal)
        @test length(revmd[ReviseTest.Internal]) == 1
        cmp = Revise.relocatable!(quote
            mult3(x) = 3*x
        end)
        @test isequal(first(revmd[ReviseTest.Internal]), collectexprs(cmp)[1])
        @test length(revmd[ReviseTest]) == 2
        cmp = Revise.relocatable!(quote
            cube(x) = x^3
            fourth(x) = x^4  # this is an addition to the file
        end)
        @test isequal(revmd[ReviseTest], OrderedSet(collectexprs(cmp)))

        Revise.eval_revised(revmd)
        @test ReviseTest.cube(2) == 8
        @test ReviseTest.Internal.mult3(2) == 6

        # Backtraces
        newmd = Revise.parse_source(fl3, Main, nothing)
        revmd = Revise.revised_statements(newmd, oldmd)
        Revise.eval_revised(revmd)
        try
            ReviseTest.cube(2)
            @test false
        catch err
            @test isa(err, ErrorException) && err.msg == "cube"
            bt = first(stacktrace(catch_backtrace()))
            @test bt.func == :cube && bt.file == Symbol(fl3) && bt.line == 7
        end
        try
            ReviseTest.Internal.mult2(2)
            @test false
        catch err
            @test isa(err, ErrorException) && err.msg == "mult2"
            bt = first(stacktrace(catch_backtrace()))
            @test bt.func == :mult2 && bt.file == Symbol(fl3) && bt.line == 13
        end
    end

    @testset "File paths" begin
        testdir = joinpath(tempdir(), randstring(10))
        mkdir(testdir)
        push!(to_remove, testdir)
        push!(LOAD_PATH, testdir)
        for (pcflag, fbase) in ((true, "pc"), (false, "npc"))  # precompiled & not
            modname = uppercase(fbase)
            # Create a package with the following structure:
            #   src/PkgName.jl   # PC.jl = precompiled, NPC.jl = nonprecompiled
            #   src/file2.jl
            #   src/subdir/file3.jl
            #   src/subdir/file4.jl
            # exploring different ways of expressing the `include` statement
            dn = joinpath(testdir, modname, "src")
            mkpath(dn)
            open(joinpath(dn, modname*".jl"), "w") do io
                println(io, """
__precompile__($pcflag)

module $modname

export $(fbase)1, $(fbase)2, $(fbase)3, $(fbase)4, $(fbase)5, using_macro_$(fbase)

$(fbase)1() = 1

include("file2.jl")
include("subdir/file3.jl")
include(joinpath(@__DIR__, "subdir", "file4.jl"))
otherfile = "file5.jl"
include(otherfile)

# Update order check: modifying `some_macro_` to return -6 doesn't change the
# return value of `using_macro_` (issue #20) unless `using_macro_` is also updated,
# *in this order*:
#   1. update the `@some_macro_` definition
#   2. update the `using_macro_` definition
macro some_macro_$(fbase)()
    return 6
end
using_macro_$(fbase)() = @some_macro_$(fbase)()

end
""")
            end
            open(joinpath(dn, "file2.jl"), "w") do io
                println(io, "$(fbase)2() = 2")
            end
            mkdir(joinpath(dn, "subdir"))
            open(joinpath(dn, "subdir", "file3.jl"), "w") do io
                println(io, "$(fbase)3() = 3")
            end
            open(joinpath(dn, "subdir", "file4.jl"), "w") do io
                println(io, "$(fbase)4() = 4")
            end
            open(joinpath(dn, "file5.jl"), "w") do io
                println(io, "$(fbase)5() = 5")
            end
            @eval using $(Symbol(modname))
            fn1, fn2 = Symbol("$(fbase)1"), Symbol("$(fbase)2")
            fn3, fn4 = Symbol("$(fbase)3"), Symbol("$(fbase)4")
            fn5 = Symbol("$(fbase)5")
            fn6 = Symbol("using_macro_$(fbase)")
            @eval @test $(fn1)() == 1
            @eval @test $(fn2)() == 2
            @eval @test $(fn3)() == 3
            @eval @test $(fn4)() == 4
            @eval @test $(fn5)() == 5
            @eval @test $(fn6)() == 6
            sleep(0.1)  # to ensure that the file watching has kicked in
            # Change the definition of function 1 (easiest to just rewrite the whole file)
            open(joinpath(dn, modname*".jl"), "w") do io
                println(io, """
__precompile__($pcflag)
module $modname
export $(fbase)1, $(fbase)2, $(fbase)3, $(fbase)4, $(fbase)5, using_mac$(fbase)
$(fbase)1() = -1
include("file2.jl")
include("subdir/file3.jl")
include(joinpath(@__DIR__, "subdir", "file4.jl"))
otherfile = "file5.jl"
include(otherfile)

macro some_macro_$(fbase)()
    return -6
end
using_macro_$(fbase)() = @some_macro_$(fbase)()

end
""")  # just for fun we skipped the whitespace
            end
            yry()
            @eval @test $(fn1)() == -1
            @eval @test $(fn2)() == 2
            @eval @test $(fn3)() == 3
            @eval @test $(fn4)() == 4
            @eval @test $(fn5)() == 5
            @eval @test $(fn6)() == 6      # because it hasn't been re-macroexpanded
            revise(eval(Symbol(modname)))
            @eval @test $(fn6)() == -6
            # Redefine function 2
            open(joinpath(dn, "file2.jl"), "w") do io
                println(io, "$(fbase)2() = -2")
            end
            yry()
            @eval @test $(fn1)() == -1
            @eval @test $(fn2)() == -2
            @eval @test $(fn3)() == 3
            @eval @test $(fn4)() == 4
            @eval @test $(fn5)() == 5
            @eval @test $(fn6)() == -6
            open(joinpath(dn, "subdir", "file3.jl"), "w") do io
                println(io, "$(fbase)3() = -3")
            end
            yry()
            @eval @test $(fn1)() == -1
            @eval @test $(fn2)() == -2
            @eval @test $(fn3)() == -3
            @eval @test $(fn4)() == 4
            @eval @test $(fn5)() == 5
            @eval @test $(fn6)() == -6
            open(joinpath(dn, "subdir", "file4.jl"), "w") do io
                println(io, "$(fbase)4() = -4")
            end
            yry()
            @eval @test $(fn1)() == -1
            @eval @test $(fn2)() == -2
            @eval @test $(fn3)() == -3
            @eval @test $(fn4)() == -4
            @eval @test $(fn5)() == 5
            @eval @test $(fn6)() == -6
            open(joinpath(dn, "file5.jl"), "w") do io
                println(io, "$(fbase)5() = -5")
            end
            yry()
            @eval @test $(fn1)() == -1
            @eval @test $(fn2)() == -2
            @eval @test $(fn3)() == -3
            @eval @test $(fn4)() == -4
            @eval @test $(fn5)() == -5
            @eval @test $(fn6)() == -6
            # Check module2files
            files = [joinpath(dn, modname*".jl"), joinpath(dn, "file2.jl"),
                     joinpath(dn, "subdir", "file3.jl"),
                     joinpath(dn, "subdir", "file4.jl"),
                     joinpath(dn, "file5.jl")]
            @test Revise.module2files[Symbol(modname)] == files
        end
        # Remove the precompiled file
        rm(joinpath(Base.LOAD_CACHE_PATH[1], "PC.ji"))

        # Test files paths that can't be statically parsed
        dn = joinpath(testdir, "LoopInclude", "src")
        mkpath(dn)
        open(joinpath(dn, "LoopInclude.jl"), "w") do io
            println(io, """
module LoopInclude

export li_f, li_g

for fn in ("file1.jl", "file2.jl")
    include(fn)
end

end
""")
        end
        open(joinpath(dn, "file1.jl"), "w") do io
            println(io, "li_f() = 1")
        end
        open(joinpath(dn, "file2.jl"), "w") do io
            println(io, "li_g() = 2")
        end
        @eval using LoopInclude
        @test li_f() == 1
        @test li_g() == 2
        sleep(0.1)  # ensure watching is set up
        open(joinpath(dn, "file1.jl"), "w") do io
            println(io, "li_f() = -1")
        end
        @test li_f() == 1  # unless the include is at toplevel it is not found

        @test isfile(Revise.sysimg_path)

        pop!(LOAD_PATH)
    end

    # issue #36
    @testset "@__FILE__" begin
        testdir = joinpath(tempdir(), randstring(10))
        mkdir(testdir)
        push!(to_remove, testdir)
        push!(LOAD_PATH, testdir)
        dn = joinpath(testdir, "ModFILE", "src")
        mkpath(dn)
        open(joinpath(dn, "ModFILE.jl"), "w") do io
            println(io, """
__precompile__()

module ModFILE

mf() = @__FILE__, 1

end
""")
        end
        @eval using ModFILE
        @test ModFILE.mf() == (joinpath(dn, "ModFILE.jl"), 1)
        sleep(0.1)
        open(joinpath(dn, "ModFILE.jl"), "w") do io
            println(io, """
__precompile__()

module ModFILE

mf() = @__FILE__, 2

end
""")
        end
        yry()
        @test ModFILE.mf() == (joinpath(dn, "ModFILE.jl"), 2)
    end

    # issue #8
    @testset "Module docstring" begin
        testdir = joinpath(tempdir(), randstring(10))
        mkdir(testdir)
        push!(to_remove, testdir)
        push!(LOAD_PATH, testdir)
        dn = joinpath(testdir, "ModDocstring", "src")
        mkpath(dn)
        open(joinpath(dn, "ModDocstring.jl"), "w") do io
            println(io, """
" Ahoy! "
module ModDocstring

include("dependency.jl")

f() = 1

end
""")
        end
        open(joinpath(dn, "dependency.jl"), "w") do io
            println(io, "")
        end
        @eval using ModDocstring
        @test ModDocstring.f() == 1
        ds = @doc ModDocstring
        @test ds.content[1].content[1].content[1] == "Ahoy! "

        sleep(0.1)  # ensure watching is set up
        open(joinpath(dn, "ModDocstring.jl"), "w") do io
            println(io, """
" Ahoy! "
module ModDocstring

include("dependency.jl")

f() = 2

end
""")
        end
        yry()
        @test ModDocstring.f() == 2
        ds = @doc ModDocstring
        @test ds.content[1].content[1].content[1] == "Ahoy! "

        open(joinpath(dn, "ModDocstring.jl"), "w") do io
            println(io, """
" Hello! "
module ModDocstring

include("dependency.jl")

f() = 3

end
""")
        end
        yry()
        @test ModDocstring.f() == 3
        ds = @doc ModDocstring
        @test ds.content[end].content[1].content[1] == "Hello! "

        pop!(LOAD_PATH)
    end

    @testset "Line numbers" begin
        # issue #27
        testdir = joinpath(tempdir(), randstring(10))
        mkdir(testdir)
        push!(to_remove, testdir)
        push!(LOAD_PATH, testdir)
        modname = "LineNumberMod"
        dn = joinpath(testdir, modname, "src")
        mkpath(dn)
        open(joinpath(dn, modname*".jl"), "w") do io
            println(io, """
module $modname
include("incl.jl")
end
""")
        end
        open(joinpath(dn, "incl.jl"), "w") do io
            println(io, """
0
0
1
2
3
4
5
6
7
8

function foo(x)
    return x+5
end

foo(y::Int) = y-51
""")
        end
        @eval using LineNumberMod
        lines = Int[]
        files = String[]
        for m in methods(LineNumberMod.foo)
            push!(files, String(m.file))
            push!(lines, m.line)
        end
        @test all(f->endswith(string(f), "incl.jl"), files)
        sleep(0.1)  # ensure watching is set up
        open(joinpath(dn, "incl.jl"), "w") do io
            println(io, """
0
0
1
2
3
4
5
6
7
8

function foo(x)
    return x+6
end

foo(y::Int) = y-51
""")
        end
        yry()
        for m in methods(LineNumberMod.foo)
            @test endswith(string(m.file), "incl.jl")
            @test m.line ∈ lines
        end
    end

    # Issue #43
    @testset "New submodules" begin
        testdir = joinpath(tempdir(), randstring(10))
        mkdir(testdir)
        push!(to_remove, testdir)
        push!(LOAD_PATH, testdir)
        dn = joinpath(testdir, "Submodules", "src")
        mkpath(dn)
        open(joinpath(dn, "Submodules.jl"), "w") do io
            println(io, """
module Submodules
f() = 1
end
""")
        end
        @eval using Submodules
        @test Submodules.f() == 1
        sleep(0.1)  # ensure watching is set up
        open(joinpath(dn, "Submodules.jl"), "w") do io
            println(io, """
module Submodules
f() = 1
module Sub
g() = 2
end
end
""")
        end
        yry()
        @test Submodules.f() == 1
        @test Submodules.Sub.g() == 2
    end

    @testset "Pkg exclusion" begin
        push!(Revise.dont_watch_pkgs, :Example)
        push!(Revise.silence_pkgs, :Example)
        @eval import Example
        for k in keys(Revise.file2modules)
            if contains(k, "Example")
                error("Should not track files in Example")
            end
        end
        # Ensure that silencing works
        sfile = Revise.silencefile[]  # remember the original
        try
            sfiletemp = tempname()
            Revise.silencefile[] = sfiletemp
            Revise.silence("GSL")
            @test isfile(sfiletemp)
            pkgs = readlines(sfiletemp)
            @test any(p->p=="GSL", pkgs)
            rm(sfiletemp)
        finally
            Revise.silencefile[] = sfile
        end
    end

    @testset "Manual track" begin
        srcfile = joinpath(tempdir(), randstring(10)*".jl")
        open(srcfile, "w") do io
            print(io, """
revise_f(x) = 1
""")
        end
        include(srcfile)
        @test revise_f(10) == 1
        Revise.track(srcfile)
        sleep(0.1)
        open(srcfile, "w") do io
            print(io, """
revise_f(x) = 2
""")
        end
        yry()
        @test revise_f(10) == 2
        push!(to_remove, srcfile)

        # Do it again with a relative path
        curdir = pwd()
        cd(tempdir())
        srcfile = randstring(10)*".jl"
        open(srcfile, "w") do io
            print(io, """
        revise_floc(x) = 1
        """)
        end
        include(joinpath(pwd(), srcfile))
        @test revise_floc(10) == 1
        Revise.track(srcfile)
        sleep(0.1)
        open(srcfile, "w") do io
            print(io, """
        revise_floc(x) = 2
        """)
        end
        yry()
        @test revise_floc(10) == 2
        push!(to_remove, joinpath(tempdir(), srcfile))
        cd(curdir)

        # Tracking Base
        Revise.track(Base)
        @test any(k->endswith(k, "number.jl"), keys(Revise.file2modules))
    end

    @testset "Cleanup" begin
        warnfile = joinpath(tempdir(), randstring(10))
        open(warnfile, "w") do io
            redirect_stderr(io) do
                for name in to_remove
                    try
                        rm(name; force=true, recursive=true)
                        deleteat!(LOAD_PATH, find(LOAD_PATH .== name))
                    end
                end
                yry()
            end
        end
        if !is_apple()
            @test contains(read(warnfile, String), "is not an existing directory")
        end
        rm(warnfile)
    end
end


text:
using Knet: cuda01

function cuda01src(f, j=f, ex="$f(s,xi)"; BLK=256, THR=256)
"""
__global__ void _$(f)_32_01(int n, float s, float *x, float *y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    float xi = x[i];
    y[i] = $ex;
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  void $(f)_32_01(int n, float s, float *x, float *y) {
    _$(f)_32_01<<<$BLK,$THR>>>(n,s,x,y);
  }
}
__global__ void _$(f)_64_01(int n, double s, double *x, double *y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    double xi = x[i];
    y[i] = $ex;
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  void $(f)_64_01(int n, double s, double *x, double *y) {
    _$(f)_64_01<<<$BLK,$THR>>>(n,s,x,y);
  }
}
"""
end

for a in cuda01
    isa(a,Tuple) || (a=(a,))
    print(cuda01src(a...))
end

ext: jl


text:
module day19

using Combinatorics

function parserule(r)
    id, spec = split(r, ": ")
    id = parse(Int32, id)
    if spec[1] == '"'
        return id => (:matchchar, spec[2])
    else
        subrules = map(x -> parse.(Int32, split(x, " ")), split(spec, " | "))
        return id => (:matchrules, subrules)
    end
end

# Try to match string s[offset:end] with a specific sub-option sequential rule list, return set of remainder offsets
function matchoption(rules, subrules, s; offset)
    offsets = Set(offset)
    for subrule in subrules
        matchsets = [match(rules, subrule, s; offset=o) for o in offsets if o <= lastindex(s)]
        if isempty(matchsets) return Set() end
        offsets = union(matchsets...)
    end
    offsets
end

# Try to match string s[offset:end] with rule ruleidx, return a set of possible remainder offsets
function match(rules, ruleidx, s; offset=1)
    if offset > lastindex(s) return Set() end
    ruletype, data = rules[ruleidx]
    if ruletype == :matchrules
        return union([matchoption(rules, ruleoption, s; offset) for ruleoption in data]...)
    elseif ruletype == :matchchar
        return s[offset] == data ? Set(offset + 1) : Set()
    end
end

rules, messages = split(read("day19_input.txt", String), "\n\n")
rules = Dict(parserule.(split(rules, "\n")))
messages = split(strip(messages), "\n")
println("Part 1: valid = $(count(map(m -> (lastindex(m) + 1) ∈ match(rules, 0, m), messages)))")

# Part 2
for r in ["8: 42 | 42 8", "11: 42 31 | 42 11 31"]
    parsed = parserule(r)
    rules[first(parsed)] = last(parsed)
end
println("Part 2: valid = $(count(map(m -> (lastindex(m) + 1) ∈ match(rules, 0, m), messages)))")

end

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
@safetestset "Responses" begin

  using Genie, HTTP, Genie.Responses

  route("/responses", method = GET) do
    setstatus(301)
    setheaders(Dict("X-Foo-Bar" => "Baz"))
    setheaders(Dict("X-A-B" => "C", "X-Moo" => "Cow"))
    setbody("Hello")
  end

  route("/broken") do
    omg!()
  end

  port = nothing
  port = rand(8500:8900)

  server = up(port)

  response = HTTP.request("GET", "http://localhost:$port/responses")
  @test response.status == 301
  @test Dict(response.headers)["X-Foo-Bar"] == "Baz"
  @test Dict(response.headers)["X-A-B"] == "C"
  @test Dict(response.headers)["X-Moo"] == "Cow"
  @test String(response.body) == "Hello"

  @test_throws HTTP.ExceptionRequest.StatusError HTTP.request("GET", "http://localhost:$port/broken", ["Content-Type"=>"text/plain"])

  down()
  sleep(1)
  server = nothing
  port = nothing
end


text:
@doc """
Check that a matrix is square and column-stochastic.

# Arguments

* `M::AbstractMatrix`: Column-stochastic square matrix.

# Returns

`true` if M is square and column-stochastic, `false` otherwise.

""" ->
function is_square_stochastic{T}(M::AbstractMatrix{T})
    r, c = size(M)
    r == c || return false
    for j = 1:c
        S = zero(T)
        @simd for i = 1:r
            @inbounds S += M[i, j]
        end
        isapprox(S, one(T)) || return false
    end
    return true
end

@doc """
Check that the matrices along the third axis are square and column-stochastic.

# Arguments

* `A::AbstractMatrix`: Column-stochastic square matrices.

# Returns

`true` if A is square and column-stochastic, `false` otherwise.

""" ->
function is_square_stochastic{_}(A::AbstractArray{_,3})
    for k = 1:size(A, 3)
        is_square_stochastic(A[:, :, k]) || return false
    end
    return true
end

function is_square_stochastic(A::Vector)
    for k = 1:length(A)
        is_square_stochastic(A[k]) || return false
    end
    return true
end

function ismdp{_}(transition::AbstractArray{_,3}, reward::AbstractVector)
    # number of rows, columns and actions according to the transition
    # probability matrix
    ST = size(transition, 1)
    # number of states according to reward vector
    SR = length(reward)
    # number of transition states must equal number of reward states
    return (ST == SR) && is_square_stochastic(transition)
end

function ismdp{_}(transition::AbstractArray{_,3}, reward::AbstractMatrix)
    ST, AT = size(transition)[[1, 3]]
    SR, AR = size(reward)
    return (ST == SR) && (AT == AR) && is_square_stochastic(transition)
end

function ismdp{P,R}(transition::AbstractArray{P,3}, reward::AbstractArray{R,3})
    S1T, S2T, AT = size(transition)
    S1R, S2R, AR = size(reward)
    return (S1T == S1R) && (S2T == S2R) && (AT == AR) && is_square_stochastic(transition)
end

_all_row_lengths_equal(len, vec::AbstractVector) = all([len == size(vec[k], 1) for k = 1:length(vec)])

function ismdp(transition::Vector, reward::AbstractVector)
    AT = length(transition)
    SR = length(reward)
    return _all_row_lengths_equal(SR, transition) && is_square_stochastic(transition)
end

function ismdp(transition::Vector, reward::AbstractMatrix)
    AT = length(transition)
    SR, AR = size(reward)
    return (AT == AR) && _all_row_lengths_equal(SR, transition) && is_square_stochastic(transition)
end

function ismdp{_}(transition::Vector, reward::AbstractArray{_,3})
    AT = length(transition)
    S1R, S2R, AR = size(reward)
    return (AT == AR) && (S1R == S2R) && _all_row_lengths_equal(S1R, transition) && is_square_stochastic(transition)
end

ismdp(transition::AbstractTransitionProbability, reward::AbstractReward) =
    ismdp(transition.array, reward.array)

ismdp(mdp::MDP) = ismdp(mdp.transition, mdp.reward)

span(a::AbstractArray) = maximum(a) - minimum(a)


text:
function GenTextures(num::Integer)
	textures = Array(Texture, num)
	ccall( (:glGenTextures, lib), Void, (GLsizei, Ptr{Texture}), num, textures)
	return textures
end

GenTexture() = GenTextures(1)[1]

function BindTexture(target::Integer, texture::Texture)
	ccall( (:glBindTexture, lib), Void, (GLenum, Texture), target, texture)
end
BindTexture(target::Integer, tex::Integer) = BindTexture(target, Texture(tex))

# BindTexture targets
const TEXTURE_1D         = 0x0DE0
const TEXTURE_2D         = 0x0DE1
const TEXTURE_3D         = 0x806F

function TexImage2D{T}(target::Integer, intformat::Integer, width::Integer,
	height::Integer, format::Integer, data::Array{T})

	ccall( (:glTexImage2D, lib), Void,
		(GLenum, GLint, GLint, GLsizei, GLsizei, GLint, GLenum, GLenum, Ptr{GLvoid}),
		target, 0, intformat, width, height, 0, format, GLtype[T], data)
	GetError()
end

# Pixel formats
const RGB  = 0x1907
const RGBA = 0x1908
const RGB8 = 0x8051

function TexParameteri(target::Integer, param::Integer, value::Integer)
	ccall( (:glTexParameteri, lib), Void, (GLenum, GLenum, GLint),
		target, param, value)
end

# Texture parameters
const TEXTURE_DEPTH      = 0x8071
const TEXTURE_MAG_FILTER = 0x2800
const TEXTURE_MIN_FILTER = 0x2801
const TEXTURE_WRAP_S     = 0x2802
const TEXTURE_WRAP_T     = 0x2803
const TEXTURE_WRAP_R     = 0x8072

const NEAREST                = 0x2600
const LINEAR                 = 0x2601
const NEAREST_MIPMAP_NEAREST = 0x2700
const LINEAR_MIPMAP_NEAREST  = 0x2701
const NEAREST_MIPMAP_LINEAR  = 0x2702
const LINEAR_MIPMAP_LINEAR   = 0x2703
const REPEAT                 = 0x2901
const CLAMP_TO_EDGE          = 0x812F

function ActiveTexture(texture::Integer)
	ccall( (:glActiveTexture, lib), Void, (GLenum,), texture)
end

const TEXTURE0  = 0x84C0
const TEXTURE1  = 0x84C1
const TEXTURE2  = 0x84C2
const TEXTURE3  = 0x84C3
const TEXTURE4  = 0x84C4
const TEXTURE5  = 0x84C5
const TEXTURE6  = 0x84C6
const TEXTURE7  = 0x84C7
const TEXTURE8  = 0x84C8
const TEXTURE9  = 0x84C9
const TEXTURE10 = 0x84CA
const TEXTURE11 = 0x84CB
const TEXTURE12 = 0x84CC
const TEXTURE13 = 0x84CD
const TEXTURE14 = 0x84CE
const TEXTURE15 = 0x84CF
const TEXTURE16 = 0x84D0
const TEXTURE17 = 0x84D1
const TEXTURE18 = 0x84D2
const TEXTURE19 = 0x84D3
const TEXTURE20 = 0x84D4
const TEXTURE21 = 0x84D5
const TEXTURE22 = 0x84D6
const TEXTURE23 = 0x84D7
const TEXTURE24 = 0x84D8
const TEXTURE25 = 0x84D9
const TEXTURE26 = 0x84DA
const TEXTURE27 = 0x84DB
const TEXTURE28 = 0x84DC
const TEXTURE29 = 0x84DD
const TEXTURE30 = 0x84DE
const TEXTURE31 = 0x84DF

function GenerateMipmap(target::Integer)
	ccall( (:glGenerateMipmap, lib), Void, (GLenum,), target)
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Composition (Modia3D/Composition/_module.jl)
#
using LinearAlgebra
EYE3() = Matrix(1.0I,3,3)

mutable struct TreeJointRevolute <: Modia3D.AbstractRevolute
   _internal::ModiaMath.ComponentInternal  # Data common to assembly component classes
   obj1::Object3D                        # Rotation of obj1 to obj2 along z-axis of obj1
   obj2::Object3D
   canCollide::Bool                        # = true, if collision detection is enabled between the two Object3Ds connected by the Revolute joint (default = false)
   isDriven::Bool                          # = true, if setAngle!(..) can be called on the joint

   # Definitions special to Revolute joints
   axis::Int            # Rotation axis (axis = 1,2,3 or -1,-2,-3)
   # phi_ref::Float64   # Value of angle phi in reference configuration
   # phi::Float64       # Actual value of relative rotation angle (obj1 and obj2 coincide, if phi = 0)

   J::Float64   # Inertia in the joint

   phi::ModiaMath.RealScalar
   w::ModiaMath.RealScalar
   a::ModiaMath.RealScalar
   tau::ModiaMath.RealScalar
   residue::ModiaMath.RealScalar

   # The following variables are only provided for checking (should be removed, once all tests successful)
   #   P_derx::ModiaMath.RealScalar  # Total power flowing in obj1, obj2, axis computed with derx-vector (must be zero; = computed with error control)
   #   P_w::ModiaMath.RealScalar     # Total power flowing in obj1, obj2, axis computed at communication points (must be zero; = computed with interpolated values of x,derx)
   #   residue_P_derx::ModiaMath.RealScalar # Residue of P_derx

   flange::RevoluteFlange

   TreeJointRevolute(_internal, obj1, obj2, canCollide, isDriven, axis, J) = new(_internal, obj1, obj2, canCollide, isDriven, axis, J)
end

function Base.show(io::IO, joint::TreeJointRevolute)
   print(io,"Revolute(", ModiaMath.fullName(joint.obj1),
                   ", ", ModiaMath.fullName(joint.obj2),
                      ", axis = ", joint.axis,
                      ", phi_start = ", Basics.radToDeg*joint.phi.start, " deg)")
end

#=
function JSON.show_json(io::JSON.StructuralContext,
                         s::JSON.CommonSerialization, joint::TreeJointRevolute)
   JSON.begin_object(io)
      JSON.show_pair(io, s, "obj_a"  , joint.obj1.name)
      JSON.show_pair(io, s, "obj_b"  , joint.obj2.name)
      JSON.show_pair(io, s, "driven"   , joint.driven)
      JSON.show_pair(io, s, "cutJoint" , joint.cutJoint)
      JSON.show_pair(io, s, "axis"     , joint.axis)
      JSON.show_pair(io, s, "phi_start", joint.phi.start)
   JSON.end_object(io)
end
=#

mutable struct CutJointRevolute <: Modia3D.AbstractRevolute
   _internal::ModiaMath.ComponentInternal  # Data common to assembly component classes
   obj1::Object3D                                   # Rotation of obj1 to obj2 along z-axis of obj1
   obj2::Object3D
   # driven::Bool                                    # = true, if setAngle!(..) can be called on the joint
   visited::Bool                                   # = true, if already visited when traversing the spanning tree

   # Definitions special to Revolute joints
   axis::Int            # Rotation axis (axis = 1,2,3 or -1,-2,-3)
   # phi_ref::Float64   # Value of angle phi in reference configuration
   # phi::Float64       # Actual value of relative rotation angle (obj1 and obj2 coincide, if phi = 0)

   lambda_1::ModiaMath.RealScalar
   lambda_2::ModiaMath.RealScalar
   mue_1::ModiaMath.RealScalar
   mue_2::ModiaMath.RealScalar
   residue_1::Union{ModiaMath.RealScalar, NOTHING}
   residue_2::Union{ModiaMath.RealScalar, NOTHING}
   residue_v1::Union{ModiaMath.RealScalar, NOTHING}
   residue_v2::Union{ModiaMath.RealScalar, NOTHING}

   CutJointRevolute(_internal, obj1, obj2, visited, axis) = new(_internal, obj1, obj2, visited, axis)
end

"""
    joint = Modia3D.Revolute(obj1, obj2; axis=3, phi_start=0, w_start=0, canCollide=false, J=NaN)

Return a Revolute `joint` that rotates `obj1::`[`Object3D`](@ref) into
`obj2::`[`Object3D`](@ref) along the axis `axis` of `obj1` (`axis = 1,2,3,-1,-2,-3`).
The initial start angle is `phi_start` and the initial angular velocity
is `w_start`. If `canCollide=false`, no collision detection will occur between `obj1` and `obj2`
(and `Object3D`s that are directly or indirectly rigidly fixed to `obj1` or `obj2`).

Optionally, an inertia can be placed on the axis of rotation (`J`).
If such an inertia is provided, the generalized coordinate is still the angle of the revolute joint and
the driving torque is still the torque acting in the revolute joint. If, for example, an electric
motor with a gearbox shall be modelled then

```julia
J = JMotor*gearRatio^2
jointTorque = motorTorque*gearRatio
motor
```
```J = JMotor*gearRatio^2; jointTorque = `.

If a `Revolute` joint *closes a kinematic loop*, then the already present objects must be consistent
to the `Revolute` joint. That is, the frames of `obj1` and `obj2` must be *parallel* to each other and
rotating `obj1` along its axis `axis` with `phi_start` results in `obj2`. If `phi_start=NaN`,
its value is computed in this case.

# Examples
```julia
using Modia3D
import ModiaMath

@assembly Pendulum(;L=1.0) begin
   world  = Modia3D.Object3D()
   obj1 = Modia3D.Object3D()
   rev    = Modia3D.Revolute(world, obj1)
   sphere = Modia3D.Object3D(obj1, Modia3D.Solid(Modia3D.SolidSphere(0.1), "Aluminium"),
                             r = [L,0,0] )
end

simulationModel = Modia3D.SimulationModel( Pendulum(L=0.8), stopTime=5.0 )
result          = ModiaMath.simulate!(simulationModel)
ModiaMath.plot(result, ("rev.phi", "rev.w"))
```
"""
function Revolute(obj1::Object3D, obj2::Object3D;
                  phi_start::Number = 0.0,
                  w_start::Number   = 0.0,
                  canCollide::Bool  = false,
                  axis::Int = 3,
                  J=0.0)
   rphi_start = Float64(phi_start)
   rw_start   = Float64(w_start)
   (obj_a,obj_b,cutJoint) = attach(obj1, obj2)
   @assert(1 <= abs(axis) <= 3)
   axis  = obj_b===obj2 ? axis : -axis
   if cutJoint
      println("... Revolute joint connecting ", ModiaMath.fullName(obj1), " with ", ModiaMath.fullName(obj2), " is a cut-joint")

      # Check that rotating from obj_a and obj_b have the same origin and
      # rotating from obj_a along z-axis with angle phi gives obj_b (or compute rphi_start if phi_start = NaN)
      diff_r_abs = norm(obj_b.r_abs - obj_a.r_abs)
      if diff_r_abs > Basics.neps
             @warn begin
                 name_obj1 = ModiaMath.fullName(obj1)
                 name_obj2 = ModiaMath.fullName(obj2)
                 "Warning from Modia3D.Revolute($name_obj1, $name_obj2, ....):\n" *
                 "The origins of the provided frames do not coincide.\n" *
                 "( norm(obj2.r_abs - obj1.r_abs) = $diff_r_abs > $(Basics.neps) )"
             end
      end
      if isnan(rphi_start)
         # Compute angle between the two axes
         # @error("brrrr!!!!!!!!!!!!!!!!!!")
         angle = acos( dot(obj_a.R_abs[abs(axis),:], obj_b.R_abs[abs(axis),:]) )
         if abs(angle) > Basics.neps
            @warn begin
                name_obj1 = ModiaMath.fullName(obj1)
                name_obj2 = ModiaMath.fullName(obj2)
                if abs(axis) == 3
                   "Warning from Modia3D.Revolute($name_obj1, $name_obj2, ....):\n" *
                   "The z-axes of the two frames do not coincide.\n" *
                   "( angle(zaxis(obj1), zaxis(obj2)) = $zaxes_angle > $(Basics.neps) )"
                elseif abs(axis) == 2
                   "Warning from Modia3D.Revolute($name_obj1, $name_obj2, ....):\n" *
                   "The y-axes of the two frames do not coincide.\n" *
                   "( angle(yaxis(obj1), yaxis(obj2)) = $yaxes_angle > $(Basics.neps) )"
                elseif abs(axis) == 1
                   "Warning from Modia3D.Revolute($name_obj1, $name_obj2, ....):\n" *
                   "The x-axes of the two frames do not coincide.\n" *
                   "( angle(xaxis(obj1), xaxis(obj2)) = $xaxes_angle > $(Basics.neps) )"
                end
            end
         end
      else
         R2_abs     = rotationMatrixAxis(axis, axis > 0 ? rphi_start : -rphi_start)*obj_a.R_abs
         diff_R_abs = norm(R2_abs*obj_b.R_abs' - ModiaMath.NullRotation)
         if diff_R_abs > Basics.neps
                 @warn begin
                     name_obj1 = ModiaMath.fullName(obj1)
                     name_obj2 = ModiaMath.fullName(obj2)
                     "Warning from Modia3D.Revolute($name_obj1, $name_obj2, ....):\n" *
                     "The orientations of the provided frames do not coincide.\n" *
                     "( norm(R2_abs*obj2.R_abs' - EYE3()) = $diff_R_abs > $(Basics.neps) )"
                 end
         end
      end

      # Define variables of cut-joint
      joint      = CutJointRevolute(ModiaMath.ComponentInternal(), obj_a, obj_b, false, axis)
      lambda_1   = ModiaMath.RealScalar("lambda_1"  , joint, unit="N", info="Constraint force in x-direction"   , numericType=ModiaMath.LAMBDA, analysis=ModiaMath.QuasiStaticAndDynamicAnalysis)
      lambda_2   = ModiaMath.RealScalar("lambda_2"  , joint, unit="N", info="Constraint force in y-direction"   , numericType=ModiaMath.LAMBDA, analysis=ModiaMath.QuasiStaticAndDynamicAnalysis)
      mue_1      = ModiaMath.RealScalar("mue_1"     , joint,           info="Stabilizing mue in x-direction"    , numericType=ModiaMath.MUE   , analysis=ModiaMath.OnlyDynamicAnalysis)
      mue_2      = ModiaMath.RealScalar("mue_2"     , joint,           info="Stabilizing mue in y-direction"    , numericType=ModiaMath.MUE   , analysis=ModiaMath.OnlyDynamicAnalysis)
      residue_1  = ModiaMath.RealScalar("residue_1" , joint,           info="Position constraint in x-direction", numericType=ModiaMath.FC)
      residue_2  = ModiaMath.RealScalar("residue_2" , joint,           info="Position constraint in y-direction", numericType=ModiaMath.FC)
      residue_v1 = ModiaMath.RealScalar("residue_v1", joint,           info="Velocity constraint in x-direction", numericType=ModiaMath.FC    , analysis=ModiaMath.OnlyDynamicAnalysis)
      residue_v2 = ModiaMath.RealScalar("residue_v2", joint,           info="Velocity constraint in y-direction", numericType=ModiaMath.FC    , analysis=ModiaMath.OnlyDynamicAnalysis)

      # Add joint to the frames
      push!(obj_a.twoObject3Dobject, joint)
      push!(obj_b.twoObject3Dobject, joint)
      obj_a.hasCutJoint = true
      obj_b.hasCutJoint = true
   else
      # Joint in the spanning tree
      if isnan(rphi_start)
         error("\nerror from Modia3D.Revolute(", ModiaMath.fullName(obj1), ", ", ModiaMath.fullName(obj2), ", phi_start=NaN):\n",
              "Argument phi_start = NaN (Not-a-Number). However, joint is not allowed if the joint is in the spanning tree")
      end
      joint   = TreeJointRevolute(ModiaMath.ComponentInternal(), obj_a, obj_b, canCollide, false, axis, J)
      phi     = ModiaMath.RealScalar("phi"    , joint, unit="rad"    , start=rphi_start, fixed=true, info="Relative rotation angle"      , numericType=ModiaMath.XD_EXP)
      w       = ModiaMath.RealScalar("w"      , joint, unit="rad/s"  , start=rw_start  , fixed=true, info="Relative angular velocity"    , numericType=ModiaMath.XD_IMP    , integral=phi, analysis=ModiaMath.OnlyDynamicAnalysis)
      a       = ModiaMath.RealScalar("a"      , joint, unit="rad/s^2", start=0.0       ,             info="Relative angular acceleration", numericType=ModiaMath.DER_XD_IMP, integral=w  , analysis=ModiaMath.OnlyDynamicAnalysis)
      tau     = ModiaMath.RealScalar("tau"    , joint, unit="Nm"     , start=0.0       ,             info="Driving torque"               , numericType=ModiaMath.WR                      , analysis=ModiaMath.QuasiStaticAndDynamicAnalysis)
      residue = ModiaMath.RealScalar("residue", joint,                 start=0.0       ,             info="Revolute residue"             , numericType=ModiaMath.FD_IMP                  , analysis=ModiaMath.OnlyDynamicAnalysis)
      #P_w     = ModiaMath.RealScalar("P_w"    , joint, unit="J", info="Total power computed with interpolated x,derx, values (should be zero)", numericType=ModiaMath.WC, analysis=ModiaMath.OnlyDynamicAnalysis)

      joint.flange = RevoluteFlange()

      obj_a.hasChildJoint = true
      obj_b.joint = joint
      obj_b.r_rel = ModiaMath.ZeroVector3D
      obj_b.r_abs = obj_a.r_abs
      obj_b.R_rel = rotationMatrixAxis(axis, axis > 0 ? rphi_start : -rphi_start)
   end
   return joint
end

function rotationMatrixAxis(axis::Int, phi::Number)
   if abs(axis) == 1
      return ModiaMath.rot1(phi)
   elseif abs(axis) == 2
      return ModiaMath.rot2(phi)
   elseif abs(axis) == 3
      return ModiaMath.rot3(phi)
   end
end

function driveJoint!(joint::TreeJointRevolute)
   joint.isDriven = true
   joint.phi.numericType = ModiaMath.WR
   joint.w.numericType   = ModiaMath.WR
   joint.a.numericType   = ModiaMath.WR
end

function setAngle!(joint::TreeJointRevolute, phi::Float64)
   #@assert( joint.isDriven )
   joint.phi.value = phi
   joint.obj2.R_rel = rotationMatrixAxis(joint.axis, joint.axis > 0 ? phi : -phi)
   return joint
end

function computeKinematics!(joint::TreeJointRevolute, obj::Object3D, analysis::ModiaMath.AnalysisType, time::Float64)::NOTHING
   parent::Object3D = obj.parent

   obj.r_abs = parent.r_abs
   obj.R_rel = rotationMatrixAxis(joint.axis, joint.axis > 0 ? joint.phi.value : -joint.phi.value)
   obj.R_abs = obj.R_rel*parent.R_abs

   if analysis == ModiaMath.DynamicAnalysis
      dynamics::Object3Ddynamics       = obj.dynamics
      parentDynamics::Object3Ddynamics = parent.dynamics

      dynamics.v0 = parentDynamics.v0
      dynamics.a0 = parentDynamics.a0

      w_rel = MVector{3,Float64}(0.0, 0.0, 0.0)
      z_rel = MVector{3,Float64}(0.0, 0.0, 0.0)

      w_rel[abs(joint.axis)] = joint.axis > 0 ? joint.w.value : -joint.w.value
      z_rel[abs(joint.axis)] = joint.axis > 0 ? joint.a.value : -joint.a.value

      dynamics.w = obj.R_rel*(parentDynamics.w + w_rel)
      dynamics.z = obj.R_rel*(parentDynamics.z + z_rel + cross(parentDynamics.w, w_rel))
   end
   return nothing
end

function computeForceTorqueAndResidue!(joint::TreeJointRevolute, obj::Object3D, analysis::ModiaMath.AnalysisType, time::Float64)::NOTHING
   parent::Object3D                 = obj.parent
   dynamics::Object3Ddynamics       = obj.dynamics
   parentDynamics::Object3Ddynamics = parent.dynamics
   parentDynamics.f    += obj.R_rel'*dynamics.f
   parentDynamics.t    += obj.R_rel'*dynamics.t
   dynamics.f          = -dynamics.f
   dynamics.t          = -dynamics.t

   # println("joint.tau.value = ", joint.tau.value)
   joint.residue.value = -joint.tau.value + joint.J*joint.a.value + (joint.axis > 0 ? dynamics.t[abs(joint.axis)] : -dynamics.t[abs(joint.axis)])
   #println("nachher joint.residue.value = ", joint.residue.value)

   #println(" ")
   # For checking: compute total power
   # P = dot(parentDynamics.f, joint.obj1.R_abs*parentDynamics.v0) +
   #    dot(dynamics.f      , joint.obj2.R_abs*dynamics.v0) +
   #    dot(parentDynamics.t, parentDynamics.w) +
   #    dot(dynamics.t      , dynamics.w) -
   #    joint.tau.value*joint.w.value
   #
   # joint.P_w.value = P

   return nothing
end

function computePositionResidues!(joint::CutJointRevolute, time::Float64)
   r_rel12 = joint.obj1.R_abs*(joint.obj2.r_abs - joint.obj1.r_abs)
   if abs(joint.axis) == 3
      joint.residue_1.value = r_rel12[1]
      joint.residue_2.value = r_rel12[2]
   elseif abs(joint.axis) == 2
      joint.residue_1.value = r_rel12[1]
      joint.residue_2.value = r_rel12[3]
   elseif abs(joint.axis) == 1
      joint.residue_1.value = r_rel12[2]
      joint.residue_2.value = r_rel12[3]
   end
   return nothing
end

function computeVelocityResidues!(joint::CutJointRevolute, time::Float64)
   v_rel12 = joint.obj1.R_abs*(joint.obj2.dynamics.v0 - joint.obj1.dynamics.v0)
   if abs(joint.axis) == 3
      joint.residue_v1.value = v_rel12[1]
      joint.residue_v2.value = v_rel12[2]
   elseif abs(joint.axis) == 2
      joint.residue_v1.value = v_rel12[1]
      joint.residue_v2.value = v_rel12[3]
   elseif abs(joint.axis) == 1
      joint.residue_v1.value = v_rel12[2]
      joint.residue_v2.value = v_rel12[3]
   end
   return nothing
end

function computeCutForcesAndToques!(joint::CutJointRevolute, time::Float64)
   if abs(joint.axis) == 3
      dynamics1.f += SVector{3,Float64}(joint.lambda_1, joint.lambda_2, 0.0)
   elseif abs(joint.axis) == 2
      dynamics1.f += SVector{3,Float64}(joint.lambda_1, 0.0, joint.lambda_2)
   elseif abs(joint.axis) == 1
      dynamics1.f += SVector{3,Float64}(0.0, joint.lambda_1, joint.lambda_2)
   end
   dynamics2.f += joint.obj2.R_abs*(joint.obj1.R_abs'*dynamics1.f)
end

function addToFlowOrPotentialArray(variable::ModiaMath.RealScalar, assembly::Modia3D.Composition.AssemblyInternal)
  if variable.flow
    if variable.causality == ModiaMath.Input
      push!(assembly.flowVarInput, variable)
    elseif variable.causality == ModiaMath.Output
      push!(assembly.flowVarOutput, variable)
    end
  else
    if variable.causality == ModiaMath.Input
      push!(assembly.potentialVarInput, variable)
    elseif variable.causality == ModiaMath.Output
      push!(assembly.potentialVarOutput, variable)
    end
  end
end

function checkLengthArrays(assembly::Modia3D.Composition.AssemblyInternal)
  @assert(length(assembly.potentialVarOutput) == length(assembly.potentialVarInput))
  @assert(length(assembly.flowVarOutput) == length(assembly.flowVarInput))
end

function addVariableToArray(flange_val::ModiaMath.RealScalar, joint_val::ModiaMath.RealScalar, assembly::Modia3D.Composition.AssemblyInternal)
  # println("flange_val.analysis = ", flange_val.analysis, " revFlange_val.analysis = ", revFlange_val.analysis, " joint_val.analysis = ", joint_val.analysis)
  if flange_val.causality != ModiaMath.Local && joint_val.causality != ModiaMath.Local
    addToFlowOrPotentialArray(flange_val, assembly)
    addToFlowOrPotentialArray(joint_val, assembly)
    checkLengthArrays(assembly)
  end
end

function assignRevoluteFlange(flange::Flange, revFlange::RevoluteFlange, joint::Modia3D.AbstractRevolute, assembly::Modia3D.Composition.AssemblyInternal)
  names = fieldnames(typeof(flange))
  for val in names # loops over phi, w, a, tau
    flange_val    = getfield(flange, val)
    revFlange_val = getfield(revFlange,val)
    joint_val     = getfield(joint,val)
    if causalitiesAreOk(flange_val, joint_val)
        setFlangeVariable!(flange_val, revFlange_val)
        setJointVariable!(joint_val, revFlange_val)
        addVariableToArray(flange_val, joint_val, assembly)
    else
      return false
    end
  end
  return true
end

#=
function addTorqueObjToJoint(torqueObj::Modia3D.AbstractForceTorque, SymbFlangeTorque::Symbol, joint::Modia3D.AbstractRevolute, SymbFlangeJoint::Symbol, assembly::Modia3D.AbstractAssembly)
  flangeTorque  = getfield(torqueObj, SymbFlangeTorque)
  flangeJoint   = getfield(joint, SymbFlangeJoint)
  if !assignRevoluteFlange(flangeTorque, flangeJoint, joint, assembly)
    error("The causalities are different. Please, check all flanges which are connected with joint = ", joint ,".")
  end
  if !isPhiInput(flangeJoint) && isTauInput(flangeJoint)
    addRevoluteTorqueObject(joint, torqueObj, SymbFlangeTorque)
    joint.isDriven = false
  elseif isPhiInput(flangeJoint) && !isTauInput(flangeJoint)
    driveJoint!(joint)
  end
end
=#

function connect(torqueObj::Modia3D.AbstractForceAdaptor, joint::Modia3D.AbstractRevolute)
  actualAssembly = torqueObj._internal.within._internal

  flangeTorque  = torqueObj.flange
  flangeJoint   = joint.flange

  if !assignRevoluteFlange(flangeTorque, flangeJoint, joint, actualAssembly) # !assignRevoluteFlangeAndSignal(signal, flangeSig, joint, flangeJoint, assembly)
    error("The causalities are different. Please, check all flanges which are connected with joint = ", joint ,".")
  end
  if !in(torqueObj.forceElement, actualAssembly.uniqueForceTorques)
    push!(actualAssembly.uniqueForceTorques, torqueObj.forceElement)
  end
  if !isPhiInput(flangeJoint) && isTauInput(flangeJoint)
    joint.isDriven = false
  elseif isPhiInput(flangeJoint) && !isTauInput(flangeJoint)
    driveJoint!(joint)
  end
end
connect(joint::Modia3D.AbstractRevolute, torqueObj::Modia3D.AbstractForceAdaptor) = connect(torqueObj::Modia3D.AbstractForceAdaptor, joint::Modia3D.AbstractRevolute)

computeTorque(obj::Any, phi::Float64, w::Float64, time::Float64) = error("Modia3D.Revolute: Missing implementation for function Modia3D.computeTorque")

function connect(signal::Modia3D.AbstractSignalAdaptor, joint::Modia3D.AbstractRevolute)
    actualAssembly = signal._internal.within._internal

    flangeSig   = signal.flange
    flangeJoint = joint.flange

    if !assignRevoluteFlange(flangeSig, flangeJoint, joint, actualAssembly) # !assignRevoluteFlangeAndSignal(signal, flangeSig, joint, flangeJoint, assembly)
      error("The causalities are different. Please, check all flanges which are connected with joint = ", joint ,".")
    end
    if !in(signal.signal, actualAssembly.uniqueSignals)
      push!(actualAssembly.uniqueSignals, signal.signal)
    end
    if !isPhiInput(flangeJoint) && isTauInput(flangeJoint)
      joint.isDriven = false
    elseif isPhiInput(flangeJoint) && !isTauInput(flangeJoint)
      driveJoint!(joint)
    end
end
connect(joint::Modia3D.AbstractRevolute, signal::Modia3D.AbstractSignalAdaptor) = connect(signal::Modia3D.AbstractSignalAdaptor, joint::Modia3D.AbstractRevolute)

computeSignal(obj::Any, time::Float64, tStart::Float64) = error("Modia3D.Revolute: Missing implementation for function Modia3D.computeSignal.")


text:
abstract type AbstractConditionalDistribution end
const ACD = AbstractConditionalDistribution

Distributions.mean(p::ACD, z::AbstractVecOrMat) = mean(condition(p,z))
Distributions.var(p::ACD, z::AbstractVecOrMat) = var(condition(p,z))
Distributions.rand(p::ACD, z::AbstractVecOrMat) = rand(condition(p,z))
Distributions.logpdf(p::ACD, x::AbstractVecOrMat, z::AbstractVecOrMat) = logpdf(condition(p,z), x)

"""
    ConditionalDistribution(T, m)

A conditional distribution that can map examples `x` to a distribution of type
`T` via the mapping `m` by calling `condition`.
There are no type restrictions on `T` so it could also be a function that
constructs and appropriate distribution (e.g. like TuringMvNormal).

# Examples
Moments, likelihoods, and samples can be obtained via e.g. `mean`, `logpdf`,
and `rand`:
```julia-repl
julia> using ConditionalDists, Flux
julia> m = ConditionalDists.SplitLayer(5,[2,2])
julia> p = ConditionalDistribution(MvNormal, m)
julia> rand(p, rand(5))
2-element Array{Float32,1}:
 0.2433714
 2.2054431
```
"""
struct ConditionalDistribution{Td,Tm} <: AbstractConditionalDistribution
    # TODO: would be nice if it was restricted to a `Distribution`
    DistConstructor::Td
    mapping::Tm
end

"""
    condition(p::ConditionalDistribution, x::AbstractVector)

Map an example `x` to the distribution type specified in `p`.
"""
function condition(p::ConditionalDistribution, x::AbstractVector)
    p.DistConstructor(p.mapping(x)...)
end

"""
    condition(p::ConditionalDistribution, x::AbstractMatrix)

Map a batch of examples `x` (examples as columns of `x`) to the distribution
type specified in `p`.
"""
function condition(p::ConditionalDistribution, xs::AbstractMatrix)
    ds = map(i -> condition(p, view(xs, :, i)), 1:size(xs,2))
    arraydist(ds)
end

function Base.show(io::IO, p::ConditionalDistribution)
    Np = nameof(typeof(p))
    Nd = nameof(p.DistConstructor)
    Nm = repr(p.mapping)
    print(io, "$Np($Nd, $Nm)")
end

@functor ConditionalDistribution


text:
"""
  logseriescdf(p::Float64)

Returns a vector{Float} of the discrete cdf of logarithmic distribution
"""
function logseriescdf(p::Float64)
  cdfs = [0.]
  for i in 1:100000000
    @inbounds push!(cdfs, cdfs[i]-(p^i)/(i*log(1-p)))
    if cdfs[i] ≈ 1.0
      return cdfs
    end
  end
  cdfs
end

"""
  logseriesquantile(p::Float64, v::Vector{Float64})

Returns a vector{Float} of the v[i] th quaintlie of logaritmic distribution with parameter p
"""
function logseriesquantile(p::Float64, v::Vector{Float64})
  w = logseriescdf(p)
  pmap(i -> findlast(w .< i), v)
end

"""
  levyel(θ::Union{Int, Float64})

An element from Levy stable distribution with parameters α = 1/θ, β = 1,
γ = (cos(pi/(2*θ)))^θ and δ = 0.
Return Float, given parameter ϴ of dostribution
"""
function levyel(θ::Union{Int, Float64})
  ϕ = pi*rand()-pi/2
  v = quantile.(Exponential(1.), rand())
  γ = (cos(pi/(2*θ)))^θ
  v = ((cos(pi/(2*θ)+(1/θ-1)*ϕ))/v)^(θ-1)
  γ*v*sin(1/θ*(pi/2+ϕ))*(cos(pi/(2*θ))*cos(ϕ))^(-θ)
end

"""
Return a Vector{Float64} of  of pseudo cdf of Levy stable distribution with parameters
α = 1/θ, β = 1, γ = (cos(pi/(2*θ)))^θ and δ = 0, given a vector of Float64 - u

"""
function levygen(θ::Union{Int, Float64}, u::Vector{Float64})
  p = invperm(sortperm(u))
  v = [levyel(θ) for a in u]
  sort(v)[p]
end

"""
  tiltedlevygen(V0::Vector{Float64}, α::Float64)

Returns a Vector{Float} genrated from the expotencialy tilted levy stable pdf
f(x; V0, α) = exp(-V0^α) g(x; α)/exp(-V0), where g(x; α) is a stable Levy pdf
with parameters α = 1/θ, β = 1, γ = (cos(pi/(2*θ)))^θ and δ = 0.

"""
function tiltedlevygen(V0::Vector{Float64}, α::Float64)
  t = length(V0)
  ret = zeros(t)
  for i in 1:t
    x = levyel(α)
    u = rand()
    while exp(-V0[i]^α*x)/(1500*exp(-V0[i])) < u
      x = levyel(α)
      u = rand()
    end
    ret[i] = x
  end
  ret.*V0.^α
end

"""
  Ginv(y::Float64, α::Float64)

Returns Float64, helper for the joe/frank nested copula generator
"""
Ginv(y::Float64, α::Float64) = ((1-y)*SpecialFunctions.gamma(1-α))^(-1/α)

"""
  InvlaJ(n::Int, α::Float64)

Returns Float64, n-th element of the inverse laplacea transform of generator of Joe nested copula
"""
InvlaJ(n::Int, α::Float64) = 1-1/(n*beta(n, 1-α))

"""
  sampleInvlaJ(α::Float64, v::Float64)

Returns Int, a sample of inverce laplacea transform of generator of Joe nested copula,
given a parameter α and a random numver v ∈ [0,1]
"""
function sampleInvlaJ(α::Float64, v::Float64)
  if v <= α
    return 1
  else
    G = Ginv(v, α)
    if G > 2^62
      return 2^62
    else
      return (InvlaJ(floor(Int, G), α) < v) ? ceil(Int, G) : floor(Int, G)
    end
  end
end

"""
  elInvlaF(θ₁::Float64, θ₀::Float64)

Returns Int, a single sample of the inverse laplacea transform of the generator
of nested Frank copula
"""
function elInvlaF(θ₁::Float64, θ₀::Float64)
  c1 = 1-exp(-θ₁)
  α = θ₀/θ₁
  if θ₀ <= 1
    v = rand()
    X = logseriesquantile(c1, rand(1))[1]
    while v > 1/((X-α)*beta(X, 1-α))
      v = rand()
      X = logseriesquantile(c1, rand(1))[1]
    end
    return X
  else
    v = rand()
    X = sampleInvlaJ(α, rand())
    while v > c1^(X-1)
      X = sampleInvlaJ(α, rand())
      v = rand()
    end
    return X
  end
end

"""
  nestedfrankgen(θ₁::Float64, θ₀::Float64, V0::Vector{Int})

Return vector of int, samples of inverse laplacea trensform of nested
Frank copula given parametes and V0 - vector of samples if invlaplace of perents copula
"""
function nestedfrankgen(θ₁::Float64, θ₀::Float64, V0::Vector{Int})
  if nprocs() == 1
    return map(k -> sum([elInvlaF(θ₁, θ₀) for j in 1:k]), V0)
  end
  u = SharedArray{Float64}(length(V0))
  @sync @distributed for i = 1:length(V0)
    u[i] = sum([elInvlaF(θ₁, θ₀) for j in 1:V0[i]])
  end
  Array{Int}(u)
end


text:
using VectorizationBase: REGISTER_SIZE

const float_to_float = [
    # "abs", "fabs",
    "exp", "exp2", "exp10", "expm1", "log", "log2", "log10", "log1p",
    "sin", "cos", "tan",  "asin", "acos", "atan",
    "sinh", "cosh", "tanh", "asinh", "acosh", "atanh",
    "erf", "erfc", "tgamma", "lgamma"
]
const floatfloat_to_float = [
    # "fmod", "remainder", "min", "max", "fmin", "fmax", "fdim",
    "pow", "hypot", "atan2"
]
const float_to_floatfloat = [
    "sincos"
]

function gen_cpp_code()
    bits = (128, 256, 512)
    native_bits = REGISTER_SIZE << 3
    cxx_code = """
    #include "xsimd/xsimd.hpp"
    #include <immintrin.h>

    namespace xs = xsimd;

    extern "C" {
    """
    for f ∈ float_to_float
        for bitcount ∈ bits
            bitcount > native_bits && continue
            cxx_code *= """
            __m$(bitcount)d xsimd_$(f)_$(bitcount)d(const __m$(bitcount)d a){
                const auto b = xs::batch<double,$(bitcount >> 6)>(a);
                auto fb = xs::$(f)(b);
                return fb;
            }
            __m$(bitcount) xsimd_$(f)_$(bitcount)(const __m$(bitcount) a){
                const auto b = xs::batch<float,$(bitcount >> 5)>(a);
                auto fb = xs::$(f)(b);
                return fb;
            }
            """
        end
    end
    for f ∈ floatfloat_to_float
        for bitcount ∈ bits
            bitcount > native_bits && continue
            cxx_code *= """
            __m$(bitcount)d xsimd_$(f)_$(bitcount)d(const __m$(bitcount)d a, const __m$(bitcount)d c){
                const auto b = xs::batch<double,$(bitcount >> 6)>(a);
                const auto d = xs::batch<double,$(bitcount >> 6)>(c);
                auto fbd = xs::$(f)(b, d);
                return fbd;
            }
            __m$(bitcount) xsimd_$(f)_$(bitcount)(const __m$(bitcount) a, const __m$(bitcount) c){
                const auto b = xs::batch<float,$(bitcount >> 5)>(a);
                const auto d = xs::batch<float,$(bitcount >> 5)>(c);
                auto fbd = xs::$(f)(b, d);
                return fbd;
            }
            """
        end
    end
    for f ∈ float_to_floatfloat
        for bitcount ∈ bits
            bitcount > native_bits && continue
            cxx_code *= """
            std::pair<__m$(bitcount)d,__m$(bitcount)d> xsimd_$(f)_$(bitcount)d(const __m$(bitcount)d a){
                const auto b = xs::batch<double,$(bitcount >> 6)>(a);
                xs::batch<double,$(bitcount >> 6)> sb, cb;
                xs::$(f)(b, sb, cb);
                __m$(bitcount)d s = sb;
                __m$(bitcount)d c = cb;
                return std::make_pair(s,c);
            }
            std::pair<__m$(bitcount),__m$(bitcount)> xsimd_$(f)_$(bitcount)(const __m$(bitcount) a){
                const auto b = xs::batch<float,$(bitcount >> 5)>(a);
                xs::batch<float,$(bitcount >> 5)> sb, cb;
                xs::$(f)(b, sb, cb);
                __m$(bitcount) s = sb;
                __m$(bitcount) c = cb;
                return std::make_pair(s,c);
            }
            """
        end
    end
    cxx_code *= "\n}"
    open("xsimdmathfunctions.cpp", "w") do f
        write(f, cxx_code)
    end
end

if isdir("xsimd")
    cd("xsimd")
    run(`git pull`)
    cd("..")
else
    run(`git clone https://github.com/QuantStack/xsimd`)
end

gen_cpp_code()
run(`$(get(ENV,"CXX","g++")) -I./xsimd/include $(split(get(ENV,"CXXFLAGS",""))) -fno-semantic-interposition -march=native -Ofast -shared -fPIC xsimdmathfunctions.cpp -o libxsimdmath.so`)


text:
export findallmax, deletefirst!, onehot,
       discounted_reward, reverse_discounted_rewards,
       importance_weight, reverse_importance_weights, huber_loss

using SparseArrays
using StatsBase:mean

"""
    findallmax(A::AbstractArray)

Like `findmax`, but all the indices of the maximum value are returned.

!!! warning
    All elements of value `NaN` in `A` will be ignored, unless all elements are `NaN`.
    In that case, the returned maximum value will be `NaN` and the returned indices will be `collect(1:length(A))`

#Examples
```julia-repl
julia> findallmax([-Inf, -Inf, -Inf])
(-Inf, [1, 2, 3])

julia> findallmax([Inf, Inf, Inf])
(Inf, [1, 2, 3])

julia> findallmax([Inf, 0, Inf])
(Inf, [1, 3])

julia> findallmax([0,1,2,1,2,1,0])
(2, [3, 5])
```
"""
function findallmax(A)
    maxval = -Inf
    idxs = Vector{Int}()
    for (i, x) in enumerate(A)
        if !isnan(x)
            if x > maxval
                maxval = x
                empty!(idxs)
                push!(idxs, i)
            elseif x == maxval
                push!(idxs, i)
            end
        end
    end
    if length(idxs) == 0
        NaN, collect(1:length(A))
    else
        maxval, idxs
    end
end

"""
    deletefirst!(A::Vector, element)

Find the first `element` in `A` and delete it.
`==` is used to compare equality.
"""
function deletefirst!(A::Vector, element)
    deleteat!(A, findfirst(x->x == element, A))
end

"""
    onehot(n::Int, x::Int, t::Type=Int; isdense::Bool=true)

If `isdense` is `false`, a `SparseArray` is returned.
"""
function onehot(n::Int, x::Int, t::Type = Int; isdense::Bool = true)
    a = isdense ? zeros(t, n) : spzeros(t, n)
    a[x] = one(t)
    a
end

discounted_reward(xs, γ) = foldr((r, g)->r + γ * g, xs)

"""
    reverse_discounted_rewards(rewards, γ)

Given the `rewards` and discount ratio `γ`, the discounted reward at each time step in the reversed order is returned.
The returned object is of type [`Reductions`](@ref).
"""
reverse_discounted_rewards(xs, γ) = Reductions((g, r)->g * γ + r, Iterators.reverse(xs))

"""
    importance_weight(π, b, states, actions)

Calculate the importance weight between the target policy `π` and behavior policy `b`
given `states` and `actions`.
"""
function importance_weight(π, b, states, actions)
    ρ = 1.
    for (s, a) in zip(states, actions)
        if π(s, a) == 0
            ρ = 0.
            break
        else
            ρ *= π(s, a) / b(s, a)
        end
    end
    ρ
end

"""
    reverse_importance_weights(π, b, states, actions)

Calculate the importance weight at each time step in the **reversed** order
between the target policy `π` and behavior policy `b` given `states` and `actions`.

The returned object is of type [`Reductions`](@ref)
"""
reverse_importance_weights(π, b, states, actions) = Reductions(
    (ρ, (s, a)) -> ρ == 0. ? 0. : ρ * π(s, a) / b(s, a),
    Iterators.reverse(zip(states, actions)),
    (init=1.,))

const is_using_gpu = false

"""
    huber_loss(labels, predictions;δ = 1.0)

See [huber_loss](https://en.m.wikipedia.org/wiki/Huber_loss)
and the [implementation](https://github.com/tensorflow/tensorflow/blob/r1.12/tensorflow/python/ops/losses/losses_impl.py#L394-L469) in TensorFlow.
"""
function huber_loss(labels, predictions;δ = 1.0)
    abs_error = abs.(predictions .- labels)
    quadratic = min.(abs_error, δ)
    linear = abs_error .- quadratic
    mean(0.5 .* quadratic .* quadratic .+ δ .* linear)
end


text:
"""
This file contains all the back-end scripts in ClimateSatellite.jl that are for the download
and extraction of precipitable water data from the Morphed Integrated Microwave Imagery at
CIMSS (MIMIC)-TPW2m product:
    * retrieval of file names
    * downloading data from ftp server
    * extraction of chosen regions from global dataset

"""

function mimiclonlat()
    lon = convert(Array,-180:0.25:180); lat = convert(Array,-90:0.25:90); pop!(lon);
    return lon,lat
end

function mimicnclist(date::TimeType)

    yr = Dates.year(date); mo = Dates.month(date); ndy = daysinmonth(date);
    fname = Array{AbstractString,2}(undef,24,ndy);

    @debug "$(Dates.now()) - Creating list of data files to download ..."
    for dy = 1 : ndy; date = Date(yr,mo,dy);
        for hh = 1 : 24

            fname[hh,dy] = "comp$(ymd2str(date)).$(@sprintf("%02d",hh-1))0000.nc";

        end
    end

    furl = "ftp://ftp.ssec.wisc.edu/pub/mtpw2/data/$(yrmo2str(date))"

    return fname,furl

end

function mimicget(
    url::AbstractString,
    file::AbstractString,
    tdir::AbstractString,
    overwrite::Bool
)
    if overwrite || !isfile(joinpath(tdir,file))
        try download("$(url)/$(file)",joinpath(tdir,file));
            @debug "$(Dates.now()) - Downloaded MIMIC tropospheric precipitable water data file $(file)"
        catch;
            @warn "$(Dates.now()) - MIMIC tropospheric precipitable water data $(file) does not exist."
        end
    end

end

function mimicretrieve(
    fname::Array{<:AbstractString,2},
    furl::AbstractString, date::TimeType,
    tdir::AbstractString, info::Dict, overwrite::Bool
)

    @info "$(Dates.now()) - Starting data download of $(info["product"]) data for $(date) ..."

    for ii = 1 : length(fname); mimicget(furl,fname[ii],tdir,overwrite); end

    @info "$(Dates.now()) - Data download of $(info["product"]) data for $(date) has been completed."

end

function mimicextract(
    fname::Array{<:AbstractString,2}, fol::AbstractString, info::Dict,
    reg::AbstractString="GLB"
)

    lon,lat = mimiclonlat(); rlon,rlat,rinfo = gregiongridvec(reg,lon,lat);
    nlon = length(lon); nrlon = length(rlon);
    nlat = length(lat); nrlat = length(rlat);
    nt = length(fname); rtime = convert(Array,1:nt); tunit = "hours";

    data = zeros(Int16,nrlon,nrlat,nt); dataii = zeros(Int16,nrlon,nrlat);
    rawi = zeros(nlon,nlat); raw = zeros(nlon,nlat); tmp = zeros(nrlon,nrlat);

    @info "$(Dates.now()) - Extracting regional $(info["product"]) data for $(rinfo["fullname"])."

    for ii = 1 : nt

        fii = joinpath(fol,fname[ii]);
        if isfile(fii)

            try
                ds = Dataset(fii); raw = ds["tpwGrid"].var[:]; close(ds);
                tmp .= regionextractgrid(raw,rinfo,lon,lat,rawi)
                real2int16!(dataii,tmp,offset=60,scale=60/32767);
                data[:,:,ii] .= dataii;
            catch
                @warn "$(Dates.now()) - Unable to extract/open $(fii).  Setting MIMIC tropospheric precipitable water data values set to 'missing'."
                data[:,:,ii] .= -32768;
            end

        else

            @info "$(Dates.now()) - $(fii) does not exist.  All values set to 'missing'."
            data[:,:,ii] .= -32768;

        end

    end

    @info "$(Dates.now()) - Extracted regional $(info["product"]) data for $(rinfo["fullname"])."
    return data,[rlon,rlat,rtime,tunit]

end

function mimicdwn(
    regions::Array{<:AbstractString,1}, date::TimeType, info::Dict; overwrite::Bool
)

    tdir = clisattmpfol(info); if !isdir(tdir) mkpath(tdir); end
    fname,furl = mimicnclist(date); mimicretrieve(fname,furl,date,tdir,info,overwrite);

    for reg in regions; isgeoregion(reg);
        data,grid = mimicextract(fname,tdir,info,reg);
        clisatrawsave(data,grid,reg,info,date)
    end

    clisatrmtmp(fname,tdir);

end


text:
import AperturePhotometry
# Prova
export dark_calibration, flat_calibration

@doc raw"""
    dark_calibration(channel::Array, path::String)

Remove the dark-current frame.
"""
function dark_calibration(channel::Array, path::String)

    dark_arr, time = load_images(path)
    dark = stack_images(dark_arr)
    dark = UInt16.(round.(dark))

    for i in 1:length(channel)
        channel[i] = UInt16.(round.(channel[i] .- dark))
    end

    return channel

end

@doc raw"""
    flat_calibration(channel::Array, path_df::String, path_flat::String)

Remove the flat's anisotrpies using a flat-referance.
"""
function flat_calibration(channel::Array, path_df::String, path_flat::String)

    flat_arr , _ = load_images(path_flat)

    dark_calibration(flat_arr, path_df)
    flat = stack_images(flat_arr)
    flat = 1 ./ flat
    flat = flat ./ maximum(flat)

    for i in 1:length(channel)
        channel[i] = UInt16.(round.(channel[i] .* flat))
    end

    return channel

end


text:
mutable struct KNNUCB <: BanditPolicy
    A::Int64
    t::Int64
    traj::Trajectory
    phi::Function
    theta::Float64
end

function KNNUCB(A::Int64, traj::Trajectory, phi::Function, theta::Float64)
    return KNNUCB(A, 0, traj, phi, theta)
end

function sort_nearest(pol::KNNUCB, x::Vector{Float64}, a::Int64)
    dist = Euclidean()
    x1 = reshape(x, (1, :))
    X = get_X(pol.traj, a, false)
    r = get_r(pol.traj, a, false)
    pwdist = pairwise(dist, X, x1; dims = 1)
    pwdist = reshape(pwdist, (:))
    inds = sortperm(pwdist)
    return pwdist[inds], r[inds]
end

function compute_upper_bound(pol::KNNUCB, x::Vector{Float64}, a::Int64)
    t = pol.t
    distances, r_nearest = sort_nearest(pol, x, a)
    uvals = Vector{Float64}(undef, length(r_nearest))
    for k = 1:length(r_nearest)
        radius = distances[k]
        uvals[k] = sqrt(pol.theta * log(t) / k) + pol.phi(t) * radius
    end
    min_uncertainty, min_k = findmin(uvals)
    return mean(r_nearest[1:min_k]) + min_uncertainty
end

function choose(x::Vector{Float64}, pol::KNNUCB)
    if pol.t + 1 <= pol.A
        return pol.t + 1
    end
    u = [compute_upper_bound(pol, x, a) for a in 1:pol.A]
    return findmax(u)[2]
end

function update!(pol::KNNUCB, x::Vector{Float64}, a::Int64, r::Float64)
    pol.t += 1
    if pol.t != pol.traj.t
        throw(ErrorException("Policy out of sync with trajectory"))
    end
end


text:
using PortNumbers
using FactCheck

facts( "Port Numbers" ) do
    p = PortNumber{UserPort}( 1234 )
    @fact Int( p ) --> 1234
    @fact string( p ) --> "1234"
    @fact ASCIIString( p ) --> "1234"
    @fact typeof( make_port_number( 22 )) --> PortNumber{EstablishedPort}
    @fact typeof( make_port_number( 512 )) --> PortNumber{RootPort}
    @fact typeof( make_port_number( 1234 )) --> PortNumber{UserPort}

end

max_stars_repo_name: sabjohnso/PortNumbers

max_issues_repo_name: sabjohnso/PortNumbers

max_forks_repo_name: sabjohnso/PortNumbers


text:
using OscillatoryIntegralsODE
using OrdinaryDiffEq
using Test

@testset "bessel" begin
    nu = 100.
    r = 100.
    f(x) = exp(-x^2/16)
    bi = BesselJIntegral{Float64}(nu, r)
    result = levintegrate(bi, f, 1.0, 5.0, Vern9(); abstol=1e-6, reltol=1e-6)
    ref = 0.006311630277650583
    @test abs( (result - ref) / ref ) < 1e-4
end

@testset "spherical bessel" begin
    nu = 100.
    r = 100.
    f(x) = exp(-x^2/16)
    bi = SphericalBesselJIntegral{Float64}(nu, r)
    result = levintegrate(bi, f, 1.0, 5.0, Vern9(); abstol=1e-6, reltol=1e-6)
    ref = 0.0008322179291456167
    @test abs( (result - ref) / ref ) < 1e-4
end

@testset "harmonic" begin
    ω = 100.
    f(x) = exp(-x^2/16)
    hi = HarmonicIntegral{Complex{Float64}}(ω)
    result = levintegrate(hi, f, 1.0, 5.0, Vern9(); abstol=1e-6, reltol=1e-6)
    ref = 0.003798145404262353+0.009953606082753418im
    @test abs(result - ref) / abs(ref) < 1e-4
end

max_stars_repo_name: xzackli/OscillatoryIntegralsODE

max_issues_repo_name: xzackli/OscillatoryIntegralsODE

max_forks_repo_name: xzackli/OscillatoryIntegralsODE


text:
using Percolation
using Base.Test

# square
## nearest neighbor
@test try
	N = 50; p = 0.6
	sq = Square(N, p, "nn")
	label_components!(sq)
	heatmap!(sq)
	true
catch
	false
end

## next nearest neighbor
@test try
	N = 100; p = 0.45
	sq = Square(N, p, "nnn")
	label_components!(sq)
	heatmap!(sq)
	true
catch
	false
end

# triangular
@test try
	N = 100; p = 0.45
	tri = Triangular(N, p, "nnn")
	label_components!(tri)
	heatmap!(tri)
	true
catch
	false
end

# honeycomb
@test try
	N = 100; p = 0.45
	h = Honeycomb(N, p, "nnn")
	label_components!(h)
	heatmap!(h)
	true
catch
	false
end

# Percolating probability
## square lattice, nearest neighbor
@test try
	linsize = 50
	ps = 0.1
	pinc = 0.025
	pf = 1
	nsample = 50
	plot_percolation_prob("square", linsize, ps, pinc, pf, nsample)
	true
catch
	false
end

## triangular lattice
@test try
	linsize = 50
	ps = 0
	pinc = 0.025
	pf = 1
	nsample = 50
	plot_percolation_prob("triangular", linsize, ps, pinc, pf, nsample)
	true
catch
	false
end

## honeycomb lattice
@test try
	linsize = 50
	ps = 0
	pinc = 0.025
	pf = 1
	nsample = 50
	plot_percolation_prob("honeycomb", linsize, ps, pinc, pf, nsample)
	true
catch
	false
end

# Cluster property
lattice_types = [:Square, :Triangular, :Honeycomb]
@test try
    for lattice_type in lattice_types
        eval(quote
    	    N = 50; p = 0.6
    	    l = $lattice_type(N, p, "nn")
    	    heatmap!(l)
    	    clustersize!(l)
    	    clustersizefreq!(l)
    	    clusternumber!(l)
    	    average_clustersize!(l)
    	    strength!(l)
    end)
    end
    true
catch
	false
end

# Lifetime of forest fire
@test try
	linsize = 50
	ps = 0
	pinc = 0.025
	pf = 1
	nsample = 50
	plot_lifetime(linsize, ps, pinc, pf, nsample)
	true
catch
	false
end


text:
#!/usr/bin/env julia

module GetPotential

using CSV, DataFrames, Serialization, HDF5, YAML
using ArgParse, Printf

using CarrierCapture

# reading potential energy surface data
get_QE_data(data::String) = names!(CSV.read(data; copycols=true), [:Q, :E])
function get_QE_data(data::Dict)
    QE_data = DataFrame()
    Q = data["Q"]; E = data["E"]
    if Q isa Number
        QE_data.Q = Q; QE_data.E = E
    else
        QE_data.Q = parse.(Float64, split(Q))
        QE_data.E = parse.(Float64, split(E))
    end
    return QE_data
end

println(raw"
      _____                _            _____            _
     / ____|              (_)          / ____|          | |
    | |     __ _ _ __ _ __ _  ___ _ __| |     __ _ _ __ | |_ _   _ _ __ ___
    | |    / _` | '__| '__| |/ _ \ '__| |    / _` | '_ \| __| | | | '__/ _ \
    | |___| (_| | |  | |  | |  __/ |  | |___| (_| | |_) | |_| |_| | | |  __/
     \_____\__,_|_|  |_|  |_|\___|_|   \_____\__,_| .__/ \__|\__,_|_|  \___|
                                                  | |
                                                  |_| v 0.1
      ____      _   ____       _             _   _       _
     / ___| ___| |_|  _ \ ___ | |_ ___ _ __ | |_(_) __ _| |
    | |  _ / _ \ __| |_) / _ \| __/ _ \ '_ \| __| |/ _` | |
    | |_| |  __/ |_|  __/ (_) | ||  __/ | | | |_| | (_| | |
     \____|\___|\__|_|   \___/ \__\___|_| |_|\__|_|\__,_|_|
")

# read arguments and input file
s = ArgParseSettings()
@add_arg_table s begin
    "--input", "-i"
      help = "Input file in a yaml format (default:input.yaml)"
      default = "input.yaml"
      arg_type = String
    "--dryrun"
    	help = "Turn off the Srödinger equation solver"
    	action = :store_true
    "--plot", "-p"
    	help = "Plot potentials"
    	action = :store_true
    "--verbose", "-v"
      help = "write verbose capture coefficient"
      action = :store_true
end

args = parse_args(ARGS, s)
input = YAML.load(open(args["input"]))
is_verbose = parse_args(ARGS, s)["verbose"]

# Set up global variables
Qi, Qf, NQ = input["Qi"], input["Qf"], input["NQ"]
Q = range(Qi, stop=Qf, length=NQ)

# solve potentials
pots = Dict{String, potential}()
for potential in input["potentials"]
	pot_cfg = potential["potential"]
	# TODO: make CSV.read compatable with a format "-.30073981E+03"
	#       consider <CSVFiles> package
    QE_data = get_QE_data(pot_cfg["data"])

	pot = pot_from_dict(QE_data, pot_cfg)
	pots[pot.name] = pot
	fit_pot!(pot, Q)
	if args["dryrun"] == false
		solve_pot!(pot)
	end
end

# print zero-phonon frequencies
if args["dryrun"] == false
	println("===========ħω0===========")
	for (name, pot) in pots
		ϵ0 = pot.ϵ[2] - pot.ϵ[1]
		ħω0 = @sprintf("%.1f", ϵ0*1000)
		println("$(name): $(ħω0) meV")
	end
	println("=========================\n")
end

# find crossing point
findcross = get(input, "findcross", nothing)
if findcross != nothing
    for crossing in findcross
        crossing = crossing["crossing"]
        Q_cross, E_cross = find_crossing(pots[crossing["pot_name_1"]],
                                         pots[crossing["pot_name_2"]])
        println("=====Crossing Points=====")
        println(": $(crossing["pot_name_1"]) and $(crossing["pot_name_2"]) ")
        @printf "Q: %0.5f amu^-1/2 Å\n" Q_cross
        @printf "E: %0.5f eV\n" E_cross
        println("=========================\n")
    end
end

# save potential structs
open("wave.jld", "w") do file
    println("====Saving wave functions====\n")
    serialize(file, pots)
end

# plot
if args["plot"] == true
    println("====Ploting potentials====\n")
	plot_cfg = get(input, "plot", Nothing)
	plot_pots(pots, plot_cfg)
end

# write potential cvs
for (name, pot) in pots
    CSV.write("pot_fit_$(name).csv", DataFrame([pot.Q, pot.E], [:Q, :E]))
end

# write eigenvalues cvs
for (name, pot) in pots
    CSV.write("eigval_$(name).csv", DataFrame([pot.ϵ .- pot.E0, pot.ϵ], [:ϵ_E0, :ϵ]))
end

# write verbose ouput in HDF5 format
if is_verbose
	for (name, pot) in pots
		h5open("wave_func_$(name).hdf5","w") do file
		    write(file, "wave_func", collect(pot.χ))
		end
	end
end

end # module


text:
"""
Using JLD macros. Danger of being unreliable.
"""
using QuNet
using LightGraphs, SimpleWeightedGraphs, GraphPlot, MetaGraphs, Plots, Colors, Statistics
using LaTeXStrings
# using DelimitedFiles
using JLD

"""
Plot the performance statistics of greedy_multi_path for 1 random end-user pair
over a range of graph sizes and for different numbers of max_path (1,2,3) (ie. the
maximum number of paths that can be purified by a given end-user pair.)
"""

generate_new_data = true
if generate_new_data == true

    # Params
    num_trials = 1000::Int64
    min_size = 2::Int64
    max_size = 20::Int64

    @assert min_size < max_size

    perf_data1 = []
    perf_data2 = []
    perf_data3 = []
    perf_data4 = []

    size_list = collect(min_size:1:max_size)
    for (j, data_array) in enumerate([perf_data1, perf_data2, perf_data3, perf_data4])
        println("Collecting for max_paths: $j")
        for i in size_list
            println("Collecting for gridsize: $i")
            # Generate ixi graph:
            net = GridNetwork(i, i)

            # Collect performance statistics
            performance, dummy, dummy, dummy = net_performance(net, num_trials, 1, max_paths=j)
            push!(data_array, performance)
        end
    end

    # Get values for x axis
    x = collect(min_size:1:max_size)

    # Save variables to julia data file
    @save "maxpaths.jdl"
        min_size
        max_size
        num_trials
        perf_data1
        perf_data2
        perf_data3
        perf_data4
end

# Load variables from data file
@load "maxpaths.jdl"
    min_size
    max_size
    num_trials
    perf_data1
    perf_data2
    perf_data3
    perf_data4

# Extract data from performance data
loss_arr1 = collect(map(x->x["loss"], perf_data1))
z_arr1 = collect(map(x->x["Z"], perf_data1))
loss_arr2 = collect(map(x->x["loss"], perf_data2))
z_arr2 = collect(map(x->x["Z"], perf_data2))
loss_arr3 = collect(map(x->x["loss"], perf_data3))
z_arr3 = collect(map(x->x["Z"], perf_data3))
loss_arr4 = collect(map(x->x["loss"], perf_data4))
z_arr4 = collect(map(x->x["Z"], perf_data4))

# Plot
plot(x, loss_arr1, linewidth=2, label=L"$\eta_1$", xlims=(0, max_size), color = :red, legend=:left)
plot!(x, z_arr1, linewidth=2, label=L"$F_1$", linestyle=:dash, color =:red)
plot!(x, loss_arr2, linewidth=2, label=L"$\eta_2$", color =:blue)
plot!(x, z_arr2, linewidth=2, label=L"$F_2$", linestyle=:dash, color =:blue)
plot!(x, loss_arr3, linewidth=2, label=L"$\eta_3$", color =:green)
plot!(x, z_arr3, linewidth=2, label=L"$F_3$", linestyle=:dash, color =:green)
plot!(x, loss_arr4, linewidth=2, label=L"$\eta_3$", color =:purple)
plot!(x, z_arr4, linewidth=2, label=L"$F_3$", linestyle=:dash, color =:purple)
xaxis!(L"$\textrm{Grid Size}$")

savefig("cost_maxpaths.png")
savefig("cost_maxpaths.pdf")

# Plot
plot(x, loss_arr1, linewidth=2, label=L"$\eta_1$", xlims=(0, max_size), yaxis=:log, color = :red, legend=:left)
plot!(x, loss_arr2, linewidth=2, label=L"$\eta_2$", yaxis=:log, color =:blue)
plot!(x[3:end], loss_arr3[3:end], linewidth=2, label=L"$\eta_3$", yaxis=:log, color =:green)
plot!(x[3:end], loss_arr4[3:end], linewidth=2, label=L"$\eta_3$", yaxis=:log, color =:purple)
xaxis!(L"$\textrm{Grid Size}$")

savefig("cost_maxpathslog.png")
savefig("cost_maxpathslog.pdf")

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
## A template for defining a new optimizer model

# Define a mutable struct for storing infinite model to optimizer model mappings
# plus other needed information
mutable struct NewReformData
    # Variable mapping (REPLACE AND REFORMAT BELOW AS NEEDED)
    infvar_mappings::Dict{GeneralVariableRef, Vector{JuMP.VariableRef}}
    finvar_mappings::Dict{GeneralVariableRef, JuMP.VariableRef}

    # Map other variable info (REPLACE AND REFORMAT BELOW AS NEEDED)
    infvar_to_supports::Dict{GeneralVariableRef, Vector{<:Tuple}}

    # Measure mappings and other info (ONLY NEEDED TO ENABLE MEASURE QUERIES)
    meas_mappings::Dict{GeneralVariableRef, Vector{JuMP.VariableRef}}
    meas_to_supports::Dict{GeneralVariableRef, Vector{<:Tuple}}

    # Constraint mapping (REPLACE AND REFORMAT BELOW AS NEEDED)
    constr_mappings::Dict{InfOptConstraintRef, Vector{JuMP.ConstraintRef}}

    # Map other constraint info (REPLACE AND REFORMAT BELOW AS NEEDED)
    constr_to_supports::Dict{InfOptConstraintRef, Vector{<:Tuple}}

    # ADD OTHER MAPPING/MODEL INFORMATION HERE

    # default constructor
    function NewReformData() # REFORMAT BELOW IN ACCORDANCE WITH ABOVE ATTRIBUTES
        return new(Dict{GeneralVariableRef, Vector{JuMP.VariableRef}}(),
                   Dict{GeneralVariableRef, JuMP.VariableRef}(),
                   Dict{GeneralVariableRef, Vector{Tuple}}(),
                   Dict{GeneralVariableRef, Vector{JuMP.VariableRef}}(),
                   Dict{GeneralVariableRef, Vector{Tuple}}(),
                   Dict{InfOptConstraintRef, Vector{JuMP.ConstraintRef}}(),
                   Dict{InfOptConstraintRef, Vector{Tuple}}()
                   )
    end
end

# Define the optimizer model key
const OptKey = :ReformData # REPLACE WITH A DESIRED UNIQUE KEY

# Make a constructor for new optimizer model type (extension of JuMP.Model)
function NewReformModel(args...; kwargs...)::JuMP.Model # ADD EXPLICT ARGS AS NEEDED
    # initialize the JuMP Model
    model = JuMP.Model(args...; kwargs...)

    # ADD ADDITIONAL OPERATIONS IF NEEDED

    # add the extension data with a chosen optimizer model key and return
    model.ext[OptKey] = NewReformData()
    return model
end

# Make function for extracting the data from the model (optional)
function reform_data(model::JuMP.Model)::NewReformData
    # UPDATE THE NOMENCLATURE AS NEEDED
    haskey(model.ext, OptKey) || error("Model is not a NewReformModel.")
    return model.ext[OptKey]
end

# Extend build_optimizer_model! (enables build_optimizer_model! and optimize!)
function InfiniteOpt.build_optimizer_model!(
    model::InfiniteModel,
    key::Val{OptKey};
    my_kwarg::Bool = true # ADD KEYWORD ARGUMENTS AS NEEDED
    )
    # clear the model for a build/rebuild
    reform_model = clear_optimizer_model_build!(model)

    # IT MAY BE USEFUL TO CALL `expand_all_measures!` TO HANDLE MEASURES FIRST
    # otherwise can extend `add_measure_variable` and `delete_semi_infinite_variable` to
    # expand in place without modifying the infinite model

    # REPLACE LINES 73-89 WITH OPERATIONS TO BUILD A `NewReformModel` BASED ON `model`
    # these lines just generate artificial data, see `TranscriptionOpt` for a thorough example of implementation
    data = reform_data(reform_model)
    for vref in all_variables(model)
        if index(vref) isa Union{InfiniteVariableIndex, SemiInfiniteVariableIndex}
            data.infvar_mappings[vref] = [@variable(reform_model) for i = 1:2]
            data.infvar_to_supports[vref] = [(0.,), (1.,)]
        else
            data.finvar_mappings[vref] = @variable(reform_model)
        end
    end
    for mref in all_measures(model)
        data.meas_mappings[mref] = [@variable(reform_model) for i = 1:2]
        data.meas_to_supports[mref] = [(-1.,), (-2.,)]
    end
    for cref in all_constraints(model)
        data.constr_mappings[cref] = [@constraint(reform_model, @variable(reform_model) >= 0) for i = 1:2]
        data.constr_to_supports[cref] = [(2.,), (3.,)]
    end
    set_objective(reform_model, objective_sense(model), @variable(reform_model))
    # update the optimizer model status
    set_optimizer_model_ready(model, true)
    return
end

# Extend optimizer_model_variable if appropriate to enable variable related queries
function InfiniteOpt.optimizer_model_variable(
    vref::GeneralVariableRef,
    key::Val{OptKey};
    my_kwarg::Bool = true # ADD KEY ARGS AS NEEDED
    )
    # REPLACE BELOW WITH ACTUAL CORRESPONDENCE TO THE OPTIMIZER MODEL VARIABLE(S)
    model = optimizer_model(JuMP.owner_model(vref))
    vindex = index(vref)
    if vindex isa Union{InfiniteVariableIndex, SemiInfiniteVariableIndex}
        map_dict = reform_data(model).infvar_mappings
    elseif vindex isa MeasureIndex
        map_dict = reform_data(model).meas_mappings
    else
        map_dict = reform_data(model).finvar_mappings
    end
    haskey(map_dict, vref) || error("Variable $vref not used in the optimizer model.")
    return map_dict[vref]
end

# Extend optimizer_model_expression if appropriate to enable expression related queries
function InfiniteOpt.optimizer_model_expression(
    expr::Union{JuMP.GenericAffExpr, JuMP.GenericQuadExpr},
    key::Val{OptKey};
    my_kwarg::Bool = true # ADD KEY ARGS AS NEEDED
    )
    # REPLACE BELOW WITH ACTUAL CORRESPONDENCE TO REFORMULATED EXPRESSIONS
    reform_expr = zero(AffExpr)
    return reform_expr
end

# Extend optimizer_model_constraint if appropriate to enable constraint related queries
function InfiniteOpt.optimizer_model_constraint(
    cref::InfOptConstraintRef,
    key::Val{OptKey}; # ADD KEY ARGS AS NEEDED
    my_kwarg::Bool = true)
    # REPLACE BELOW WITH ACTUAL CORRESPONDENCE TO THE OPTIMIZER MODEL CONSTRAINT(S)
    model = optimizer_model(JuMP.owner_model(cref))
    map_dict = reform_data(model).constr_mappings
    haskey(map_dict, cref) || error("Constraint $cref not used in the optimizer model.")
    return map_dict[cref]
end

# If appropriate extend variable_supports (enables support queries of infinite variables)
function InfiniteOpt.variable_supports(
    model::JuMP.Model,
    vref::Union{InfiniteVariableRef, SemiInfiniteVariableRef},
    key::Val{OptKey};
    my_kwarg::Bool = true # ADD KEY ARGS AS NEEDED
    )
    # REPLACE BELOW WITH ACTUAL CORRESPONDENCE TO THE INFINITE VARIABLE SUPPORT VALUES
    map_dict = reform_data(model).infvar_to_supports
    gvref = InfiniteOpt._make_variable_ref(owner_model(vref), index(vref))
    haskey(map_dict, gvref) || error("Variable $gvref not used in the optimizer model.")
    return map_dict[gvref]
end

# If appropriate extend variable_supports for measures (enables support queries of measures)
function InfiniteOpt.variable_supports(
    model::JuMP.Model,
    mref::MeasureRef,
    key::Val{OptKey};
    my_kwarg::Bool = true # ADD KEY ARGS AS NEEDED
    )
    # REPLACE BELOW WITH ACTUAL CORRESPONDENCE TO THE INFINITE VARIABLE SUPPORT VALUES
    map_dict = reform_data(model).meas_to_supports
    gvref = InfiniteOpt._make_variable_ref(owner_model(mref), index(mref))
    haskey(map_dict, gvref) || error("Variable $gvref not used in the optimizer model.")
    return map_dict[gvref]
end

# If appropriate extend expression_supports (enables support queries of expressions)
function InfiniteOpt.expression_supports(
    model::JuMP.Model,
    expr::Union{JuMP.GenericAffExpr, JuMP.GenericQuadExpr},
    key::Val{OptKey};
    my_kwarg::Bool = true # ADD KEY ARGS AS NEEDED
    )
    # REPLACE BELOW WITH ACTUAL CORRESPONDENCE TO SUPPORT(S) OF THE EXPRESSION(S)
    supps = [(-42.,), (1.,)]
    return supps
end

# If appropriate extend constraint_supports (enables support queries of constraints)
function InfiniteOpt.constraint_supports(
    model::JuMP.Model,
    cref::InfOptConstraintRef,
    key::Val{OptKey}; # ADD KEY ARGS AS NEEDED
    my_kwarg::Bool = true)
    # REPLACE BELOW WITH ACTUAL CORRESPONDENCE TO THE INFINITE VARIABLE SUPPORT VALUES
    map_dict = reform_data(model).constr_to_supports
    haskey(map_dict, cref) || error("Constraint $cref does not have supports in the optimizer model.")
    return length(map_dict[cref]) == 1 ? first(map_dict[cref]) : map_dict[cref]
end

# If it doesn't make sense to extend optimizer_model_variable,
# optimizer_model_expression and/or optimizer_model_constraint then you'll need
# to extend the following:
# - InfiniteOpt.map_value to enable JuMP.value
# - InfiniteOpt.map_optimizer_index to enable JuMP.optimizer_index
# - InfiniteOpt.map_dual to enable JuMP.dual
# - InfiniteOpt.map_shadow_price to enable JuMP.shadow_price


text:
using Oceananigans.Operators
using Oceananigans.Solvers: FFTBasedPoissonSolver, FourierTridiagonalPoissonSolver, solve!
using Oceananigans.Architectures: device_event
using Oceananigans.Distributed: DistributedFFTBasedPoissonSolver

#####
##### Calculate the right-hand-side of the non-hydrostatic pressure Poisson equation.
#####

@kernel function calculate_pressure_source_term_fft_based_solver!(rhs, grid, Δt, U★)
    i, j, k = @index(Global, NTuple)
    @inbounds rhs[i, j, k] = divᶜᶜᶜ(i, j, k, grid, U★.u, U★.v, U★.w) / Δt
end

@kernel function calculate_pressure_source_term_fourier_tridiagonal_solver!(rhs, grid, Δt, U★)
    i, j, k = @index(Global, NTuple)
    @inbounds rhs[i, j, k] = Δzᶜᶜᶜ(i, j, k, grid) * divᶜᶜᶜ(i, j, k, grid, U★.u, U★.v, U★.w) / Δt
end

#####
##### Solve for pressure
#####

function solve_for_pressure!(pressure, solver::DistributedFFTBasedPoissonSolver, Δt, U★)
    rhs = first(solver.storage)
    arch = architecture(solver)
    grid = solver.local_grid

    rhs_event = launch!(arch, grid, :xyz, calculate_pressure_source_term_fft_based_solver!,
                        rhs, grid, Δt, U★, dependencies = device_event(arch))

    wait(device(arch), rhs_event)

    # Solve pressure Poisson equation for pressure, given rhs
    solve!(pressure, solver)

    return pressure
end

function solve_for_pressure!(pressure, solver::FFTBasedPoissonSolver, Δt, U★)

    # Calculate right hand side:
    rhs = solver.storage
    arch = architecture(solver)
    grid = solver.grid

    rhs_event = launch!(arch, grid, :xyz, calculate_pressure_source_term_fft_based_solver!,
                        rhs, grid, Δt, U★, dependencies = device_event(arch))

    wait(device(arch), rhs_event)

    # Solve pressure Poisson given for pressure, given rhs
    solve!(pressure, solver, rhs)

    return nothing
end

function solve_for_pressure!(pressure, solver::FourierTridiagonalPoissonSolver, Δt, U★)

    # Calculate right hand side:
    rhs = solver.source_term
    arch = architecture(solver)
    grid = solver.grid

    rhs_event = launch!(arch, grid, :xyz, calculate_pressure_source_term_fourier_tridiagonal_solver!,
                        rhs, grid, Δt, U★, dependencies = device_event(arch))

    wait(device(arch), rhs_event)

    # Pressure Poisson rhs, scaled by Δzᶜᶜᶜ, is stored in solver.source_term:
    solve!(pressure, solver)

    return nothing
end


text:
@testset "Essentials" begin
    @testset "Spatials" begin include("Spatials.jl") end
    @testset "DegreesOfFreedom" begin include("DegreesOfFreedom.jl") end
    @testset "Terms" begin include("Terms.jl") end
    @testset "FockPackage" begin include("FockPackage.jl") end
    @testset "SpinPackage" begin include("SpinPackage.jl") end
    @testset "Frameworks" begin include("Frameworks.jl") end
end


text:
using DrWatson
using DataFrames
using Statistics
using CairoMakie
using LaTeXStrings
using CSVFiles
using Dates
using Interpolations

sims = DataFrame(load(projectdir("figures", "figure_chemotaxis", "data", "simulation_2021-10-12.csv")))
sims_wo_methylation = DataFrame(load(projectdir("figures", "figure_chemotaxis", "data", "simulation_2022-01-27.csv")))
imi_data = DataFrame(load(projectdir("figures", "figure_chemotaxis", "data", "imi.csv")))
imi_data_no_methylation = DataFrame(load(projectdir("figures", "figure_chemotaxis", "data", "imi_no_methylation.csv")))
gaus_data = DataFrame(load(projectdir("figures", "figure_chemotaxis", "data", "gaussian_rate.csv")))

normalfont = projectdir("figures", "fonts", "NotoSans-Condensed.ttf")
italicfont = projectdir("figures", "fonts", "NotoSans-Italic.ttf")
boldfont = projectdir("figures", "fonts", "NotoSans-CondensedSemiBold.ttf")
subfigurelabelfont = projectdir("figures", "fonts", "NotoSans-Bold.ttf")

plot_theme = Theme(
    font = normalfont,
    fontsize = 8,
    colgap = 4,
    rowgap = 0,
    linewidth = 2,
    markersize = 5,
    spinewidth = 0.5,
    figure_padding = (2, 2, 2, 2),
    Axis = (
        spinewidth = 1.0,
        xgridwidth = 1.0,
        xminorgridwidth = 1.0,
        yminorgridwidth = 1.0,
        ygridwidth = 1.0,
        xtickwidth = 1.0,
        ytickwidth = 1.0,
        xticksize = 3,
        yticksize = 3,
        xlabelpadding = 0,
        ylabelpadding = 0,
        xticklabelpad = 0,
        yticklabelpad = 1,
    ),
    Legend = (
        framevisible = false,
        framewidth = 1.0,
        rowgap = 2,
        titlegap = 2,
        titlefont = boldfont,
        titlehalign = :left,
        halign = :center,
        gridshalign = :left,
        titlefontsize = 11,
        groupgap = 5,
        patchsize = (10, 5),
        padding = (1, 1, 1, 1),
        linewidth = 1.5,
    )
)
set_theme!(plot_theme)

fig = Figure(resolution = (246, 240))
ax1 = Axis(fig[1, 1], xlabel = "traj. duration (s)", title = "mutual\ninformation (bits)")

groups = filter(g -> g.TauL[1] ∈ [0.01, 0.1, 1.0, 10.0], groupby(sims, [:M, :TauL], sort = true))
cmap = cgrad(:OrRd_9)
colors = cmap[3:2:9]

for (g, c) in zip(groups, colors)
    lines!(ax1, g.DiscreteTimes, g.Mean ./ log(2), label = "$(g.TauL[1]) s", color = c)
end
Legend(fig[1, 3], ax1, "input\ncorrelation\ntime", orientation = :vertical, tellwidth = true)
# colgap!(fig.layout, 2, Fixed(-35))

ax2 = Axis(fig[1, 2], xlabel = "traj. duration (s)", title = "information\n rate (bits/s)")
for (g, c) in zip(groups, colors)
    nodes = g.DiscreteTimes[1]:10*diff(g.DiscreteTimes)[1]:g.DiscreteTimes[end]
    itp = CubicSplineInterpolation(nodes, g.Mean[begin:10:end] ./ log(2))
    us = range(extrema(g.DiscreteTimes)...; step = 2)
    vs = [Interpolations.gradient(itp, u)[1] for u in us]
    lines!(ax2, us, vs, color = c)
end
ylims!(ax2, 0, nothing)

subgl = GridLayout()
fig.layout[2, 1:2] = subgl
rowsize!(fig.layout, 2, Relative(0.7))

ax3 = Axis(
    subgl[1, 1],
    title = "optimal input timescale",
    ylabel = "bits/s",
    xscale = log10,
    xminorticks = IntervalsBetween(9),
    xminorgridvisible = true,
    titlevisible = false
)

lines!(ax3, gaus_data.TauL, gaus_data.InfoRate, color = "#2b8cbe", label = "Gaussian")

function get_info_rates(sims)
    xvals = Float64[]
    yvals = Float64[]
    for g in groupby(sims, [:M, :TauL], sort = true)
        push!(xvals, g.TauL[1])
        indexfrom = searchsortedfirst(g.DiscreteTimes, 150)
        rates = diff(g.Mean) ./ diff(g.DiscreteTimes) ./ log(2)
        asymptotic_rate = mean(rates[indexfrom:end])
        push!(yvals, asymptotic_rate)
    end
    xvals, yvals
end

xvals, yvals = get_info_rates(sims)
lines!(ax3, xvals, yvals, color = :black, label = "PWS")

xvals, yvals = get_info_rates(sims_wo_methylation)
lines!(ax3, xvals, yvals, color = :black, linestyle = :dot, label = "push-pull")
xlims!(ax3, 1e-2, 10)
ylims!(ax3, 0, 25)
text!(ax3, "information transmission rate R", position = Point(1.1e-2, 20), textsize = 8, align = (:left, :baseline), font = boldfont)

c4 = :black
ax4 = Axis(
    subgl[2, 1],
    xminorticks = IntervalsBetween(9),
    xminorgridvisible = true,
    ytickcolor = c4,
    yticklabelcolor = c4,
    yaxisposition = :left,
    xscale = log10,
    ylabel = "bits",
    xlabel = L"input correlation time $\tau_L$ (s)",
    ylabelcolor = c4)
text!(ax4, "instantaneous mutual information", position = Point(1.1e-2, 1.15), textsize = 8, align = (:left, :baseline), font = boldfont)

# ax4.xticks = [0.01, 0.1, 1.0, 10.0]
# ax4.xtickformat = "{:.2f}"
hidexdecorations!(ax3, grid = false, minorgrid = false)
hidespines!(ax4, :t)

imi_line = lines!(ax4, imi_data.TauL, imi_data.IMI, color = c4, label = "chemotaxis")
translate!(imi_line, 0, 0, 9) # bring to foreground
lines!(ax4, gaus_data.TauL, gaus_data.InstantaneousInfo, color = "#2b8cbe", label = "Gaussian\n(LNA)")
lines!(ax4, imi_data_no_methylation.TauL, imi_data_no_methylation.IMI, color = c4, linestyle = :dot, label = "push-pull")
xlims!(ax4, 1e-2, 10)
ylims!(ax4, 0, 1.49)
# annotations!(ax4, ["no methylation"], [Point2(4, 1.7)], font = boldfont, align = (:right, :center), textsize = 7, rotation = 0.3)

Legend(fig[2, 3], ax4, "systems", valign=0.2, orientation = :vertical, tellwidth = true)

for (label, layout) in zip(["b", "c", "d"], [fig[1, 1], fig[1, 2], fig[2, 1:3]])
    Label(layout[1, 1, TopLeft()], label,
        textsize = 14,
        font = subfigurelabelfont,
        padding = (3, 5, 4, -10),
        halign = :left)
end

# create inset

inset_ax = Axis(fig[2, 1:2],
    width = Relative(0.25),
    height = Relative(0.25),
    titlevisible = false,
    # xminorgridvisible = true,
    # xminorticksvisible = true,
    yminorgridvisible = true,
    xticklabelsize = 7,
    yticklabelsize = 7,
    xlabelpadding = -2,
    ylabelpadding = -2,
    halign = 1.0,
    valign = 1.0,
    xtickalign = 1.0,
    ytickalign = 1.0,
    backgroundcolor = :white,
    xticklabelsvisible = true,
    xgridvisible = true,
    yticklabelsvisible = true,
    xscale = log10,
    xlabel = L"\tau_L",
    ylabel = L"R \tau_L",
    xticks = ([0.01, 10.0], ["0.01", "10"]),
    # xminorticks = 0.01:0.1:10.0,
    yminorticks = IntervalsBetween(4),
    topspinevisible = false,
    rightspinevisible = false,
    yticks = [0, 20],
    xticklabelspace = 0,
    yticklabelspace = 0,
    xlabelfont = italicfont,
    ylabelfont = italicfont
)

xlims!(inset_ax, 1e-2, 1e1)
ylims!(inset_ax, 0, 20)
xvals, yvals = get_info_rates(sims)
lines!(inset_ax, xvals, yvals .* xvals, color = :black, label = "PWS")

translate!(inset_ax.scene, 0, 0, 10)
translate!(inset_ax.elements[:background], 0, 0, 9)
translate!(inset_ax.elements[:xgridlines], 0, 0, 9)
translate!(inset_ax.elements[:ygridlines], 0, 0, 9)
translate!(inset_ax.elements[:xminorgridlines], 0, 0, 9)
translate!(inset_ax.elements[:yminorgridlines], 0, 0, 9)

save(projectdir("figures", "figure_chemotaxis", "chemotaxis_mi.svg"), fig, pt_per_unit = 1)
save(projectdir("figures", "figure_chemotaxis", "chemotaxis_mi.pdf"), fig, pt_per_unit = 1)

fig

# presentation
# f = Figure(resolution = (300, 150))
# ax1 = Axis(f[1, 1], xscale = log10)
# lines!(ax1, 1 ./ imi_data_no_methylation.TauL, imi_data_no_methylation.IMI * 10, color = c4, linestyle = :dot)
# xvals, yvals = get_info_rates(sims_wo_methylation)
# lines!(ax1, 1 ./ xvals, yvals, color = :black)
# save(projectdir("figures", "figure_chemotaxis", "chemotaxis_pres.pdf"), f, pt_per_unit = 1)

# f

# fnew = Figure(resolution = (300, 200))
# ax = Axis(
#     fnew[1, 1],
#     xscale = log10,
#     xminorticks = IntervalsBetween(9),
#     xminorgridvisible = true,
#     xlabel = "input timescale",
#     ylabel = "information accumulated per run (bits)"
# )
# ax.xticks = [0.01, 0.1, 1.0, 10.0]
# ax.xtickformat = "{:.2f}"
# xvals, yvals = get_info_rates(sims)
# lines!(ax, xvals, yvals .* xvals, color = :black, label = "PWS")
# xvals, yvals = get_info_rates(sims_wo_methylation)
# lines!(ax, xvals, yvals .* xvals, color = :black, linestyle = :dot, label = "push-pull")
# lines!(ax, gaus_data.TauL, gaus_data.InfoRate .* gaus_data.TauL, color = "#2b8cbe", label = "Gaussian")

# save(projectdir("figures", "figure_chemotaxis", "info_per_run.pdf"), fnew, pt_per_unit = 1)
# fnew


text:
using Sigma
## From probmods.org / Kahnenman and Tversky
## ==========================================

# The probability of breast cancer is 1% for a woman at 40 who
# participates in a routine screening. If a woman has breast cancer, the probability is
# 80% that she will have a positive mammography. If a woman does not have breast cancer,
# the probability is 9.6% that she will also have a positive mammography.
# A woman in this age group had a positive mammography in a routine screening.
function ground_truth(cancer_rate)
  (cancer_rate * 0.008) / ((0.008 * cancer_rate) + (0.00096 * 0.99))
end

# Model
cancer_base_rate = 0.01
breast_cancer = flip(cancer_base_rate)
positive_mammogram = ifelse(breast_cancer,flip(0.8),flip(0.096))

# Queries
prob(positive_mammogram)
prob(breast_cancer, positive_mammogram)