text:
const BAREHTML = "<meta http-equiv=\"Content-Type\" content=\"text/html; charset=UTF-8\">
 <title>Empty.html</title>"
import Sockets
const BODY =  "
API:

/api/runIBox

/api/getDataTree/
                "
"""
структура для хранения константного состояния сервера
"""
struct ServerState
    path::String
    obj::Dict #храним все объекты на основе файла
end
ServerState() = ServerState("")
using HTTP
"""
Формирование стандартной шапки HTTP-ответа
"""
@inline function addResponseHeader(res)
    push!(res.headers, "Connection" => "keep-alive")
    push!(res.headers, "Access-Control-Allow-Origin" => "*")
    push!(res.headers, "Content-Type" => "text/html; charset=utf-8")
    push!(res.headers, "Access-Control-Allow-Methods" => "POST, GET")
    res
end

"""
корневой запрос - показать страницу с описанием API
"""
handle(reg::HTTP.Request) = replace(BAREHTML, "" => BODY) |> HTTP.Response

"""
`filename, filepath, datapath, param = parse_uri(reqstr, srvpath)`

Входные аргументы:
`reqstr` - строка запрсоа
`srvpath` - путь на сервере

Выходные аргументы:
`filename` - имя запрошенной записи
`filepath` - полный путь hdf5-файла на сервере
`datapath` - путь к данным внутри hdf-файла
`param` - параметры запроса (словарь)
"""
function parse_uri(reqstr::AbstractString, srv::ServerState)
    uri = parse(HTTP.URI, reqstr)
    param = HTTP.URIs.queryparams(uri)
    args = HTTP.URIs.splitpath(uri.path)

    if !haskey(param,"res")
        #старый вариант работы
        filename = "" # omit /api/command/...
        filepath = srv.obj["filepath"] #
    else
        filepath = joinpath(srv.path, param["res"])
        filename = basename(filepath)

        srv.obj["filename"]  = filename
        srv.obj["filepath"]  = filepath
    end

    if haskey(param,"dataName")
        datapath = string(param["dataName"])
        # собираем или ищем датагруппу
    else
        datapath = ""
    end
    # manualpath = joinpath(srvpath, filename, "manual_mark.hdf")
    return filename, filepath, datapath, param
end

#стартуем бокс на файле данных
function runIBox(srv::ServerState, req::HTTP.Request)
    filename, filepath, datapath, param = parse_uri(req.target, srv)

    # @info filename, filepath, datapath, param

    IBox_path = param["IBox_path"]
    port = param["IBox_port"]

    if haskey(param, "config") #какую конфигурацию бокса брать
        conf = param["config"]
    else
        conf = "ConfigClsWebApi"
    end
    if haskey(param, "resName") #в какую разметку писать или из какой читать
        resName = param["resName"]
    else
        resName = "000"
    end
    if haskey(param, "arg") #если надо открыть, то пишем -open
        arg = param["arg"]
    else
        arg = ""
    end

    out = "Тут тип загрузился бокс"
    # args = `-config:IBTestWebApi -WebAPISrc[port=$port|debug=true] -finalize -res:000`

    args = `-config:$conf -WebAPISrc[port=$port|debug=true] -finalize -res:$resName $arg`

    # args = `-config:IBOpen -WebAPISrc[port=$port|debug=true] -finalize -res:000 -open`

    command = `$IBox_path $filepath $args`
    @show command
    @async run(command)
    out = "Бокс запущен..."
    srv.obj["IBox_port"] = port
    srv.obj["IBox_path"] = IBox_path
    srv.obj["IBox_host"] = IPv4(param["IBox_host"])

    res = out |> HTTP.Response |> addResponseHeader
end

#закрываем бокс на файле данных
function Close(srv::ServerState, req::HTTP.Request)
    filename, filepath, datapath, param = parse_uri(req.target, srv)

    # @info filename, filepath, datapath, param

    localIP = srv.obj["IBox_host"]
    port = srv.obj["IBox_port"]

    r = HTTP.request("GET", "http://$localIP:$port/api/Close")
    out = "Бокс закрыт..."

    res = out |> HTTP.Response |> addResponseHeader
end

function redirectRequest(srv::ServerState, req::HTTP.Request)
    filename, filepath, datapath, param = parse_uri(req.target, srv)
    regstr = req.target #само содержание запроса
    @info filename, filepath, datapath, param
    if haskey(srv.obj["dataStorage"],datapath) #был объект, значит были ручные правки - читаем из него

        DG = srv.obj["dataStorage"][datapath]
        if haskey(param,"from")
            from = parse(Int32,param["from"])
        end
        if haskey(param,"to")
            to = parse(Int32,param["to"])
        end

        if occursin("getData",regstr)
            data = getData(DG,from,to,param["fields"])
        elseif occursin("getStructData",regstr)
            data = getStructData(DG,from,to,param["fields"])
        end
        if isa(data,Array{Any})
            t_data = (data...,)
        else
            t_data = data
        end
        println(data)
        if !isempty(data)
            out = pack_vec(t_data) |> base64encode
        else
            out = "null"
        end
        println("Объект имеется")
    else #объекта нет, значит читаем напрямую
        port = srv.obj["IBox_port"]
        localIP = srv.obj["IBox_host"]
        println("http://$localIP:$port$regstr")
        r = HTTP.request("GET", "http://$localIP:$port$regstr")
        # @info r
        out = r.body
    end

    res = out |> HTTP.Response |> addResponseHeader
end

#бокс уже должен быть запущен! читаем данные из бокса по getData
function getData(srv::ServerState, req::HTTP.Request)
    filename, filepath, datapath, param = parse_uri(req.target, srv)

    @info filename, filepath, datapath, param
    port = srv.obj["IBox_port"]
    localIP = srv.obj["IBox_host"]

    data = getData(localIP,port,param)
    if !isempty(data)
        out = pack_vec(data) |> base64encode
    else
        out = "null" # gvg: это стандартный вывод?
    end

    res = out |> HTTP.Response  |> addResponseHeader
end

"""
запись ручных правок
{
    "chName": "",
    "targetData": "",
    "command": {
        "id": "",
        "args": {
            "ibeg": [],
            "iend": [],
            "type": [""],
            "mode": ""
        }
    }
}
"""
function manualChange(srv::ServerState, req::HTTP.Request)
    #информация о добавляемом сегменте лежит в JSONe, он может передаваться вместе с данными о файле
    #пока файл подчитывается из папки
    # println("я зашел")

    filename, filepath, datapath, param = parse_uri(req.target, srv)

    infoEvent = String(req.body)
    # println(infoEvent)
    #парсим команду
    port = srv.obj["IBox_port"]
    localIP = srv.obj["IBox_host"]
    #парсим команду
    segInRange = parseCommand(localIP, port, srv.obj,infoEvent)

    if isa(segInRange,StructArray)
        out = pack_vec(segInRange) |> base64encode
    else
        out = "Command was parsed"
    end

    res = out |> HTTP.Response |> addResponseHeader

    return res
end

"""
получение дерева классов QRS-комплекса в json-формате
работает ТОЛЬКО С ИСХОДНЫМИ ДАННЫМИ БОКСА - правки не подтягивает
"""
function getQRStree(srv::ServerState, req::HTTP.Request)

    filename, filepath, datapath, param = parse_uri(req.target, srv)

    port = srv.obj["IBox_port"]
    localIP = srv.obj["IBox_host"]

    #запрашиваем из Бокса данные для дерева
    r = HTTP.request("GET", "http://$localIP:$port/api/getData?dataName=QPoint&all")
    QPoint = reinterpret(Int32, base64decode(r.body)) |> collect

    r = HTTP.request("GET", "http://$localIP:$port/api/getData?dataName=ClassQRS&all")
    Class = reinterpret(Int32, base64decode(r.body)) |> collect

    r = HTTP.request("GET", "http://$localIP:$port/api/getData?dataName=SubClassQRS&all")
    SubClass = reinterpret(Int16, base64decode(r.body)) |> collect

    #строим дерево
    result = buildQRStree(QPoint, Class, SubClass)

    res = result |> HTTP.Response |> addResponseHeader

end

function undo!(srv::ServerState, req::HTTP.Request)
    port = srv.obj["IBox_port"]
    localIP = srv.obj["IBox_host"]

    undo!(srv.obj,localIP, port)
    @info srv.obj["history"]
    @info srv.obj["state"]
    out = "State "*string(srv.obj["state"])
    res = out|> HTTP.Response |> addResponseHeader

    return res
end

function redo!(srv::ServerState, req::HTTP.Request)
    port = srv.obj["IBox_port"]
    localIP = srv.obj["IBox_host"]

    redo!(srv.obj,localIP, port)

    out = "State "*string(srv.obj["state"])
    res = out|> HTTP.Response |> addResponseHeader

    return res
end

"""
Запуск сервера (в асинхронном режиме):

`server, task = start_server()`
`using Sockets
server, task = start_server(dir; localIP = Sockets.localhost, port = 8080)`

dir - локальная директория на сервере для работы с файлами, по умолчанию test\\files\n
localIP - локальный IP сервера, по умолчанию - IP машины в локальной сети, но можно указать просто "localhost"\n
port - порт, по умолчанию 8080
"""
start_server() = start_server("")
function start_server(dir::AbstractString; localIP = Sockets.getipaddr(), port = 8080)
    if isempty(dir)
        dir = joinpath(Base.@__DIR__, "..", "test", "files")
    end
    AllObj = Dict()
    AllObj["history"] = [] #тут храним действия над ВСЕМИ датагруппами подряд
    AllObj["state"] = 0
    AllObj["dataStorage"] = Dict{String,Any}()
    srv = ServerState(dir,AllObj)

    # define REST endpoints to dispatch to "service" functions
    #=const=# H5_ROUTER = HTTP.Router()

    # note the use of `*` to capture the path segment "variables"
    HTTP.@register(H5_ROUTER, "GET", "", handle)
    HTTP.@register(H5_ROUTER, "GET", "/api/runIBox", x->runIBox(srv, x))
    HTTP.@register(H5_ROUTER, "GET", "/api/Close", x->Close(srv, x))
    HTTP.@register(H5_ROUTER, "GET", "/api/getDataTree", x->redirectRequest(srv, x))
    HTTP.@register(H5_ROUTER, "GET", "/api/getData", x->redirectRequest(srv, x))
    HTTP.@register(H5_ROUTER, "GET", "/api/getStructData", x->redirectRequest(srv, x))
    HTTP.@register(H5_ROUTER, "GET", "/api/getTag", x->redirectRequest(srv, x))
    HTTP.@register(H5_ROUTER, "GET", "/api/getType", x->redirectRequest(srv, x))
    HTTP.@register(H5_ROUTER, "GET", "/api/getAttributes", x->redirectRequest(srv, x))
    HTTP.@register(H5_ROUTER, "POST", "/api/manualChange", x->manualChange(srv, x))
    HTTP.@register(H5_ROUTER, "GET", "/api/getFileInfo", x->redirectRequest(srv, x))
    HTTP.@register(H5_ROUTER, "GET", "/api/getQRStree", x->getQRStree(srv, x))
    HTTP.@register(H5_ROUTER, "GET", "/api/undo", x->undo!(srv, x))
    HTTP.@register(H5_ROUTER, "GET", "/api/redo", x->redo!(srv, x))

    # HTTP.@register(H5_ROUTER, "GET", "/api/getAttributes", x->getAttributes(srv, x))
    #
    server = Sockets.listen(localIP, port)
    HTTP.@register(H5_ROUTER, "GET", "/api/closeServer", req->(close(server); HTTP.Response("server closed")))

    task = @async HTTP.serve(H5_ROUTER, localIP, port; verbose = false, server = server)

    sleep(1) # пауза, чтобы вывод текста в консоль не перебивался потоком сервера

    @info "In browser:"
    @info "$localIP:$port/api/runIBox?res=oxy115829.dat&IBox_port=8888&IBox_path='C:/Temp/IBox/IBoxLauncher.exe'&IBox_host=127.0.0.1"
    @info "$localIP:$port/api/getDataTree"
    @info "$localIP:$port/api/getData?dataName=QPoint&index=0&from=0&to=10&count=100"

    @info "$localIP:$port/api/Close"
    @info "$localIP:$port/api/closeServer"

    server, task
end


text:
using Attrs
using Test

mutable struct Foo{Tag, T}
    x::T
end

@inline Foo{Tag}(x::T) where {Tag, T} = Foo{Tag, T}(x)

@defattrs Foo

@inline Attrs.getattr(x::Foo{Tag}, ::Attr{:tag}) where Tag = Tag

@literalattrs Attrs.getattr(x::Foo, ::Attr{:y}) = 2 * x.x

@literalattrs Attrs.getattr(x::Foo, ::Attr{:z}) = 2 * x.x + 1

@literalattrs Attrs.getattr(x::Foo{:even}, ::Attr{:z}) =
    2 * x.x + 2

