text:
module Nonlinear

import StatsBase
import Statistics

#=
This function calculates the approximate entropy
:param n: the array that contains the NN-inetrvals
:param m: the embedding dimension, default=2
:param r: the tolerance, default=6
:return: the approximate entropy
=#
function apen(n,m,r)
    c1=get_apen_dist(n,m,r)
    c2=get_apen_dist(n,m+1,r)
    return log(c1/c2)
end # apen

#=
This function calculates the sample entropy
:param n: the array that contains the NN-inetrvals
:param m: the embedding dimension, default=2
:param r: the tolerance, default=6
:return: the sample entropy
=#
function sampen(n,m,r)
    c1=get_sampen_dist(n,m,r,1)
    c2=get_sampen_dist(n,m+1,r,0)
    return -log(c2/c1)
end # sampen

#=
This function creates a template of a given array over an embedding dimension
:param n: the array that contains the NN-inetrvals
:param m: the embedding dimension, default=2
:return template: the created template
=#
function get_template(n,m)
    template=[]
    for i in 1:length(n)-m+1
        push!(template,n[i:i+m-1])
    end
    return template
end # get_template

#=
This function calculates the distances for the approximate entropy
:param n: the array that contains the NN-inetrvals
:param m: the embedding dimension, default=2
:param r: the tolerance, default=6
:return: the distance for the approximate entropy
=#
function get_apen_dist(n,m,r)
    template=get_template(n,m)
    count=zeros(length(template))
    for i in 1:length(template)
        for j in i+1:length(template)
            if maximum(abs.(template[i].-template[j]))<=r
                count[i]+=1
                count[j]+=1
            end
        end
    end
    return sum(count./(length(n)-m+1))/(length(n)-m+1)
end # get_apen_dist

#=
This function calculates the distances for the sample entropy
:param n: the array that contains the NN-inetrvals
:param m: the embedding dimension, default=2
:param r: the tolerance, default=6
:param l: a value to limit the for-loops
:return: the distance for the sample entropy
=#
function get_sampen_dist(n,m,r,l)
    template=get_template(n,m)
    counts=[]
    count=0
    for i in 1:length(template)-l
        for j in 1:length(template)-l
            if maximum(abs.(template[i].-template[j]))>=r || i==j
                push!(counts,count)
                count=0
            else
                count+=1
            end
        end
    end
    return sum(counts)
end # get_sampen_dist

#=
This function calculates the renyi entropy of a given order
:param n: the array that contains the NN-inetrvals
:param a: the order of the renyi entropy
:return: the calculated renyi entropy
=#
function renyi(n,a)
    return StatsBase.renyientropy(n,a)
end # renyi

#=
This function calculates the hurst coefficient
It was inspired by the python hurst package by Dmitry A. Mottl (https://github.com/Mottl/hurst)
:param n: the array that contains the NN-inetrvals
:return H: the hurst coefficient
=#
function hurst(n)
    ws=Array(range(log10(10),stop=log10(length(n)),step=0.25))
    window=[]
    for x in ws
        push!(window,round(Int64,exp10(x),RoundDown))
    end
    if !(length(n) in window)
        push!(window,length(n))
        push!(ws,log10(length(n)))
    end
    RS=[]
    for w in window
        rs=[]
        for start in (range(0,stop=length(n),step=w))
            if (start+w)>length(n)
                break
            end
            RS_part= get_rs(n[start+1:start+w])
            if RS_part != 0
                push!(rs,RS_part)
            end
        end
        if length(rs)>0
            push!(RS,Statistics.mean(rs))
        end
    end
    A=Array{Float64}([ws ones(length(RS))])
    RSlog=[]
    for r in RS
        push!(RSlog,log10(r))
    end
    B=Array{Float64}(RSlog)
    H,c=A\B
    c=exp10(c)
    return H
end # hurst

#=
This function calculates the rescaled range of a time series
It was inspired by the python hurst package by Dmitry A. Mottl (https://github.com/Mottl/hurst)
:param n: the array that contains the NN-inetrvals
:return: the rescaled range
=#
function get_rs(n)
    incs=n[2:end].-n[1:end-1]
    mean_inc=(n[end]-n[1])/length(incs)
    deviations=incs.-mean_inc
    Z=cumsum(deviations)
    R=maximum(Z)-minimum(Z)
    S=Statistics.std(incs)
    if R==0 || S==0
        return 0
    else
        return R/S
    end
end # get_rs

end # module


text:
function rot(ch::Char, key::Integer)
    if key < 1 || key > 25 end
    if isalpha(ch)
        shft = ifelse(islower(ch), 'a', 'A')
        ch = (ch - shft + key) % 26 + shft
    end
    return ch
end
rot(str::AbstractString, key::Integer) = map(x -> rot(x, key), str)

msg = "The five boxing wizards jump quickly"
key = 3
invkey = 26 - 3

println("# original: $msg\n  encrypted: $(rot(msg, key))\n  decrypted: $(rot(rot(msg, key), invkey))")

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@testset "GeometricBM" begin

using DiffEqNoiseProcess, DiffEqBase, DiffEqMonteCarlo, Test, Statistics

μ = 1.0
σ = 2.0

W = GeometricBrownianMotionProcess(μ,σ,0.0,1.0,1.0)
dt = 0.1
calculate_step!(W,dt)
for i in 1:10
  accept_step!(W,dt)
end

W = GeometricBrownianMotionProcess(μ,σ,0.0,1.0,1.0)
prob = NoiseProblem(W,(0.0,1.0))
sol = solve(prob;dt=0.1)

W = GeometricBrownianMotionProcess(μ,σ,0//1,1.0,1.0)
prob_rational = NoiseProblem(W,(0,1))
dt = 1//10
sol = solve(prob_rational;dt=dt)

dt = dt/100
for t in dt:dt:1-dt
  sol(t)
end

t = 1.0
u0 = 1.0
expected_mean = u0*exp(μ*t)
expected_variance = u0^2*exp(2μ*t)*(exp(σ^2*t)-1)
monte_prob = MonteCarloProblem(prob;output_func = (sol,i)-> (sol[end],false))
sol = solve(monte_prob;dt=0.1,num_monte=40000)
@test abs(mean(sol) - expected_mean) < 0.4
#abs(var(sol) - expected_variance) < 0.4 # Converges slowly

W = GeometricBrownianMotionProcess!(μ,σ,0.0,ones(2),ones(2))
prob = NoiseProblem(W,(0.0,1.0))
sol = solve(prob;dt=0.1)

end


text:
module Queues
    include("Consumers.jl")
    import .Consumers

    include("../../utils/Utils.jl")
    import .Utils

    using AMQPClient

    export delete, declare, bind, unbind, purge, unsubscribe, subscribe

    """Функция подписки на сообщения очереди.

    channel: канал связи с сервисом очередей
    queue_name: имя очереди в сервисе очередей
    user_consumer_reaction: пользовательская функция консьюмера
    """
    function subscribe(channel::AMQPClient.MessageChannel, queue_name::String, user_consumer_reaction::Function)::Bool
        @info Utils.Chronometer.message_with_time("Вызвана операция подписки на очередь $queue_name")
        # Логирование для режима отладки
        @debug "Параметры вызова функции:"
        @debug "> канал соединения с серверов RabbitMQ " channel
        @debug "> имя очереди " queue_name
        @debug "> функция консьюмера пользователя " user_consumer_reaction
        try
            consumer(message) = Consumers.base_consumer(channel, message, user_consumer_reaction)
            @debug "Итоговый консьюмер " consumer
            success, consumer_tag = AMQPClient.basic_consume(channel, queue_name, consumer)
            @debug "Результат выполненя подписки " success
            @debug "Тег консьюмера " consumer_tag
            @assert success
            @info Utils.Chronometer.message_with_time("Зарегистрирован консьюмер с тегом $consumer_tag")
            return success
        catch error
            @error Utils.Chronometer.message_with_time("Не удалось подписаться на сообщения очереди ") error
        end
    end

    """ Функция, которая позволяет отписаться от прослушивания очереди.

    channel: канал связи с сервисом очередей,
    consume_tag: тег консьюмера
    """
    function unsubscribe(channel::AMQPClient.MessageChannel, consume_tag::String)::Bool
        @info Utils.Chronometer.message_with_time("Вызвана функция отписки консьюмера от очереди")
        @debug "Параметры вызова функции:"
        @debug "> канал соединения с серверов RabbitMQ " channel
        @debug "> тег консьюмера " consume_tag
        try
            return AMQPClient.basic_cancel(channel, consume_tag)
        catch error
            @error Utils.Chronometer.message_with_time("Не удалось отписаться от очереди, возникла ошибка ") error
            throw(error)
        end
    end

    """Функция очистки очереди.

    channel: канал связи с сервисом очередей,
    queue_name: имя очереди
    """
    function purge(channel::AMQPClient.MessageChannel, queue_name::String)::Dict{String, Any}
        # Вызываем функция логирования
        @info Utils.Chronometer.message_with_time("Вызвана функция очистки очериди...")
        # Логирование для режима отладки
        @debug "Параметры вызова функции:"
        @debug "> канал соединения с серверов RabbitMQ" channel
        @debug "> имя очереди" queue_name
        # Пробуем очистить очередь в канале по имени
        try
            # Если очередь очистится получим два значения - результат и
            # количество сообщений
            result, message_count = AMQPClient.queue_purge(channel, queue_name)
            # Дополняем логирование для режима отладки
            @debug ">> результат выполнения операции очистки очереди" result
            @debug ">> количество сообщений в очереди" message_count
            # вычисляем сообщение для логирования уровня info и выводим его
            message::String = result ? "Очередь успешно очищена" : "Не удалось совершить очистку очереди"
            @info string("=> ", message)
            # Возвращаем словарь с результатами работы функци, содержащий ключи:
            # STATUS - результат выполнения операции очистки очереди
            # QUEUE - наименование очереди
            # MSG_COUNT - количество сообщений
            return Dict{String, Any}(
                "STATUS" => result,
                "QUEUE" => queue_name,
                "MSG_COUNT" => message_count
            )
        catch error
            # Обрабатываем исключение. Выводим его в терминал и кидаем ошибку
            @error Utils.Chronometer.message_with_time("Очистка очереди завершилась с ошибкой ") error
            throw(error)
        end
    end

    """Функция, которая позволяет отвязать очередь от маршрута.

    channel: канал связи с сервисом очередей,
    queue_name: имя очереди,
    exchanger_name: имя обменника,
    route_name: имя маршрута
    """
    function unbind(channel::AMQPClient.MessageChannel, queue_name::String, exchanger_name::String, route_name)::Bool
        # Вызываем функця логирования
        @info Utils.Chronometer.message_with_time("Вызвана функция отвызывания от очереди...")
        # Логирование для режима отладки
        @debug "Параметры вызова функции:"
        @debug "> канал соединения с серверов RabbitMQ" channel
        @debug "> имя очереди" queue_name
        @debug "> имя обменника" exchanger_name
        @debug "> имя маршрута" route_name
        # Пытаемся отвязать маршрут от череди
        try
            # Получаем результат отвязывания маршрута от очереди
            result = AMQPClient.queue_unbind(channel, queue_name, exchanger_name, route_name)
            # Дополняем логи для уровня отладки
            @debug ">> результат отвязывания очереди, обменника и маршрута" result
            # Собираем сообщение из результат
            message::String = result ? "Отвязывание от очереди и обменника прошло успешно" : "Не удалось совершить отвязывание"
            @info string("=> ", message)
            # Возвращаем результат выполнения операции отвязывания
            return result
        catch error
            # Обрабатываем ошибку: выводим в терминал и кидаем исключение
            @error Utils.Chronometer.message_with_time("Возникла ошибка при отвязывание")
            throw(error)
        end
    end

    """Функция связывание очереди и маршрута.

    channel: канал связи с сервисом очередей,
    queue_name: имя очереди,
    exchanger_name: имя обменнника,
    route_name: имя маршрута
    """
    function bind(channel::AMQPClient.MessageChannel, queue_name::Any, exchanger_name::Any, route_name::String)::Dict{String, Any}
        # Выводим информационное сообщение в логер
        @info Utils.Chronometer.message_with_time("Вызвана функция связывания с очередью")
        # Логирование для режима отладки
        @debug "Параметры вызова функции:"
        @debug "> канал соединения с серверов RabbitMQ" channel
        @debug "> имя очереди" queue_name
        @debug "> имя обменника" exchanger_name
        @debug "> имя маршрута" route_name
        # Пытаемся связать очередь, обменник и маршрут в канале
        try
            # Запускаем операцию связывания. Её реультат всегда типа Bool
            result::Bool = AMQPClient.queue_bind(channel, queue_name, exchanger_name, route_name)
            # Дополняем логирование для режима отладки. Выводим результат операции
            @debug ">> результат связывания очереди, обменника и маршрута" result
            # Вычисляем какое сообщение показывать в информационных логах
            message::String = result ? "Связывание очереди и обменника прошло успешно" : "Не удалось совершить связывание"
            @info string("=> ", message)
            # Возвращаем результат работы функции.
            # В качестве возвращаемого значеня используем словарь (маршрут),
            # в котором содержатся следующие ключи:
            # NAME - наименование роута
            # QUEUE - имя очереди
            # ECHANGE - имя обменника
            return Dict{String, Any}(
                "NAME" => route_name,
                "QUEUE" => queue_name,
                "EXCHANGE" => exchanger_name
            )
        catch error
            # Обрабатываем ошибки: выводим в терминал и кидаем исключение
            @error Utils.Chronometer.message_with_time("Возникла ошибка при связывании")
            throw(error)
        end
    end

    """Функция регистрации очереди.

    channel: канал связи с сервисом очередей
    queue_id: идетификатор очереди
    """
    function declare(channel::AMQPClient.MessageChannel, queue_id::String)::Dict{String, Any}
        # Выводим информационное сообщение в терминал
        @info Utils.Chronometer.message_with_time("Вызвана функция регистрации очереди...")
        # Логирование для режима отладки
        @debug "Параметры вызова функции:"
        @debug "> канал соединения с серверов RabbitMQ" channel
        @debug "> идетификатор очереди" queue_id
        # Пытаемся объявить очередь
        try
            # Вызываем функцию объявления очереди. В качестве результата
            # возвращаются значения
            # result - результат выполнеиня операции
            # queue_name - имя очереди
            # message_count - количество сообщений
            # consumer_count - количество потребителей
            result, queue_name, message_count, consumer_count = AMQPClient.queue_declare(channel, queue_id)
            # Дополняем логирование для режима отладки
            @debug ">> результат операции объявления очереди" result
            @debug ">> имя очереди" queue_name
            @debug ">> количество сообщений в очереди" message_count
            @debug ">> число слушателей/потребителей" consumer_count
            # Исходя из результата вычисляем сообщение для информационного лога
            message::String = result ? "Очередь зарегистрирована удачно" : "Не удалось зарегистрировать очередь"
            @info string("=> ", message)
            # Возвращаем результат (словарь), который содержит ключи
            # STATUS - результат выполнения операции объявления
            # NAME - наименование очереди
            # MSG_COUNT - количество сообщений
            # CONS_COUNT - количество потребителей
            return Dict{String, Any}(
                "STATUS" => result,
                "NAME" => queue_name,
                "MSG_COUNT" => message_count,
                "CONS_COUNT" => consumer_count
            )
        catch error
            # Обрабатываем ошибки: выводим в терминал и кидаем исключение
            @error Utils.Chronometer.message_with_time("Регистрация очереди завершилась ошибкой ") error
            throw(error)
        end
    end

    """Функция удаления очереди.

    channel: канал связи с сервисом очередей
    queue_id: идетификатор очереди
    """
    function delete(channel::AMQPClient.MessageChannel, queue_id::Any)::Dict{String, Any}
        # Выводим информационное сообщение
        @info Utils.Chronometer.message_with_time("Запущена функуция удаления очереди...")
        # Логирование для режима отладки
        @debug "Параметры вызова функции:"
        @debug "> канал соединения с серверов RabbitMQ" channel
        @debug "> идетификатор очереди" queue_id
        #
        try
            # Выполняем функцию удаления очереди, возвращается два результата
            # result - результат операции удаления
            # message_count - количество сообщений в очереди
            result, message_count = queue_delete(channel, queue_id)
            # Дополняем логирование в режиме отладки. Выводим результаты
            @debug ">> результат выполненя операции удаления очереди" result
            @debug ">> количество сообщений в очереди" message_count
            # Собираем сообщение в заивисимости от результата
            message::String = result ? "Удаление очереди прошло успешно" : "Не удалось удалить очередь"
            @info string("=> ", message)
            # Возвращаем словарь с результатами, содержащий поля:
            # STATUS - статус выполнения операции
            # MSG_COUNT - количество сообщений
            return Dict{String, Any}(
                "STATUS" => result,
                "MSG_COUNT" => message_count
            )
        catch error
            # Обрабатываем ошибки: выводим в терминал и кидаем исключение
            @error Utils.Chronometer.message_with_time("Удаление очереди завершилось ошибкой ") error
            throw(error)
        end
    end
end


text:
using Documenter, JuliaRunClient

#push!(LOAD_PATH,"../src/")
makedocs()
#deploydocs(
#    repo = "github.com/JuliaCOmputing/JuliaRunClient.jl.git"
#)


text:
@testset "Technical Indicators" begin

    @testset "Interval, Time Period, Series Type" begin
        for ti in Symbol.(AlphaVantage.interval_timeperiod_seriestype_indicators[1:MAX_TESTS])
            @eval begin
                tiname = string($ti)
                @testset "technical_indicator: $tiname" begin
                    @testset "JSON" begin
                        data = $ti("MSFT", "weekly", 10, "open")
                        @test typeof(data) === Dict{String, Any}
                        @test length(data) === 2
                    end
                    sleep(TEST_SLEEP_TIME + 2*rand())
                    @testset "CSV" begin
                        data = $ti("MSFT", "weekly", 10, "open", datatype="csv")
                        @test typeof(data) === Tuple{Array{Any, 2}, Array{AbstractString, 2}}
                        @test length(data) === 2
                    end
                end
            end
            sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
        end
    end

    @testset "Interval, Time Period" begin
        for ti in Symbol.(AlphaVantage.interval_timeperiod_indicators[1:MAX_TESTS])
            @eval begin
                tiname = string($ti)
                @testset "technical_indicator: $tiname" begin
                    @testset "JSON" begin
                        data = $ti("MSFT", "weekly", 10)
                        @test typeof(data) === Dict{String, Any}
                        @test length(data) === 2
                    end
                    sleep(TEST_SLEEP_TIME + 2*rand())
                    @testset "CSV" begin
                        data = $ti("MSFT", "weekly", 10, datatype="csv")
                        @test typeof(data) === Tuple{Array{Any, 2}, Array{AbstractString, 2}}
                        @test length(data) === 2
                    end
                end
            end
            sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
        end
    end

    @testset "Interval, Series Type" begin
        for ti in Symbol.(AlphaVantage.interval_seriestype_indicators[1:MAX_TESTS])
            @eval begin
                tiname = string($ti)
                @testset "technical_indicator: $tiname" begin
                    @testset "JSON" begin
                        data = $ti("MSFT", "weekly", "open")
                        @test typeof(data) === Dict{String, Any}
                        @test length(data) === 2
                    end
                    sleep(TEST_SLEEP_TIME + 2*rand())
                    @testset "CSV" begin
                        data = $ti("MSFT", "weekly", "open", datatype="csv")
                        @test typeof(data) === Tuple{Array{Any, 2}, Array{AbstractString, 2}}
                        @test length(data) === 2
                    end

                end
            end
            sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
        end
    end

    @testset "Interval" begin
        for ti in Symbol.(AlphaVantage.interval_indicators[1:MAX_TESTS])
            @eval begin
                tiname = string($ti)
                @testset "technical_indicator: $tiname" begin
                    @testset "JSON" begin
                        data = $ti("MSFT", "15min")
                        @test typeof(data) === Dict{String, Any}
                        @test length(data) === 2
                    end
                    sleep(TEST_SLEEP_TIME + 2*rand())
                    @testset "CSV" begin
                        data = $ti("MSFT", "15min", datatype="csv")
                        @test typeof(data) === Tuple{Array{Any, 2}, Array{AbstractString, 2}}
                        @test length(data) === 2
                    end
                end
            end
            sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
        end
    end

    @testset "Optional Arguments" begin

        data = MACD("MSFT", "5min", "high", fastperiod = 13, slowperiod = 25)
        @test typeof(data) === Dict{String, Any}
        @test length(data) === 2
        @test data["Meta Data"]["5.2: Slow Period"] == 25
        @test data["Meta Data"]["5.1: Fast Period"] == 13

    end

end


text:
using YaoClifford
using Yao
using Test

@testset "GHZ states" begin
    n = 3
    reg1 = zero_stabilizer_state(n)
    @test measure!(reg1) == 0
    circ = chain(put(n, 1 => Yao.H), chain(cnot(i, i+1) for i = 1:n-1))
    Yao.apply!(reg1, circ)
    reg2 = zero_stabilizer_state(n)
    fidelity(reg1, reg2) ≈ sqrt(1/2)
    @test measure!(reg1) in (0b000, 0b111)
    @test fidelity(reg1, reg2) in (0, 1)
end

@testset "Reorder" begin
    reg = zero_stabilizer_state(2)
    reg |> put(2, 1 => Yao.X)
    @test measure!(invorder!(reg)) == 0b10
end

@testset "Partial trace" begin
    reg = zero_stabilizer_state(2)
    reg |> chain(put(2, 1 => H), cnot(1, 2))
    reg2 = partial_tr(reg, [1])
    @test expect(PauliString(X), reg2) == expect(PauliString(Y), reg2) == expect(PauliString(Z), reg2) == 0
end

@testset "Expectation of Hamiltonian" begin
    reg = zero_stabilizer_state(2)
    reg |> chain(put(2, 1 => H), cnot(1, 2))
    Hmts = [PauliString(X, X), PauliString(Z, Z), PauliString(Y, Y), -2im*PauliString(I2, I2)]
    @test [expect(Hmt, reg) for Hmt in Hmts] == [1, 1, -1, -2im]
    rand_reg = rand_stabilizer_state(2)
    @test sum(expect(Hmt, rand_reg) for Hmt in Hmts) == expect(sum(Hmts), rand_reg)
end


text:
# IMPORT CONTROL SYSTEMS
using ControlSystems
# LAPLACE OPERATOR
s = tf("s")
# SYSTEM DEFINITION
sys = 2/(s+1)
# IMPORT PLOTS
using Plots
# PLOT BODE FREQUENCY RESPONSE
bodeplot(sys)
# PLOT NICHOLS FREQUENCY RESPONSE
nicholsplot(sys)
# PLOT NYQUIST FREQUENCY RESPONSE
nyquistplot(sys)

max_stars_repo_name: a-mhamdi/Eg-Softs

max_issues_repo_name: a-mhamdi/Eg-Softs

max_forks_repo_name: a-mhamdi/Eg-Softs


text:
"""
#/* compute Laplace coefficients and Leverrier derivative/
#          j
#     j   d     i
#    a   ---   b (a)
#          j    s
#        da
#
#   by series summation */
#/* Code due to Jack Wisdom */
"""
function laplace_wisdom(s::Rational,i::Integer,j::Integer,a::T) where T<:Real
# function laplace_wisdom,s,i,j,a  IDL
# double laplace(double s, int i, int j, double a);  c

##define LAPLACE_EPS 1.0e-12
  LAPLACE_EPS = convert(T,1.0e-12)

  #if (i lt 0) then i = -i
  i=abs(i)

  if j <= i     #/* compute first term in sum */
    factor4 = one(T)
    for k=0:j-1
      factor4 *= i - k
    end
    lap_coef_sum = factor4
    q0 = 0
  else
    q0 = fld(j + 1 - i,2) #largest integer less than or equal to x/y
    lap_coef_sum = zero(T)
    factor4 = one(T)
  end

  #  /* compute factors for terms in lap_coef_sum */

  factor1 = s
  factor2 = s + i
  factor3 = i + 1
  for q=1:q0-1       #/* no contribution for q = 0 */
    factor1 *= s + q
    factor2 *= s + i + q
    factor3 *= i + q + 1
  end
  if q0 > 1
    q=q0
  else
    q=1
  end
  #println(j+1-i,q0)
  term = a*a * factor1 * factor2 / (factor3 * q)

  #  /* sum series */

  while  (term*factor4) > LAPLACE_EPS
    factor4 = one(T)
    for k=0:j-1
      factor4 *= (2*q + i - k)
    end
    lap_coef_sum += term * factor4
    factor1 += 1
    factor2 += 1
    factor3 += 1
    q = q+1
    term *= a*a * factor1 * factor2 / (factor3 * q)
  end

  #  /* fix coefficient */

  for k=0:i-1
    lap_coef_sum *= (s+k)/(k+1)
  end

  apower = (q0 <= 0) ?  i : 2*q0 + i - 2
  lap_coef_sum *= 2 * a^apower
  # Return the Laplace Coefficient:
  return lap_coef_sum
end

max_stars_repo_name: bmlindor/TTVFaster

max_issues_repo_name: bmlindor/TTVFaster

max_forks_repo_name: bmlindor/TTVFaster


text:
############################################################################################
"""
$(TYPEDEF)

Default arguments for MCMC constructor.

# Fields
$(TYPEDFIELDS)
"""
struct MCMCDefault{K<:NamedTuple, S<:ConfigStepsize, P<:ConfigProposal, G, U<:BaytesCore.UpdateBool}
    "Individual keyword arguments for tuning different MCMC kernels."
    kernel::K
    "Stepsize default configuration for adaption."
    stepsize::S
    "Proposal distribution default configuration for adaption."
    proposal::P
    "Gradient backend used in MCMC step. Not used if Metropolis sampler is chosen."
    GradientBackend::G
    "Boolean if initial parameter are fixed or resampled."
    TunedModel::Bool
    "Boolean if generate(_rng, objective) for corresponding model is stored in MCMC Diagnostics."
    generated::U
    function MCMCDefault(;
        kernel=(;),
        stepsize = ConfigStepsize(),
        proposal = ConfigProposal(),
        GradientBackend=:ForwardDiff,
        TunedModel=true,
        generated=BaytesCore.UpdateFalse()
    )
        ArgCheck.@argcheck (
            isa(GradientBackend, Symbol) || isa(GradientBackend, AnalyticalDiffTune)
        ) "GradientBackend keywords has to be either an AD symbol (:ForwardDiff, :ReverseDiff, :ReverseDiffUntaped, :Zyogte), or an AnalyticalDiffTune object."
        return new{
            typeof(kernel), typeof(stepsize), typeof(proposal), typeof(GradientBackend), typeof(generated)
        }(
            kernel, stepsize, proposal, GradientBackend, TunedModel, generated
        )
    end
end

############################################################################################
"""
$(TYPEDEF)

Stores information for proposal step.

# Fields
$(TYPEDFIELDS)
"""
struct MCMC{M<:MCMCKernel,N<:MCMCTune} <: AbstractAlgorithm
    "MCMC sampler"
    kernel::M
    "Tuning configuration for kernel."
    tune::N
    function MCMC(kernel::M, tune::N) where {M<:MCMCKernel,N<:MCMCTune}
        return new{M,N}(kernel, tune)
    end
end

############################################################################################
function MCMC(
    _rng::Random.AbstractRNG,
    kernel::Type{M},
    objective::Objective,
    default::MCMCDefault=MCMCDefault(),
    info::BaytesCore.SampleDefault = BaytesCore.SampleDefault()
) where {M<:MCMCKernel}
    ## Checks before algorithm is initiated
    @unpack output = objective.model.info.reconstruct.default
    @unpack GradientBackend, TunedModel, generated = default
    ## Sample from prior if TunedModel == false
    if !TunedModel
        sample!(objective.model, objective.tagged)
    end
    @unpack model, data, tagged = objective
    ## Initialize default configuration for chosen algorithm
    config = init(AbstractConfiguration, kernel, objective, default.proposal; default.kernel...)
    ##	If a valid AD backend is provided, change it to an AutomaticDifftune Object
    if isa(GradientBackend, Symbol)
        GradientBackend = AutomaticDiffTune(GradientBackend, objective)
    end
    ## Initiate MCMC Algorithm
    mcmc = init(kernel, config, objective, GradientBackend)
    ## Initial Phase tune
    phasetune = PhaseTune(;
        window=config.window.window, buffer=config.window.buffer, phasenames=config.window.phasenames
    )
    ## Initial Proposal tune
    proposal = Proposal(
        output,
        default.proposal.proposaladaption,
        MatrixTune(default.proposal.metric, default.proposal.shrinkage),
        length(objective.tagged),
        phasetune.iterations[phasetune.counter.current],
    )
    ## Tune stepsize in case it is not assigned via local acceptance rate bisection
    if default.stepsize.initialstepsize isa BaytesCore.UpdateTrue
        acceptrateᵗᵉᵐᵖ = get_acceptrate(_rng, mcmc, objective, proposal.Σ) # Functor to calculate local accpetance rate
        ϵ = find_initial_stepsize(
            InitialStepsizeSearch(output), acceptrateᵗᵉᵐᵖ
        )
    else
        ϵ = output(default.stepsize.ϵ)
    end
    stepsize = StepSizeTune(default.stepsize.stepsizeadaption, init(DualAverage, config.δ, ϵ), ϵ)
    ## Initial MCMC Tune struct, and return MCMC container
    mcmctune = MCMCTune(objective, phasetune, stepsize, proposal, generated)
    ## Return MCMC container
    return MCMC(mcmc, mcmctune)
end

############################################################################################
"""
$(SIGNATURES)
Propose new parameter with mcmc sampler. If update=true, objective function will be updated with input model and data.

# Examples
```julia
```

"""
function propose(_rng::Random.AbstractRNG, mcmc::MCMC, objective::Objective)
    #!NOTE: Temperature is fixed for propose() step and will not be adjusted
    ## Make MCMC Proposal step
    resultᵖ, divergent, accept, kernel_diagnostics = propagate(
        _rng, mcmc.kernel, mcmc.tune, objective
    )
    ## If accepted, update kernel and Model parameter
    if accept.accepted
        mcmc.kernel.result = resultᵖ
        ModelWrappers.unflatten_constrain!(objective.model, mcmc.tune.tagged, resultᵖ.θᵤ)
    end
    ## Upate tuning container - includes storing current parameter, updating sampling phase, discretization steps, and proposal distribution (in that order)
    update!(mcmc.tune, mcmc.kernel.result, accept.rate)
    diagnostics = MCMCDiagnostics(
        BaytesCore.BaseDiagnostics(
            mcmc.kernel.result.ℓθᵤ,
            objective.temperature,
            ModelWrappers.predict(_rng, mcmc, objective),
            mcmc.tune.iter.current
        ),
        kernel_diagnostics,
        divergent,
        accept,
        ModelWrappers.generate(_rng, objective, mcmc.tune.generated)
    )
    return objective.model.val, diagnostics
end

############################################################################################
"""
$(SIGNATURES)
Inplace version of propose.

# Examples
```julia
```

"""
function propose!(
    _rng::Random.AbstractRNG,
    mcmc::MCMC,
    model::ModelWrapper,
    data::D,
    temperature::F = model.info.reconstruct.default.output(1.0),
    update::U=BaytesCore.UpdateTrue(),
) where {D,F<:AbstractFloat, U<:BaytesCore.UpdateBool}
    ## Update Objective with new model parameter from other MCMC samplers and/or new/latent data
    objective = Objective(model, data, mcmc.tune.tagged, temperature)
    update!(mcmc.kernel, objective, update) #Update Kernel with current objective/configs
    ## Compute MCMC step
    val, diagnostics = propose(_rng, mcmc, objective)
    ## If accepted, update kernel and Model parameter
    if diagnostics.accept.accepted
        model.val = val
    end
    return val, diagnostics
end

############################################################################################
# Export
export MCMC, MCMCDefault, propose, propose!


text:
using Random
using Plots
using Distances
using ProgressMeter

function time_driven(nstep)

    np = 50
    rng = MersenneTwister(1234)
    q = -0.5 .+ rand(2,np)
    v = randn(2,np)

    dt = 0.01

    bounds = [-2, 2, -2, 2]
    ϵ = 0.04

    bar = Progress(nstep)

    anim = @animate for step in 1:nstep

        next!(bar)

        q .+= dt .* v

        d = pairwise(Euclidean(), q)
        p = findall(x -> x < 2ϵ, d)

        ind1 = getindex.(p,1)
        ind2 = getindex.(p,2)

        unique = ind1 .< ind2

        ind1 = ind1[unique]
        ind2 = ind2[unique]

        for (i1, i2) in zip(ind1, ind2)

            q1 = q[:, i1]
            q2 = q[:, i2]

            v1 = v[:, i1]
            v2 = v[:, i2]

            q_rel = q1 - q2
            v_rel = v1 - v2

            v_cm = 0.5 * (v1+v2)

            qq_rel = q_rel'q_rel
            vq_rel = v_rel'q_rel
            v_rel = 2 * q_rel * vq_rel / qq_rel - v_rel

            v[:, i1] .= v_cm + v_rel * 0.5
            v[:, i2] .= v_cm - v_rel * 0.5

        end

        crossed_x1 = (q[1,:] .< (bounds[1] + ϵ))
        crossed_x2 = (q[1,:] .> (bounds[2] - ϵ))
        crossed_y1 = (q[2,:] .< (bounds[3] + ϵ))
        crossed_y2 = (q[2,:] .> (bounds[4] - ϵ))

        q[1, crossed_x1] .= bounds[1] + ϵ
        q[1, crossed_x2] .= bounds[2] - ϵ

        q[2, crossed_y1] .= bounds[3] + ϵ
        q[2, crossed_y2] .= bounds[4] - ϵ

        v[1, crossed_x1] .*= -1
        v[1, crossed_x2] .*= -1
        v[2, crossed_y1] .*= -1
        v[2, crossed_y2] .*= -1

        plt = plot(size  = (200,200),
                   xlims = (-2,2),
                   ylims = (-2,2),
                   grid  = false,
                   axis  = nothing,
                   framestyle= :box,
                   legend=false,
                   widen = false)

        c = trunc(Int, 100ϵ)
        scatter!( view(q,1,:), view(q,2,:), markershape  = :circle,
              markersize   = c ,
              aspect_ratio = :equal)

    end every 10

    gif(anim, joinpath(@__DIR__, "time_driven.gif"), fps = 10)

end


text:
module SuperGLUE
using Pkg.Artifacts

function download()
    artifact_toml = joinpath(@__DIR__, "Artifacts.toml")
    superglue_hash = artifact_hash("superglue", artifact_toml)
    if superglue_hash === nothing || !artifact_exists(superglue_hash)
        superglue_hash = create_artifact() do artifact_dir
            path = tempdir()*"/combined.zip"
            Base.download("https://dl.fbaipublicfiles.com/glue/superglue/data/v2/combined.zip", path)
            run(`unzip -x $path -d $artifact_dir`)
            rm(path)
        end
        bind_artifact!(artifact_toml, "superglue", superglue_hash)
    end
    return artifact_path(superglue_hash)
end

end


text:
function make_mask(
    Nside, θpix;
    edge_padding_deg = 2,
    edge_rounding_deg = 1,
    apodization_deg = 1,
    ptsrc_radius_arcmin = 7, # roughly similar to Story et al. 2015
    num_ptsrcs = round(Int,((Nside*θpix)/60)^2 * 120/100) # SPT-like density of sources
    )

    deg2npix(x) = round(Int, x/θpix*60)
    arcmin2npix(x) = round(Int, x/θpix)

    ptsrc = .!bleed(sim_ptsrcs(Nside, num_ptsrcs), arcmin2npix(ptsrc_radius_arcmin))
    boundary = boundarymask(Nside, deg2npix(edge_padding_deg))
    mask_array = if apodization_deg in (false, 0)
        boundary .& ptsrc
    else
        cos_apod(boundary, deg2npix(apodization_deg), deg2npix(edge_rounding_deg)) .* cos_apod(ptsrc, arcmin2npix(ptsrc_radius_arcmin));
    end

    FlatMap(mask_array, θpix=θpix)
end

make_mask(::FlatField{<:Flat{Nside,θpix}}; kwargs...) where {Nside,θpix} = make_mask(Nside, θpix; kwargs...)

# all padding/smoothing/etc... quantities below are in units of numbers of pixels

function boundarymask(Nside, pad)
    m = fill(true,Nside,Nside)
    m[1:pad,:]          .= false
    m[:,1:pad]          .= false
    m[end-pad+1:end,:]  .= false
    m[:,end-pad+1:end]  .= false
    m
end

function bleed(img, w)
    Nside,Nside = size(img)
    nearest = getfield.(@ondemand(Images.feature_transform)(img),:I)
    [norm(nearest[i,j] .- [i,j]) < w for i=1:Nside,j=1:Nside]
end

function cos_apod(img, w, smooth_distance=false)
    Nside,Nside = size(img)
    nearest = getfield.(@ondemand(Images.feature_transform)(.!img),:I)
    distance = [norm(nearest[i,j] .- [i,j]) for i=1:Nside,j=1:Nside]
    if smooth_distance!=false
        distance = @ondemand(Images.imfilter)(distance, @ondemand(Images.Kernel.gaussian)(smooth_distance))
    end
    @. (1-cos(min(distance,w)/w*pi))/2
end

function round_edges(img, w)
    .!(@ondemand(Images.imfilter)(img, @ondemand(Images.Kernel.gaussian)(w)) .< 0.5)
end

function sim_ptsrcs(Nside,nsources)
    m = fill(false,Nside,Nside);
    for i=1:nsources
        m[rand(1:Nside),rand(1:Nside)] = true
    end
    m
end


text:
using JuMP, Ipopt

@testset "Rosenbrock Ipopt" begin
    # min (1-x)^2 + 100*(y-x^2)^2
    # solution: (x,y) = (1,1)
    # optimal objective 0
    rosenbrock = Problem()

    x = MPC.Variable("x", Real)
    y = MPC.Variable("y", Real)
    set_objective!(rosenbrock, Base.power_by_squaring((1-x),2) + 100*Base.power_by_squaring((y-Base.power_by_squaring(x,2)),2))

    mysolver = Ipopt.Optimizer
    m, variables_jump, ctr_jump, ctr_exp = get_JuMP_cartesian_model(rosenbrock, mysolver)
    optimize!(m)

    sol = get_JuMP_solution(m, variables_jump, rosenbrock)

    @test sol[x] ≈ 1.0 atol=1e-5
    @test sol[y] ≈ 1.0 atol=1e-5
    @test get_objective(rosenbrock, sol) ≈ 0.0 atol=1e-5
end

@testset "QCQP Ipopt" begin
    # min x - y
    # st  x + x^2 + x*y + y^2 <= 1
    #     -2 <= x, y <= 2
    # solution: x+y = -1/3
    # optimal objective -1-4/sqrt(3)

    qcqp = Problem()

    x = MPC.Variable("x", Real)
    y = MPC.Variable("y", Real)

    set_objective!(qcqp, x - y)

    add_constraint!(qcqp, "ctr", (x + Base.power_by_squaring(x, 2) + x*y + Base.power_by_squaring(y,2)) << 1)
    add_constraint!(qcqp, "x bounds", -2 << x << 2)
    add_constraint!(qcqp, "y bounds", -2 << y << 2)

    mysolver = Ipopt.Optimizer
    m, variables_jump, ctr_jump, ctr_exp = get_JuMP_cartesian_model(qcqp, mysolver)
    optimize!(m)

    sol = get_JuMP_solution(m, variables_jump, qcqp)

    ## Test for optimal solution
    @test sol[x] + sol[y] ≈ -1/3 atol=1e-5
    @test get_objective(qcqp, sol) ≈ -1-4/sqrt(3) atol=1e-5

    ## Test slack functions
    slacks = get_slacks(qcqp, sol)
    @test real(slacks.coords[MPC.Variable("x bounds", Complex)]) > 0
    @test real(slacks.coords[MPC.Variable("y bounds", Complex)]) > 0

    @test imag(slacks.coords[MPC.Variable("x bounds", Complex)]) == 0
    @test imag(slacks.coords[MPC.Variable("y bounds", Complex)]) == 0

    @test get_minslack(qcqp, sol)[1] < 1e-6
end

@testset "Matpower case9 Ipopt" begin
    instancepath = joinpath(pwd(), "test", "instances")
    # instancepath = joinpath(Pkg.dir("MathProgComplex"), "test", "instances")
    pb_real, ~ = import_from_dat(joinpath(instancepath,"case9real.dat"))

    mysolver = Ipopt.Optimizer
    m, variables_jump = get_JuMP_cartesian_model(pb_real, mysolver)
    optimize!(m)
    @test JuMP.objective_value(m) ≈ 1.45883471040144e+003 atol=1e-6
end


text:
type CuDNNSoftmaxState
  inputs_desc  :: Vector{CuDNN.Tensor4dDescriptor}
  outputs_desc :: Vector{CuDNN.Tensor4dDescriptor}
end

function setup_etc(backend::GPUBackend, layer::SoftmaxLayer, dims::Vector{Int}, data_type, inputs)
  inputs_desc = Array{CuDNN.Tensor4dDescriptor}(length(inputs))
  outputs_desc = Array{CuDNN.Tensor4dDescriptor}(length(inputs))
  for i = 1:length(inputs)
    dim_sp, dim_prob, dim_num = split_dims(inputs[i], dims[i])
    inputs_desc[i] = CuDNN.create_tensor4d_descriptor(data_type, (1,dim_sp,dim_prob,dim_num))
    outputs_desc[i] = CuDNN.create_tensor4d_descriptor(data_type, (1,dim_sp,dim_prob,dim_num))
  end
  etc = CuDNNSoftmaxState(inputs_desc, outputs_desc)
  return etc
end
function shutdown(backend::GPUBackend, state::SoftmaxLayerState)
  map(destroy, state.blobs)
  map(CuDNN.destroy_tensor4d_descriptor, state.etc.inputs_desc)
  map(CuDNN.destroy_tensor4d_descriptor, state.etc.outputs_desc)
end

function forward(backend::GPUBackend, state::SoftmaxLayerState, inputs::Vector{Blob})
  alpha = one(eltype(inputs[1]))
  beta = zero(eltype(inputs[1]))
  for i = 1:length(inputs)
    CuDNN.softmax_forward(backend.cudnn_ctx, CuDNN.CUDNN_SOFTMAX_ACCURATE,
        CuDNN.CUDNN_SOFTMAX_MODE_CHANNEL, alpha, state.etc.inputs_desc[i], inputs[i].ptr,
        beta, state.etc.outputs_desc[i], state.blobs[i].ptr)
  end
end

function backward(backend::GPUBackend, state::SoftmaxLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
  alpha = one(eltype(inputs[1]))
  beta = zero(eltype(inputs[1]))

  for i = 1:length(inputs)
    diff = diffs[i]
    if !isa(diff, NullBlob)
      CuDNN.softmax_backward(backend.cudnn_ctx, CuDNN.CUDNN_SOFTMAX_ACCURATE,
          CuDNN.CUDNN_SOFTMAX_MODE_CHANNEL, alpha, state.etc.outputs_desc[i], state.blobs[i].ptr,
          state.etc.outputs_desc[i], state.blobs_diff[i].ptr, beta, state.etc.inputs_desc[i], diff.ptr)
    end
  end
end


text:
# Density test for manipulation in the running variable.
# This work is based on McCrary, J. (2008). Manipulation of the running variable in the regression discontinuity design: A density test. Journal of econometrics, 142(2), 698-714.

"""
    McCraryTest(bin::Union{Real,Nothing} = nothing, bw::Union{Real,Nothing} = nothing)

A McCrary (2008) test for manipulation in the running variable for `SharpRD` estimator.

The `bin` object is the width of each bin. The default is 2*sd(runvar)*length(runvar)^(-0.5), following McCrary 2008.

The `bw` object is the bandwidth. The default uses the calculation from McCrary (2008), pg. 705.
"""
Base.@kwdef struct McCraryTest{K}
    bin::Union{Real,Nothing} = nothing
    bw::Union{Real,Nothing} = nothing
    kernel::K = SymTriangularDist()
end

Base.@kwdef struct FittedMcCraryTest{T,M,V,C}
    b::T
    h::T
    θhat::T
    σθ::T
    pval::T
    Js::M
    fitted_res::V
    coeftable::C
end

# 1. Descriptive Statistics

function descriptive_stats(RunVar::RunningVariable)
    n = Base.length(RunVar)
    sd = StatsBase.std(RunVar)
    rmax = Base.maximum(RunVar)
    rmin = Base.minimum(RunVar)
    return n, sd, rmax, rmin
end

# 2. Binsize (uses descriptive_stats)
# function Binsize(n, sd, bin::Union{Real,Nothing}=nothing)
function bin_width(n, sd)
    b = 2 * sd * n^(-.5)
end

function bin_width(Zs::RunningVariable)
    n, sd, rmax, rmin = descriptive_stats(Zs)
    return b = 2 * sd * (n)^(-.5)
end

# 3. Number of bins
function bin_number(c, b, rmin, rmax)
    Jl = convert(Int64, round(floor((c - rmin) / b) + 1))
    J = convert(Int64, floor((rmax - rmin) / b) + 2)
    Jr = J - Jl
    return J, Jl, Jr
end

function bin_number(Zs::RunningVariable)
    n, sd, rmax, rmin = descriptive_stats(Zs)
    b = bin_width(Zs)
    c = Zs.cutoff
    Jl = convert(Int64, round(floor((c - rmin) / b) + 1))
    J = convert(Int64, floor((rmax - rmin) / b) + 2)
    Jr = J - Jl
    return J, Jl, Jr
end

function bin_number(method::McCraryTest, Zs::RunningVariable)
    n, sd, rmax, rmin = descriptive_stats(Zs)

    if method.bin === nothing
        b = bin_width(Zs)
    else
        b = method.bin
    end
    c = Zs.cutoff
    Jl = convert(Int64, round(floor((c - rmin) / b) + 1))
    J = convert(Int64, floor((rmax - rmin) / b) + 2)
    Jr = J - Jl
    return J, Jl, Jr
end

# 4. Kernel Constant for McCrary
### According to footnote 18 of McCrary (2008), this constant cannot be computed by the standard formula in the book (Eq 4.3) because is a boundary case. What I understand is that the constant is calculted from 3.20 and 3.22 of Fan and Gijbels.

function kernel_constant(method::McCraryTest)
    # Here I need to see how to calculate this. For now its output is just for the Triangular Kernel
    # kernel = EquivalentKernel(method.kernel)

    κ = 3.348
end

# 5. Bandwidth, needs (4), (3), (2), (1)
# function bandwidth(mccrary::McCrary, kernel::SymTriangularDist, X, Y, Jl, c)

function bandwidth(method::McCraryTest, Zs::RunningVariable)
        # I will proceed to estimate the optimal bandwidth using McCrary (2008) optimal method.
    κ = kernel_constant(method)
    c = Zs.cutoff

    # Discretized Running Variable
    n, sd, rmin, rmax = descriptive_stats(Zs)

    #### NEED TO ADAPT THIS
    if method.bin === nothing
        b = bin_width(Zs)
    else
        b = method.bin
    end

    J, Jl, Jr = bin_number(method, Zs)

    DRV = DiscretizedRunningVariable(Zs, J + 1, b)
    X = collect(DRV)
    Y = DRV.weights / (n * b)

    # Left part
    ldf = DataFrame(Y=Y[1:Jl], X=X[1:Jl])
    lols = lm(@formula(Y ~ X + X^2 + X^3 + X^4 ), ldf)
    lcoef = coef(lols)
    lmse =  deviance(lols) / dof_residual(lols)

    lf = 2 * lcoef[3] .+ 6 * lcoef[4] * ldf.X + 12 * lcoef[5] * ldf.X.^2
    lh = κ * (lmse * (c - ldf.X[1]) / sum(lf'lf))^(0.2)

    # Right part
    rdf = DataFrame(Y=Y[Jl + 1:end], X=X[Jl + 1:end])
    rols = lm(@formula(Y ~ X + X^2 + X^3 + X^4 ), rdf)
    rcoef = coef(rols)
    rmse =  deviance(rols) / dof_residual(rols)

    rf = 2 * rcoef[3] .+ 6 * rcoef[4] * rdf.X + 12 * rcoef[5] * rdf.X.^2
    rh = κ * (rmse * (rdf.X[end] - c) / sum(rf'rf))^(0.2)

    # Taking the average
    h = 0.5 * (rh + lh)
end

# 6. McCraryTest
"""
    fit(method::McCraryTest, data::RunningVariable)

Test for manipulation in the running variable following McCrary (2008).

# Arguments

- `data::RunningVariable`: The running variable.

# Returns
Returns  a FittedMcCraryTest object with fields `b` and `h` which are the bin size and the bandwidth used, respectively. Also, it returns `θhat`, `σθ`, `pval`, `Js`, and `fitted_res` which are the estimated difference in the intercepts, the standard error of the estimate, the pval of the test H0: θhat = 0, the number of bins used, and the fitted values for the support of the running variable (as a DataFrame). Finally, `coeftable` provides an object of type CoefTable for the estimator.
"""
function fit(method::McCraryTest, Zs::RunningVariable)

    # Cutoff selection
    c = Zs.cutoff
    R = Zs
    kernel(t) = pdf(method.kernel, t)

    # Call the functions
    # 1.
    n, sd, rmax, rmin = descriptive_stats(Zs)

    # 2.
    if method.bin === nothing
        b = bin_width(Zs)
    else
        b = method.bin
    end

    # 3.
    J, Jl, Jr = bin_number(method, Zs)
    Js = [J; Jl; Jr]

    # 4.
    DRV = DiscretizedRunningVariable(Zs, J + 1, b)
    X = collect(DRV)
    Y = DRV.weights / (n * b)

    # 5.
    if method.bw === nothing
        h = bandwidth(method, Zs)
    else
        h = method.bw
    end

    ##############################################################
    # Smooth histograms to the left and to the right of the cutoff
    ldf = DataFrame(Y=Y[1:Jl], X=X[1:Jl], pred=zeros(Jl))
    rdf = DataFrame(Y=Y[Jl + 1:end], X=X[Jl + 1:end], pred=zeros(Jr))

    # Left side
    Sminus(k) = sum(kernel.((ldf.X .- c) ./ h) .* (ldf.X .- c).^k)
    lfhat = sum(kernel.((ldf.X .- c) ./ h) .* ((Sminus(2) .- Sminus(1) .* (ldf.X .- c)) ./ (Sminus(2) * Sminus(0) - (Sminus(1))^2 )) .* ldf.Y)

    # Right side
    Splus(k) = sum(kernel.((rdf.X .- c) ./ h) .* (rdf.X .- c).^k)
    rfhat = sum(kernel.((rdf.X .- c) ./ h) .* ((Splus(2) .- Splus(1) .* (rdf.X .- c)) ./ (Splus(2) * Splus(0) - (Splus(1))^2)) .* rdf.Y)

    # The estimates:
    θhat = log(rfhat) - log(lfhat)
    σθ = sqrt((1 / (n * h)) * (24 / 5) * ( (1 / rfhat) + (1 / lfhat)))
    z = θhat / σθ
    pval = 2(1 -  cdf(Normal(), abs(z)))

    # Creating the table
    res = [θhat σθ z pval]
    colnms = ["θ̂"; "σ̂"; "z"; "p-val"]
    rownms = ["McCrary Test"]
    coeftbl = CoefTable(res, colnms, rownms, 4, 3)

    # Fitted Values:
        # Left side
    for i in eachindex(ldf.X)
        ldf[!,:dist] = ldf.X  .-  ldf.X[i]
        ldf.pred[i] = predict(lm(@formula(Y ~ dist), ldf, wts=kernel.(ldf.dist / h)))[i]
    end

    m = min.((c .- ldf.X) / h, 1)

    lVarf = (12 / (5 * n * h)) .* ldf.Y .* (2.0 .+ 3 * m.^11 .- 24 * m.^10 .+ 83 * m.^9 .- 72 * m.^8 .- 42 * m.^7 .+ 18 * m.^6 .+ 18 * m.^5 .+ 18 * m.^4 .+ 3 * m.^3 .+ 18 * m.^2 .+ 15 * m) ./ ((1.0 .+ m.^6 .- 6 * m.^5 .- 3 * m.^4 .+ 4  * m.^3 .+ 9 * m.^2 .+ 6 * m))

    ldf[!, :se_pred] = sqrt.(lVarf)

        # Right side
    for i in eachindex(rdf.X)
        rdf[!,:dist] = rdf.X  .- rdf.X[i]
        rdf.pred[i] = predict(lm(@formula(Y ~ dist), rdf, wts=kernel.(rdf.dist / h)))[i]
    end

    m = max.(-1, (-rdf.X .+ c) / h)

    rVarf = (12 / (5 * n * h)) .* rdf.Y .* (2.0 .- 3 * m.^11 .- 24 * m.^10 .- 83 * m.^9 .- 72 * m.^8 .+ 42 * m.^7 .+ 18 * m.^6 .- 18 * m.^5 .+ 18 * m.^4 .- 3 * m.^3 .+ 18 * m.^2 .- 15 * m) ./ ((1.0 .+ m.^6 .+ 6 * m.^5 .- 3 * m.^4 .- 4  * m.^3 .+ 9 * m.^2 .- 6 * m).^2 )

    rdf[!, :se_pred] = sqrt.(rVarf)

fitted_res = DataFrame(Y=Y, X=X, fhat=vcat(ldf.pred, rdf.pred), se_fhat=vcat(ldf.se_pred, rdf.se_pred))

FittedMcCraryTest(
        b=b,
        h=h,
        θhat=θhat,
        σθ=σθ,
        pval=pval,
        Js=Js,
        fitted_res=fitted_res,
        coeftable=coeftbl
    )

end
function Base.show(io::IO, rdd_fit::FittedMcCraryTest)
    println(io, "The McCrary (2008) test for manipulation in the")
    println(io, "running variable for RDD.")
    println(io, "          ⋅⋅⋅⋅ " * "Bin size: ", string(round(rdd_fit.b, digits=4)))
    println(io, "          ⋅⋅⋅⋅ " * "Bandwidth size: ", string(round(rdd_fit.h, digits=4)))
    println(io, "          ⋅⋅⋅⋅ " * "Number of iterations: ", string(rdd_fit.Js[1]))
    Base.show(io, rdd_fit.coeftable)
end

# 7. Plot

@recipe function f(rdd_fit::FittedMcCraryTest)
    # Recovering the data
    X = rdd_fit.fitted_res[!,:X]
    Y = rdd_fit.fitted_res[!,:Y]
    f̂ = rdd_fit.fitted_res[!,:fhat]
    Jl = rdd_fit.Js[2]
    lX = X[1:Jl]
    rX = X[Jl + 1:end]
    lf̂ = f̂[1:Jl]
    rf̂ = f̂[Jl + 1:end]

    # set up
    legend := :false
    grid := false

    # scatter
    @series begin
        seriestype := :scatter
        seriescolor := :lightgray
        linecolor --> :lightgray
        seriesalpha := 0.99
        X, Y
    end

    # common to both series
    seriestype := :line
    linecolor --> :firebrick
    linewidth --> 2
    label --> "McCrary Test"

    # left side
    @series begin
        lX, lf̂
    end

    # left side
    @series begin
        rX, rf̂
end
end


text:
#struct Vertex { union { struct { _float x, y, z; }; _float pos[3]; }; };

type ShapeCoordinates
	#std::string& name
	positions::AVec3f
	normals::AVec3f
	texcoords::AVec2f
	origin::Vec3f
	scale::Vec3f

	ShapeCoordinates() = new(AVec3f(),AVec3f(),AVec1f(),Vec3f(),Vec3f())
	ShapeCoordinates(positions,normals,texcoords,origin,scale) = new(positions,normals,texcoords,origin,scale)
end

type VertexIndex
	v::Integer
	vt::Integer
	vn::Integer

	VertexIndex() = new(0,0,0)
	VertexIndex(a::Integer) = new(a,a,a)
	VertexIndex(v::Integer, vn::Integer, vt::Integer) = new(v,vt,vn)
end

const AVertexIndex = Array{VertexIndex,1}
const AAVertexIndex = Array{AVertexIndex,1}

function clear(shape::ShapeCoordinates)
	shape.positions = AVec3f()
	shape.normals = AVec3f()
	shape.texcoords = AVec1f()
	shape.origin = Vec3f()
	shape.scale = Vec3f()
end

# for std::map
function operator<(a::VertexIndex, b::VertexIndex)
	if a.v != b.v return (a.v < b.v) end
	if a.vn != b.vn return (a.vn < b.vn) end
	if a.vt != b.vt return (a.vt < b.vt) end
	false
end

isSpace(_char_ c) = (c == ' ') || (c == '\t')
isNewLine(_char_ c) = (c == '\r') || (c == '\n') || (c == '\0')

function getOrigin(min_BBOX::Vec3f, max_BBOX::Vec3f, scale::Float32)
	(max_BBOX + min_BBOX) * Vec3f(scale * 0.5f)
end

function getScale(min_BBOX::Vec3f, max_BBOX::Vec3f)
	size = max_BBOX - min_BBOX
	scale_value = size.x > size.z ? size.x : size.z
	scale_value = scale_value > size.y ? scale_value : size.y
	if scale_value > 0 scale_value = 1 / scale_value end
	scale_value
end

function parseFloat(token::String)
	token += strspn(token, " \t")
	f = (Float32)atof(token)
	token += strcspn(token, " \t\r")
	f
end

function parsevec(token::String, count::Number)
	a = zeros(count)
	for i=1:count a[i] = parseFloat(token) end
	Vec{count, Float32}(a)
end

# Make index zero-base, and also support relative index.
function fixIndex(idx::Integer, n::Integer)
	if idx > 0 i = idx - 1
	else if idx == 0  i = 0
	else i = n + idx # negative value = relative
	end
	i
end

# Parse triples: i, i/j/k, i//k, i/j
function parseTriple(token::String, vsize::Integer, vnsize::Integer, vtsize::Integer)
	vi = VertexIndex(-1)

	vi.v = fixIndex(atoi(token), vsize)
	token += strcspn(token, "/ \t\r");
	if token[0] != '/' return vi end
	++token;

	# i//k
	if token[0] == '/'
		++token
		vi.vn = fixIndex(atoi(token), vnsize)
		token += strcspn(token, "/ \t\r")
		return vi
	end

	# i/j/k or i/j
	vi.vt = fixIndex(atoi(token), vtsize)
	token += strcspn(token, "/ \t\r")
	if (token[0] != '/') return vi; else ++token;  # skip '/'
	vi.vn = fixIndex(atoi(token), vnsize)
	token += strcspn(token, "/ \t\r")
	vi
end

function updateVertex(vi::VertexIndex, vCache::Dict{VertexIndex, _uint},
	coord::ShapeCoordinates, shape::Mesh)

	idx = 0

	it = vCache.find(vi)

	if it != vCache.end() return it->second end

	if vi.v >= length(coord.position)
		error("_DEBUG_BREAK_IF(((_uint)(vi.v) >= (_uint)(coord.positions.size())))")
		return
	end

	shape->addPosition(coord.positions[vi.v] * coord.scale - coord.origin)

	if vi.vn >= 0
		shape->addNormal(coord.normals[vi.vn])
	end

	if vi.vt >= 0
		shape->addTexcoord(coord.texcoords[vi.vt])
	end

	idx = shape->getPositionSize() - 1
	vCache[vi] = idx

	idx;
end

function parseIndicies(
	std::map<VertexIndex, _uint> vCache,
	const std::vector<std::vector<VertexIndex> >& group,
	#const int material_id,
	#const std::string &name,
	#bool clearCache,
	ShapeCoordinates in,
	shape::Mesh)

	clear(shape) # remove old

	if empty(group) return false end

	# Flatten vertices and indices
	for i=0:(group.size()-1)
		face = group[i]

		i0 = face[0]
		i1 = VertexIndex(-1)
		i2 = face[1]

		npolys = length(face)

		# Polygon -> triangle fan conversion
		for k=2:(npolys-1)
			i1 = i2
			i2 = face[k]

			v0 = updateVertex(i0, vCache, in, shape)
			v1 = updateVertex(i1, vCache, in, shape)
			v2 = updateVertex(i2, vCache, in, shape)

			addIndex(shape,v0)
			addIndex(shape,v1)
			addIndex(shape,v2)

			#shape.material_ids.push_back(material_id);
		end
	end

	#shape.name = name;
	#if (clearCache)	vCache.clear();

	true;

end

function parseIndicies(tmp_vi::Array{Array{VertexIndex,1},1},tmp::ShapeCoordinates,shape::Mesh)

	clear(shape) # remove old

	ig = 0
	iCount = length(tmp_vi)
	vCount = length(tmp.positions)
	vnCount = length(tmp.normals)
	vtCount = length(tmp.texcoords)

	if !iCount
		error("No index coordinates! (tmp)")
		return false
	end

	if !vCount
		error("No vertex coordinates! (tmp)")
		return false
	end

	# For each vertex of each triangle
	for (AAVertexIndex::iterator it = tmp_vi.begin(); it < tmp_vi.end(); ++it, ++ig)
		for (AVertexIndex::iterator jt = (*it).begin(); jt < (*it).end(); ++jt)
			vi = jt
			iv = vi.v
			ivn = vi.vn
			ivt = vi.vt

			if iv == -1
				error("Invalid index for vertex coordinates! %u (%d = -1).", ig, iv)
				return false
			else
				if iv < vCount
					addPosition(shape,*(tmp.positions.begin() + iv))
				else
					error("Missmatch count of vertex coordinates! %u (%d >= %u).", ig, iv, vCount)
					return false
				end
			end

			if ivn > -1
				if ivn < vnCount
					addNormal(shape,*(tmp.normals.begin() + ivn))
				else
					error("Missmatch count of normal coordinates! %u (%d = -1 or >= %u).", ig, ivn, vnCount)
					return false
				end
			end

			if ivt > -1
				if ivt < vtCount
					addTexcoord(shape,*(tmp.texcoords.begin() + ivt))
				else
					error("Missmatch count of texture coordinates! %u (%d = -1 or >= %u).", ig, ivt, vtCount)
					return false
				end
			end
		end
	end

	vCount = length(shape.positions)
	vnCount = length(shape.normals)
	vtCount = length(shape.texcoords)

	if !vCount
		error("No vertex coordinates!")
		return false
	end

	if vnCount > 0 && vCount != vnCount
		error("Missmatch count of vertex and normal coordinates! ($vCount != $vnCount).")
		return false
	end

	if vtCount > 0 && vCount != vtCount
		error("Missmatch count of vertex and texture coordinates! ($vCount != $vtCount).")
		return false
	end

	true
end

function loadOBJ(path::String, shape::Mesh)
	println("Load OBJ file $path...");

	tmp_v = AVec3f()
	tmp_vn = AVec3f()
	tmp_vt = AVec2f()
	tmp_g = AAVertexIndex()

	std::filebuf fb
	std::filebuf* fb_result = fb.open(path, std::ios::in)

	if fb_result == C_NULL
		error("Cannot open file $path!")
		return false
	end

	std::istream is(&fb)
	maxchars = 8192
	buf = zeros(Arrray{Char,1},maxchars)
	vertexCache = Dict{VertexIndex,UInt32}()

	min_BBOX = zeros(Vec3f)
	max_BBOX = zeros(Vec3f)

	while is.peek() != -1
		is.getline(&buf[0], maxchars)
		std::string linebuf(&buf[0])
		size = linebuf.size()
		id = size - 1

		# trim '\r\n', '\n'
		if size > 0 && linebuf[id] == '\n' linebuf.erase(id)
			size = linebuf.size()
			id = size - 1
		end

		if size > 0 && linebuf[id] == '\r' linebuf.erase(id)
			size = linebuf.size()
			id = size - 1
		end

		if empty(linebuf) continue end

		token = linebuf.c_str()
		token += strspn(token, " \t")
		if !token error("!token") end

		first = linebuf[0]
		second = linebuf[1]
		third = linebuf[2]

		if first == '\0' continue # empty line
		elseif first == '#' continue  # comment line

		elseif first == 'v' && isSpace(second)
			token += 2
			v = parsevec3(token)

			# set bounding box
			if v.x < min_BBOX.x min_BBOX.x = v.x end
			if v.y < min_BBOX.y min_BBOX.y = v.y end
			if v.z < min_BBOX.z min_BBOX.z = v.z end
			if v.x > max_BBOX.x max_BBOX.x = v.x end
			if v.y > max_BBOX.y max_BBOX.y = v.y end
			if v.z > max_BBOX.z max_BBOX.z = v.z end

			push!(tmp_v,v)
			continue

		elseif first == 'v' && second == 'n' && isSpace(third)
			token += 3
			push!(tmp_vn,parsevec3(token))
			continue

		elseif first == 'v' && second == 't' && isSpace(third)
			token += 3
			push!(tmp_vt,parsevec2(token))
			continue

		elseif first == 'f' && isSpace(second)
			token += 2
			token += strspn(token, " \t")

			group = Array{VertexIndex,1}()
			while !isNewLine(token[0])
				vi = parseTriple(token, tmp_v.size() / 3, tmp_vn.size() / 3, tmp_vt.size() / 2)
				push!(group,vi)
				n = strspn(token, " \t\r")
				token += n
			end

			push!(tmp_g,group)

			continue
		end
	end
	close(fb)

	scale_value = getScale(min_BBOX, max_BBOX)
	coord = ShapeCoordinates( tmp_v, tmp_vn, tmp_vt, getOrigin(min_BBOX, max_BBOX, scale_value), Vec3f(scale_value) )

	result = parseIndicies(vertexCache, tmp_g, coord, shape) # parseIndicies(tmp_g, in, shape);
	clear(tmp_g)

	if result println("Finsihed loading OBJ file $path.") end

	result
end

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
using Documenter, OceanDistributions
import PlutoSliderServer

makedocs(;
    modules=[OceanDistributions],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/gaelforget/OceanDistributions.jl/blob/{commit}{path}#L{line}",
    sitename="OceanDistributions.jl",
    authors="gaelforget <gforget@mit.edu>",
    assets=String[],
)

lst=("one_dim_diffusion.jl",)
for i in lst
    fil_in=joinpath(@__DIR__,"..", "examples",i)
    fil_out=joinpath(@__DIR__,"build", i[1:end-2]*"html")
    PlutoSliderServer.export_notebook(fil_in)
    mv(fil_in[1:end-2]*"html",fil_out)
    cp(fil_in,fil_out[1:end-4]*"jl")
end

deploydocs(;
    repo="github.com/gaelforget/OceanDistributions.jl",
)


text:
function Relationship(e::EzXML.Node) :: Relationship
    @assert EzXML.nodename(e) == "Relationship" "Unexpected XMLElement: $(EzXML.nodename(e)). Expected: \"Relationship\"."

    return Relationship(
        e["Id"],
        e["Type"],
        e["Target"]
    )
end

function get_relationship_target_by_id(wb::Workbook, Id::String) :: String
    for r in wb.relationships
        if Id == r.Id
            return r.Target
        end
    end
    error("Relationship Id=$(Id) not found")
end

function get_relationship_target_by_type(wb::Workbook, _type_::String) :: String
    for r in wb.relationships
        if _type_ == r.Type
            return r.Target
        end
    end
    error("Relationship Type=$(_type_) not found")
end

function has_relationship_by_type(wb::Workbook, _type_::String) :: Bool
    for r in wb.relationships
        if _type_ == r.Type
            return true
        end
    end
    false
end


text:
function reset_shadow!(cur_player::AbstractPlayer)
  cur_piece = cur_player.piece
  cur_shadow = cur_player.shadow

  cur_shadow.color = cur_piece.color
  cur_shadow.width = cur_piece.width

  cur_shadow.i = cur_piece.i
  cur_shadow.j = cur_piece.j

  for (piece_block, shadow_block) in zip(cur_piece.blocks, cur_shadow.blocks)
    shadow_block.x = piece_block.x
    shadow_block.y = piece_block.y
  end

  hard_drop(cur_player, cur_shadow)
end


text:
# """
# Visualisation methods, Bloch sphere plotting + other things
# """

# using LinearAlgebra
# using Plots

# theme(:juno)

# """
# Visualise the expectation value of op under the evoltuion of the ctrl_arr, at the moment we will ask for both the problem and the solution.
# """
# function visualise_expt_val(H_drift, H_ctrl, ctrl_arr, duration, n_timeslices, X_init::T, op; label = "") where T
#     # dealing with the evolution of the system is a bit of a problem
#     n_pulses = length(H_ctrl)
#     timestep = duration / n_timeslices
#     P_list = pw_evolve_save(H_drift, H_ctrl, ctrl_arr, n_pulses, timestep, n_timeslices)

#     states = similar(P_list)
#     states[1] = X_init
#     for i = 2:n_timeslices
#         states[i] = P_list[i - 1] * states[i - 1]
#     end

#     exp_op = [real(tr(op * i)) for i in states]
#     time = 0:timestep:duration - timestep

#     lab = label

#     plot(time, exp_op, label = lab, xlabel = "Time", ylabel = "Expt. value")

# end

# """
# Visualise the expectation value of an array of ops under the evoltuion of the ctrl_arr, at the moment we will ask for both the problem and the solution.
# """
# function visualise_expt_vals(H_drift::K, H_ctrl, ctrl_arr, duration, n_timeslices, X_init::T, ops::Array{K,1}, labels) where {T,K}
#     N = length(ops)
#     p = []
#     for i = 1:N
#         append!(p, [visualise_pulse(H_drift, H_ctrl, ctrl_arr, duration, n_timeslices, X_init, ops[i], label = labels[i])])
#     end

#     if sqrt(N) - floor(sqrt(N)) == 0
#         layout = (Int(sqrt(N)), Int(sqrt(N)))
#     else
#         layout = (N, 1)
#     end

#     plot(p..., layout = layout)
# end

# """
# Visualise a single pulse
# """
# function visualise_pulse(ctrl_arr::Array{T,1}, duration) where T
#     # we plot a single array
#     n_timeslices = length(ctrl_arr)
#     time = range(0, duration, length = n_timeslices)
#     bar(time, ctrl_arr, xlabel = "Time", ylabel = "Pulse Amplitude", label = "1")
# end

# """
# Visualise pulses that are stacked as a 2D array
# """
# function visualise_pulse(ctrl_arr::Array{T,2}, duration) where T
#     n_pulses = size(ctrl_arr)[1]
#     n_timeslices = size(ctrl_arr)[2]
#     time = range(0, duration, length = n_timeslices)
#     p = bar(time, ctrl_arr[1, :], xlabel = "Time", ylabel = "Pulse Amplitude", label = "1")
#     for i = 2:n_pulses
#         bar!(time, ctrl_arr[i, :], label = i)
#     end
#     current()
#     # show(p)
# end


text:
@testset "Angles" begin
  # 2D points
  @test ∠(P2(0,1),P2(0,0),P2(1,0)) ≈ T(-π/2)
  @test ∠(P2(1,0),P2(0,0),P2(0,1)) ≈ T(π/2)
  @test ∠(P2(-1,0),P2(0,0),P2(0,1)) ≈ T(-π/2)
  @test ∠(P2(0,1),P2(0,0),P2(-1,0)) ≈ T(π/2)
  @test ∠(P2(0,-1),P2(0,0),P2(1,0)) ≈ T(π/2)
  @test ∠(P2(1,0),P2(0,0),P2(0,-1)) ≈ T(-π/2)
  @test ∠(P2(0,-1),P2(0,0),P2(-1,0)) ≈ T(-π/2)
  @test ∠(P2(-1,0),P2(0,0),P2(0,-1)) ≈ T(π/2)

  # 3D points
  @test ∠(P3(1,0,0),P3(0,0,0),P3(0,1,0)) ≈ T(π/2)
  @test ∠(P3(1,0,0),P3(0,0,0),P3(0,0,1)) ≈ T(π/2)
  @test ∠(P3(0,1,0),P3(0,0,0),P3(1,0,0)) ≈ T(π/2)
  @test ∠(P3(0,1,0),P3(0,0,0),P3(0,0,1)) ≈ T(π/2)
  @test ∠(P3(0,0,1),P3(0,0,0),P3(1,0,0)) ≈ T(π/2)
  @test ∠(P3(0,0,1),P3(0,0,0),P3(0,1,0)) ≈ T(π/2)

  # 2D vectors
  @test ∠(V2(1, 0), V2(0, 1)) ≈ T(π / 2)
  @test ∠(V2(1, 0), V2(0, -1)) ≈ T(-π / 2)
  @test ∠(V2(1, 0), V2(-1, 0)) ≈ T(π)
  @test ∠(V2(0, 1), V2(-1, 0)) ≈ T(π / 2)
  @test ∠(V2(0, 1), V2(0, -1)) ≈ T(-π)
  @test ∠(V2(0, 1), V2(1, 1)) ≈ T(-π / 4)
  @test ∠(V2(0, -1), V2(1, 1)) ≈ T(π * 3 / 4)
  @test ∠(V2(-1, -1), V2(1, 1)) ≈ T(π)

  # 3D vectors
  @test ∠(V3(0, 0, 1), V3(1, 1, 0)) ≈ T(π / 2)
  @test ∠(V3(1, 0, 1), V3(1, 1, 0)) ≈ T(π / 3)
  @test ∠(V3(-1, -1, 0), V3(1, 1, 0)) ≈ T(π)
  @test ∠(V3(0, -1, -1), V3(0, 1, 1)) ≈ T(π)
  @test ∠(V3(0, -1, -1), V3(0, 1, 0)) ≈ T(π * 3 / 4)
  @test ∠(V3(0, 1, 1), V3(1, 1, 0)) ≈ T(π / 3)
end


text: read_docx(fn; flags="--self-contained") = read(`pandoc -f docx -t html $flags $fn`, String)


text:
@testset "GaussianMixtureModel" begin
    mu = randn(3)
    sigma = 1+rand(3)
    weights = [0.2,0.35,0.45]
    dm = GaussianMixtureModel(mu,sigma, weights)

    g = DataModel.getgrad(dm)
    l = DataModel.getllik(dm)

    for i in 1:10
        @test checkgrad(randn(1),l,g, eps = 1e-8) < 1e-5
    end
end


text:
using SparseArrays
using LinearAlgebra
using Logging


text:
module Jui

using MLStyle

using CImGui
using CImGui.CSyntax
using CImGui.CSyntax.CStatic
using CImGui.GLFWBackend
using CImGui.OpenGLBackend
using CImGui.GLFWBackend.GLFW
using CImGui.OpenGLBackend.ModernGL
using Printf

@static if Sys.isapple()
    # OpenGL 3.2 + GLSL 150
    const glsl_version = 150
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
    GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
    GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # required on Mac
else
    # OpenGL 3.0 + GLSL 130
    const glsl_version = 130
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 0)
end

# @data JWidget begin
#     JButton(String)
#     JText(String)
# end

macro jwindow(name, expr)
    # gen_window(name, expr) |> esc
    @when :(begin $(es...) end) = expr begin
        map(es) do e
            @match e begin
                ln::LineNumberNode => ln
                :(button($bname)) => :(CImGui.Button($bname))
                :(text($s)) => :(CImGui.Text($s))
                :(image($w, $h)) => :(CImGui.Image(Ptr{Cvoid}(image_id), ($w, $h)))
                node => node
            end
        end |> exprs ->
            Expr(:block,
                :(CImGui.Begin($name)),
                exprs...,
                :(CImGui.End())
            )
    end
end

# function gen_window(name :: String, expr)
#     @when :(begin $(es...) end) = expr begin
#         map(es) do e
#             @match e begin
#                 ln::LineNumberNode => ln
#                 JButton(bname::String) => :(CImGui.Button($bname))
#                 JText(s::String) => :(CImGui.Text($s))
#                 node => node
#             end
#         end |> exprs ->
#             Expr(:block,
#                 :(CImGui.Begin($name)),
#                 exprs...,
#                 :(CImGui.End())
#             )
#     end
# end

# function addWindow(name::String)
#     CImGui.Begin(name)  # create a window called "Hello, world!" and append into it.
#     CImGui.Text("This is some useful text.")  # display some text
#
#     CImGui.Button("Button")
#
#     CImGui.SameLine()
#     CImGui.Text("counter")
#
#     CImGui.End()
# end

# setup GLFW error callback
error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
GLFW.SetErrorCallback(error_callback)

# create window
window = GLFW.CreateWindow(1280, 720, "Demo")
@assert window != C_NULL
GLFW.MakeContextCurrent(window)
GLFW.SwapInterval(1)  # enable vsync

# setup Dear ImGui context
ctx = CImGui.CreateContext()

# setup Dear ImGui style
CImGui.StyleColorsDark()
# CImGui.StyleColorsClassic()
# CImGui.StyleColorsLight()

# setup Platform/Renderer bindings
ImGui_ImplGlfw_InitForOpenGL(window, true)
ImGui_ImplOpenGL3_Init(glsl_version)

img_width, img_height = 256, 256
image_id = ImGui_ImplOpenGL3_CreateImageTexture(img_width, img_height)
clear_color = Cfloat[0.45, 0.55, 0.60, 1.00]
while !GLFW.WindowShouldClose(window)
    # oh my global scope
    # TODO: remove all global variables
    global clear_color
    global image_id, img_width, img_height

    GLFW.PollEvents()
    # start the Dear ImGui frame
    ImGui_ImplOpenGL3_NewFrame()
    ImGui_ImplGlfw_NewFrame()
    CImGui.NewFrame()

    # show a simple window that we create ourselves.
    # we use a Begin/End pair to created a named window.

    image = rand(GLubyte, 4, img_width, img_height)
    ImGui_ImplOpenGL3_UpdateImageTexture(image_id, image, img_width, img_height)
    @jwindow "my window" begin
        button("Button111")
        text("some text")
        image(256, 256)
    end

    @cstatic f=Cfloat(0.0) counter=Cint(0) begin
        CImGui.Begin("Hello, world!")  # create a window called "Hello, world!" and append into it.
        CImGui.Text("This is some useful text.")  # display some text

        @c CImGui.SliderFloat("float", &f, 0, 1)  # edit 1 float using a slider from 0 to 1
        CImGui.ColorEdit3("clear color", clear_color)  # edit 3 floats representing a color
        CImGui.Button("Button") && (counter += 1)

        # CImGui.Button("Button-addWindow") && addWindow("new windows")

        CImGui.SameLine()
        CImGui.Text("counter = $counter")
        CImGui.Text(@sprintf("Application average %.3f ms/frame (%.1f FPS)", 1000 / CImGui.GetIO().Framerate, CImGui.GetIO().Framerate))

        CImGui.End()
    end

    # rendering
    CImGui.Render()
    GLFW.MakeContextCurrent(window)
    display_w, display_h = GLFW.GetFramebufferSize(window)
    glViewport(0, 0, display_w, display_h)
    glClearColor(clear_color...)
    glClear(GL_COLOR_BUFFER_BIT)
    ImGui_ImplOpenGL3_RenderDrawData(CImGui.GetDrawData())

    GLFW.MakeContextCurrent(window)
    GLFW.SwapBuffers(window)
end

# cleanup
ImGui_ImplOpenGL3_Shutdown()
ImGui_ImplGlfw_Shutdown()
CImGui.DestroyContext(ctx)

GLFW.DestroyWindow(window)

end # module Jui


text:
function miyak02()
#
#  test polynomial suggested by Goedecker
#  square of Fibocacci polynomial
#
    p = poly([(1.1+1.1im)*ones(4); (3.2+2.3im)*[1;1]; 2.1+1.5im])
    p = reverse(conv(p.a, p.a))
    z = [[(1.1+1.1im);(3.2+2.3im);2.1+1.5im] 2*[4;2;1]]
    p, PolyZeros(z)
end


text:
@doc raw"""
    costIntrICTV12(M, f, u, v, α, β)

Compute the intrinsic infimal convolution model, where the addition is replaced
by a mid point approach and the two functions involved are [`costTV2`](@ref)
and [`costTV`](@ref). The model reads

```math
E(u,v) =
  \frac{1}{2}\sum_{i ∈ \mathcal G}
    d_{\mathcal M}\bigl(g(\frac{1}{2},v_i,w_i),f_i\bigr)
  +\alpha\bigl( \beta\mathrm{TV}(v) + (1-\beta)\mathrm{TV}_2(w) \bigr).
```
"""
function costIntrICTV12(M::Manifold, f, u, v, α, β)
    IC = 1/2*distance(M, shortest_geodesic(M, u, v, 0.5), f)^2
    TV12 = β * costTV(M, u) + (1-β) * costTV2(M, v)
    return  IC + α*TV12
end

@doc raw"""
    costL2TV(M, f, α, x)

compute the $\ell^2$-TV functional on the `PowerManifold manifold `M` for given
(fixed) data `f` (on `M`), a nonnegative weight `α`, and evaluated at `x` (on `M`),
i.e.

```math
E(x) = d_{\mathcal M}^2(f,x) + \alpha \operatorname{TV}(x)
```

# See also
[`costTV`](@ref)
"""
costL2TV(M, f, α, x) = 1/2 * distance(M, f, x)^2  +  α*costTV(M, x)

@doc raw"""
    costL2TVTV2(M, f, α, β, x)

compute the $\ell^2$-TV-TV2 functional on the `PowerManifold` manifold `M` for
given (fixed) data `f` (on `M`), nonnegative weight `α`, `β`, and evaluated
at `x` (on `M`), i.e.

```math
E(x) = d_{\mathcal M}^2(f,x) + \alpha\operatorname{TV}(x)
  + \beta\operatorname{TV}_2(x)
```

# See also
[`costTV`](@ref), [`costTV2`](@ref)
"""
costL2TVTV2(M::PowerManifold, f, α, β, x) = 1/2*distance(M,f,x)^2 + α*costTV(M,x) + β*costTV2(M,x)

@doc raw"""
    costL2TV2(M, f, β, x)

compute the $\ell^2$-TV2 functional on the `PowerManifold` manifold `M`
for given data `f`, nonnegative parameter `β`, and evaluated at `x`, i.e.

```math
E(x) = d_{\mathcal M}^2(f,x) + \beta\operatorname{TV}_2(x)
```

# See also
[`costTV2`](@ref)
"""
function costL2TV2(M::PowerManifold, f, β, x)
    return 1/2*distance(M,f,x)^2 + β*costTV2(M,x)
end

@doc raw"""
    costTV(M, x, p)

Compute the $\operatorname{TV}^p$ functional for a tuple `pT` of pointss
on a [Manifold](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.Manifold) `M`, i.e.

```math
E(x_1,x_2) = d_{\mathcal M}^p(x_1,x_2), \quad x_1,x_2 ∈ \mathcal M
```

# See also

[`∇TV`](@ref), [`prox_TV`](@ref)
"""
function costTV(M::Manifold, x::Tuple{T,T}, p::Int=1) where {T}
  return distance(M,x[1],x[2])^p
end
@doc raw"""
    costTV(M,x [,p=2,q=1])

Compute the $\operatorname{TV}^p$ functional for data `x`on the `PowerManifold`
manifold `M`, i.e. $\mathcal M = \mathcal N^n$, where $n ∈ \mathbb N^k$ denotes
the dimensions of the data `x`.
Let $\mathcal I_i$ denote the forward neighbors, i.e. with $\mathcal G$ as all
indices from $\mathbf{1} ∈ \mathbb N^k$ to $n$ we have
$\mathcal I_i = \{i+e_j, j=1,\ldots,k\}\cap \mathcal G$.
The formula reads

```math
E^q(x) = \sum_{i ∈ \mathcal G}
  \bigl( \sum_{j ∈  \mathcal I_i} d^p_{\mathcal M}(x_i,x_j) \bigr)^{q/p}.
```

# See also
[`∇TV`](@ref), [`prox_TV`](@ref)
"""
function costTV(M::PowerManifold, x, p=1, q=1)
    power_size = power_dimensions(M)
    R = CartesianIndices(Tuple(power_size))
    d = length(power_size)
    maxInd = last(R)
    cost = fill(0.,Tuple(power_size))
    for k in 1:d # for all directions
        ek = CartesianIndex(ntuple(i  ->  (i==k) ? 1 : 0, d) ) #k th unit vector
        for i in R # iterate over all pixel
            j = i+ek # compute neighbor
            if all( map(<=, j.I, maxInd.I)) # is this neighbor in range?
                cost[i] += costTV(
                    M.manifold,
                    (get_component(M,x,i),get_component(M,x,j)),
                    p
                )
            end
        end
    end
    cost = (cost).^(1/p)
    if q > 0
        return sum(cost.^q)^(1/q)
    else
        return cost
    end
end
@doc raw"""
    costTV2(M,(x1,x2,x3) [,p=1])

Compute the $\operatorname{TV}_2^p$ functional for the 3-tuple of points
`(x1,x2,x3)`on the [Manifold](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.Manifold) `M`. Denote by

```math
  \mathcal C = \bigl\{ c ∈  \mathcal M \ |\ g(\tfrac{1}{2};x_1,x_3) \text{ for some geodesic }g\bigr\}
```

the set of mid points between $x_1$ and $x_3$. Then the functionr reads

$d_2^p(x_1,x_2,x_3) = \min_{c ∈ \mathcal C} d_{\mathcal M}(c,x_2).$

# See also
[`∇TV2`](@ref), [`prox_TV2`](@ref)
"""
function costTV2(M::MT, x::Tuple{T,T,T}, p=1) where {MT <: Manifold, T}
  # note that here mid_point returns the closest to x2 from the e midpoints between x1 x3
  return 1/p*distance(M,mid_point(M,x[1],x[3]),x[2])^p
end
@doc raw"""
    costTV2(M,x [,p=1])

compute the $\operatorname{TV}_2^p$ functional for data `x` on the
`PowerManifold` manifoldmanifold `M`, i.e. $\mathcal M = \mathcal N^n$,
where $n ∈ \mathbb N^k$ denotes the dimensions of the data `x`.
Let $\mathcal I_i^{\pm}$ denote the forward and backward neighbors, respectively,
i.e. with $\mathcal G$ as all indices from $\mathbf{1} ∈ \mathbb N^k$ to $n$ we
have $\mathcal I^\pm_i = \{i\pm e_j, j=1,\ldots,k\}\cap \mathcal I$.
The formula then reads

```math
E(x) = \sum_{i ∈ \mathcal I,\ j_1 ∈  \mathcal I^+_i,\ j_2 ∈  \mathcal I^-_i}
d^p_{\mathcal M}(c_i(x_{j_1},x_{j_2}), x_i),
```

where $c_i(\cdot,\cdot)$ denotes the mid point between its two arguments that is
nearest to $x_i$.

# See also
[`∇TV2`](@ref), [`prox_TV2`](@ref)
"""
function costTV2(M::PowerManifold, x, p::Int=1, Sum::Bool=true)
    Tt = Tuple( power_dimensions(M) )
    R = CartesianIndices( Tt )
    d = length(Tt)
    minInd, maxInd = first(R), last(R)
    cost = fill(0., Tt)
    for k in 1:d # for all directions
        ek = CartesianIndex(ntuple(i  ->  (i==k) ? 1 : 0, d) ) #k th unit vector
        for i in R # iterate over all pixel
            jF = i+ek # compute forward neighbor
            jB = i-ek # compute backward neighbor
            if all( map(<=, jF.I, maxInd.I) ) && all( map(>=, jB.I, minInd.I)) # are neighbors in range?
                cost[i] += costTV2(
                    M.manifold,
                    (get_component(M,x,jB),get_component(M,x,i),get_component(M,x,jF)),
                    p,
                )

            end
        end # i in R
  end # directions
  if p != 1
    cost = (cost).^(1/p)
  end
  if Sum
    return sum(cost)
  else
    return cost
  end
end


text:
#=
This document implements a WrightSampler mutable struct type.
This object is used to generate Gibbs-sampler output on Wright equilibria
Also included are methods of the current state to generate stochastic gradients
And fisher information estimates
 =#

# This lists the neccesary parameters
# for sampling from the Wright-equilibrium distribution. The partition function is the most costly to compute and so it cahched
# θ_mut, θ_wt, σ are forward and backward diversity parameters respectively. At the moment, there is only code for  σ>0.
mutable struct WrightSampler
    θ_wt::Float64
    θ_mut::Float64
    σ::Float64
    k::Int64   # auxiliary variable for gibbs sampling
    x::Float64 # The freqeuncy of wt
end

function WrightSampler(θ_wt, θ_mut, σ)
#    z = gamma(θ_wt) * HypergeometricFunctions.pFq([θ_wt],[θ_mut + θ_wt],σ) /gamma(θ_wt +θ_mut)
    ws = WrightSampler(θ_wt, θ_mut, σ, 0, .5)
    for ii in 1:10
    incrementk!(ws)
    incrementx!(ws)
    end
    return ws
end

# This sampler is robuset for small-exponent gamma distributions
# Code was copied from somewhere unknown (citation needed)
function randlogGamma(a;scale = 1)
    if a > .2
        return log.(rand(Gamma(a)))
    else
        L = 1 / a - 1
        w = a / MathConstants.e / (1 - a)
        ww = 1 / (1 + w)
        η = z -> ( z >= 0 ? exp(-z) : w * L * exp(L * z))
        h = z -> exp(-z - exp(-z / a))
        function rh(a)
            while true
                U = rand()
                z = (U <= ww) ? -log(U / ww) : log(rand()) / L
                h(z) / η(z) > rand() && return(z)
            end
        end
        return log(scale) - rh(a) / a
    end
end

function sigmoid(x) # numerically stable sigmoid
    if x >= 0
        1 / ( 1 + exp(-x) )
    else
        exp(x) / ( 1 + exp(x) )
    end
end

function incrementk!(ws::WrightSampler)
    if ws.σ >= 0
        ws.k = rand(Poisson(abs(ws.σ * ws.x)))
    else
        ws.k = rand(Poisson(abs(ws.σ * (1-ws.x))))
    end
end

function incrementx!(ws::WrightSampler)
    if ws.σ >= 0
        a = randlogGamma(ws.θ_wt + ws.k;scale = 1)
        b = randlogGamma(ws.θ_mut;scale = 1)
    else
        a = randlogGamma(ws.θ_wt;scale = 1)
        b = randlogGamma(ws.θ_mut + ws.k;scale = 1)
    end
    ws.x = sigmoid(a-b)
end

function sample!(ws::WrightSampler)
    incrementk!(ws)
    incrementx!(ws)
    return ws.x
end

function rejection_sample_k(θ_wt, θ_mut,σ)
    k = 0
    while true
        k = rand(Poisson(σ))
        k == 0 ? break : nothing
        p = exp( loggamma(θ_wt + θ_mut) + loggamma(k + θ_wt) - loggamma(k + θ_wt + θ_mut) -loggamma(θ_wt))
        rand(Bernoulli(p)) ? break : nothing
    end
    return k
end

## Generate rebound values

struct MutantState
	n::Int # number of individuals
	t::Float64
end

function reboundtime_step_closure(;theta = 0.02, mut_per_gen = 10^5, extra_mut_flux = 1.0, delta_per_mut = 100, pop_size::Int = 10^4, decay_rate = .5, x0 = 10^-3)
	# x0 is the inital mutant fraction
	let f = 1/3,
	µ = mut_per_gen * f,
	θ = .02, K = pop_size,
	λ = μ*2*K/θ, r =decay_rate,
	delta = delta_per_mut*μ # paramters in absolute units (running on f = 1/3 days)
	if λ<f
		error("pop_size too small to support growth rate")
	end
	function step_state(state::MutantState)
		mut = extra_mut_flux*K*μ*exp(-r*state.t)
		rate = mut + λ*state.n
		delta_t = randexp()/rate
		if mut/rate > rand()
			return MutantState(state.n+1,state.t+delta_t)
		else
			birth = (λ + (f - delta) -
					((1-x0)*exp(-r*state.t)*f + state.n*(f - delta)/K)
					)/2
			if birth/λ > rand()
				return MutantState(state.n+1,state.t+delta_t)
			else
				return MutantState(state.n-1,state.t+delta_t)
			end
		end
	end

	end

end

function run_reboundtimes(start::Int ; theta = 0.02, nsample = 1000, mut_per_gen = 10^(-5.), extra_mut_flux = 1.0, delta_per_mut = 100, pop_size=10^4, decay_rate = .5)
	steppy = reboundtime_step_closure(;theta, mut_per_gen, extra_mut_flux, delta_per_mut, pop_size, decay_rate, x0 = start/pop_size)
    results = Float64[]
    q = 0.5
	for _ in 1:nsample
	state = MutantState(start,0);
	while (&)(state.t < 56, state.n <  pop_size*q)
		state = steppy(state)
	end
	if state.t < 56
		t_reb = 3*log(1 + exp(state.t / 3)*(1-q)/q) # Interpolate to deterministic rebound time
	else
		t_reb = 57
	end
	push!(results,t_reb)
	end
	return results
end

# function ab_escapes(theta::Array{Float64,1}, abprofiles::Array{Array{Float64,1},1}...; samples = 1)
#     out = Float64[]
#     for θ in theta
#         samplingarray = [map(x->WrightSampler( (θ.*x)... ), ab) for ab in abprofiles] # array of array of wright samplers
#         for _ in 1:samples
#         push!(out, prod( sample!(ab)) for ab in abprofiles)) # escape probability is the probability of being escaped to all of the abs
#         #each ab consists of an array of arrays of site profiles
#         end
#     end
# end

max_stars_repo_name: StatPhysBio/Tomoko

max_issues_repo_name: StatPhysBio/Tomoko

max_forks_repo_name: StatPhysBio/Tomoko


text:
struct EvaluationEnvironment{T}
    molecule_indices::Vector{Int}
    system_size::Tuple{Int,Int}
    other_atoms::Matrix{T}
    offset::Vector{T}
    orthogonal_vectors::Matrix{T}
end

Base.broadcastable(p::EvaluationEnvironment) = Ref(p)

function EvaluationEnvironment(molecule_indices, system_size, other_atoms, height, surface_normal)
    if size(other_atoms, 2) != 0
        height += maximum(other_atoms[3,:])
    end
    offset = surface_normal .* height
    orthogonal_vectors = nullspace(reshape(surface_normal, 1, :)) # Plane parallel to surface
    return EvaluationEnvironment(molecule_indices, system_size, other_atoms, offset, orthogonal_vectors)
end

"""
    separate_slab_and_molecule(atom_indices, r)

Get the coordinates of the molecule and slab separately.
"""
function separate_slab_and_molecule(atom_indices, r)
    molecule = r[:,atom_indices]
    slab_indices = [i for i in axes(r, 2) if i ∉ atom_indices]
    slab = r[:,slab_indices]
    return (molecule, slab)
end

"""
    combine_slab_and_molecule(atom_indices, molecule, slab)

Revert the transformation `separate_slab_and_molecule`
"""
function combine_slab_and_molecule(atom_indices, molecule, slab)
    r = zeros(size(molecule, 1), size(molecule, 2) + size(slab, 2))
    slab_indices = [i for i in axes(r, 2) if i ∉ atom_indices]
    r[:,atom_indices] .= molecule
    r[:,slab_indices] .= slab
    return r
end

"""
    calculate_diatomic_energy(model::AdiabaticModel, bond_length::Real;
        height=10, normal_vector=[0, 0, 1])

Returns potential energy of diatomic with `bond_length` at `height` from surface.

Orients molecule parallel to the surface at the specified height, the surface is assumed
to intersect the origin.
This requires that the model implicitly provides the surface, or works fine without one.
"""
function calculate_diatomic_energy(
    bond_length::Real, model::AdiabaticModel, environment::EvaluationEnvironment
)
    r = assemble_evaluation_geometry(bond_length, environment)
    return NQCModels.potential(model, r)
end

function assemble_evaluation_geometry(bond_length::Real, environment::EvaluationEnvironment)
    molecule = build_molecule(bond_length, environment)
    r = combine_slab_and_molecule(environment.molecule_indices, molecule, environment.other_atoms)
    return r
end

function build_molecule(bond_length::Real, environment::EvaluationEnvironment)
    if environment.system_size == (1, 1)
        r = [bond_length;;]
    else
        r = environment.offset .+ environment.orthogonal_vectors .* bond_length ./ sqrt(2)
    end
    return r
end


text:
module LatentDiffEq

using OrdinaryDiffEq
using DiffEqFlux
using DiffEqSensitivity
using Flux
using Statistics
using CUDA

## Types definitions
abstract type LatentDE end

## Models definitions
include("./models/LatentDiffEqModel.jl")
include("./models/GOKU.jl")
include("./models/LatentODE.jl")
export LatentDiffEqModel, GOKU, LatentODE
export default_layers

include("./utils/utils.jl")
export vector_mse, kl, vector_kl, frange_cycle_linear
export normalize_to_unit_segment, time_loader, rand_time

end # end LatentDiffEq module

max_issues_repo_name: gabrevaya/LatentDiffEq

max_forks_repo_name: gabrevaya/LatentDiffEq


text:
# v1.0 only
# Import dependency.
using BenchmarkTools, Roots, Statistics

#Create benchmark group and benchmarks
benchmarks = BenchmarkGroup()

#Put in specific benchmarks

bracket_methods = (bisection=Roots.Bisection(), a42=Roots.A42(),
                   aps=Roots.AlefeldPotraShi())
derivative_free_methods =  (o0=Roots.Order0(), o1=Roots.Order1(), o1b=Roots.Order1B(),
                            o2=Roots.Order2(), o2b=Roots.Order2B(),
                            o5=Roots.Order5(), o8=Roots.Order8(), o16=Roots.Order16()
)

# collection of doable problems
problems = Dict("f1" => (x -> sin(x), 3.0,  (3.0, 4.0)),
                "f2" => (x -> x^5 - x - 1, 1.0, (0.5, 5.0)),
                "f3" => (x -> exp(x) - x^4, 7.0, (5.0, 20.0)),
                "f4" => (x -> cos(x) - x/2, pi/4, (0.0, pi/2)),
                "f5" => (x -> x^2 - exp(x) - 3x + 2, -0.5, (-1.0, 1.0)),
                "f6" => (x -> x^2 - exp(x) - 3x + 2, 2.0, (0.0, 3.0)),
                "f7" => (x -> tanh(x) - tan(x), 7.6, (4.0, 8.0)),
"f8" => (x -> exp(-x^2 + x + 2) - cos(x) + x^3 + 1, -0.5, (-2.0, 1.0)),
"f9" => (x -> log(x) + sqrt(x) - 5, 7, (7.0, 10.0)),
"f10" => (x -> log(x) + sqrt(x) - 5, 20, (7.0, 10.0))
                )

function run_bracket(problems, Ms)

    for (nm, prob) in problems
        fn, x0, ab = prob
        for (mnm, M) in zip(fieldnames(typeof(Ms)), Ms)
            solve(ZeroProblem(fn,ab),M)
            #find_zero(fn, ab, M)
        end
    end

end

function run_bracketing(problems, Ms)
    rts = Float64[]
    for (nm, prob) in problems
        fn, x0, ab = prob
        for M in Ms
            rt  = solve(ZeroProblem(fn,ab),M)
            #rt = find_zero(fn, ab, M)
            push!(rts, rt)
        end
    end
    rts

end

function run_derivative_free(problems, Ms)
    rts = Float64[]
    for (nm, prob) in problems
        fn, x0, ab = prob
        for M in Ms
            if M == Order0()
                rt = find_zero(fn, x0, M)
            else
                rt  = solve(ZeroProblem(fn,ab),M)
            end

            push!(rts, rt)
        end
    end
    rts

end

function run_simple(problems)
    rts = Float64[]
    for (nm, prob) in problems
        fn, x0, ab = prob
        push!(rts, Roots.bisection(fn, ab[1], ab[2]))
        push!(rts, Roots.bisection(fn, ab[1], ab[2], xatol=1e-6))
        push!(rts, Roots.secant_method(fn, x0))
    end
    rts
end

benchmarks = BenchmarkGroup()

benchmarks["bracketing"] = @benchmarkable run_bracketing($problems, $bracket_methods)
benchmarks["derivative_free"] = @benchmarkable run_derivative_free($problems, $derivative_free_methods)
benchmarks["simple"] = @benchmarkable run_simple($problems)

for (nm, prob) in problems
    fn, x0, ab = prob
    @assert fn(ab[1]) * fn(ab[2]) < 0

    Ms =  bracket_methods
    for (mnm, M) in zip(fieldnames(typeof(Ms)), Ms)
        benchmarks[nm * "-" * string(mnm)] = @benchmarkable find_zero($fn, $ab, $M)
    end

    Ms =  derivative_free_methods
    for (mnm, M) in zip(fieldnames(typeof(Ms)), Ms)
        benchmarks[nm * "-" * string(mnm)] = @benchmarkable find_zero($fn, $x0, $M)
    end

    # simple methods
    u,v = ab
    benchmarks[nm * "-bisection"] = @benchmarkable Roots.bisection($fn, $u, $v)
    benchmarks[nm * "-bisection-atol"] = @benchmarkable Roots.bisection($fn, $u, $v, xatol=1e-6)
        benchmarks[nm * "-secant"] = @benchmarkable Roots.secant_method($fn, $x0)
end

results = run(benchmarks) # Get results.
results = median(results) # Condense to median.

nm = "benchmarks.json"
fname = joinpath(@__DIR__,  nm)
if isinteractive()
    println("""
To save results, manually call in the REPL: BenchmarkTools.save("benchmarks.json", results)
""")
end

#Compare to old results
try
    oldresults= BenchmarkTools.load(fname)[1]
    judge(oldresults, results)
catch err
    error("Couldn't load file- make sure that you've previously saved results.", err.prefix)
end


text:
using Random, Plots, DataFrames, ColorSchemes

"""
    PlotBasic(df::DataFrame,nn::Integer,dMax::Float64=0.)

Plot random subset of size nn trajectories.
"""
function PlotBasic(df::DataFrame,nn::Integer,dMax::Float64=0.)
   IDs = randperm(maximum(df.ID))
   COs=["w" "y" "g" "k"]

   plt=plot(leg=false)
   df_by_ID = groupby(df, :ID)
   ID_list=[df_by_ID[i][1,:ID] for i in 1:length(df_by_ID)]
   for ii=1:nn
      jj=findall(ID_list.==IDs[ii])[1]
      tmp=df_by_ID[jj]
      if dMax > 0.
         d=abs.(diff(tmp[!,:lon]))
         jj=findall(d .> dMax)
         tmp[jj,:lon].=NaN; tmp[jj,:lat].=NaN
         d=abs.(diff(tmp[!,:lat]))
         jj=findall(d .> dMax)
         tmp[jj,:lon].=NaN; tmp[jj,:lat].=NaN
      end
      CO=COs[mod(ii,4)+1]

      hasproperty(df,:z) ? plot!(tmp[!,:lon],tmp[!,:lat],tmp[!,:z],linewidth=0.3) : nothing
      !hasproperty(df,:z) ? plot!(tmp[!,:lon],tmp[!,:lat],linewidth=0.3) : nothing
      #plot!(tmp[!,:lon],tmp[!,:lat],tmp[!,:z],linewidth=0.3)
   end
   return plt
end

"""
    scatter_subset(df,t)

```
nf=size(u0,2)
t=[ceil(i/nf) for i in 1:367*nf]
df[!,:t]=2000 .+ 10/365 * t

@gif for t in 2000:0.1:2016
   scatter_subset(df,t)
end
```
"""
function scatter_subset(df,t)
    dt=0.25
    df_t = df[ (df.t.>t-dt).&(df.t.<=t) , :]
    scatter(df_t.lon,df_t.lat,markersize=2,
    xlims=(-180.0,180.0),ylims=(-90.0,90.0))
end

"""
    scatter_zcolor(df,t,zc; plt=plot(),cam=(0, 90))

```
t=maximum(df[!,:t])
scatter_zcolor(df,t,df.z)
```
"""
function scatter_zcolor(df,t,zc; plt=plot(),cam=(0, 90), dt=1.0)
    lo=extrema(df.lon); lo=round.(lo).+(-5.0,5.0)
    la=extrema(df.lat); la=round.(la).+(-5.0,5.0)
    de=extrema(zc); de=round.(de).+(-5.0,5.0)

    df_t = df[ (df.t.>t-dt).&(df.t.<=t) , :]
    zc_t = Float64.(zc[ (df.t.>t-dt).&(df.t.<=t)])

    #fig=deepcopy(plt)
    scatter(df_t.lon,df_t.lat,zc_t,zcolor = zc_t,
    markersize=2,markerstrokewidth=0.1,camera = cam,
    xlims=lo,ylims=la,zlims=de,clims=de)
end

"""
    scatter_movie(𝐼; cam=(0, 90))

Animation using `scatter_zcolor()
```
𝐼,Γ=example3("OCCA", lon_rng=(-165.0,-155.0),lat_rng=(25.0,35.0), z_init=5.5,)
scatter_movie(𝐼,cam=(70, 70))
```
"""
function scatter_movie(𝐼; cam=(0, 90))
   df=𝐼.🔴
   nf=maximum(df.ID)
   nt=min(size(df,1)/nf,100)
   dt=maximum(df.t)/(nt-1)
   #println("nt="*"$nt"*"dt="*"$dt")
   return @gif for t in 0:nt-1
        scatter_zcolor(df,t*dt,df.z;cam=cam,dt=dt)
   end
end

"""
    phi_scatter(ϕ,df)

```
phi_scatter(ϕ,df)
```
"""
function phi_scatter(ϕ,df)
    nx,ny=size(ϕ)
    contourf(-0.5 .+ (1:nx),-0.5 .+ (1:ny),
             transpose(ϕ),c = :blues,linewidth = 0.1)
    scatter!(df.x,df.y,markersize=4.0,c=:red,marker = (:circle, stroke(0)),
             xlims=(0,nx),ylims=(0,ny),leg=:none)
end

"""
    DL()

Compute Ocean depth logarithm.
"""
function DL(Γ)
    lon=[i for i=19.5:1.0:379.5, j=-78.5:1.0:78.5]
    lat=[j for i=19.5:1.0:379.5, j=-78.5:1.0:78.5]
    DL=interp_to_lonlat(Γ["Depth"],Γ,lon,lat)
    DL[findall(DL.<0)].=0
    DL=transpose(log10.(DL))
    DL[findall((!isfinite).(DL))].=NaN
    return lon[:,1],lat[1,:],DL
end

"""
    plot_end_points(𝐼::Individuals,Γ)

Plot initial and final positions, superimposed on a map of ocean depth log.
"""
function plot_end_points(𝐼::Individuals,Γ)
    lo,la,dl=DL(Γ)
    xlims=extrema(lo)
    ylims=extrema(la)
    plt=contourf(lo,la,dl,clims=(1.5,5),c = :ice, colorbar=false, xlims=xlims,ylims=ylims)

    t=𝑃.𝑇[2]
    df = 𝐼.🔴[ (𝐼.🔴.t.>t-1.0).&(𝐼.🔴.t.<=t) , :]
    lo=deepcopy(df.lon); lo[findall(lo.<xlims[1])]=lo[findall(lo.<xlims[1])].+360
    scatter!(lo,df.lat,markersize=1.5,c=:red,leg=:none,marker = (:circle, stroke(0)))

    t=0.0
    df = 𝐼.🔴[ (𝐼.🔴.t.>t-1.0).&(𝐼.🔴.t.<=t) , :]
    lo=deepcopy(df.lon); lo[findall(lo.<xlims[1])]=lo[findall(lo.<xlims[1])].+360
    scatter!(lo,df.lat,markersize=1.5,c=:yellow,leg=:none,marker = (:dot, stroke(0)))

    return plt
end