@literalattrs Attrs.setattr!(x::Foo, ::Attr{:z}, y) = x.x = div(y - 1, 2)

@literalattrs Attrs.setattr!(x::Foo{:even}, ::Attr{:z}, y) =
    x.x = div(y, 2) - 1

h(x) = x.x
@literalattrs h1(x) = x.x

g(x) = x.y
@literalattrs g1(x) = x.y

k(x) = x.z
@literalattrs k1(x) = x.z

@testset "Attrs" begin

    x = Foo{:foo}(42)

    @test h(x) == 42
    @test g(x) == 84
    @test k(x) == 85
    @test h1(x) == 42
    @test g1(x) == 84
    @test k1(x) == 85

    y = Foo{:even}(42)

    @test h(y) == h(x)
    @test g(y) == g(x)
    @test k(y) == 86

    @test begin
        x2 = deepcopy(x)
        x2.x = -1
        x2.x == -1
    end

    @test begin
        x2 = deepcopy(x)
        x2.z = 85
        x2.x == 42
    end

    @test begin
        y2 = deepcopy(y)
        y2.z = 86
        y2.x == 42
    end

    @test_throws MethodError x.a
    @test_throws MethodError x.a = 42
    @test_throws MethodError x.y = 1

end

ext: jl


text:
"""
    This implementation is based on the port from GSL to Python by
    J. Michael Burgess for his pynchrotron package:
    https://github.com/grburgess/pynchrotron
"""

GSL_SQRT_DBL_EPSILON = 1.4901161193847656e-08
GSL_LOG_DBL_MIN = -7.0839641853226408e02

global const synch_c0 = π / sqrt(3.0)

global const c01 = 0.2257913526447274323630976
global const cond1 = 2 * sqrt(2.0) * GSL_SQRT_DBL_EPSILON
global const cond3 = -8.0 * GSL_LOG_DBL_MIN / 7.0

global const synchrotron1_data = [  30.364682982501076273,
                                    17.079395277408394574,
                                    4.560132133545072889,
                                    0.549281246730419979,
                                    0.372976075069301172e-01,
                                    0.161362430201041242e-02,
                                    0.481916772120371e-04,
                                    0.10512425288938e-05,
                                    0.174638504670e-07,
                                    0.22815486544e-09,
                                    0.240443082e-11,
                                    0.2086588e-13,
                                    0.15167e-15
                                  ]

global const synchrotron1_cs_order = 12
global const synchrotron1_cs_a = -1.0
global const synchrotron1_cs_b = 1.0

global const synchrotron1a_data = [ 2.1329305161355000985,
                                    0.741352864954200240e-01,
                                    0.86968099909964198e-02,
                                    0.11703826248775692e-02,
                                    0.1645105798619192e-03,
                                    0.240201021420640e-04,
                                    0.35827756389389e-05,
                                    0.5447747626984e-06,
                                    0.838802856196e-07,
                                    0.13069882684e-07,
                                    0.2053099071e-08,
                                    0.325187537e-09,
                                    0.517914041e-10,
                                    0.83002988e-11,
                                    0.13352728e-11,
                                    0.2159150e-12,
                                    0.349967e-13,
                                    0.56994e-14,
                                    0.9291e-15,
                                    0.152e-15,
                                    0.249e-16,
                                    0.41e-17,
                                    0.7e-18
                                ]

global const synchrotron1a_cs_order = 22
global const synchrotron1a_cs_a = -1.0
global const synchrotron1a_cs_b = 1.0

global const synchrotron2_data = [  0.4490721623532660844,
                                    0.898353677994187218e-01,
                                    0.81044573772151290e-02,
                                    0.4261716991089162e-03,
                                    0.147609631270746e-04,
                                    0.3628633615300e-06,
                                    0.66634807498e-08,
                                    0.949077166e-10,
                                    0.1079125e-11,
                                    0.10022e-13,
                                    0.77e-16,
                                    0.5e-18,
                                ]

global const synchrotron2_cs_order = 11
global const synchrotron2_cs_a = -1.0
global const synchrotron2_cs_b = 1.0

@inline function cheb_eval(coeff::typeof([1.0]), order::Integer, a::Real, b::Real, x::Real)

    d = 0.0
    dd = 0.0

    y = (2x - a - b) / (b - a)
    y2 = 2y

    @inbounds for j = (order+1):-1:2
        temp = d
        d = y2 * d - dd + coeff[j]
        dd = temp
    end

    temp = d
    d = y * d - dd + 0.5 * coeff[1]

    return d
end

@inline function synchrotron_kernel(x::Real)

    if x < cond1

        z = x^(1/3)
        cf = 1 - 8.43812762813205e-01 * z * z
        return 2.14952824153447863671 * z * cf

    elseif x <= 4.0

        px = x^(1/3)
        px11 = px^11
        t = x^2 / 8.0 - 1.0
        result_c1 = cheb_eval( synchrotron1_data,
                               synchrotron1_cs_order,
                               synchrotron1_cs_a,
                               synchrotron1_cs_b,
                               t )

        result_c2 = cheb_eval( synchrotron2_data,
                               synchrotron2_cs_order,
                               synchrotron2_cs_a,
                               synchrotron2_cs_b,
                               t )

        return px * result_c1 - px11 * result_c2 - synch_c0 * x

    elseif x < cond3

        t = (12.0 - x) / (x + 4.0)

        result_c1 = cheb_eval( synchrotron1a_data,
                               synchrotron1a_cs_order,
                               synchrotron1a_cs_a,
                               synchrotron1a_cs_b,
                               t )

        return sqrt(x) * result_c1 * exp(c01 - x)

    else
        return 0.0
    end
end

@inline function ndensity_integral(bound_low::Real, bound_up::Real, norm::Real,
                                  slope::Real, density_in::Real)

    nb = 4π * norm * bound_low^3 / density_in
    density = nb * ( (bound_up/bound_low)^(3.0-slope) - 1.0 ) / ( 3.0 - slope )

	Δq = 1.e-6 # "softening" to avoid divergence in denominator
    if ( (3.0 - Δq) < slope < (3.0 + Δq) )

	    slope_var = ( slope - 3.0 ) / Δq
	    density2 = nb * log(bound_up/bound_low)
	    if ( slope_var != 0.0 )
	      density = density * slope_var + density2 * (1.0 - slope_var)
	    else
	      density = density2
	  	end
	end

	density
end

function calculate_synch_intensity(CReNorm, CReSlope, bounds, bin_width::Real,
                                   B::Real, density::Real, ν0::Real=1.4e9)

    m_e = 9.10953e-28
    c_l = 2.9979e10
    qe  = 4.8032e-10
    erg2eV = 6.242e+11

    E_cntr_prefac = m_e*c_l*c_l * sqrt(2.0/3.0 * (2π * m_e*c_l)/qe /0.29)
    nu_c_prefac = 3.0 * qe / (4π * m_e^3 * c_l^5)
    j_nu_prefac = qe * qe * qe * sqrt(3.0) / (m_e * c_l * c_l)

    LMB_SPECTRAL_CRs = length(CReNorm)

    F     = zeros(LMB_SPECTRAL_CRs)
    F_mid = zeros(LMB_SPECTRAL_CRs)
    E     = zeros(LMB_SPECTRAL_CRs)
    dE    = zeros(LMB_SPECTRAL_CRs)
    J     = zeros(LMB_SPECTRAL_CRs)
    N     = zeros(LMB_SPECTRAL_CRs)
    N_mid = zeros(LMB_SPECTRAL_CRs)
    x     = zeros(LMB_SPECTRAL_CRs)

    E_min = bounds[1] * m_e * c_l^2
    E_max = bounds[end] * m_e * c_l^2
    di = log(E_max/E_min)/LMB_SPECTRAL_CRs

    for i = 1:LMB_SPECTRAL_CRs
        E[i]     = bounds[i] * m_e * c_l^2
        x[i]     = ν0 / (nu_c_prefac * E[i]^2 * B)

		# old version
        # N[i]     = CReNorm[i]
        # bound_up = bounds[i] * 10^(bin_width*0.5)
        # N_mid[i] = CReNorm[i] * (bound_up/bounds[i])^(-CReSlope[i])

		# integrate over half a bin
		bound_up = bounds[i] * 10^(bin_width*0.5)
		N[i]     = ndensity_integral(bounds[i], bound_up, CReNorm[i],
									CReSlope[i], density)

		Norm_mid = CReNorm[i] * (bound_up/bounds[i])^(-CReSlope[i])
        N_mid[i] = ndensity_integral(bound_up, bounds[i+1], Norm_mid,
									CReSlope[i], density)
    end

    dE[1] = E[1] - E_min * exp(-1*di)
    for i = 2:LMB_SPECTRAL_CRs
        dE[i] = E[i] - E[i-1]
    end

    for i = 2:LMB_SPECTRAL_CRs

        K = synchrotron_kernel(x[i])
        F[i] = N[i] * K

        x_mid = 0.5 * (x[i-1] + x[i])
        K_mid = synchrotron_kernel(x_mid)

        F_mid[i] = N_mid[i] * K_mid

        J[i] = dE[i] / 6.0 * (F[i] + F[i-1] + 4*F_mid[i])
    end

    return J
end

max_stars_repo_name: LudwigBoess/GadJet

max_issues_repo_name: LudwigBoess/GadJet

max_forks_repo_name: LudwigBoess/GadJet


text:
# test with the real data set
# Test Bayesian Detctor
include("./lib/metadata.jl")
include("./2bayes_detector_replica.jl")
include("./lib/SignalGenerator.jl")
include("./lib/roc.jl")

using SignalGenerator
using BayesChpOrig
using Metadata
using Roc

function test_bd_octave_vs_julia()
    print(projectdir * "common_data/sig-activ.dat")
    signal = readdlm(projectdir * "common_data/sig-activ.dat", ',')
    x = signal[:,1]
    y = signal[:,2]
    # START to compare with Octave/Matlab
    x = x[1:10:end]
    # END  to compare with Octave/Matlab
    matrix_max_vals, matrix = bayes_detector(x, lambda = 50)
    detections = extract_changes(matrix_max_vals)
    file_name = "test_bd_octave_vs_julia_detections.txt"

    println("... test_bd_octave_vs_julia > save detections into ", outdir * file_name)
    writedlm(outdir * file_name, detections)

    file_name2 = "test_bd_octave_vs_julia_signal.txt"
    writedlm(outdir * file_name2, x)
end
test_bd_octave_vs_julia()

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
module SecondQuantization

using SymbolicUtils
using SymbolicUtils: Sym, Symbolic, Term, Add, Mul, Div, Pow, @ordered_acrule
using SymbolicUtils: isliteral, is_literal_number, has_trig_exp, _name_type, sym_isa
using TermInterface
using Latexify
import Base: (+), (-), (*), (//), (/), (\), (^)

# SymbolicUtils tools
export simplify, arguments, operation

hilbert_space_iterator = 0
include("qoperators.jl")
export @boson, @parameter, commutator

include("utils.jl")

include("basic_operations.jl")

include("normal_ordering.jl")
export normal_order

include("simplify_rules.jl")
export default_simplifier, serial_simplifier, threaded_simplifier, serial_expand_simplifier
export normal_order_simplifier, serial_normal_order_simplifier, threaded_normal_order_simplifier

include("latexify_recipes.jl")

end


text:
using LinearAlgebra, OrdinaryDiffEq, Test
f = (du,u,p,t) -> du .= u ./ t
jac = (J,u,p,t) -> (J[1,1] = 1/t; J[2,2] = 1/t; J[1,2] = 0; J[2,1] = 0)

jp_diag = Diagonal(zeros(2))
fun = ODEFunction(f; jac=jac, jac_prototype=jp_diag)
prob = ODEProblem(fun,ones(2),(1.0,10.0))
sol = solve(prob,Rosenbrock23())
@test sol[end] ≈ [10.0,10.0]
@test length(sol) < 60

sol = solve(prob,Rosenbrock23(autodiff=false))
@test sol[end] ≈ [10.0,10.0]
@test length(sol) < 60

jp = Tridiagonal(jp_diag)
fun = ODEFunction(f; jac=jac, jac_prototype=jp)
prob = ODEProblem(fun,ones(2),(1.0,10.0))

sol = solve(prob,Rosenbrock23())
@test sol[end] ≈ [10.0,10.0]
@test length(sol) < 60

sol = solve(prob,Rosenbrock23(autodiff=false))
@test sol[end] ≈ [10.0,10.0]
@test length(sol) < 60

#=
jp = SymTridiagonal(jp_diag)
fun = ODEFunction(f; jac=jac, jac_prototype=jp)
prob = ODEProblem(fun,ones(2),(1.0,10.0))
sol = solve(prob,Rosenbrock23())
@test sol[end] ≈ [10.0,10.0]
@test length(sol) < 60
=#

# Don't test the autodiff=false version here because it's not as numerically stable,
# so lack of optimizations would lead to unsymmetric which causes an error:
# LoadError: ArgumentError: broadcasted assignment breaks symmetry between locations (1, 2) and (2, 1)

@test_broken begin
    jp = Hermitian(jp_diag)
    fun = ODEFunction(f; jac=jac, jac_prototype=jp)
    prob = ODEProblem(fun,ones(2),(1.0,10.0))
    sol = solve(prob,Rosenbrock23(autodiff=false))
    @test sol[end] ≈ [10.0,10.0]
    @test length(sol) < 60
end

@test_broken begin
    jp = Symmetric(jp_diag)
    fun = ODEFunction(f; jac=jac, jac_prototype=jp)
    prob = ODEProblem(fun,ones(2),(1.0,10.0))
    sol = solve(prob,Rosenbrock23(autodiff=false))
    @test sol[end] ≈ [10.0,10.0]
    @test length(sol) < 60
end


text:
using Plasmo
using Plots, PlasmoPlots
gr()

graph1 = OptiGraph()

@optinode(graph1,n1)
@variable(n1, x >= 0)
@variable(n1, y >= 2)
@constraint(n1,x + y >= 3)
@objective(n1, Min, y)

@optinode(graph1,n2)
@variable(n2,x >= 0)
@variable(n2, y >= 2)
@constraint(n2,x + y >= 3)
@objective(n2, Min, y)

@optinode(graph1,n3)
@variable(n3,x >= 0)
@variable(n3, y >= 2)
@constraint(n3,x + y >= 3)
@objective(n3, Min, y)

@linkconstraint(graph1, n1[:x] + n2[:x] + n3[:x] == 3)

graph2 = OptiGraph()

@optinode(graph2,n4)
@variable(n4, x >= 0)
@variable(n4, y >= 2)
@constraint(n4,x + y >= 5)
@objective(n4, Min, y)

@optinode(graph2,n5)
@variable(n5, x >= 0)
@variable(n5, y >= 2)
@constraint(n5,x + y >= 5)
@objective(n5, Min, y)

@optinode(graph2,n6)
@variable(n6, x >= 0)
@variable(n6, y >= 2 )
@constraint(n6,x + y >= 5)
@objective(n6, Min, y)

@linkconstraint(graph2, n4[:x] + n5[:x] + n6[:x] == 5)

graph3 = OptiGraph()

@optinode(graph3,n7)
@variable(n7, x >= 0)
@variable(n7, y >= 2)
@constraint(n7,x + y >= 7)
@objective(n7, Min, y)

@optinode(graph3,n8)
@variable(n8, x >= 0)
@variable(n8, y >= 2)
@constraint(n8,x + y >= 7)
@objective(n8, Min, y)

@optinode(graph3,n9)
@variable(n9,x >= 0)
@variable(n9, y >= 2)
@constraint(n9,x + y >= 7)
@objective(n9, Min, y)

@linkconstraint(graph3, n7[:x] + n8[:x] + n9[:x] == 7)

graph0 = OptiGraph()
@optinode(graph0,n0)
@variable(n0,x)
@constraint(n0,x >= 0)

add_subgraph!(graph0,graph1)
add_subgraph!(graph0,graph2)
add_subgraph!(graph0,graph3)
@linkconstraint(graph0,n3[:x] + n5[:x] + n7[:x] == 10)

@linkconstraint(graph0,n0[:x] + n3[:x] == 3)
@linkconstraint(graph0,n0[:x] + n5[:x] == 5)
@linkconstraint(graph0,n0[:x] + n7[:x] == 7)

for (i,node) in enumerate(all_nodes(graph0))
    node.label = "n$i"
end

plt_graph2 = layout_plot(graph0,node_labels = true,markersize = 60,labelsize = 30,linewidth = 4,subgraph_colors = true,
layout_options = Dict(:tol => 0.001,:C => 2, :K => 4, :iterations => 5));

plt_matrix2 = matrix_plot(graph0,node_labels = true,subgraph_colors = true,markersize = 16);


text:
export stask_run, start_stasks, count_stasks, stop_stasks, pmap_stasks

typealias STasks Array{RemoteRef, 1}

function stask_run(wrr::RemoteRef)
    msg = take(wrr)

    while(true)
        # execute the function
#        println(msg)
        (msg_type, resp_rr, f, args) = msg
        local response
        if (msg_type == :execfun)
            try
                response = apply(f, args)
                put(resp_rr, (:ok, response))
            catch e
                println(e)
                put(resp_rr, (:err, "exception : " * string(e)))
            end
        elseif (msg_type == :exit)
            put(resp_rr, :ok)
            return
        end

        # wait for the next request
        msg = take(wrr)
    end
end

function start_stasks(shmcfg=false, shmpfx=false)
    # Creates one persistent task per Worker
    # First create remote refs for each of them...this will be used to comunicate with the tasks...
    # If nprocs > 1 then only the workers will be used for task processing.

    np = nprocs()

    wp_start = (np > 1) ? 2 : 1
    pid_add = (np > 1) ? 1 : 0

    stasks = [RemoteRef(i)::RemoteRef for i in wp_start:np]

    # start the tasks
    [remotecall(i+pid_add, stask_run, stasks[i]) for i in 1:length(stasks)]

    if shmcfg != false
        unlink_shm(shmcfg, shmpfx)
        # local
        setup_shm(shmcfg, shmpfx)

        # Map each symbol into the stasks task_local_storage()
        pmap_stasks(stasks, () -> setup_shm(shmcfg, shmpfx, false))
    end

    stasks::STasks
end

function stop_stasks(stasks::STasks, shmcfg=false, shmpfx=false)
    srr_arr = cell(length(stasks))
    i = 1
    for rr in stasks
        srr_arr[i] = RemoteRef()
        put(rr, (:exit, srr_arr[i], nothing, ()))
        i += 1
    end

    # Just unlink shmcfg is present
    if (shmcfg != false)
        unlink_shm(shmcfg, shmpfx)
    end
end

function pmap_stasks(stasks::STasks, f::Function, lsts...)
    nw = length(stasks)
    if lsts != ()
        if nw != length(lsts[1]) error ("dimension of persistent tasks must be same as jobs") end
    end

    srr_arr = cell(nw)
    i = 1
    for rr in stasks
        srr_arr[i] = RemoteRef()
        args = map(L -> L[i], lsts)
        put(rr, (:execfun, srr_arr[i], f, args))
        i += 1
    end

    resplist = [take(srr) for srr in srr_arr]
    statuses = [e[1] for e in resplist]
    responses = [e[2] for e in resplist]

    if all(x -> x == :ok, statuses) == true
        return (:ok, responses)
    else
        return (statuses, responses)
    end
end

count_stasks(stasks::STasks) = length(stasks)


text:
#=
  Beta comparison

  From infer.net src/tests/Test/BlogTests.cs
  """
  a = Beta(11,500)[mean=0.02153], b = Beta(16,695)[mean=0.0225]
  aGreaterThanConstant = Bernoulli(0.9849)
  P(A > B) =  Bernoulli(0.4467)
  """

  (Yet another A/B test.)

  Summary Statistics
         parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
             Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

              aRate    0.0215    0.0065     0.0001    0.0002   1208.1553    1.0015       35.7887
              bRate    0.0238    0.0056     0.0001    0.0002   1433.9797    1.0000       42.4782
               aGtB    0.3834    0.4862     0.0049    0.0174    709.5007    0.9999       21.0173
            aGr0_01    0.9853    0.1204     0.0012    0.0031   1829.2015    1.0001       54.1857
  aSuccessCountPost   10.7714    4.5923     0.0459    0.0939   2137.7235    1.0006       63.3249
  bSuccessCountPost   16.5616    5.5983     0.0560    0.1184   2440.4334    1.0001       72.2920

  Cf ~/webppl/beta_comparison.wppl

=#

using Turing
include("jl_utils.jl")

@model function beta_comparison(aTrialCount=400,aSuccessCount=10,bTrialCount=700,bSuccessCount=15)
    aRate ~ Beta(1,10)
    bRate ~ Beta(1,10)

    aSuccessCount ~ Binomial(aTrialCount,aRate)
    bSuccessCount ~ Binomial(bTrialCount,bRate)

    aGtB ~ Dirac(aRate > bRate)
    aGt0_01 ~ Dirac(aRate > 0.01)

    aSuccessCountPost ~ Binomial(aTrialCount,aRate)
    bSuccessCountPost ~ Binomial(bTrialCount,bRate)
end

aTrialCount=500
aSuccessCount=10
bTrialCount=700
bSuccessCount=16
model = beta_comparison(aTrialCount,aSuccessCount,bTrialCount,bSuccessCount)

num_chns = 4
# chs = sample(model, Prior(), 1000)
# chs = sample(model, MH(), 10_000)
# chs = sample(model, PG(15), 10_000)
# chs = sample(model, IS(), 10_000)
chs = sample(model, SMC(), 10_000)
# chs = sample(model, SMC(), MCMCThreads(), 10_000, num_chs)

# chs = sample(model, SGLD(), 10_000)
# chs = sample(model,NUTS(), 10_000)
# chs = sample(model,HMC(0.01,5), 10_000)

display(chs)
# display(plot(chs))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# script for arbitrary computation of EE value

"""
    computevalue(session::EESession, value::EE.AbstractEEObject)

Fuction to request data from any arbitrary EarthEngine computation and return as
the appropriate Julia object.
See https://developers.google.com/earth-engine/reference/rest/v1beta/projects.value/compute
"""
function computevalue(session::EESession, value::EE.AbstractEEObject)
    endpoint = "value:compute"

    serialized = ee.serializer.encode(value, for_cloud_api = true)

    payload = Dict(:expression => serialized)

    response = _sendrequest(session, endpoint, payload)

    return response.json()["result"]
end


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
using DelimitedFiles
Random.seed!(233)

function halo_exchange_two_d(u,fill_value,m,n)
    halo_exchange_two_d_ = load_op_and_grad("./build/libHaloExchangeTwoD","halo_exchange_two_d")
    u,fill_value,m,n = convert_to_tensor(Any[u,fill_value,m,n], [Float64,Float64,Int64,Int64])
    halo_exchange_two_d_(u,fill_value,m,n)
end

mpi_init()
U = reshape(1:24, 4, 6)'|>Array

m = 3
n = 2
fill_value = 1.0

M = mpi_rank()÷n+1
N = mpi_rank()%n+1
ulocal = U[(M-1)*2 + 1: M * 2, (N-1)*2+1:N*2]
# TODO: specify your input parameters
u = halo_exchange_two_d(ulocal,fill_value,m,n)
sess = Session(); init(sess)
Uval = run(sess, u)

sleep(mpi_rank())
println("=========================== rank = $(mpi_rank()) ======================")
writedlm(stdout, Uval)
mpi_finalize()


text:
### A Pluto.jl notebook ###
# v0.14.2

using Markdown
using InteractiveUtils

# ╔═╡ 17ada06e-a49a-11eb-29f8-35891d4afe6c
begin
	using Random
	using Plots
	using Base.Threads
	using Distributions
	using DataFrames
end

# ╔═╡ d79fcb85-28ca-4e32-bc4d-8c442b771cd6
md"""
# Random Election Upsets Riddler

This notebook solves the [Riddler Classic problem from April 23rd 2021](https://fivethirtyeight.com/features/can-you-cut-the-perfect-pancake/), which asks us to consider a large nation of voters who independently submit their vote between two candidates (A and B) by flipping identical fair coins. We're asked: if 80% of votes are tallied on election day (regular votes), while the remaining 20% (early votes) are tallied later, what is the probability of an upset -- i.e. what is the probability that the election day outcome doesn't match the final outcome?

As usual, I'll solve this problem in two ways. First I'll compute the result using some gratuitous math, and then I'll run a simulation of the random process to cross-check my analytical results. Since this problem isn't substantially harder if we assume that 20% of votes are submitted early, or that an arbitrary fraction of votes are submitted early, I'll solve the prompt in general and then report the solution for the prompt as a special case.
"""

# ╔═╡ 7c4713bb-7e75-4c0b-815c-e817b71e2f8b
md"""
## Analytical Solution

### Notation
First let's introduce a bit of notation. Throughout, we will use the subscripts $E$ and $R$ to indicate quantities computed respectively for the early voter and regular voter subpopulations. Similarly, we will use a subscript $T$ to refer to the total population. We will use $v$ to denote the number of votes in favor of candidate A in a population and $n$ will denote the number of voters in a population. E.g. $v_E$ is the number of voters who voted for candidate A in the early voter subpopulation, and $n_E$ is the number of voters. Finally, we will use $f$ to denote the fracction of voters who belong to a subpopulation, e.g. $f_E = n_E / n_T$.

### General Solution
Since each voter places their vote by independently flipping a fair coin, the total number of votes for candidate A is distributed as $v_\text{T} \sim \text{Binomial}(n_T,~0.5)$, and this remains true when we sub-divide into regular and early voter sub-populations. That is,

$$\begin{array}{ccl}
v_R &\sim& \text{Binomial}(n_R,~0.5). \\
v_E &\sim& \text{Binomial}(n_E,~0.5), \\
\end{array}$$

The prompt tells us that the population is "large," which is a hint that we should  think about using the Normal approximation to the Binomial. To that end, we recall that if $m \sim \text{Binomial}(n, p)$, then $m$ is approximately $\text{Normal}(\mu_m,\sigma_m)$, where $\mu_m = E[m] = pn$, and $\sigma_m = \sqrt{np(1-p)}$. Hence,

$$\begin{array}{ccl}
v_R &\approx& \text{Normal}(\mu_R,\sigma_R). \\
v_E &\approx& \text{Normal}(\mu_E,\sigma_E), \\
\end{array}$$

where $\mu_R = 0.5n_R$, $\sigma_R = 0.5 \sqrt{n_R}$ and an analogous expressions hold for $\mu_E$ and $\sigma_E$.

An election upset will occur if candidate A wins the regular vote but loses the total vote, or if candidate A loses the regular vote and wins the total vote. Since both events occur with equal probability (since candidate A and B are interchangeable), it follows that

$$P(\text{upset}) = 2 \cdot P(\text{A loses the total vote} \cap \text{A wins the regular vote}).$$

We can then re-write this as,

$$\begin{eqnarray}
P(\text{upset}) &=& 2 \cdot P(\text{A wins the regular vote} \cap \text{A loses the total vote}) \\
&=& 2 \cdot P(v_R > 0.5 \cdot n_R \text{ and }  v_T < 0.5 \cdot n_T)\\
&=& 2 \cdot P(v_R > 0.5 \cdot n_R \text{ and } v_E + v_R < 0.5 \cdot (n_E + n_R)) \\
&=& 2 \cdot P(v_R > \mu_R \text{ and } v_E + v_R < \mu_E+\mu_R) \\
&=& 2 \cdot P(v_R > \mu_R \text{ and } v_E < \mu_E + \mu_R - v_R). \\
\end{eqnarray}$$

Here we have used the fact that candidate A wins a round of votes if $v > 0.5\cdot n$, and used the expressions for $\mu_R, \mu_E$ from above.

Next, we note that $v_E$ and $v_R$ are independent so their joint density is just the product of their respective densities. Using this obsevation and the fact that both variables are normally distributed we obtain:

$$P(\text{upset}) = 2 \cdot \int_{\mu_R}^\infty \,dv_R \int_{-\infty}^{\mu_R + \mu_E - v_R}\,dv_E
\frac{e^{-\frac{1}{2}(\frac{v_R-\mu_R}{\sigma_R})^2}}{\sigma_R\sqrt{2\pi}}
\frac{e^{-\frac{1}{2}(\frac{v_E-\mu_E}{\sigma_E})^2}}{\sigma_E\sqrt{2\pi}}.$$

We will now manipulate this expression so that it matches an expression that I managed to find in a table of integrals 🙂. To that end, we make the substitution $w_R = \frac{v_R-\mu_R}{\sigma_R\sqrt{2}}$, $w_E = \frac{v_E -\mu_E}{\sigma_E\sqrt{2}}$, which (after some book-keeping) yields:

$$P(\text{upset}) = 2 \cdot \int_{0}^\infty \,dw_R \frac{e^{-w_R^2}}{\sqrt{\pi}} \int_{-\infty}^{-\frac{\sigma_R}{\sigma_E}w_R}\,dw_E \frac{e^{-w_E^2}}{\sqrt{\pi}}$$

Next, we recall that $$\text{erfc}(x) = \frac{2}{\sqrt{\pi}} \int_x^\infty e^{-t^2}\,dt$$, and that $$\text{erf}(x) = 1 - \text{erfc}(x)$$ so the above can be written:

$$\begin{eqnarray}
P(\text{upset}) &=&\int_{0}^\infty \,dw_R \frac{e^{-w_R^2}}{\sqrt{\pi}} \left(1 - \text{erf}(\frac{\sigma_R}{\sigma_E}w_R)\right)\\
&=&\int_{0}^\infty \,dw_R \frac{e^{-w_R^2}}{\sqrt{\pi}} - \frac{1}{\sqrt{\pi}} \cdot \int_{0}^\infty\,dw_R~e^{-w_R^2} \cdot \text{erf}(\frac{\sigma_R}{\sigma_E}w_R)
\end{eqnarray}$$

The first integral is equal to $\frac{1}{2}\text{erfc}(0) = \frac{1}{2}$. For the second integral, we apply formula 2 from section 4.3 of [this table of integrals involving error functions](https://nvlpubs.nist.gov/nistpubs/jres/73b/jresv73bn1p1_a1b.pdf) with $a = \sigma_R/\sigma_E$ and $b = 1$, to obtain:

$$\begin{eqnarray}
P(\text{upset})&=&\frac{1}{2} - \frac{1}{\sqrt{\pi}} \left[\frac{\sqrt{\pi}}{2} - \frac{1}{\sqrt{\pi}}\tan^{-1}\left(\frac{\sigma_E}{\sigma_R}\right)\right]\\
&=&\frac{1}{\pi}\tan^{-1}\left(\frac{\sigma_E}{\sigma_R}\right)
\end{eqnarray}$$

Finally, since $\sigma_E = \sqrt{\frac{1}{2}f_E\cdot n_\text{total}}$, $\sigma_R = \sqrt{\frac{1}{2}f_R\cdot n_\text{total}}$, and $f_E + f_R = 1$, we conclude that:

$$\boxed{P(\text{upset})=\frac{1}{\pi}\tan^{-1}\left(\sqrt{\frac{f_E}{1-f_E}}\right)}$$

### Solution for the Prompt
In particular, for the case outlined in prompt with $f_E = 0.2$, we find:

$$\boxed{P\left(\text{upset}_\text{prompt}\right) = \frac{1}{\pi} \tan^{-1}\left(\frac{1}{2}\right) \approx 0.14758}$$

"""

# ╔═╡ 70ef11be-81b9-4e4d-b83a-f7a225a9001b
function exact_prob_upset(frac_early)
	return (1/π) * atan(sqrt(frac_early / (1 - frac_early)))
end

# ╔═╡ 24e2efe2-f6f2-4648-8c28-bffedf98095b
exact_prob_upset(0.2)

# ╔═╡ 2383b70b-0693-4adf-99d4-9bcedec5d6cb
md"""
## Simulation Solution

As discussed above, the number of votes for candidate A in the early and regular vote subpopulations are both Binomial random variables. So, we  simulate the election tally by simulating appropriate Binomial random variables. From the simulated results, we next determine if A wins on election night (i.e. if more than 50% of the regular vote is for A) and determine the overall winner (i.e. if more than 50% of the total vote is for A). An election upset results if A wins in either the regular or total tally, but loses in the other tally. We simulate this process 100K times and tally up the fraction of times that an election upset occurred.

This process is repeated for 1001 evenly spaced early voter share fractions in [0,1], and the result is then compared (graphically) against the analytical solution.
"""

# ╔═╡ 562b2637-ec4f-4be2-a525-09e83c6458fa
"""
    simulate_election(frac_early, n_simulations, [n_country, p_A])

	Inputs:
	frac_early - percentage of early voters, 20% in the problem prompt
	n_simulations - number of simulations to run

	Optional:
	n_country - Has to be big-enough™. defaults to 100_001 votes
	p_A - probability that a voter selects candidate A. defaults to 0.5

	Output:
	estimated probability of an upset to the election night results
"""
function simulate_election(frac_early, n_simulations, n_country = 100_001, p_A=0.5)

	# number of early and regular voters
	n_early = floor(Int, frac_early * n_country)
	n_regular = n_country - n_early

	# distributions of votes for A for early and regular voters
	d_early = Binomial(n_early, p_A)
	d_regular = Binomial(n_regular, p_A)

	# simulate votes for A for early, regular, and total
	n_votes_A_early = rand(d_early, n_simulations)
	n_votes_A_regular = rand(d_regular, n_simulations)
	n_votes_A = n_votes_A_early + n_votes_A_regular

	# fraction of vote for A between regular and all votes
	frac_votes_A_regular = n_votes_A_regular / n_regular
	frac_votes_A = n_votes_A / n_country

	# determine if A wins on election night and if A wins overall
	a_wins_regular = frac_votes_A_regular .> 0.5
	a_wins = frac_votes_A .> 0.5

	# an upset occurs if
	#  1. A wins  regular, but loses overall (i.e. True, False)
	#  2. A loses regular, but wins overall (i.e. False, True)
	upset = a_wins_regular .!= a_wins
	prob_upset = sum(upset) / n_simulations

	return prob_upset
end

# ╔═╡ 5f82ec16-38a7-4c44-aa5a-f76d18d521dd
begin
	n_sim = 100_000
	n_pop = 10_00_001
	frac_earlys = LinRange(0, 1, 1_001)

	prob_upset_sim = map(fe->simulate_election(fe, n_sim, n_pop), frac_earlys)
	prob_upset_exact = map(exact_prob_upset, frac_earlys)
end

# ╔═╡ b2346518-43b5-4c4d-9536-244514a4d5a3
plot(
	frac_earlys,
	hcat(prob_upset_exact, prob_upset_sim),
	label=["Exact" "Simulated"],
	legend=:bottomright,
	linewidth=[5 2],
	xlabel="Fraction of Early Voters",
	ylabel="Probability of Election Upset",
	title="Probability of Election Upset \nvs Fraction of Early Voters",
	size=(800,500)
)

# ╔═╡ 06cdba4f-426d-47e5-a4c0-c35f9a41d547
md"""
The above figure shows great agreement between the exact solution discussed above and the simulation results -- so I'd wager that my exact solution is probably correct! 🙂
"""

# ╔═╡ 4a8020c4-b9f4-4c22-8919-eb0f135b1725
savefig("election_upset_probability.png")

# ╔═╡ e1bfd335-b818-4931-8e71-03b46c90c9ec
md"""
### Simulating the Prompt Solution
Finally, we'll simulate the prompt problem 10 million times to get a bit more accuracy in our simulation.
"""

# ╔═╡ e55c8934-7268-4c5f-af6c-2098a48cc094
begin
	Random.seed!(2021)
	simulate_election(0.2, 100_000_000, n_pop)
end

# ╔═╡ 537e2e5d-378b-4676-8ba4-3a50f54eb84c
md"""
Setting the seed for reproducibility (but not selecting it to cherry-pick!), the simulation result agrees with the theoretical result to around four decimal places, which is pretty good!
"""

# ╔═╡ Cell order:
# ╟─d79fcb85-28ca-4e32-bc4d-8c442b771cd6
# ╟─7c4713bb-7e75-4c0b-815c-e817b71e2f8b
# ╠═70ef11be-81b9-4e4d-b83a-f7a225a9001b
# ╠═24e2efe2-f6f2-4648-8c28-bffedf98095b
# ╟─2383b70b-0693-4adf-99d4-9bcedec5d6cb
# ╠═562b2637-ec4f-4be2-a525-09e83c6458fa
# ╠═5f82ec16-38a7-4c44-aa5a-f76d18d521dd
# ╠═b2346518-43b5-4c4d-9536-244514a4d5a3
# ╟─06cdba4f-426d-47e5-a4c0-c35f9a41d547
# ╟─4a8020c4-b9f4-4c22-8919-eb0f135b1725
# ╟─e1bfd335-b818-4931-8e71-03b46c90c9ec
# ╠═e55c8934-7268-4c5f-af6c-2098a48cc094
# ╟─537e2e5d-378b-4676-8ba4-3a50f54eb84c
# ╠═17ada06e-a49a-11eb-29f8-35891d4afe6c


text:
using HSA
include("helpers.jl")
using FactCheck

facts("The Agents") do
    rt = Runtime()
    agents = Array(HSA.Agent,0)

    context("Can be iterated over") do
        context("with a callback returning true to continue") do
            HSA.iterate_agents(a -> begin
                push!(agents, a)
                true # continue
            end)

            for a in agents
                @fact isa(a, HSA.Agent) --> true
            end
        end

        context("with a callback returning not(nothing) to continue") do
            count = 0

            HSA.iterate_agents(a -> begin
                count += 1
            end)

            @fact count --> size(agents,1)
        end

        context("with a callback returning false to terminate early") do
            count = 0

            HSA.iterate_agents(a -> begin count += 1; return false end)

            @with_agents @fact count --> 1
        end

        context("rethrowing exceptions from the callback") do
            @with_agents @fact_throws HSA.iterate_agents(rt, a -> error("testerror"))
        end
    end

    context("can all be retrieved as an array") do
        count = 0

        HSA.iterate_agents(a -> begin count += 1 end)

        agents = HSA.all_agents()

        @fact size(agents,1) --> count
    end

    context("can be retrieved by device type") do
        count = 0
        t = HSA.HSA_DEVICE_TYPE_CPU

        HSA.iterate_agents(a -> begin
            if HSA.agent_info_device(a) & t == t
                count += 1
            end
        end)

        agents = HSA.all_agents(dev = t)

        @fact size(agents,1) --> count
    end

    context("can be queried for information") do

        @with_agents context("property by property") do
            a = agents[1]

            name = HSA.agent_info_name(a)
            @fact length(name) --> greater_than(1)
            @fact HSA.agent_info_name(a) --> not(nothing)
            @fact HSA.agent_info_vendor_name(a) --> not(nothing)
            @fact HSA.agent_info_feature(a) --> not(nothing)
            @fact HSA.agent_info_wavefront_size(a) --> not(nothing)

            wg_dim = HSA.agent_info_workgroup_max_dim(a)
            @fact wg_dim --> not(nothing)
            @fact typeof(wg_dim) --> Tuple{UInt16, UInt16, UInt16}

            @fact HSA.agent_info_workgroup_max_size(a) --> not(nothing)
            @fact HSA.agent_info_grid_max_dim(a) --> not(nothing)
            @fact HSA.agent_info_grid_max_size(a) --> not(nothing)
            @fact HSA.agent_info_fbarrier_max_size(a) --> not(nothing)
            @fact HSA.agent_info_queues_max(a) --> not(nothing)
            @fact HSA.agent_info_queue_max_size(a) --> not(nothing)
            @fact HSA.agent_info_queue_type(a) --> not(nothing)
            @fact HSA.agent_info_node(a) --> not(nothing)
            @fact HSA.agent_info_device(a) --> not(nothing)
            @fact HSA.agent_info_cache_size(a) --> not(nothing)
        end

        @with_agents context("collectively for all basic agent info") do
            a = agents[1]

            @fact typeof(HSA.agent_info(a)) --> HSA.AgentInfo
        end
    end

    finalize(rt)
end


text:
"""
A file with structs and functions for setting up and using the topology of the spatial system (e.g. a graph or a Cartesian grid)
"""

################### LightGraph ########################
num_sites(graph::AbstractGraph) = LightGraphs.nv(graph)
# neighbors(graph::AbstractGraph, site) = LightGraphs.neighbors(graph, site)
# outdegree(graph::AbstractGraph, site) = LightGraphs.outdegree(graph, site)
rand_nbr(rng, graph::AbstractGraph, site) = rand(rng, neighbors(graph, site))
nth_nbr(graph::AbstractGraph, site, n) = @inbounds neighbors(graph, site)[n]

################### CartesianGrid ########################
const offsets_1D = [CartesianIndex(-1),CartesianIndex(1)]
const offsets_2D = [CartesianIndex(0,-1),CartesianIndex(-1,0),CartesianIndex(1,0),CartesianIndex(0,1)]
const offsets_3D = [CartesianIndex(0,0,-1), CartesianIndex(0,-1,0),CartesianIndex(-1,0,0),CartesianIndex(1,0,0),CartesianIndex(0,1,0),CartesianIndex(0,0,1)]

"""
    potential_offsets(dimension::Int)

NOTE: dimension is assumed to be 1, 2, or 3
"""
function potential_offsets(dimension::Int)
    if dimension==1
        return offsets_1D
    elseif dimension==2
        return offsets_2D
    else # else dimension == 3
        return offsets_3D
    end
end

dimension(grid) = length(grid.dims)
num_sites(grid) = prod(grid.dims)
outdegree(grid, site) = grid.nums_neighbors[site]

"""
    nth_potential_nbr(grid, site, n)

return nth Cartesian neighbor ignoring boundaries
"""
nth_potential_nbr(grid, site, n) = grid.LI[grid.offsets[n]+grid.CI[site]]

"""
    nth_nbr(grid, site, n)

return the nth neighbor of site in grid, in ascending order
"""
function nth_nbr(grid, site, n)
    CI = grid.CI; offsets = grid.offsets
    @inbounds I = CI[site]
    @inbounds for off in offsets
        nb = I + off
        if nb in CI
            n -= 1
            @inbounds n == 0 && return grid.LI[nb]
        end
    end
end

"""
    neighbors(grid, site)

return an iterator over neighbors of site in ascending order. Do not use in hot loops
"""
function neighbors(grid, site)
    CI = grid.CI
    LI = grid.LI
    I = CI[site]
    Iterators.map(off -> LI[off+I], Iterators.filter(off -> off+I in CI, grid.offsets))
end

"""
given a vector of hopping constants of length num_neighbors(grid, site), form the vector of size 2*(dimension of grid) with zeros at indices where the neighbor is out of bounds. Store it in to_pad
"""
function pad_hop_vec!(to_pad::AbstractVector{F}, grid, site, hop_vec::Vector{F}) where F <: Number
    CI = grid.CI
    I = CI[site]
    nbr_counter = 1
    for (i,off) in enumerate(grid.offsets)
        if I + off in CI
            to_pad[i] = hop_vec[nbr_counter]
            nbr_counter += 1
        else
            to_pad[i] = zero(F)
        end
    end
    to_pad
end

CartesianGrid(dims) = CartesianGridRej(dims) # use CartesianGridRej by default

# neighbor sampling is rejection-based
struct CartesianGridRej{N,T}
    "dimensions (side lengths) of the grid"
    dims::NTuple{N, Int}

    "number of neighbor for each site"
    nums_neighbors::Vector{Int}
    CI::CartesianIndices{N, T}
    LI::LinearIndices{N, T}

    "offsets, e.g. [-1, 1] for 1D"
    offsets::Vector{CartesianIndex{N}}
end

"""
    CartesianGridRej(dims::Tuple)

initialze CartesianGridRej
"""
function CartesianGridRej(dims::Tuple)
    dim = length(dims)
    CI = CartesianIndices(dims)
    LI = LinearIndices(dims)
    offsets = potential_offsets(dim)
    nums_neighbors = [count(x -> x+CI[site] in CI, offsets) for site in 1:prod(dims)]
    CartesianGridRej(dims, nums_neighbors, CI, LI, offsets)
end
CartesianGridRej(dims) = CartesianGridRej(Tuple(dims))
CartesianGridRej(dimension, linear_size::Int) = CartesianGridRej([linear_size for i in 1:dimension])
function rand_nbr(rng, grid::CartesianGridRej, site::Int)
    CI = grid.CI; offsets = grid.offsets
    @inbounds I = CI[site]
    while true
        @inbounds nb = rand(rng, offsets) + I
        @inbounds nb in CI && return grid.LI[nb]
    end
end

# neighbor sampling is iterator-based
struct CartesianGridIter{N,T}
    dims::NTuple{N, Int}
    nums_neighbors::Vector{Int}
    CI::CartesianIndices{N, T}
    LI::LinearIndices{N, T}
    offsets::Vector{CartesianIndex{N}}
end
function CartesianGridIter(dims::Tuple)
    dim = length(dims)
    CI = CartesianIndices(dims)
    LI = LinearIndices(dims)
    offsets = potential_offsets(dim)
    nums_neighbors = [count(x -> x+CI[site] in CI, offsets) for site in 1:prod(dims)]
    CartesianGridIter(dims, nums_neighbors, CI, LI, offsets)
end
CartesianGridIter(dims) = CartesianGridIter(Tuple(dims))
CartesianGridIter(dimension, linear_size::Int) = CartesianGridIter([linear_size for i in 1:dimension])
function rand_nbr(rng, grid::CartesianGridIter, site::Int)
    nth_nbr(grid, site, rand(rng, 1:outdegree(grid,site)))
end

function Base.show(io::IO, ::MIME"text/plain", grid::Union{CartesianGridRej, CartesianGridIter})
    println(io, "A Cartesian grid with dimensions $(grid.dims)")
end


text:
export Circle

##  Circle

# T Must be in an Algebra
"""
    Circle(c,r,o)

represents the circle centred at `c` with radius `r` which is positively (`o=true`)
or negatively (`o=false`) oriented.
"""
struct Circle{T,V<:Real,TT} <: PeriodicDomain{TT}
	center::T
	radius::V
	orientation::Bool

end

Circle(c::Number,r::Real,o::Bool) = Circle{typeof(c),typeof(r),Complex{typeof(r)}}(c,r,o)
Circle(c::Vec,r::Real,o::Bool) = Circle{typeof(c),typeof(r),typeof(c)}(c,r,o)

Circle(::Type{T1},c::T2,r::V,o::Bool) where {T1,T2,V<:Real} = Circle(convert(promote_type(T1,T2,V),c),
															  convert(promote_type(real(T1),real(T2),V),r),o)
Circle(::Type{T1},c,r::Bool) where {T1<:Number} = Circle(T1,c,r)
Circle(::Type{T1},c,r::Real) where {T1<:Number} = Circle(T1,c,r)
Circle(c,r::Real) = Circle(c,r,true)
Circle(r::Real) = Circle(zero(r),r)
Circle(r::Int) = Circle(Float64,0.,r)
Circle(a::Tuple,r::Real) = Circle(Vec(a...),r)

Circle(::Type{V}) where {V<:Real} = Circle(one(V))
Circle() = Circle(1.0)

isambiguous(d::Circle{T}) where {T<:Number} = isnan(d.center) && isnan(d.radius)
isambiguous(d::Circle{T}) where {T<:Vec} = all(isnan,d.center) && isnan(d.radius)
convert(::Type{Circle{T,V}},::AnyDomain) where {T<:Number,V<:Real} = Circle{T,V}(NaN,NaN)
convert(::Type{IT},::AnyDomain) where {IT<:Circle} = Circle(NaN,NaN)

function tocanonical(d::Circle{T},ζ) where T<:Number
    v=mappoint(d,Circle(),ζ)
    atan(imag(v)-0.0,real(v))  # -0.0 to get branch cut right
end

function tocanonical(d::Circle{T},ζ) where T<:Vec
    v=mappoint(d,Circle((0.0,0.0),1.0),ζ)
    atan(v[2]-0.0,v[1])  # -0.0 to get branch cut right
end

fromcanonical(d::Circle{T,V,Complex{V}},θ) where {T<:Number,V<:Real} =
	d.radius*exp((d.orientation ? 1 : -1)*1.0im*θ) + d.center
fromcanonicalD(d::Circle{T},θ) where {T<:Number} =
	(d.orientation ? 1 : -1)*d.radius*1.0im*exp((d.orientation ? 1 : -1)*1.0im*θ)

fromcanonical(d::Circle{T},θ::Number) where {T<:Vec} =
	d.radius*Vec(cos((d.orientation ? 1 : -1)*θ),sin((d.orientation ? 1 : -1)*θ)) + d.center
fromcanonicalD(d::Circle{T},θ::Number) where {T<:Vec} =
	d.radius*(d.orientation ? 1 : -1)*Vec(-sin((d.orientation ? 1 : -1)*θ),cos((d.orientation ? 1 : -1)*θ))

indomain(z,d::Circle) = norm(z-d.center) ≈ d.radius

arclength(d::Circle) = 2π*d.radius
complexlength(d::Circle) = (d.orientation ? 1 : -1)*im*arclength(d)  #TODO: why?

==(d::Circle,m::Circle) = d.center == m.center && d.radius == m.radius && d.orientation == m.orientation

mappoint(d1::Circle{T},d2::Circle{V},z) where {T<:Vec,V<:Number} =
	mappoint(Circle(complex(d1.center...),d1.radius),d2,z[1]+im*z[2])

mappoint(d1::Circle{T},d2::Circle{V},z) where {T<:Number,V<:Vec} =
	mappoint(Circle(Vec(d1.center...),d1.radius),d2,Vec(real(z),imag(z)))

function mappoint(d1::Circle,d2::Circle,z)
   v=(z-d1.center)/d1.radius
   d1.orientation != d2.orientation && (v=1/v)
   v*d2.radius+d2.center
end

reverseorientation(d::Circle) = Circle(d.center,d.radius,!d.orientation)
conj(d::Circle) = Circle(conj(d.center),d.radius,!d.orientation)

for op in (:+,:-)
    @eval begin
        $op(c::Number,d::Circle) = Circle($op(c,d.center),d.radius,d.orientation)
        $op(d::Circle,c::Number) = Circle($op(d.center,c),d.radius,d.orientation)
    end
end

*(c::Real,d::Circle) = Circle(*(c,d.center),*(abs(c),d.radius),sign(c)<0 ? !d.orientation : d.orientation)
*(d::Circle,c::Real) = Circle(*(c,d.center),*(abs(c),d.radius),sign(c)<0 ? !d.orientation : d.orientation)

/(c::Number,d::Circle) =
	c==1 ? (d.center==0 ? Circle(d.center,1/d.radius,!d.orientation) :
				Circle(1/d.center,abs(1/(d.center+d.radius)-1/(d.center)),!d.orientation)) :
				c*(1/d)


text:
using CPIDataBase, Dates
using CPIDataBase.TestHelpers
using DataFrames
using Test

countrydata = getrandomcountryst()
infl_cpi = InflationTotalCPI()(countrydata)

@testset "InflationCombination" begin

    @testset "Pruebas sin name y sin tag" begin
        combfn = InflationCombination(
            InflationTotalCPI(),
            InflationTotalCPI(),
            Float32[0.5, 0.5]
        )

        @test contains(lowercase(measure_name(combfn)), "promedio")
        @test measure_tag(combfn) == "COMBFN"
        @test all(combfn(countrydata) .≈ infl_cpi)
    end

    @testset "Pruebas con name y no tag" begin
        combfn = InflationCombination(
            InflationTotalCPI(),
            InflationTotalCPI(),
            Float32[0.5, 0.5],
            "Promedio InflationTotalCPI"
            )

        @test contains(measure_name(combfn), "InflationTotalCPI")
        @test measure_tag(combfn) == "COMBFN"
        @test all(combfn(countrydata) .≈ infl_cpi)
    end

    @testset "Pruebas con name y tag" begin
        combfn = InflationCombination(
            InflationTotalCPI(),
            InflationTotalCPI(),
            Float32[0.5, 0.5],
            "Promedio InflationTotalCPI",
            "COMBFNTOTAL"
        )

        @test contains(measure_name(combfn), "InflationTotalCPI")
        @test measure_tag(combfn) == "COMBFNTOTAL"
        @test all(combfn(countrydata) .≈ infl_cpi)
    end

    @testset "Pruebas con ponderadores con suma = 1.5" begin
        combfn = InflationCombination(
            InflationTotalCPI(),
            InflationTotalCPI(),
            Float32[1, 0.5],
            "Promedio InflationTotalCPI mayor que 1",
            "COMBFNTOTAL1"
        )

        # Promedio interanual
        @test all(combfn(countrydata) .≈ 1.5*infl_cpi)
    end

end


text:
import FoldsThreads
let actual = realpath(pathof(FoldsThreads)),
    expected = realpath(joinpath(@__DIR__, "..", "src", "FoldsThreads.jl"))

    if actual != expected
        msg = ("FoldsThreads.jl loaded from an unexpected path. This may be due to" *
            " misconfigured load-path.")
        @warn msg actual expected
        if get(ENV, "CHECK_LOAD_PATH", "false") == "true"
            error(msg)
        end
    end
end

using BenchmarkTools
SUITE = BenchmarkGroup()
for file in readdir(@__DIR__)
    if startswith(file, "bench_") && endswith(file, ".jl")
        SUITE[file[length("bench_")+1:end-length(".jl")]] = include(file)
    end
end


text:
@testset "generic API sanity checks" begin
    @test !is_function_basis("a fish")
end


text:
"""
    ZeroRange(n)

Defines an `AbstractUnitRange` corresponding to the half-open interval
[0,n), equivalent to `0:n-1` but with the lower bound guaranteed to be
zero by Julia's type system.
"""
struct ZeroRange{T<:Integer} <: AbstractUnitRange{T}
    len::T
    ZeroRange{T}(len) where {T} = new{T}(max(zero(T), len))
end
ZeroRange(len::T) where {T<:Integer} = ZeroRange{T}(len)

Base.length(r::ZeroRange) = r.len

Base.length(r::ZeroRange{T}) where {T<:Union{Int,Int64}} = T(r.len)

let smallint = (Int === Int64 ?
                Union{Int8,UInt8,Int16,UInt16,Int32,UInt32} :
                Union{Int8,UInt8,Int16,UInt16})
    Base.length(r::ZeroRange{T}) where {T<:smallint} = Int(r.len)
end

Base.first(r::ZeroRange{T}) where {T} = zero(T)
Base.last(r::ZeroRange{T})  where {T} = r.len-one(T)

function Base.iterate(r::ZeroRange{T}) where T
    r.len <= 0 && return nothing
    x = oftype(one(T) + one(T), zero(T))
    return (x,x)
end
function Base.iterate(r::ZeroRange{T}, i) where T
    x = i + one(T)
    return (x >= oftype(x, r.len) ? nothing : (x,x))
end

@inline function Base.getindex(v::ZeroRange{T}, i::Integer) where T
    @boundscheck ((i > 0) & (i <= length(v))) || Base.throw_boundserror(v, i)
    convert(T, i-1)
end

@inline function Base.getindex(r::ZeroRange{R}, s::AbstractUnitRange{S}) where {R,S<:Integer}
    @boundscheck checkbounds(r, s)
    R(first(s)-1):R(last(s)-1)
end

Base.intersect(r::ZeroRange, s::ZeroRange) = ZeroRange(min(r.len,s.len))

Base.promote_rule(::Type{ZeroRange{T1}},::Type{ZeroRange{T2}}) where {T1,T2} =
    ZeroRange{promote_type(T1,T2)}
Base.convert(::Type{ZeroRange{T}}, r::ZeroRange{T}) where {T<:Real} = r
Base.convert(::Type{ZeroRange{T}}, r::ZeroRange) where {T<:Real} = ZeroRange{T}(r.len)
ZeroRange{T}(r::ZeroRange) where {T} = convert(ZeroRange{T}, r)

Base.show(io::IO, r::ZeroRange) = print(io, typeof(r).name, "(", r.len, ")")


text:
function get_Yᵏ(R,k,n,ℓ,n′,ℓ′,Z,orbital_mode,coulomb_mode,V,poissons,args...)
    u = get_orbitals(R, ℓ, Z, n-ℓ, orbital_mode)[end][2]
    v = get_orbitals(R, ℓ′, Z, n′-ℓ′, orbital_mode)[end][2]
    r = CoulombIntegrals.locs(R)
    if coulomb_mode == :poisson
        Π = if k ∈ keys(poissons)
            poissons[k]
        else
            poissons[k] = PoissonProblem(R, k)
        end
        Π.Y .= false
        ρ = Density(u, v)
        solve!(Π, ρ)
        copy(Π.Y)
    elseif coulomb_mode == :direct
        LCk = LazyCoulomb(R,k)
        M = materialize(applied(*,u',LCk,v))
        w = M .* r
    else
        throw(ArgumentError("Unknown coulomb_mode $(coulomb_mode)"))
    end
end

function get_Fᵏ(R,k,n,ℓ,n′,ℓ′,Z,orbital_mode,coulomb_mode, V, S, r⁻¹,poissons,args...)
    u = get_orbitals(R, ℓ, Z, n-ℓ, orbital_mode)[end][2]
    v = get_orbitals(R, ℓ′, Z, n′-ℓ′, orbital_mode)[end][2]
    r = CoulombIntegrals.locs(R)

    # The F integrals are symmetric, and thus which orbital is used to
    # form the Yᵏ potential should not matter.
    Yᵏ = get_Yᵏ(R,k,n,ℓ,n,ℓ,Z,orbital_mode,coulomb_mode,V,poissons,args...)
    Yᵏ′ = get_Yᵏ(R,k,n′,ℓ′,n′,ℓ′,Z,orbital_mode,coulomb_mode,V,poissons,args...)

    ρv = Density(v, v)
    ρu = Density(u, u)

    # Perform the integral over the other coordinate.
    Fᵏ = dot(ρv.ρ, S, r⁻¹*Yᵏ)
    Fᵏ′ = dot(ρu.ρ, S, r⁻¹*Yᵏ′)

    Fᵏ,Fᵏ′
end

function get_Gᵏ(R,k,n,ℓ,n′,ℓ′,Z,orbital_mode,coulomb_mode, V, S, r⁻¹,poissons,args...)
    u = get_orbitals(R, ℓ, Z, n-ℓ, orbital_mode)[end][2]
    v = get_orbitals(R, ℓ′, Z, n′-ℓ′, orbital_mode)[end][2]
    r = CoulombIntegrals.locs(R)

    Yᵏ = get_Yᵏ(R,k,n,ℓ,n′,ℓ′,Z,orbital_mode,coulomb_mode,V,poissons,args...)

    ρ = Density(u, v)

    # Perform the integral over the other coordinate.
    Gᵏ = dot(ρ.ρ, S, r⁻¹*Yᵏ)
end

function Yᵏ_error(cache, k, n, ℓ, n′, ℓ′, Z, orbital_mode, coulomb_mode, args...)
    @unpack R,ρ,V,poissons = cache
    N = size(R,2)
    # r = CoulombIntegrals.locs(R)
    r = axes(R,1)

    t = time()
    Yᵏ = get_Yᵏ(R,k,n,ℓ,n′,ℓ′,Z,orbital_mode,coulomb_mode,V,poissons,args...)
    el = time()-t
    Ỹᵏ = R \ exact_Yᵏs[((n,ℓ),(n′,ℓ′),k)].(Z*r)

    dYᵏ = Yᵏ - Ỹᵏ
    lδ = maximum(abs, dYᵏ)
    # The way we estimate the global error, by simply summing the
    # magnitudes of the local errors, decreases the order of
    # convergence by 1.
    gδ = sum(abs, dYᵏ)

    a = (n,ℓ)
    b = (n′,ℓ′)

    [N ρ orb_label(a) orb_label(b) k lδ gδ el]
end

function Fᵏ_error(cache, k, n, ℓ, n′, ℓ′, Z, orbital_mode, coulomb_mode, args...)
    @unpack R, ρ, V, S, r⁻¹, poissons = cache
    N = size(R,2)
    r = CoulombIntegrals.locs(R)

    t = time()
    Fᵏ,Fᵏ′ = get_Fᵏ(R,k,n,ℓ,n′,ℓ′,Z,orbital_mode,coulomb_mode, V, S, r⁻¹,poissons,args...)
    el = time()-t

    a = (n,ℓ)
    b = (n′,ℓ′)
    ev = exact_Fᵏs[(a,b,k)]*Z
    evf = convert(Float64, ev)
    [N ρ orb_label(a) orb_label(b) orb_label(a) orb_label(b) k ev evf Fᵏ Fᵏ-evf Fᵏ′ Fᵏ′-evf Fᵏ-Fᵏ′ el]
end

function Gᵏ_error(cache, k, n, ℓ, n′, ℓ′, Z, orbital_mode, coulomb_mode, args...)
    @unpack R, ρ, V, S, r⁻¹, poissons = cache
    N = size(R,2)
    r = CoulombIntegrals.locs(R)

    t = time()
    Gᵏ = get_Gᵏ(R,k,n,ℓ,n′,ℓ′,Z,orbital_mode,coulomb_mode, V, S, r⁻¹,poissons,args...)
    el = time()-t

    a = (n,ℓ)
    b = (n′,ℓ′)
    ev = exact_Gᵏs[(a,b,k)]*Z
    evf = convert(Float64, ev)
    [N ρ orb_label(a) orb_label(b) orb_label(b) orb_label(a) k ev evf Gᵏ Gᵏ-evf el]
end


text:
# This file was generated, do not modify it. # hide
@load OneHotEncoder
@load MultinomialClassifier pkg="MLJLinearModels"


text:
function get_acceptance(acceptance_dist)
    acceptance(states::StatePair) = min(1, acceptance_dist(states))
    return acceptance
end

function update(states::StatePair, alpha::Float64)::State
    u = rand()
    if u <= alpha
        return states.proposal
    else
        return states.current
    end
end

function MH_init(init::Function, x=State[])
    push!(x, init())
    return x
end

function MH_body!(x::Vector{State}, next_state::Function, acceptance::Function, markov_step::Int64)
    for i in 1:markov_step
        states = StatePair(x[i], next_state())
        alpha = acceptance(states)
        push!(x, update(states, alpha))
    end
    return x
end

doc"""
Metropolis-Hastings algorithm
"""
function MH_algorithm(initiate::Function, next_state::Function, acceptance::Function, markov_step::Int64=10000)
    x = MH_init(initiate)
    MH_body!(x, next_state, acceptance, markov_step)
    return x
end


text:
using DiffRules
using SpecialFunctions

mybesselix(ν, x) = besselix.(ν, x)

∇mybesselix(ν, x::Flux.Tracker.TrackedMatrix) = ∇mybesselix(Flux.Tracker.data(ν), Flux.Tracker.data(x))
∇mybesselix(ν, x::Flux.Tracker.TrackedReal) = ∇mybesselix(Flux.Tracker.data(ν), Flux.Tracker.data(x))
∇mybesselix(ν, x::Flux.Tracker.TrackedArray) = ∇mybesselix(Flux.Tracker.data(ν), Flux.Tracker.data(x))
∇mybesselix(ν, x) = @. mybesselix(ν - 1, x) - mybesselix(ν, x) * (ν + x) / x

SpecialFunctions.besselix(ν, x::Flux.Tracker.TrackedMatrix) = Flux.Tracker.track(mybesselix, ν, x)
SpecialFunctions.besselix(ν::Real, x::Flux.Tracker.TrackedReal) = Flux.Tracker.track(mybesselix, ν, x)
SpecialFunctions.besselix(ν, x::Flux.Tracker.TrackedArray) = Flux.Tracker.track(mybesselix, ν, x)
SpecialFunctions.besselix(ν, x::AbstractArray) = besselix.(ν, x)

Flux.Tracker.@grad function mybesselix(ν, x)
    return mybesselix(Flux.Tracker.data(ν), Flux.Tracker.data(x)), Δ -> (nothing, ∇mybesselix(ν, x) .* Δ)
end

# mybesseli(ν, x) = besseli.(ν, x)
#
# ∇mybesseli(ν, x::Flux.Tracker.TrackedMatrix) = ∇mybesseli(Flux.Tracker.data(ν), Flux.Tracker.data(x))
# ∇mybesseli(ν, x::Flux.Tracker.TrackedReal) = ∇mybesseli(Flux.Tracker.data(ν), Flux.Tracker.data(x))
# ∇mybesseli(ν, x::Flux.Tracker.TrackedArray) = ∇mybesseli(Flux.Tracker.data(ν), Flux.Tracker.data(x))
# ∇mybesseli(ν, x) = @. (besseli(ν - 1, x) + besseli(ν + 1, x)) / 2
#
# SpecialFunctions.besseli(ν, x::Flux.Tracker.TrackedMatrix) = Flux.Tracker.track(mybesseli, ν, x)
# SpecialFunctions.besseli(ν::Real, x::Flux.Tracker.TrackedReal) = Flux.Tracker.track(mybesseli, ν, x)
# SpecialFunctions.besseli(ν, x::Flux.Tracker.TrackedArray) = Flux.Tracker.track(mybesseli, ν, x)
# SpecialFunctions.besseli(ν, x::AbstractArray) = besseli.(ν, x)
#
# Flux.Tracker.@grad function mybesseli(ν, x)
#     return mybesseli(Flux.Tracker.data(ν), Flux.Tracker.data(x)), Δ -> (nothing, ∇mybesseli(ν, x) .* Δ)
# end

max_stars_repo_name: honzabim/FewShotAnomalyDetection

max_issues_repo_name: honzabim/FewShotAnomalyDetection

max_forks_repo_name: honzabim/FewShotAnomalyDetection


text:
using Distributed
# Add a process with two threads
addprocs(1; exeflags = "-t 2")
# addprocs(1)

import ConstraintDomains
import CompositionalNetworks
using ConstraintModels
using Constraints
using Dictionaries
@everywhere using LocalSearchSolvers
using Test

# @testset "Distributed" begin
#     @test workers() == [2]
# end

const LS = LocalSearchSolvers

@testset "LocalSearchSolvers.jl" begin
    include("internal.jl")
    include("raw_solver.jl")
end


text:
using Manopt, Documenter, Literate

# generate examples using Literate
tutorialsInputPath = joinpath(@__DIR__, "..", "src/tutorials")
tutorialsRelativePath = "tutorials/"
tutorialsOutputPath = joinpath(@__DIR__, "src/" * tutorialsRelativePath)
tutorials = [
    "MeanAndMedian",
    "StochasticGradientDescent",
    "BezierCurves",
    "GradientOfSecondOrderDifference",
    "JacobiFields",
]
menuEntries = [
    "get Started: Optimize!",
    "do stochastic gradient descent",
    "work with Bézier curves",
    "see the gradient of \$d_2\$",
    "use Jacobi Fields",
]
TutorialMenu = Array{Pair{String,String},1}()
for (i, tutorial) in enumerate(tutorials)
    global TutorialMenu
    sourceFile = joinpath(tutorialsInputPath, tutorial * ".jl")
    targetFile = joinpath(tutorialsOutputPath, tutorial * "md")
    Literate.markdown(sourceFile, tutorialsOutputPath; name=tutorial, credit=false)
    push!(TutorialMenu, menuEntries[i] => joinpath(tutorialsRelativePath, tutorial * ".md"))
end
makedocs(;
    format=Documenter.HTML(; prettyurls=false),
    modules=[Manopt],
    sitename="Manopt.jl",
    pages=[
        "Home" => "index.md",
        "About" => "about.md",
        "How to..." => TutorialMenu,
        "Plans" => "plans/index.md",
        "Solvers" => [
            "Introduction" => "solvers/index.md",
            "Chambolle-Pock" => "solvers/ChambollePock.md",
            "Conjugate gradient descent" => "solvers/conjugate_gradient_descent.md",
            "Cyclic Proximal Point" => "solvers/cyclic_proximal_point.md",
            "Douglas–Rachford" => "solvers/DouglasRachford.md",
            "Gradient Descent" => "solvers/gradient_descent.md",
            "Nelder–Mead" => "solvers/NelderMead.md",
            "Particle Swarm Optimization" => "solvers/particle_swarm.md",
            "Quasi-Newton" => "solvers/quasi_Newton.md",
            "Stochastic Gradient Descent" => "solvers/stochastic_gradient_descent.md",
            "Subgradient method" => "solvers/subgradient.md",
            "Steihaug-Toint TCG Method" =>
                "solvers/truncated_conjugate_gradient_descent.md",
            "Trust-Regions Solver" => "solvers/trust_regions.md",
        ],
        "Functions" => [
            "Introduction" => "functions/index.md",
            "Bézier curves" => "functions/bezier.md",
            "Cost functions" => "functions/costs.md",
            "Differentials" => "functions/differentials.md",
            "Adjoint Differentials" => "functions/adjointdifferentials.md",
            "Gradients" => "functions/gradients.md",
            "Jacobi Fields" => "functions/Jacobi_fields.md",
            "Proximal Maps" => "functions/proximal_maps.md",
            "Specific Manifold Functions" => "functions/manifold.md",
        ],
        "Helpers" => [
            "Data" => "helpers/data.md",
            "Error Measures" => "helpers/errorMeasures.md",
            "Exports" => "helpers/exports.md",
        ],
        "Function Index" => "list.md",
    ],
)
deploydocs(; repo="github.com/JuliaManifolds/Manopt.jl", push_preview=true)


text:
@testset "Document" begin

    sample_text1 = "This is a string"
    sample_text2 = "This is also a string"
    sample_file = joinpath(dirname(@__FILE__), "data", "poem.txt")

    sd = StringDocument(sample_text1)
    fd = FileDocument(sample_file)
    td = TokenDocument(sample_text1)
    ngd = NGramDocument(sample_text1)

    @test isequal(text(sd), sample_text1)
    text!(sd, sample_text2)
    @test isequal(text(sd), sample_text2)
    text!(sd, sample_text1)
    @test isequal(text(sd), sample_text1)

    @test all(tokens(sd) .== ["This", "is", "a", "string"])
    @test "This" in keys(ngrams(sd, 1))
    @test "is" in keys(ngrams(sd, 1))
    @test "a" in keys(ngrams(sd, 1))
    @test "string" in keys(ngrams(sd, 1))

    @test length(sd) == 16

    hamlet_text = "To be or not to be..."
    sd = StringDocument(hamlet_text)
    @test isa(sd, StringDocument)
    @test isequal(text(sd), hamlet_text)

    @test isa(fd, FileDocument)
    @test length(text(fd)) > 0

    my_tokens = ["To", "be", "or", "not", "to", "be..."]
    td = TokenDocument(my_tokens)
    @test isa(td, TokenDocument)
    @test all(tokens(td) .== my_tokens)

    my_ngrams = Dict{String,Int}()
    my_ngrams["To"] = 1
    my_ngrams["be"] = 2
    my_ngrams["or"] = 1
    my_ngrams["not"] = 1
    my_ngrams["to"] = 1
    my_ngrams["be..."] = 1
    ngd = NGramDocument(my_ngrams)
    @test isa(ngd, NGramDocument)
    @test "To" in keys(ngrams(ngd))

    sd = StringDocument(hamlet_text)
    td = TokenDocument(hamlet_text)
    ngd = NGramDocument(hamlet_text)

    d = Document("To be or not to be...")
    @test isa(d, StringDocument)
    d = Document(joinpath(dirname(@__FILE__), "data", "poem.txt"))
    @test isa(d, FileDocument)
    d = Document(["To", "be", "or", "not", "to", "be..."])
    @test isa(d, TokenDocument)
    ng = Dict{String,Int}()
    ng["a"] = 1
    ng["b"] = 3
    d = Document(ng)
    @test isa(d, NGramDocument)

    @test isequal(length(Document("this is text")), 12)
end


text:
using ExtensibleUnions
using Test
using Traceur

abstract type AbstractCar end
abstract type AbstractFireEngine end

struct RedCar <: AbstractCar
end
struct BlueCar <: AbstractCar
    x
end

struct LadderTruck{T} <: AbstractFireEngine
    x::T
end
mutable struct WaterTender{T} <: AbstractFireEngine
    x::T
    y::T
end

struct RedColorTrait end
struct BlueColorTrait end

extensibleunion!(RedColorTrait)
extensibleunion!(BlueColorTrait)

describe(x) = "I don't know anything about this object"

methods(describe)
@test length(methods(describe)) == 1
@test @inferred describe(RedCar()) == "I don't know anything about this object"
@test @inferred describe(BlueCar(1)) == "I don't know anything about this object"
@test @inferred describe(LadderTruck{Int}(2)) == "I don't know anything about this object"
@test @inferred describe(WaterTender{Int}(3,4)) == "I don't know anything about this object"
@check describe(RedCar()) nowarn=[describe] maxdepth=typemax(Int)
@check describe(BlueCar(1)) nowarn=[describe] maxdepth=typemax(Int)
@check describe(LadderTruck{Int}(2)) nowarn=[describe] maxdepth=typemax(Int)
@check describe(WaterTender{Int}(3,4)) nowarn=[describe] maxdepth=typemax(Int)

describe(x::RedColorTrait) = "The color of this object is red"
extensiblefunction!(describe, RedColorTrait)
@test length(methods(describe)) == 2

describe(x::BlueColorTrait) = "The color of this object is blue"
extensiblefunction!(describe, BlueColorTrait)
@test length(methods(describe)) == 3

methods(describe)
@test length(methods(describe)) == 3
@test @inferred describe(RedCar()) == "I don't know anything about this object"
@test @inferred describe(BlueCar(1)) == "I don't know anything about this object"
@test @inferred describe(LadderTruck{Int}(2)) == "I don't know anything about this object"
@test @inferred describe(WaterTender{Int}(3,4)) == "I don't know anything about this object"
@check describe(RedCar()) nowarn=[describe] maxdepth=typemax(Int)
@check describe(BlueCar(1)) nowarn=[describe] maxdepth=typemax(Int)
@check describe(LadderTruck{Int}(2)) nowarn=[describe] maxdepth=typemax(Int)
@check describe(WaterTender{Int}(3,4)) nowarn=[describe] maxdepth=typemax(Int)

addtounion!(RedColorTrait, RedCar)

methods(describe)
@test length(methods(describe)) == 3
@test @inferred describe(RedCar()) == "The color of this object is red"
@test @inferred describe(BlueCar(1)) == "I don't know anything about this object"
@test @inferred describe(LadderTruck{Int}(2)) == "I don't know anything about this object"
@test @inferred describe(WaterTender{Int}(3,4)) == "I don't know anything about this object"
@check describe(RedCar()) nowarn=[describe] maxdepth=typemax(Int)
@check describe(BlueCar(1)) nowarn=[describe] maxdepth=typemax(Int)
@check describe(LadderTruck{Int}(2)) nowarn=[describe] maxdepth=typemax(Int)
@check describe(WaterTender{Int}(3,4)) nowarn=[describe] maxdepth=typemax(Int)

addtounion!(BlueColorTrait, BlueCar)

methods(describe)
@test length(methods(describe)) == 3
@test @inferred describe(RedCar()) == "The color of this object is red"
@test @inferred describe(BlueCar(1)) == "The color of this object is blue"
@test @inferred describe(LadderTruck{Int}(2)) == "I don't know anything about this object"
@test @inferred describe(WaterTender{Int}(3,4)) == "I don't know anything about this object"
@check describe(RedCar()) nowarn=[describe] maxdepth=typemax(Int)
@check describe(BlueCar(1)) nowarn=[describe] maxdepth=typemax(Int)
@check describe(LadderTruck{Int}(2)) nowarn=[describe] maxdepth=typemax(Int)
@check describe(WaterTender{Int}(3,4)) nowarn=[describe] maxdepth=typemax(Int)

addtounion!(RedColorTrait, AbstractFireEngine)

methods(describe)
@test length(methods(describe)) == 3
@test @inferred describe(RedCar()) == "The color of this object is red"
@test @inferred describe(BlueCar(1)) == "The color of this object is blue"
@test @inferred describe(LadderTruck{Int}(2)) == "The color of this object is red"
@test @inferred describe(WaterTender{Int}(3,4)) == "The color of this object is red"
@check describe(RedCar()) nowarn=[describe] maxdepth=typemax(Int)
@check describe(BlueCar(1)) nowarn=[describe] maxdepth=typemax(Int)
@check describe(LadderTruck{Int}(2)) nowarn=[describe] maxdepth=typemax(Int)
@check describe(WaterTender{Int}(3,4)) nowarn=[describe] maxdepth=typemax(Int)


text:
module BasisFunctionsModule
export PrimitiveGaussianBasisFunction, ContractedGaussianBasisFunction, identityPGF, identityCGF
import Base.display
import QuantumLab.BaseModule.evaluateFunction
using ..BaseModule

"""
A PrimitiveGaussianBasisFunction is a function of the form x^a y^b z^c exp(-ζ r^2) centered at a position in space. (a,b,c) here are by definition the cartesian MQuantumNumber tuple.
"""
type PrimitiveGaussianBasisFunction
  center::Position
  exponent::Float64
  mqn::MQuantumNumber
end

"""
A ContractedGaussianBasisFunction is defined as a sum of PrimitiveGaussianBasisFunctions with fixed relative coefficients. Different from other Quantum Chemistry codes we do not require
an additional scaling for the total ContractedGaussianBasisFunction - scaling the whole ContractedGaussianBasisFunction happens by modifying the contraction coefficients accordingly.
"""
type ContractedGaussianBasisFunction
  coefficients::Array{Float64,1}
  primitiveBFs::Array{PrimitiveGaussianBasisFunction,1}
end

function display(cgbf::ContractedGaussianBasisFunction,indent="")
  print(indent); dump(typeof(cgbf))
  for idx in 1:length(cgbf.coefficients)
    print(indent * "  $(cgbf.coefficients[idx]) × "); dump(typeof(cgbf.primitiveBFs[idx]))
    print(indent * "    center:   "); println(cgbf.primitiveBFs[idx].center)
    print(indent * "    exponent: "); println(cgbf.primitiveBFs[idx].exponent)
    print(indent * "    mqn:      "); println(cgbf.primitiveBFs[idx].mqn)
  end
end

function evaluateFunction(x::Position,pgbf::PrimitiveGaussianBasisFunction)
  α = pgbf.exponent
  O = pgbf.center
  r = distance(x,O)
  return (x.x - O.x)^(pgbf.mqn.x) * (x.y - O.y)^(pgbf.mqn.y) * (x.z - O.z)^(pgbf.mqn.z) * exp(-α*r^2)
end

function evaluateFunction(x::Position,cgbf::ContractedGaussianBasisFunction)
  result = 0.::Float64
  for (coeff,pgbf) in zip(cgbf.coefficients,cgbf.primitiveBFs)
    result += coeff * evaluateFunction(x,pgbf)
  end
  return result
end

const identityPGF = PrimitiveGaussianBasisFunction(origin, 0., MQuantumNumber(0,0,0))
const identityCGF = ContractedGaussianBasisFunction([1.],[identityPGF])

end


text:
@testset "599.minimum-index-sum-of-two-lists.jl" begin
    list1 = ["Shogun", "Tapioca Express", "Burger King", "KFC"]
    list2 = ["Piatti", "The Grill at Torrey Pines", "Hungry Hunter Steakhouse", "Shogun"]
    list3 = ["KFC", "Shogun", "Burger King"]
    list4 = ["KFC", "Burger King", "Tapioca Express", "Shogun"]
    list5 = ["KNN", "KFC", "Burger King", "Tapioca Express", "Shogun"]
    @test find_restaurant(list1, list2) == ["Shogun"]
    @test find_restaurant(list1, list3) == ["Shogun"]
    @test find_restaurant(list1, list4) == list4
    @test find_restaurant(list1, list5) ==
        ["KFC", "Burger King", "Tapioca Express", "Shogun"]
    @test find_restaurant(["KFC"], ["KFC"]) == ["KFC"]
end


text:
using SparseArrays

is_int(x::Float64) = x == floor(x)

"""
    extract_implied_integer(sf)

Detect implied integer variables.

Given a problem
```
min     c'x
s.t.    A * x = b
        x ∈ K
```
detects all components of `x` that must be integer.

A variable `xj` is implied integer if it appears in a constraint
```
± x_j0 = b_i - Σ a_i,j x_j
```
with ``x_j``, ``j ≠ j0`` (implied) integer and ``b_i, a_i,j`` integers.
"""
function extract_implied_integer(sf::StandardProblem)
    m, n = size(sf.A)
    A = sf.A
    b = sf.b

    At = sparse(sf.A')

    rhs_integer_flag = is_int.(sf.b)

    nz_per_row = zeros(Int, m)
    n_intvar_per_row = zeros(Int, m)
    integer_row_flags = copy(rhs_integer_flag)

    arows, acols, avals = findnz(sf.A)

    # Initialize number of integer
    for (i, j, v) in zip(arows, acols, avals)
        nz_per_row[i] += 1
        n_intvar_per_row[i] += sf.vartypes_implied[j]
        integer_row_flags[i] &= is_int(v)
    end

    has_changed = true
    while has_changed
        has_changed = false

        # Go through rows
        for i in 1:m
            integer_row_flags[i] || continue  # Skip non-integer rows

            # Go through coefficients
            row = At[:, i]

            j0 = 0
            has_implied_integer = true
            for (j, v) in zip(row.nzind, row.nzval)
                # Check if row writes ± x0 = b - Σ aj xj
                # where b, aj, are integers, and xj, j ≠ j0 are integer variables

                if sf.vartypes_implied[j]
                    # Check if coefficient is integer
                    if !is_int(v)
                        # Not good
                        has_implied_integer = false
                        break
                    end
                else
                    # This variable is fractional
                    if j0 == 0 && (abs(v) == 1.0)
                        j0 = j
                    else
                        # Two fractional variables, or wrong coefficient
                        has_implied_integer = false
                        break
                    end
                end

            end

            has_implied_integer && (j0 != 0) || continue

            # Flag variable as implied integer
            sf.vartypes_implied[j0] = true
            # @info "Variable $(j0) flagged as integer"
            has_changed = true
        end
    end

    @info "Model pre-processing" m n sum(integer_row_flags) sum(sf.vartypes_implied)

    return sf
end

max_stars_repo_name: mtanneau/CLaP

max_issues_repo_name: mtanneau/CLaP

max_forks_repo_name: mtanneau/CLaP


text:
# using DelimitedFiles
# const UCR_DATASETS = """
# 50words,
# Adiac,
# ArrowHead,
# Beef,
# BeetleFly,
# BirdChicken,
# CBF,
# Car,
# ChlorineConcentration,
# CinC_ECG_torso,
# Coffee,
# Computers,
# Cricket_X,
# Cricket_Y,
# Cricket_Z,
# DiatomSizeReduction,
# DistalPhalanxOutlineAgeGroup,
# DistalPhalanxOutlineCorrect,
# DistalPhalanxTW,
# ECG200,
# ECG5000,
# ECGFiveDays,
# Earthquakes,
# ElectricDevices,
# FISH,
# FaceAll,
# FaceFour,
# FacesUCR,
# FordA,
# FordB,
# Gun_Point,
# Ham,
# HandOutlines,
# Haptics,
# Herring,
# InlineSkate,
# InsectWingbeatSound,
# ItalyPowerDemand,
# LargeKitchenAppliances,
# Lighting2,
# Lighting7,
# MALLAT,
# Meat,
# MedicalImages,
# MiddlePhalanxOutlineAgeGroup,
# MiddlePhalanxOutlineCorrect,
# MiddlePhalanxTW,
# MoteStrain,
# NonInvasiveFatalECG_Thorax1,
# NonInvasiveFatalECG_Thorax2,
# OSULeaf,
# OliveOil,
# PhalangesOutlinesCorrect,
# Phoneme,
# Plane,
# ProximalPhalanxOutlineAgeGroup,
# ProximalPhalanxOutlineCorrect,
# ProximalPhalanxTW,
# RefrigerationDevices,
# ScreenType,
# ShapeletSim,
# ShapesAll,
# SmallKitchenAppliances,
# SonyAIBORobotSurface,
# SonyAIBORobotSurfaceII,
# StarLightCurves,
# Strawberry,
# SwedishLeaf,
# Symbols,
# ToeSegmentation1,
# ToeSegmentation2,
# Trace,
# TwoLeadECG,
# Two_Patterns,
# UWaveGestureLibraryAll,
# Wine,
# WordsSynonyms,
# Worms,
# WormsTwoClass,
# synthetic_control,
# uWaveGestureLibrary_X,
# uWaveGestureLibrary_Y,
# uWaveGestureLibrary_Z,
# wafer,
# yoga
# """
#
# """
#     DynamicAxisWarping.download_data([path="/DynamicAxisWarping/data/"])
#
# Downloads the UC Riverside Time Series Classification Archive to
# the specified path.
# """
# function download_ucr()
#
#     UCR_URL = "http://www.cs.ucr.edu/~eamonn/time_series_data/UCR_TS_Archive_2015.zip"
#
#     # the password is "attempttoclassify"
#     @info(
#         """
#         A simple password is required to extract the UCR Time Series Classification. While the
#         data is downloading, please read the following paper
#
#         Bing Hu, Yanping Chen, Eamonn J. Keogh: Time Series Classification under More Realistic
#         Assumptions. SDM 2013: 578-586.
#         http://dx.doi.org/10.1137/1.9781611972832.64
#
#         Find the following sentence (it's near the end of the first page):
#
#         “Every item that we ******* ## @@@@@@@ belongs to exactly one of our well defined classes”
#
#         Enter the three redacted words (without spaces) to extract the data. The purpose of this
#         exercise is to encourage you to read the associated paper, see:
#
#         http://www.cs.ucr.edu/~eamonn/time_series_data/
#
#         """
#     )
#
#     download(UCR_URL, "$DATAPATH/UCR.zip")
#
#     # should unzip the file correctly in platform-specific manner
#     run(unpack_cmd("UCR", DATAPATH , ".zip", ""))
#
#     # delete the zip file when we're done
#     rm("$DATAPATH/UCR.zip")
#
#     @info("Download and extraction successful!")
#
# end
#
# """
#     data,labels = DynamicAxisWarping.traindata(name)
#
# Loads the training set of the specified dataset. Returns a matrix `data` where each column
# holds a 1-dimensional time series. The class label for each column is held `labels`
# which is a vector of length `size(data,2)`.
#
# Available datasets:
#
# $UCR_DATASETS
# """
# function ucr_traindata(name::AbstractString)
#     try
#         Y = readdlm(DATAPATH*"/UCR_TS_Archive_2015/"*name*"/"*name*"_TRAIN", ',')
#         labels = round(Int,vec(Y[:,1]))
#         data = transpose(Y[:,2:end])
#         return data,labels
#     catch err
#         showerror(stdout, err, backtrace());println()
#         @info("You may have recieved this error because you haven't downloaded the database yet.")
#         @info("Try running download_ucr() first.")
#         @info("You may also have mispelled the name of the dataset.")
#     end
# end
#
# """
#     data,labels = DynamicAxisWarping.traindata(name)
#
# Loads the test set of the specified dataset. Returns a matrix `data` where each column
# holds a 1-dimensional time series. The class label for each column is held `labels`
# which is a vector of length `size(data,2)`.
#
# Available datasets:
#
# $UCR_DATASETS
# """
# function ucr_testdata(name::AbstractString)
#     try
#         Y = return readcsv(DATAPATH*"/UCR_TS_Archive_2015/"*name*"/"*name*"_TEST")
#         labels = round(Int,vec(Y[:,1]))
#         data = transpose(Y[:,2:end])
#         return data,labels
#     catch err
#         showerror(stdout, err, backtrace());println()
#         @info("You may have recieved this error because you haven't downloaded the database yet.")
#         @info("Try running download_ucr() first.")
#         @info("You may also have mispelled the name of the dataset.")
#     end
# end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "MozillaCACerts"
# Info and new versions here: https://curl.haxx.se/docs/caextract.html
cacert_version = "2022-02-01"
version = VersionNumber(replace(cacert_version, '-'=>'.'))

# Collection of sources required to build MozillaCACerts
sources = [
    FileSource("https://curl.haxx.se/ca/cacert-$cacert_version.pem",
               "1d9195b76d2ea25c2b5ae9bee52d05075244d78fcd9c58ee0b6fac47d395a5eb",
               filename="cacert.pem"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/
mkdir -p $prefix/share
cp cacert.pem $prefix/share/cacert.pem
install_license LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [AnyPlatform()]

# The products that we will ensure are always built
products = [
    FileProduct("share/cacert.pem", :cacert)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil