text:
# -*- coding: utf-8 -*-
"""
Created on Mon Sep 21 21:00:39 2021

@author: jaewoo.joung
"""

function makeadeck(조커=2)
    if (조커>2 || 조커<0) #혹시 조커를 두장이상 또는 0이하로 썻다면 그냥 두장으로 합니다.
        조커 = 2
    end
    #Tuple suit{무늬(카드뜻/함축의미)}: ♥(성직자/방어),♦(상인/돈),♣(농민/지혜),♠(기사/공격)
    Suit = ('♥','♦','♣','♠')
    DeckRank = ('A',2,3,4,5,6,7,8,9,10,'잭','퀸','킹')
    oneDeck =[]
    for i in 1:length(Suit)
        for j in 1:length(DeckRank)
            push!(oneDeck,(DeckRank[j],Suit[i],j))
        end
    end
    while 조커 > 0
        push!(oneDeck,('조','커',14))
        조커 -= 1
    end
    println("카드를 만듭니다.")
    return oneDeck
end

function shuffle(deck)
    maxnum = length(deck)
    for i in 1:1000
        mtemp = popat!(deck,rand(1:maxnum)) #하나를 골라서 뽑습니다.
        push!(deck,mtemp)                   #고른것을 다시 넣습니다.
    end
    println("카드를 섞습니다.")
    return deck
end

function F_Ace(m_sumVal,m_ACE)
    while m_ACE >0
        if ((m_sumVal >= 7)&&(m_sumVal <= 10))
             m_sumVal += 11
        else
             m_sumVal += 1
        end
        m_ACE -=1
    end
    return m_sumVal
end

function F_Jocker(m_sumVal,m_Jocker)
    while m_Jocker >0
        for i in 1:11
            m_sumVal += 1
            (m_sumVal >= 21) && break
        end
        m_Jocker -=1
    end
    return m_sumVal
end

function F_CalculateBest(m_deck,myDeck)
    m_ACE = 0 ; m_Jocker = 0
    m_size = length(m_deck)
    m_sumVal = 0
    local m_value
    while m_size >0
        m_card = popfirst!(m_deck)
        if (m_card[3]>1 && m_card[3]<11)
            m_sumVal += m_card[3]
        elseif m_card[3] == 1
            m_ACE +=1
        elseif ((m_card[3] > 10) && (m_card[3] < 14))
            m_sumVal += 10
        elseif m_card[3] == 14
            m_Jocker +=1
        end
        m_size -= 1
        push!(m_deck,m_card)
    end
    m_sumVal = F_Ace(m_sumVal,m_ACE)
    m_sumVal = F_Jocker(m_sumVal,m_Jocker)
    return m_sumVal
end

function F_display(myDeck)
    for i in 1:length(myDeck)
        m_card = popfirst!(myDeck)
        print("""( '$(m_card[1])❕':'$(m_card[2])' )""")
        push!(myDeck,m_card)
    end
end

function F_DealersTurn(dealerdeck,myDeck)
    sumofdealer = F_CalculateBest(dealerdeck,myDeck)
    while sumofdealer < 17
        push!(dealerdeck,deal(myDeck))
        sumofdealer = F_CalculateBest(dealerdeck,myDeck)
    end
    println("딜러")
    F_display(dealerdeck)
    sumofdealer > 21 ? (println(""" $sumofdealer 은 꽝(버스트)입니다 \n  제가 졌습니다.""")) : (println("컴퓨터는 $sumofdealer <-- 이숫자입니다. 도전하세요!"))
    return sumofdealer
end

function F_PlayerChoice(m_deck)
    m_size = length(m_deck)
    m_sumVal = 0
    local m_value
    while m_size >0
        m_card = popfirst!(m_deck)
        if (m_card[3]>1 && m_card[3]<11)
            m_sumVal += m_card[3]
        elseif m_card[3] == 1
            @label start1
            print("ACE를 1 로 할까요? 아니면 11로 할까요?: ")
            try
                m_value = parse(Int, readline())
            catch
                println("죄송합니다, 숫자가 아닙니다.")
                @goto start1
            end
            if ((m_value == 1) || (m_value == 11))
                println("감사합니다,",m_value, " 로 정했습니다.")
                m_sumVal +=m_value
            else
                println("죄송합니다, 1또는 11이 아닙니다.")
                @goto start1
            end
        elseif ((m_card[3] > 10) && (m_card[3] < 14))
            m_sumVal += 10
        elseif m_card[3] == 14
            @label start2
            print("조커를 1 부터 11중  무었으로 할까요?: ")
            try
                m_value = parse(Int, readline())
            catch
                println("죄송합니다, 숫자가 아닙니다.")
                @goto start2
            end
            if ((m_value >= 1) && (m_value <= 11))
                println("감사합니다,",m_value, " 로 정했습니다.")
                m_sumVal +=m_value
            else
                println("죄송합니다, 1또는 11이 아닙니다.")
                @goto start2
            end
        end
        m_size -= 1
        push!(m_deck,m_card)
    end
    return m_sumVal
end

function F_PlayerTurn(playerdeck,myDeck)
    Sumofuser = F_PlayerChoice(playerdeck)
    local m_recieve
    while Sumofuser < 21
        println("더 받으시겠습니까? Y/N  ")
        m_recieve = readline()
        if !((m_recieve == "Y") || (m_recieve == "y") || (m_recieve == "N")|| (m_recieve == "n"))
            println("죄송합니다, (원했던)문자가 아닙니다.")
            F_PlayerTurn(playerdeck,myDeck)
        elseif (m_recieve == "Y")|| (m_recieve == "y")
            println("카드를 뽑겠습니다.")
            push!(playerdeck,deal(myDeck))
            Sumofuser = F_PlayerChoice(playerdeck)
            F_display(playerdeck)
        elseif (m_recieve == "N")|| (m_recieve == "n")
            println("여기서 멈추겠습니다.")
            break
        else
            println("죄송합니다, Y또는 N이 아닙니다.")
            F_PlayerTurn(playerdeck,myDeck)
        end
    end
    return  Sumofuser
end

deal = deck -> return pop!(deck) #function

firstdeal(deck) = return [(deal(deck)) for i in 1:2]
########### Main program ###############
@time myDeck = makeadeck()
#println(myDeck)
@time myDeck = shuffle(myDeck)
#println(myDeck)
dealerdeck = firstdeal(myDeck)
playerdeck = firstdeal(myDeck)
Dealer_SUM = F_DealersTurn(dealerdeck,myDeck)
println("나의덱: ")
F_display(playerdeck)

if (Dealer_SUM < 22)
    Sumofuser = F_PlayerTurn(playerdeck,myDeck)
    if (Sumofuser > 21)
        println("당신이 졌습니다. TT","합산하면 $Sumofuser 입니다.")
    elseif (Dealer_SUM == Sumofuser)
        println("우리는 비겼습니다. :)","합산하면 $Sumofuser 입니다.")
    elseif (Dealer_SUM < Sumofuser)
        println("당신이 이겼습니다. :)","합산하면 $Sumofuser 입니다.")
    else
        println("당신이 졌습니다. TT","합산하면 $Sumofuser 입니다.")
    end
else
    println("당신이 이겼습니다. :)","합산하면 $Sumofuser 입니다.")
end

max_stars_repo_name: JaewooJoung/julia

max_issues_repo_name: JaewooJoung/julia

max_forks_repo_name: JaewooJoung/julia


text:
test_integrity_status_flag = false
test_TOW = 34945
test_alert_flag = false
test_anti_spoof_flag = true
test_trans_week = 67
test_codeonl2 = 1
test_ura =  2.0
test_svhealth = "000000"
test_IODC = "0001001000"
test_l2pcode = false
test_T_GD = -1.0710209608078003e-8
test_t_oc =  216000
test_a_f2 =  0.0
test_a_f1 =  -4.774847184307873e-12
test_a_f0 =  -0.00018549291417002678
test_IODE_Sub_2 =  "01001000"
test_C_rs =  70.65625
test_Δn =  3.930878022562108e-9
test_M_0 =  2.4393048719362045
test_C_uc =  3.604218363761902e-6
test_e =  0.01144192845094949
test_C_us =  1.3023614883422852e-5
test_sqrt_A =  5153.7995529174805
test_t_oe =  216000
test_fit_interval =  false
test_AODO = 31
test_C_ic =  -1.73225998878479e-7
test_Ω_0 =  0.0600607702978756
test_C_is =  -2.2351741790771484e-7
test_i_0 =  0.9781895349147778
test_C_rc =  136.34375
test_ω =  0.635978551768012
test_Ω_dot =  -7.383521839035659e-9
test_IODE_Sub_3 =  "01001000"
test_IDOT =  -3.4465721349922174e-10

@testset "Test Decoding true data" begin
    dc = deepcopy(DECODER)
    dc.prev_30 = PREV_30
    dc.prev_29 = PREV_29
    buff = deepcopy(BUFFER_TEST_TRUE)
    GNSSDecoder.decode_words(dc, buff)

    @test dc.data.integrity_status_flag == test_integrity_status_flag #1
    @test dc.data.TOW == test_TOW #2
    @test dc.data.alert_flag == test_alert_flag #3
    @test dc.data.anti_spoof_flag == test_anti_spoof_flag #4
    @test dc.data.trans_week == test_trans_week #5
    @test dc.data.codeonl2 == test_codeonl2 #6
    @test dc.data.ura == test_ura #7
    @test dc.data.svhealth == test_svhealth #8
    @test dc.data.IODC == test_IODC#9
    @test dc.data.l2pcode == test_l2pcode #10
    @test dc.data.T_GD == test_T_GD #11
    @test dc.data.t_oc == test_t_oc #12
    @test dc.data.a_f2 == test_a_f2 #13
    @test dc.data.a_f1 == test_a_f1 #14
    @test dc.data.a_f0 == test_a_f0 #15
    @test dc.data.IODE_Sub_2 == test_IODE_Sub_2 #16
    @test dc.data.C_rs == test_C_rs #17
    @test dc.data.Δn == test_Δn #18
    @test dc.data.M_0 == test_M_0 #19
    @test dc.data.C_uc == test_C_uc #20
    @test dc.data.e == test_e #21
    @test dc.data.C_us == test_C_us #22
    @test dc.data.sqrt_A == test_sqrt_A #23
    @test dc.data.t_oe == test_t_oe #24
    @test dc.data.fit_interval == test_fit_interval #25
    @test dc.data.AODO == test_AODO #26
    @test dc.data.C_ic == test_C_ic #27
    @test dc.data.Ω_0 == test_Ω_0 #28
    @test dc.data.C_is == test_C_is #29
    @test dc.data.i_0 == test_i_0 #30
    @test dc.data.C_rc == test_C_rc #31
    @test dc.data.ω == test_ω #32
    @test dc.data.Ω_dot == test_Ω_dot #33
    @test dc.data.IODE_Sub_3 == test_IODE_Sub_3 #34
    @test dc.data.IDOT == test_IDOT #35
end

@testset "Test parity errors" begin

    PREV_29_FALSE = 1
    PREV_30_FALSE = 1
    dc = deepcopy(DECODER)
    @suppress begin
    dc.prev_30 = PREV_30_FALSE
    dc.prev_29 = PREV_29_FALSE
    buff = deepcopy(BUFFER_TEST_TRUE)
    GNSSDecoder.decode_words(dc, buff)
    end
    @test dc.data_integrity == false

    dc = deepcopy(DECODER)
    @suppress begin
        BAD_LINE = BitArray([1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1])
        dc.prev_30 = PREV_30
        dc.prev_29 = PREV_29
        buff = deepcopy(BUFFER_TEST_TRUE)
        buff[3][3] = BAD_LINE
        GNSSDecoder.decode_words(dc, buff)
    end

    @test dc.data_integrity == false
    @test dc.data.TOW != test_TOW
end

@testset "Test Decoding CEI-cutovers" begin

    IODC_CHANGE_LSB = BitArray([1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1])
    dc = deepcopy(DECODER)
    dc.prev_30 = PREV_30
    dc.prev_29 = PREV_29
    buff = deepcopy(BUFFER_TEST_TRUE)
    buff[1][8] = IODC_CHANGE_LSB
    out = true
    @suppress begin
        TLM_HOW_data,sub_1_data = GNSSDecoder.decode_subframe_1(buff[1])
        TLM_HOW_data,sub_2_data = GNSSDecoder.decode_subframe_2(buff[2])
        TLM_HOW_data,sub_3_data = GNSSDecoder.decode_subframe_3(buff[3])

        out = GNSSDecoder.control_data(TLM_HOW_data, sub_1_data, sub_2_data, sub_3_data)
    end

    @test out == false

    out = true

    IODE_SUB2_CHANGE = BitArray([1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0])
    dc = deepcopy(DECODER)
    dc.prev_30 = PREV_30
    dc.prev_29 = PREV_29
    buff = deepcopy(BUFFER_TEST_TRUE)
    buff[2][3] = IODE_SUB2_CHANGE
    @suppress begin
        TLM_HOW_data,sub_1_data = GNSSDecoder.decode_subframe_1(buff[1])
        TLM_HOW_data,sub_2_data = GNSSDecoder.decode_subframe_2(buff[2])
        TLM_HOW_data,sub_3_data = GNSSDecoder.decode_subframe_3(buff[3])

    out = GNSSDecoder.control_data(TLM_HOW_data, sub_1_data, sub_2_data, sub_3_data)
    end
    @test out == false

    out = true
    IODE_SUB3_CHANGE = BitArray([1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0])
    dc = deepcopy(DECODER)
    dc.prev_30 = PREV_30
    dc.prev_29 = PREV_29
    buff = deepcopy(BUFFER_TEST_TRUE)
    buff[3][10] = IODE_SUB3_CHANGE

    @suppress begin
        TLM_HOW_data,sub_1_data = GNSSDecoder.decode_subframe_1(buff[1])
        TLM_HOW_data,sub_2_data = GNSSDecoder.decode_subframe_2(buff[2])
        TLM_HOW_data,sub_3_data = GNSSDecoder.decode_subframe_3(buff[3])
        out = GNSSDecoder.control_data(TLM_HOW_data, sub_1_data, sub_2_data, sub_3_data)
    end
    @test out == false

end

ext: jl


text:
using OpenCV
using BenchmarkTools
using ProgressMeter

using OpenCV: Mat

# Stuff to add to OpenCV:
Base.materialize!(dest::OpenCV.Mat, src::Broadcast.Broadcasted{Style,Axes,typeof(+),Tuple{<:OpenCV.Mat,<:OpenCV.Mat}}) where {Style,Axes} =
    OpenCV.add(src.args..., dest)

run_complement(img::OpenCV.Mat{T}) where T<:Integer = OpenCV.bitwise_not(img)
# run_complement(img::OpenCV.Mat{T}) where T<:AbstractFloat = 1 .- img

# Not quite equivalent to Julia & Python as it takes a tuple-of-images as input
function run_mean(imgr)
    n = length(imgr)
    n < 2 && return first(imgr)
    img = OpenCV.addWeighted(imgr[1], 1/n, imgr[2], 1/n, 0.0; dtype=Cint(OpenCV.CV_64F))
    for j = 3:n
        OpenCV.addWeighted(img, 1.0, imgr[j], 1/n, 0.0, img, Cint(OpenCV.CV_64F))
    end
    return img
end

run_gradient(img) = (OpenCV.Sobel(img, Cint(OpenCV.CV_64F), Cint(0), Cint(1)), OpenCV.Sobel(img, Cint(OpenCV.CV_64F), Cint(1), Cint(0)))

run_blur(img) = OpenCV.GaussianBlur(img, OpenCV.Size{Int32}(21, 21), 5.0)

function run_histeq(img::OpenCV.Mat)
    if size(img, 1) == 1
        return OpenCV.equalizeHist(img)
    end
    ym = OpenCV.cvtColor(img, OpenCV.COLOR_BGR2LAB)
    vm = OpenCV.split(ym)
    OpenCV.equalizeHist(first(vm), first(vm))
    ym = OpenCV.merge(vm)
    return OpenCV.cvtColor(ym, OpenCV.COLOR_LAB2BGR)
end

const tasknames = ["complement", "mean", "gradient", "blur", "histeq"]
function time_generics(workdir)
    tasks = Dict(task=>Dict{String,Float64}() for task in tasknames)
    @showprogress "Timing generic operations: " for fn in readdir(workdir)
        img = OpenCV.imread(joinpath(workdir, fn), OpenCV.IMREAD_UNCHANGED)
        fn, _ = splitext(fn)
        if isempty(img)
            for tn in tasknames
                tasks[tn][fn] = NaN
            end
            continue
        end
        tasks["complement"][fn] = try
            @belapsed run_complement($img) samples=16 evals=1
        catch
            NaN
        end
        tasks["mean"][fn] = @belapsed run_mean(imgr) samples=16 evals=1 setup=(imgr = ($img, $img, $img, $img, $img))
        tasks["gradient"][fn] = @belapsed run_gradient($img) samples=16 evals=1
        tasks["blur"][fn] = @belapsed run_blur($img) samples=16 evals=1
        tasks["histeq"][fn] = try
            @belapsed run_histeq($img) samples=16 evals=1
        catch
            NaN
        end
    end
    return tasks
end

# tdata = time_generics("/tmp/imgs")

# For things that OpenCV just can't do
return_err(img) = error("not supported")

run_label(img) = OpenCV.connectedComponents(img; connectivity=Cint(8))

run_disttform(img) = OpenCV.distanceTransform(img, OpenCV.DIST_L2, OpenCV.DIST_MASK_3)

function run_flood(img)
    middle = map(Base.tail(axes(img))) do ax
        (first(ax) + last(ax)) ÷ 2
    end
    val = img[:, middle...]
    @show val
    loDiff = eltype(img).(val .- half(img))
    upDiff = white(img) .- val
    masksz = ntuple(d -> d == 1 ? size(img, 1) : size(img, d) + 2, ndims(img))
    mask = falses(masksz)
    newval = (fill(white(img), size(img, 1))...,)
    @show loDiff upDiff newval middle
    OpenCV.floodFill(img, Mat(UInt8.(mask)), OpenCV.Point{Cint}(reverse(middle)...), newval, (loDiff...,), (upDiff...,), Cint(4) | OpenCV.FLOODFILL_MASK_ONLY)
end

white(img::Mat{T}) where T<:Integer = typemax(T)
white(img::Mat{T}) where T<:AbstractFloat = oneunit(T)
half(img::Mat{T}) where T<:Integer = white(img) ÷ 2
half(img::Mat{T}) where T<:AbstractFloat = white(img) / 2

function time_special(workdir)
    imgs = ["components2d" => run_label, "components3d" => return_err,
            "dblcone2d" => run_disttform, "dblcone3d" => return_err,
            "spiral2d" => run_flood, "spiral3d" => return_err]
    fls = readdir(workdir)
    tdata = Dict{String,Float64}()
    @showprogress "Timing special operations: " for (name, f) in imgs
        @show name
        idx = findfirst(str->startswith(str, name), fls)
        img = OpenCV.imread(joinpath(workdir, fls[idx]), OpenCV.IMREAD_UNCHANGED)
        if f === run_label
            img = Mat(UInt8.(img .> 0))
        end
        tdata[name] = try @belapsed $f($img) evals=1 catch NaN end
    end
    return tdata
end

# tdata = time_special("/tmp/imgspecial")


text:
const OVERDUB_ARGS_SYMBOL = gensym("cassette_overdub_arguments")

const BEGIN_OVERDUB_REGION = gensym("cassette_begin_overdub_region")

# Return the CodeInfo method body for signature `S` and `world`,
# if it exists in the method table. Otherwise, return `nothing`.
function lookup_method_body(::Type{S};
                            world::UInt = typemax(UInt),
                            pass::DataType = Unused) where {S<:Tuple}
    @safe_debug "looking up method" signature=S world=world
    S.parameters[1].name.module === Core.Compiler && return nothing

    # retrieve initial Method + CodeInfo
    _methods = Base._methods_by_ftype(S, -1, world)
    length(_methods) == 1 || return nothing
    type_signature, raw_static_params, method = first(_methods)
    method_instance = Core.Compiler.code_for_method(method, type_signature, raw_static_params, world, false)
    method_signature = method.sig
    static_params = Any[raw_static_params...]
    code_info = Core.Compiler.retrieve_code_info(method_instance)
    isa(code_info, CodeInfo) || return nothing
    code_info = Core.Compiler.copy_code_info(code_info)
    @safe_debug "retrieved initial method" method code_info

    # execute user-provided pass if present
    if !(pass <: Unused)
        code_info = pass(S, code_info)
    end

    # substitute static parameters
    body = Expr(:block)
    body.args = code_info.code
    Core.Compiler.substitute!(body, 0, Any[], method_signature, static_params, 0, :propagate)

    # construct new slotnames/slotflags, offset non-self slotnumbers
    code_info.slotnames = Any[code_info.slotnames[1], OVERDUB_ARGS_SYMBOL, code_info.slotnames[2:end]...]
    code_info.slotflags = Any[code_info.slotflags[1], 0x00,                code_info.slotflags[2:end]...]
    replace_match!(x -> isa(x, SlotNumber) || isa(x, NewvarNode), code_info.code) do x
        if isa(x, NewvarNode) && x.slot.id > 1
            return NewvarNode(SlotNumber(x.slot.id + 1))
        elseif x.id > 1
            return SlotNumber(x.id + 1)
        end
        return x
    end

    # construct new `code_info.code` in which original arguments are properly destructured
    new_code = copy_prelude_code(code_info.code)
    prelude_end = length(new_code)
    n_actual_args = fieldcount(S) - 1
    n_method_args = Int64(method.nargs) - 1
    for i in 1:n_method_args
        slot = i + 2
        actual_argument = Expr(:call, GlobalRef(Core, :getfield), SlotNumber(2), i)
        push!(new_code, :($(SlotNumber(slot)) = $actual_argument))
        code_info.slotflags[slot] |= 0x01 << 0x01 # make sure the "assigned" bitflag is set
    end
    if method.isva
        isempty(new_code) || pop!(new_code) # remove the slot reassignment that we're replacing
        final_arguments = Expr(:call, GlobalRef(Core, :tuple))
        for i in n_method_args:n_actual_args
            ssaval = SSAValue(code_info.ssavaluetypes)
            actual_argument = Expr(:call, GlobalRef(Core, :getfield), SlotNumber(2), i)
            push!(new_code, :($ssaval = $actual_argument))
            push!(final_arguments.args, ssaval)
            code_info.ssavaluetypes += 1
        end
        push!(new_code, :($(SlotNumber(n_method_args + 2)) = $final_arguments))
    end
    push!(new_code, BEGIN_OVERDUB_REGION)
    append!(new_code, code_info.code[(prelude_end + 1):end])
    code_info.code = fix_labels_and_gotos!(new_code)
    return code_info
end


text:
#=
   Same rank problem

   From
   Gunnar Blom, Lars Holst, Dennis Sandell:
   "Problems and Snapshots from the World of Probability"
   Page 22.

   100 students are seating (randomly) on seats numbered 0 to 99.
   The lectures says: assume that I rank you according to your age.
   How many will have the same rank number as the seat number?

   Here are three different approaches to this problem which all
   give the same expectation, i.e. E(same_rank) = 1.

   1) Implementing an "all_different" constraint, i.e.
      ensure that it's an permutation.

      n: 6 all_diff: true
      mean(numSameRank: 0.9653

      Distributions of variable numSameRank (num:0)
      1.00000 =>    3942  (0.394200)
      0.00000 =>    3638  (0.363800)
      2.00000 =>    1796  (0.179600)
      3.00000 =>     377  (0.037700)
      4.00000 =>     247  (0.024700)

   2) Generating randomInteger(n) _without_ the all_different
      constraint.

      n: 6 all_diff: false
      mean(numSameRank: 1.0145

      Distributions of variable numSameRank (num:0)
      1.00000 =>    4024  (0.402400)
      0.00000 =>    3312  (0.331200)
      2.00000 =>    1982  (0.198200)
      3.00000 =>     586  (0.058600)
      4.00000 =>      82  (0.008200)
      5.00000 =>      13  (0.001300)
      6.00000 =>       1  (0.000100)

      n: 100 all_diff: false
      mean(numSameRank: 1.0009

   3) Generate n 0/1 numbers with probability of 1 is 1/n.
      Summing this array gives E(same_rank).
      This is the method mentioned in the book, and is implemented
      as same_rank2.

      same_rank2 n: 6
      mean(numSameRank: 1.0064
      Distributions of variable numSameRank (num:0)
      1.00000 =>    4017  (0.401700)
      0.00000 =>    3321  (0.332100)
      2.00000 =>    2035  (0.203500)
      3.00000 =>     539  (0.053900)
      4.00000 =>      80  (0.008000)
      5.00000 =>       8  (0.000800)

      same_rank2 n: 100
      mean(numSameRank: 1.0219

      same_rank2 n: 1000
      mean(numSameRank: 0.9891

     I.e. all these has - also - the same expectation: E(same_rank) = (about) 1.

   Cf ~/webppl/same_rank.wppl

=#
using Turing, StatsPlots
include("jl_utils.jl")

@model function same_rank1(n=6,alldiff=true)
    x ~ filldist(DiscreteUniform(1,n),n)
    if (alldiff)
        # Ensure permutation (see all_different.jl)
        for i in 1:n
            for j in 1:n
                if i != j
                    true ~ Dirac(x[i] != x[j])
                end
            end
        end
    end
    numSameRank ~ Dirac(sum([x[i] == i ? 1 : 0 for i in 1:n]))
end

function run_same_rank1(n, all_diff)
    println("\nsame rank1 n: $n all_diff: $all_diff")
    model = same_rank1(n,all_diff)

    # chns = sample(model, Prior(), 10_000)
    # chns = sample(model, MH(), 10_000)
    # chns = sample(model, PG(5), 10_000)
    chns = sample(model, SMC(), 10_000)
    # chns = sample(model, IS(), 10_000)

    # display(chns)
    # display(plot(chns))
    println("mean(numSameRank: ", mean_val(chns,:numSameRank))

    if n < 10
        show_var_dist_pct(chns,:numSameRank)
    end
end

run_same_rank1(6,true)
run_same_rank1(6,false)
run_same_rank1(100,false)

# Generate n 0/1 values with probability 1/n of getting a 1
@model function same_rank2(n=6)
    x ~ filldist(flip(1/n),n)
    numSameRank ~ Dirac(sum(x))
end

function run_same_rank2(n)
    println("\nsame_rank2 n: $n")
    model = same_rank2(n)

    # chns = sample(model, Prior(), 10_000)
    # chns = sample(model, MH(), 10_000)
    # chns = sample(model, PG(5), 10_000)
    chns = sample(model, SMC(), 10_000)
    # chns = sample(model, IS(), 10_000)

    # display(chns)
    # display(plot(chns))
    println("mean(numSameRank: ", mean_val(chns,:numSameRank))
    if n < 10
        show_var_dist_pct(chns,:numSameRank)
    end
end

run_same_rank2(6)
run_same_rank2(100)
run_same_rank2(1000)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
import StateSpaceReconstruction:
    cembed,
    assign_bin_labels,
    groupslices, groupinds,
    marginal_visitation_freq

import PerronFrobenius:
        InvariantDistribution

"""
    transferentropy_visitfreq(
        E::AbstractEmbedding,
        ϵ::Union{Int, Float64, Vector{Float64}},
        v::TransferEntropy.TEVars) -> Float64

Using the traditional method of estimation probability
distribution by visitation frequencies [1], calculate
transfer entropy from the embedding `E`, given a
discretization scheme controlled by `\epsilon` and
information `v::TEVars`about which columns of the embedding to
consider for each of the marginal distributions. From these
marginal distributions, we calculate marginal entropies and
insert these into the transfer entropy expression.
"""

function transferentropy_visitfreq(
                    E::Embeddings.AbstractEmbedding,
                    ϵ::Union{Int, Float64, Vector{Float64}, Vector{Int}},
                    v::TEVars; b = 2)

    all_inds = unique(vcat(v.target_future,
                    v.target_presentpast,
                    v.source_presentpast,
                    v.conditioned_presentpast))

    visited_bin_inds = assign_bin_labels(E, ϵ)
    slices = groupslices(visited_bin_inds, 2)
    n_visited_bins = length(groupinds(slices))

    # The total number of variables in the embedding
    nvars = size(visited_bin_inds, 1)
    C = v.conditioned_presentpast
    XY = [v.target_future;      v.target_presentpast; C]
    YZ = [v.target_presentpast; v.source_presentpast; C]
    Y =  [v.target_presentpast;                       C]

    p_Y  = marginal_visitation_freq(Y, visited_bin_inds, n_visited_bins)
    p_XY = marginal_visitation_freq(XY, visited_bin_inds, n_visited_bins)
    p_YZ = marginal_visitation_freq(YZ, visited_bin_inds, n_visited_bins)
    P_joint = marginal_visitation_freq(all_inds, visited_bin_inds, n_visited_bins)

    te = StatsBase.entropy(p_YZ, b) +
        StatsBase.entropy(p_XY, b) -
        StatsBase.entropy(p_Y, b) -
        StatsBase.entropy(P_joint, b)

    return te
end

"""
Compute transfer entropy over a range of bin sizes.

Using the traditional method of estimation probability
distribution by visitation frequencies [1], calculate
transfer entropy from the embedding `E`, given a
discretization scheme controlled by the `\epsilon`s and
information `v::TEVars`about which columns of the embedding to
consider for each of the marginal distributions. From these
marginal distributions, we calculate marginal entropies and
insert these into the transfer entropy expression.
"""
function transferentropy_visitfreq(E::Embeddings.AbstractEmbedding,
        ϵ::Vector{Union{Int, Float64, Vector{Float64}, Vector{Int}}},
        v::TEVars)
    map(ϵᵢ -> transferentropy_visitfreq(E, ϵᵢ, v), ϵ)
end

"""
    transferentropy_visitfreq(pts, ϵ, vars::TEVars; b = 2)

Compute transfer entropy for a set of ordered points representing
an appropriate embedding of some time series. See documentation for
`TEVars` for info on how to specify the marginals (i.e. which variables
of the embedding are treated as what).

`b` sets the base of the logarithm (e.g `b = 2` gives the transfer
entropy in bits).
"""
function transferentropy_visitfreq(pts, ϵ, vars::TEVars; b = 2)

    # Collect variables for the marginals
    C = vars.conditioned_presentpast
    XY = [vars.target_future;      vars.target_presentpast; C]
    YZ = [vars.target_presentpast; vars.source_presentpast; C]
    Y =  [vars.target_presentpast;                          C]

    # Find the bins visited by the joint system (and then get
    # the marginal visits from that, so we don't have to encode
    # bins multiple times).
    joint_bin_visits = joint_visits(pts, ϵ)

    # Compute visitation frequencies for nonempty bi
    p_Y = non0hist(marginal_visits(joint_bin_visits, Y))
    p_XY = non0hist(marginal_visits(joint_bin_visits, XY))
    p_YZ = non0hist(marginal_visits(joint_bin_visits, YZ))
    p_joint = non0hist(joint_bin_visits)

    te = StatsBase.entropy(p_YZ, b) +
            StatsBase.entropy(p_XY, b) -
            StatsBase.entropy(p_Y, b) -
            StatsBase.entropy(p_joint, b)
end

tefreq = transferentropy_visitfreq

export
transferentropy_visitfreq,
tefreq


text:
using DataStructures

"""
    PrioritizedSweepingSampleModel <: AbstractSampleModel
    PrioritizedSweepingSampleModel(θ::Float64=1e-4)

See more details at Section (8.4) on Page 168 of the book *Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.*
"""
mutable struct PrioritizedSweepingSampleModel <: AbstractSampleModel
    experiences::Dict{Tuple{Any,Any},Tuple{Float64,Bool,Any}}
    PQueue::PriorityQueue{Tuple{Any,Any},Float64}
    predecessors::Dict{Any,Set{Tuple{Any,Any,Float64,Bool}}}
    θ::Float64
    sample_count::Int
    PrioritizedSweepingSampleModel(θ::Float64=1e-4) = new(Dict{Tuple{Any,Any},Tuple{Float64,Bool,Any}}(),
                                                     PriorityQueue{Tuple{Any,Any},Float64}(Base.Order.Reverse),
                                                     Dict{Any,Set{Tuple{Any,Any,Float64,Bool}}}(),
                                                     θ, 0)
end

function update!(m::PrioritizedSweepingSampleModel, buffer::EpisodeSARDBuffer, learner::AbstractLearner)
    s, a, r, d, s′ = buffer.state[end - 1], buffer.action[end - 1], buffer.reward[end], buffer.isdone[end], buffer.state[end]
    m.experiences[(s, a)] = (r, d, s′)
    P = priority(learner, s, a, r, d, s′)
    if P >= m.θ
        m.PQueue[(s, a)] = P
    end
    if !haskey(m.predecessors, s′)
        m.predecessors[s′] = Set{Tuple{Any,Any,Float64,Bool}}()
    end
    push!(m.predecessors[s′], (s, a, r, d))
end

function sample(m::PrioritizedSweepingSampleModel)
    if length(m.PQueue) > 0
        m.sample_count += 1
        s, a = dequeue!(m.PQueue)
        r, d, s′ = m.experiences[(s, a)]
        s, a, r, d, s′
    else
        nothing
    end
end


text:
function sampleExact(nodepot, edgepot, es::EdgeStruct; Z = -1, nsamples=-1)

  nnodes = size(nodepot, 1)
  nedges = size(edgepot, 3)

  if nsamples < 0
    nsamples = es.maxiter
  end

  samples = zeros(int(nsamples),int(nnodes));

  if Z < 0
    nothing, nothing, nothing, Z = inferExact(nodepot, edgepot, es)
  end

  for s = 1:nsamples
     samples[s,:] = sampleY(nodepot,edgepot,es,Z)
  end
  return samples
end

function sampleY(nodepot, edgepot, es::EdgeStruct, Z)

  nnodes = size(nodepot,1)
  nedges = size(edgepot,3)

  y = ones(1,int(nnodes));
  cumpot = 0
  U = rand()
  i = 1
  while true

      pot = configPotential(y, nodepot, edgepot, es::EdgeStruct)
      cumpot = cumpot + pot;

      if cumpot/Z >= U
          break;
      end

      # next y
      for i = 1:nnodes
          y[i] = y[i] + 1;
          if y[i] <= es.nstates[i]
              break
          else
              y[i] = 1
          end
      end
  end
  return y
end


text:
function velToConductMod(v,mid,a,b)
d = (b-a)./2.0;
dinv = 10;
tt = dinv.*(mid - v);
t = (d.*(tanh.(tt)+1) + a);
dt = -(dinv*d)*(sech.(tt)).^2;
dt = (2.0-v./mid).*dt + (-1./mid).*t;
t = t.*(2.0-v./mid);
return vec(t),spdiagm(vec(dt))
end


text:
"""
    DirectMCEstimate(M::Int)

This is the simplest marginalization strategy, estimating the marginal path
probabilities using a direct Monte Carlo simulation.

`M` specifies the number of samples to be used for the Monte Carlo average.
Larger `M` improves the accuracy of the marginalization integral, but increases
the computational cost.

# Mathematical Description

The marginal probability
```math
\\mathrm{P}[\\bm{x}] = \\int\\mathrm{d}\\bm{s} \\mathrm{P}[\\bm{x}|\\bm{s}] \\mathrm{P}[\\bm{s}]
```
can be computed via a Monte Carlo estimate by sampling `M` trajectories from
``\\mathrm{P}[\\bm{s}]`` and taking the average of the likelihoods:
```math
\\mathrm{P}[\\bm{x}] = \\langle \\mathrm{P}[\\bm{x}|\\bm{s}] \\rangle_{\\mathrm{P}[\\bm{s}]}\\,.
```
"""
struct DirectMCEstimate
    num_samples::Int
end

name(x::DirectMCEstimate) = "Direct MC"

struct DirectMCResult{Samples} <: SimulationResult
    samples::Samples
end

log_marginal(result::DirectMCResult{<:AbstractMatrix}) = vec(logmeanexp(result.samples, dims=2))
log_marginal(result::DirectMCResult{<:AbstractVector}) = logmeanexp(result.samples)

function simulate(algorithm::DirectMCEstimate, initial, system; kwargs...)
    samples = zeros(Float64, algorithm.num_samples)
    for i in 1:algorithm.num_samples
        signal = sample(initial, system)
        samples[i] = -energy_difference(signal, system)
    end
    DirectMCResult(samples)
end


text:
using Documenter, ScenTrees

makedocs(
	sitename =  "ScenTrees.jl",
	authors = "Kipngeno Kirui",
	clean = true,
	doctest = true,
	format = Documenter.HTML(
		assets = ["exampleTree1.png",
		          "Tree402.png",
		          "Tree4022.png",
			  "Tree40221.png",
			  "Tree40222.png",
			  "example1.png",
			  "example21.png",
			  "example22.png",
		          "treeapprox1.png",
		          "treeapprox2D.png",
		          "ExampleLattice2.png",
		          "LatticeApprox.png",
			  "diffHeights.png",
			  "KernLattice.png"
		          ],
		prettyurls = get(ENV, "CI", nothing) == "true"),
	pages = ["Home" => "index.md",
		"Tutorials" => Any["tutorial/tutorial1.md",
				    "tutorial/tutorial2.md",
				    "tutorial/tutorial3.md",
				    "tutorial/tutorial31.md",
				    "tutorial/tutorial4.md",
				    "tutorial/tutorial41.md",
				    "tutorial/tutorial5.md"
				]

		]
)

deploydocs(repo = "github.com/kirui93/ScenTrees.jl.git")


text:
using Test, FinancialMonteCarlo, DualNumbers, BenchmarkTools;
S0 = 100.0;
K = 100.0;
r = 0.02;
T = 1.0;
d = 0.01;
D = 90.0;

underlying_name = "ENJ"
underlying_name_2 = "ABPL"

underlying_ = underlying_name * "_" * underlying_name_2 * "_TESL"

Nsim = 10000;
Nstep = 30;
sigma = 0.2;
mc = MonteCarloConfiguration(Nsim, Nstep);
mc1 = MonteCarloConfiguration(Nsim, Nstep, FinancialMonteCarlo.AntitheticMC());
toll = 0.8

rfCurve = ZeroRate(r);

FwdData = Forward(T)
EUData = EuropeanOptionND(T, K)
Model_enj = BlackScholesProcess(sigma, Underlying(S0, d));
Model_enj_2 = BlackScholesProcess(dual(sigma, 1.0), Underlying(S0, d));
Model_abpl = BlackScholesProcess(sigma, Underlying(S0, d));
Model_tesl = BlackScholesProcess(sigma, Underlying(S0, d));
rho_1 = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0];
Model = GaussianCopulaNVariateProcess(rho_1, Model_enj, Model_abpl, Model_tesl)
ModelDual = GaussianCopulaNVariateProcess(rho_1, Model_enj_2, Model_abpl, Model_tesl)
Model2 = GaussianCopulaNVariateProcess([1.0 0.3 0.0; 0.3 1.0 0.0; 0.0 0.0 1.0], Model_enj_2, Model_abpl, Model_tesl)

# display(Model)

mktdataset = underlying_ → Model
mktdataset_dual = underlying_ → ModelDual
mktdataset_aug = underlying_ → Model2

portfolio_ = [EUData];
portfolio = underlying_ → EUData

@btime price_mkt = pricer(mktdataset, rfCurve, mc, portfolio)
@btime price_mkt = pricer(mktdataset_dual, rfCurve, mc, portfolio)
@btime price_mkt = pricer(mktdataset_aug, rfCurve, mc, portfolio)


text:
function unassigned_test_unit()
    test_is_unassigned_shift_success()
    test_is_unassigned_shift_failure()
    test_meet_base_coverage_sufficient()
    test_meet_base_coverage_generate()
    test_meet_base_coverage_empty_employees()
    test_meet_unassigned_base_coverage()
end

function test_is_unassigned_shift_success()
    # Run it a few times - it generates randomly
    for i in 1:20
        name, employee = generate_unassigned_shift({
            "shift_time_min" => 4,
            "shift_time_max" => 8,
            "coverage" => Array[[2,3,4,],[1,2,3],[1,0,0]],
        })

        @test is_unassigned_shift(name)
    end
end

function test_is_unassigned_shift_failure()
    not_unassigned_shifts = Any[
        "Lenny",
        "Lenny Euler",
        "1",
        1,
        "StaffJoy.com rocks",
        "-1",
    ]

    for s in not_unassigned_shifts
        @test !is_unassigned_shift(s)
    end
end

function test_meet_base_coverage_sufficient()
    # Enough coverage - don't need to generate any more shifts
    # We'll use tests we know pass in functional testing
    new_employees = meet_base_coverage(test_bike_week_1_env, test_bike_week_1_employees)

    @test length(new_employees) == length(test_bike_week_1_employees)
end

function test_meet_base_coverage_generate()
    employees = {
        "Lenny" => {
            "hours_min" => 18,
            "hours_max" => 32,
            "shift_count_max" => 5,
            "availability" => Array[
                [ones(Int, 12)],
                [ones(Int, 12)],
                [ones(Int, 12)],
                [ones(Int, 12)],
                [ones(Int, 12)],
                [zeros(Int, 4), ones(Int, 8)],
                [ones(Int, 12)],
            ],
        },
    }
    new_employees = meet_base_coverage(test_bike_week_1_env, employees)
    @test length(new_employees) > length(employees)
    for e in keys(new_employees)
        if !(e in keys(employees))
            @test is_unassigned_shift(e)
        end
    end
end

function test_meet_base_coverage_empty_employees()
    # Test again with null employees
    employees = Dict()
    new_employees = meet_base_coverage(test_bike_week_1_env, employees)
    @test length(new_employees) > 0
    for e in keys(new_employees)
        @test is_unassigned_shift(e)
    end
end

function test_meet_unassigned_base_coverage()
    env = {
        "coverage" => Array[
            [1, 1, 1, 1, 1],
            [2, 2, 2, 2, 2],
            [3, 3, 3, 3, 3],
        ],
        "shift_time_min" => 1,
        "shift_time_max" => 3,
    }
    day = 2
    expected_max = 10
    expected_min = 4

    employees = meet_unassigned_base_coverage(env, day)
    shift_count = length(employees)
    @test shift_count > expected_min
    @test shift_count < expected_max
    for name in keys(employees)
        @test is_unassigned_shift(name)
    end
end

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
function acceleration(obj::ForceModel, arg0::FieldSpacecraftState, arg1::Vector{RealFieldElement})
    return jcall(obj, "acceleration", FieldVector3D, (FieldSpacecraftState, Vector{RealFieldElement}), arg0, arg1)
end

function acceleration(obj::ForceModel, arg0::SpacecraftState, arg1::Vector{jdouble})
    return jcall(obj, "acceleration", Vector3D, (SpacecraftState, Vector{jdouble}), arg0, arg1)
end

function add_contribution(obj::ForceModel, arg0::FieldSpacecraftState, arg1::FieldTimeDerivativesEquations)
    return jcall(obj, "addContribution", void, (FieldSpacecraftState, FieldTimeDerivativesEquations), arg0, arg1)
end

function add_contribution(obj::ForceModel, arg0::SpacecraftState, arg1::TimeDerivativesEquations)
    return jcall(obj, "addContribution", void, (SpacecraftState, TimeDerivativesEquations), arg0, arg1)
end

function depends_on_position_only(obj::AbstractRadiationForceModel)
    return jcall(obj, "dependsOnPositionOnly", jboolean, ())
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_events_detectors(obj::AbstractRadiationForceModel)
    return jcall(obj, "getEventsDetectors", Stream, ())
end

function get_field_events_detectors(obj::AbstractRadiationForceModel, arg0::Field)
    return jcall(obj, "getFieldEventsDetectors", Stream, (Field,), arg0)
end

function get_parameter_driver(obj::AbstractForceModel, arg0::JString)
    return jcall(obj, "getParameterDriver", ParameterDriver, (JString,), arg0)
end

function get_parameters(obj::ForceModel)
    return jcall(obj, "getParameters", Vector{jdouble}, ())
end

function get_parameters(obj::ForceModel, arg0::Field)
    return jcall(obj, "getParameters", Vector{RealFieldElement}, (Field,), arg0)
end

function get_parameters_drivers(obj::ForceModel)
    return jcall(obj, "getParametersDrivers", Vector{ParameterDriver}, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function init(obj::ForceModel, arg0::SpacecraftState, arg1::AbsoluteDate)
    return jcall(obj, "init", void, (SpacecraftState, AbsoluteDate), arg0, arg1)
end

function is_supported(obj::AbstractForceModel, arg0::JString)
    return jcall(obj, "isSupported", jboolean, (JString,), arg0)
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
# this assumes we have already installed LogicCircuits
using LogicCircuits
# Loads our local package
# We need to "activate" our local version first
import Pkg
# This assumes you are running the file from inside the Project's folder
Pkg.activate(".")
# If this is not the case, you need to give the full path address, e.g.
# Pkg.activate("/Users/denis/ProbabilisticCircuits.jl/")
# Alternatively, we can import our packages from github (but then local changes are not included)
#Pkg.add("https://github.com/denismaua/ProbabilisticCircuits.jl")
using ProbabilisticCircuits

# psdd = learn_probabilistic_circuit(WXData(data)); #using clt
# Loads a PSDD file
pc = zoo_psdd("little_4var.psdd")
# Set some query. 1 assignts True, 0 assigns False, -1 missing
data = XData(Int8.([1 1 -1 -1]))
prob = exp.(log_proba(pc, data))
# println("Prob: $(prob)")

# loads a .vtree file
# vtree_lines = parse_vtree_file("examples/4vars.vtree");
# vtree = compile_vtree_format_lines(vtree_lines);
# println(vtree)

#loads formula and vtree from  the repo https://github.com/UCLA-StarAI/Circuit-Model-Zoo
#
# cnf = zoo_cnf("easy/C17_mince.cnf")
# vtree = zoo_vtree("easy/C17_mince.min.vtree");
# mgr = SddMgr(TrimSddMgr, vtree)
# cnfΔ = node2dag(compile_cnf(mgr, cnf))

#loads a .sdd  file
# sdd = load_logical_circuit("examples/random.sdd")

# Learns CSDD from data based on learn_probabilistic_circuit
data = train(dataset(twenty_datasets("nltcs"); do_shuffle=false, batch_size=-1));
# learn_credal_circuit(data, s_idm) using CSDD/CredalCircuits.jl:estimate_credal_parameters(CredalΔ, data, s_idm)
csdd = learn_credal_circuit(WXData(data), 40.0); #using clt

# println(csdd)

# for node in csdd
#     if node isa Credal⋁  # typeof(node) == Credal⋁{UnstLogicalΔNode}
#         println(node)
#         for j=1:length(node.children)
#             println(typeof(node.children[j]))
#             println(length(node.children[j].children))
#             println(length(node.children[j].children[1]))
#         end
#     end
# end

#         println("lower: ",exp.(node.log_thetas))
#         println("upper: ",exp.(node.log_thetas_u))
#     end
# end

# Testing marginal upper and lower flows

obs = XData(Int8.([1 -1 -1 0 0 -1 -1 1 1 -1 -1 1 -1 -1 -1 -1 ; -1 -1 -1 0 0 -1 -1 1 1 -1 -1 1 -1 -1 -1 -1]))
lower_marg = exp.(log_marginal_lower(csdd, obs))
upper_marg = exp.(log_marginal_upper(csdd, obs))
println("Lower marginal: $(lower_marg)")
println("Upper marginal: $(upper_marg)")

# Testing complete evidence likelihood upper and lower flows

complete_obs = XData(Bool.([1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1; 1 1 1 0 0 1 1 1 1 1 0 0 0 0 0 0]))
lower_prob = exp.(log_prob_lower(csdd, complete_obs))
upper_prob = exp.(log_prob_upper(csdd, complete_obs))
println("Lower prob: $(lower_prob)")
println("Upper prob: $(upper_prob)")

# Testing conditional inference upper and lower bounds
## 1,0 observed/evidence variables values; 2,3 query varables values for 0,1;  -1 marginalize variables.
## cond_queries:    X_1=0,X_2=0|X_3=0,X_4=0
##             X_15=1,X_16=1|X_1=1,X_2=1
cond_queries= XData(Int8.([2 2 0 0 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1; 1 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 3 3; 1 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 3 3 ]))
lower_cond = conditional_lower(csdd, cond_queries, [1.0, 1.0, 1.0])
#upper_cond = conditional_upper(csdd, cond_queries)
println("Lower cond. prob: $(lower_cond)")
#println("Upper cond. prob: $(upper_cond)")


text:
using AODictionary
using Test

@testset "AODictionary.jl" begin
    @test AODict(:a => 2, :b => 4, :c => 6)[:b] == 4
    @test AODict(:a => 2, :b => 4, :c => 6).seq[2] == (:b => 4)
    @test AODict(:a => 2, :b => 4, :c => 6).seqvals[2] == 4
    @test AODict(:a => 2, :b => 4, :c => 6).seqvals[2:end] == [4,6]
    @test AODict(:a => 2, :b => 4, :c => 6).index[:b] == 2
    @test (AODict{Int,Int}()[1]=2) == 2
    @test :a in AODict(:a => 2, :b => 4, :c => 6)
end


text:
module factade

using ForwardDiff
```# use this line and the next few lines to define all external packages required
one after the other
```

``` Use this line below to include the base files as required```
include("extra_file.jl")

export my_f, der_of_my_f

end # module


text:
"""
    pbc(index::Int, N)

Convert the site index to one that is consistent to the periodic boundary condition.

# Examples
```julia-repl
julia> pbc(0, 3)
3

julia> pbc(1, 3)
1

julia> pbc(2, 3)
2

julia> pbc(3, 3)
3

julia> pbc(4, 3)
1
```
"""
function pbc(index::Int, N)
    i = index % N
    i > 0 ? i : i+N
end

"""
    pbc(indeces::Vector{Int}, N)

Convert the site index to one that is consistent to the periodic boundary condition.

# Examples
```julia-repl
julia> pbc([0, 1, 2, 3, 4], 3)
[3, 1, 2, 3, 1]
```
"""
function pbc(indeces::Vector{Int}, N)
    [pbc(i, N) for i in indeces]
end

"""
    pbc(indeces::UnitRange{Int}, N)

Convert the site index to one that is consistent to the periodic boundary condition.

# Examples
```julia-repl
julia> pbc(-7:7, 3)
[2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1]
```
"""
function pbc(indeces::UnitRange{Int}, N)
    [pbc(i, N) for i in indeces]
end


text:
using YaoBase, YaoArrayRegister
import StaticArrays: SMatrix
export RotationGate, Rx, Ry, Rz, rot

"""
    RotationGate{N, T, GT <: AbstractBlock{N, Complex{T}}} <: PrimitiveBlock{N, Complex{T}}

RotationGate, with GT both hermitian and isreflexive.

# Definition

```math
\\mathbf{I} cos(θ / 2) - im sin(θ / 2) * mat(U)
```
"""
mutable struct RotationGate{N,T,GT<:AbstractBlock{N}} <: PrimitiveBlock{N}
    block::GT
    theta::T
    function RotationGate{N,T,GT}(block::GT, theta) where {N,T,GT<:AbstractBlock{N}}
        ishermitian(block) && isreflexive(block) ||
            throw(ArgumentError("Gate type $GT is not hermitian or not isreflexive."))
        new{N,T,GT}(block, T(theta))
    end
end

RotationGate(block::GT, theta::T) where {N,T,GT<:AbstractBlock{N}} =
    RotationGate{N,T,GT}(block, theta)
RotationGate(block::GT, theta::Integer) where {N,GT<:AbstractBlock{N}} =
    RotationGate{N,Float64,GT}(block, Float64(theta))

# bindings
"""
    Rx(theta)

Return a [`RotationGate`](@ref) on X axis.

# Example

```jldoctest; setup=:(using YaoBlocks)
julia> Rx(0.1)
rot(X, 0.1)
```
"""
Rx(theta) = RotationGate(X, theta)

"""
    Ry(theta)

Return a [`RotationGate`](@ref) on Y axis.

# Example

```jldoctest; setup=:(using YaoBlocks)
julia> Ry(0.1)
rot(Y, 0.1)
```
"""
Ry(theta) = RotationGate(Y, theta)

"""
    Rz(theta)

Return a [`RotationGate`](@ref) on Z axis.

# Example

```jldoctest; setup=:(using YaoBlocks)
julia> Rz(0.1)
rot(Z, 0.1)
```
"""
Rz(theta) = RotationGate(Z, theta)

"""
    rot(U, theta)

Return a [`RotationGate`](@ref) on U axis.
"""
rot(axis::AbstractBlock, theta) = RotationGate(axis, theta)

content(x::RotationGate) = x.block
chcontent(x::RotationGate, gen::AbstractBlock) = RotationGate(gen, x.theta)
# General definition
function mat(::Type{T}, R::RotationGate{N}) where {N,T}
    I = IMatrix{1 << N,T}()
    return I * cos(T(R.theta) / 2) - im * sin(T(R.theta) / 2) * mat(T, R.block)
end

# Specialized
mat(::Type{T}, R::RotationGate{1,<:Any,<:XGate}) where {T} =
    T[cos(R.theta / 2) -im * sin(R.theta / 2); -im * sin(R.theta / 2) cos(R.theta / 2)]
mat(::Type{T}, R::RotationGate{1,<:Any,<:YGate}) where {T} =
    T[cos(R.theta / 2) -sin(R.theta / 2); sin(R.theta / 2) cos(R.theta / 2)]
# mat(R::RotationGate{1, T, ZGate{Complex{T}}}) where T =
#     SMatrix{2, 2, Complex{T}}(cos(R.theta/2)-im*sin(R.theta/2), 0, 0, cos(R.theta/2)+im*sin(R.theta/2))

function apply!(r::ArrayReg, rb::RotationGate)
    v0 = copy(r.state)
    apply!(r, rb.block)
    # NOTE: we should not change register's memory address,
    # or batch operations may fail
    r.state .= -im * sin(rb.theta / 2) * r.state + cos(rb.theta / 2) * v0
    return r
end

# parametric interface
niparams(::Type{<:RotationGate}) = 1
getiparams(x::RotationGate) = x.theta
# no need to specify the type of param, Julia will try to do the conversion
setiparams!(r::RotationGate, param::Number) where {N,T} = (r.theta = param; r)

# fallback to matrix methods if it is not real
YaoBase.isunitary(r::RotationGate{N,<:Real}) where {N} = true

function YaoBase.isunitary(r::RotationGate)
    isreal(r.theta) && return true
    @warn "θ in RotationGate is not real, got θ=$(r.theta), fallback to matrix-based method"
    return isunitary(mat(r))
end

Base.adjoint(blk::RotationGate) = RotationGate(blk.block, -blk.theta)
Base.copy(R::RotationGate) = RotationGate(R.block, R.theta)
Base.:(==)(lhs::RotationGate{TA,GTA}, rhs::RotationGate{TB,GTB}) where {TA,TB,GTA,GTB} = false
Base.:(==)(lhs::RotationGate{TA,GT}, rhs::RotationGate{TB,GT}) where {TA,TB,GT} =
    lhs.theta == rhs.theta

cache_key(R::RotationGate) = R.theta


text:
using Base.Test
# spec.jl
import WaveletScattering: tune_motherfrequency, default_motherfrequency
# meta.jl
import WaveletScattering: get_γ
# spec1d.jl
import WaveletScattering: Spec1D

# tune_motherfrequency
nfos = [1, 2, 4, 8, 12, 16, 24]
pitchforks = [392, 415, 422, 430, 435, 440, 442, 444, 466]
for nfo in nfos, pitchfork in pitchforks
    tuningfrequency = pitchfork / 44100.0
    spec = Spec1D()
    ξ = tune_motherfrequency(tuningfrequency, spec.class, nfo)
    γs = map(get_γ, spec.ψmetas)
    ωs = ξ * exp2.(-γs / nfo)
    @test any(abs.(ωs - ξ) .< 1e-4)
    max_ξ = default_motherfrequency(spec.class, nfo)
    @test ξ < max_ξ
    @test ξ * exp2(inv(nfo)) > max_ξ
end


text:
using DataStructures

#mutable struct Node
#    ip::Int
#    neighbours::Array{Tuple{Node,Int}}
#    risk::Union{Float64, Nothing}
#    predecessor::Union{Node, Nothing}
#    probability::Float64
#    Node() = (obj = new(0);)
#end

function create_network()
    nodes = []
    for i = 1:4
        push!(nodes, Node() )
        nodes[i].ip = i
    end

    nodes[1].probability = 1.0
    nodes[1].risk = 0
    nodes[1].neighbours = [(nodes[2], 6), (nodes[3], 3)]

    nodes[2].probability = 0.5
    nodes[2].risk = 12
    nodes[2].neighbours = [(nodes[4], 3)]

    nodes[3].probability = 0.75
    nodes[3].risk = 4
    nodes[3].neighbours = [(nodes[4], 4)]

    nodes[4].probability = 0.8
    nodes[4].risk = 9

    return nodes
end

function initialize_single_source!(graph, start)
    for v in graph
        v.risk = typemax(Float64)
        v.predecessor = nothing
    end
    start.risk = 0
end

function relax!(from_node,to_node,cost)
    tot_risk = from_node.risk + (cost/to_node.probability)
    if (tot_risk < to_node.risk)
        to_node.risk = tot_risk
        to_node.predecessor = from_node
    end
end

function extract_min(list)
    min = typemax(Float64)
    ret_node = list[1]
    delete_idx = 1
    for i = 1:length(list)
        if list[i].risk < min
            ret_node = list[i]
            delete_idx = i
        end
    end
    deleteat!(list, delete_idx)
    return ret_node
end

function dijkstra!(graph,start)
    initialize_single_source!(graph, start)
    pq = PriorityQueue()                  # create a priority que
    enqueue!(pq, start, start.risk)
    while !isempty(pq)
        from_node = dequeue!(pq)
        for to_node in from_node.neighbours
            relax!(from_node, to_node[1], to_node[2])
            for node in pq # check if node already exists in que
                if node[1] == to_node[1] # if so change it
                    pq[to_node[1]] = to_node[1].risk # priority changed to new risk
                end
                println(node[to_node[1]])
            end
            enqueue!(pq, to_node[1], to_node[1].risk)
        end
    end
end

test = PriorityQueue()
network = create_network()
enqueue!(test, network[1].neighbours[1], network[1].neighbours[1][1].risk)
enqueue!(test, network[1].neighbours[2], network[1].neighbours[2][1].risk)

for i in test
    println(i[1])
    if i[1] == network[1].neighbours[2]
        #println(i[1])
    end
end

max_stars_repo_name: evenlwanvik-student/algdat

max_issues_repo_name: evenlwanvik-student/algdat

max_forks_repo_name: evenlwanvik-student/algdat


text:
#
# The functions in this file look at the status of the iteration at the
# end and set idid and errcode.
#
# Nothing exciting in here, but it must be done.
#
"""
NewtonOK: Figure out idid and errcode for Newton's method
"""
function NewtonOK(resnorm, iline, tol, toosoon, itc, ItRules)
    maxit = ItRules.maxit
    armmax = ItRules.armmax
    printerr = ItRules.printerr
    resfail = (resnorm > tol)
    idid = ~(resfail || toosoon)
    errcode = 0
    if ~idid
        errcode =
            NewtonError(resfail, iline, resnorm, toosoon, tol, itc, maxit, armmax, printerr)
    end
    return (idid, errcode)
end

"""
PTCOK: Figure out idid and errcode
"""
function PTCOK(resnorm, tol, toosoon, ItRules, printerr)
    maxit = ItRules.maxit
    delta0 = ItRules.delta0
    errcode = 0
    resfail = (resnorm > tol)
    idid = ~(resfail || toosoon)
    errcode = 0
    if ~idid
        (errcode = PTCError(resnorm, maxit, delta0, toosoon, tol, printerr))
    end
    return (idid, errcode)
end


text:
################################################################################
## Loading packages that are used.                                            ##
################################################################################
using CSV, DataFrames, DataFramesMeta
using Gadfly
using ZipFile
#################################################################################

#################################################################################
## Loading data from CDPR FTP server.                                          ##
#################################################################################
"""
    choose_years(year1[, year2])

Use this function to start analysis, by choosing the years you want to use for your study.
This will create the proper filepaths.
"""
function choose_years(year1::Int, year2=Union{Nothing,Int})
    if year2 != nothing
        first_year = minimum([year1, year2])
        second_year = maximum([year1, year2])
        years = collect(UnitRange(first_year, second_year))
    else
        years = [year1]
    end
end
#################################################################################

#################################################################################
## Unzipping and loading files into memory                                      ##
#################################################################################
function uncompress_file(year::Int)
    path = "Pesticide_Use/pur" * string(year) *".zip"
    io = open(path, "r")
    io = IOBuffer(read(io))
    fl = ZipFile.Reader(io)
    return fl
end
"""
    get_files(:kind, year)

Unizips and loads the desires files from the zip file from the year and save it in
a Dataframe.

kind = :chemicals loads the file "chemical.txt", which contains a list with
all approved chemicals for each year and their chem_codes.
kind = :sites loads the file "site.txt", which contains a list with all sites
(commodities) that can be reported in CA for that year and their site_codes.
kind = :counties loads the the file "county.txt", which contains a list with
all counties in CA and their county_cd (codes).
kind = :products loads the files "product.txt", which contains a list with all approved
chemical products in CA for any year. The resulting dataframe contains
"""
function get_files(kind::Symbol, year::Int)
    reading = uncompress_file(year)
    # the zile file for 2016 has and extra "pur2016/" on the file name that has to be removed
    if year == 2016
        for i in 1:length(reading.files)
            (Y, N) = split(reading.files[i].name, "/")
            reading.files[i].name = N
        end
    end

    for i in 1:length(reading.files)
        if kind == :chemicals && reading.files[i].name == "chemical.txt"
            fl = CSV.read(IOBuffer(read(reading.files[i])))
            return fl
        elseif kind == :sites && reading.files[i].name == "site.txt"
            fl = CSV.read(IOBuffer(read(reading.files[i])))
            return name
        elseif kind == :counties && reading.files[i].name == "county.txt"
            fl = CSV.read(IOBuffer(read(reading.files[i])))
            return fl
        elseif kind == :products && reading.files[i].name == "product.txt"
            fl = CSV.read(IOBuffer(read(reading.files[i])))
            return fl
        end
    end
    throw("Wrong symbol for kind")
end

# experimental function. need to time with BenchmarkTools
function get_files1(kind::Symbol, year::Int)
    reading = uncompress_file(year)
    # the zile file for 2016 has and extra "pur2016/" on the file name that has to be removed
    if year == 2016
        for i in 1:length(reading.files)
            (Y, N) = split(reading.files[i].name, "/")
            reading.files[i].name = N
        end
    end

    if kind == :chemicals
        for i in 1:length(reading.files)
            if reading.files[i].name == "chemical.txt"
                name = CSV.read(IOBuffer(read(reading.files[i])))
                return name
            end
        end
    elseif kind == :sites
        for i in 1:length(reading.files)
            if reading.files[i].name == "site.txt"
                name = CSV.read(IOBuffer(read(reading.files[i])))
                return name
            end
        end
    elseif kind == :counties
        for i in 1:length(reading.files)
            if reading.files[i].name == "county.txt"
                name = CSV.read(IOBuffer(read(reading.files[i])))
                return name
            end
        end
    elseif kind == :products
        for i in 1:length(reading.files)
            if reading.files[i].name == "product.txt"
                name = CSV.read(IOBuffer(read(reading.files[i])))
                return name
            end
        end
    end
    throw("Wrong symbol for kind")
end

# Function to extract the differences between the approved chemicals between years
function compare_info_years(kind::Symbol, year1::Int, year2::Int)
    if kind == :chemicals
        join(get_files(maximum([year1, year2]), :chemicals),
            get_files(minimum([year1, year2]), :chemicals),
            on =:chem_code, kind = :anti)
    elseif kind == :sites
        join(get_files(maximum([year1, year2]), :sites),
             get_files(minimum([year1, year2]), :sites),
             on =:site_code, kind = :anti)
    elseif kind == :products
        join(get_files(maximum([year1, year2]), :products),
             get_files(minimum([year1, year2]), :products),
             on =:prodno, kind = :anti)
    else throw("Wrong symbol for kind")
    end
end

# find sites using any word (site) as a string.
function find_product_code(product::AbstractString, year::Int)
    products = get_files(:products, year)
    filtered = DataFrame(prodno=[], product_name=[])
    for i in 1:size(products,1)
        matching = match(Regex(uppercase(product)), products[i,8])
        if matching != nothing
            push!(filtered, (prodno = products[i,1], product_name=products[i,8]))
        end
    end
    if nrow(filtered) != 0
        return filtered
    else throw("Non-existent product")
    end
end

# find sites using any word (site) as a string.
function find_site_code(site::AbstractString, year::Int)
    sites = get_files(:sites, year)
    filtered = DataFrame(site_code=[], site_name=[])
    for i in 1:size(sites,1)
        matching = match(Regex(uppercase(site)), sites[i,2])
        if matching != nothing
            push!(filtered, (site_code = sites[i,1], site_name=sites[i,2]))
        end
    end
    if nrow(filtered) != 0
        return filtered
    else throw("Non-existent commodity")
    end
end

function find_chem_code(chemical::AbstractString, year::Int)
    chemicals = get_files(:chemicals, year)
    filtered = DataFrame(chem_code=[], chem_name=[])
    for i in 1:size(chemicals,1)
        matching = match(Regex(uppercase(chemical)), chemicals[i,3])
        if matching != nothing
            push!(filtered, (chem_code = chemicals[i,1], chem_name=chemicals[i,3]))
        end
    end
    if nrow(filtered) != 0
        return filtered
    else throw("Non-existent chemical")
    end
end

# Extract the county code
function find_county_code(county_name::AbstractString, year::Int)
    counties = get_files(year, :counties)
    for i in 1:length(counties[:couty_name])
        if counties[i,2] == uppercase(county_name)
            county_code = counties[i,1]
            return county_code
        end
    end
    throw("Wrong county name")
end

################################################################################
## Dealing with the PUR.txt file                                              ##
################################################################################
# Function to load pur file from disk, using the county name as a string and the year
function load_pur(county_name::AbstractString, year::Int, dropmissing::Bool = true)
    county_code = find_county_code(county_name, year)
    string_year = string(year)

    reading = uncompress_file(year)
    # the zip file for 2016 has and extra "pur2016/" on the file name that has to be removed
    if year == 2016
        for i in 1:length(reading.files)
            (Y, N) = split(reading.files[i].name, "/")
            reading.files[i].name = N
        end
    end

    for i in 1:length(reading.files)
        if reading.files[i].name == "udc" * string_year[3:4] * "_"* string(county_code) * ".txt"
            pur = CSV.File(IOBuffer(read(reading.files[i])), escapechar='#',
            types=Dict(18 => Union{Int64, Char, Missing},20 => Union{Int64, Char, Missing}))
            pur_df = DataFrame(pur)
        end
    end

    pur_df = @select(pur_df, :use_no, :prodno, :chem_code, :prodchem_pct, :lbs_chm_used,
    :lbs_prd_used, :amt_prd_used, :unit_of_meas, :acre_planted, :unit_planted,
    :acre_treated, :unit_treated, :applic_cnt, :applic_dt, :grower_id, :site_code)

    # add column for year
    pur_df[:year] = repeat([year], outer = size(pur_df, 1))

    # add chem name
    pur_df[:chem_name] = repeat(["0"],outer = size(pur_df, 1))

    if dropmissing == true
        pur_df = dropmissing!(pur_df)
    end

    chem_codes = get_files(year, :chemicals)
    chemcodes = Dict()
    for i in 1:size(chem_codes,1)
        push!(chemcodes, (chem_codes[i,1] => chem_codes[i,3]))
    end

    for i in 1:size(pur_df,1)
        chemical = pur_df[i,:chem_code]
        if typeof(chemical) != Missing
            pur_df[i,:chem_name] = chemcodes[chemical]
        end
    end
    return pur_df
end

# Select a site code using the find_site_code function
function select_site_code(df::DataFrame, site_code::Int)
    @where(df, :site_code .== site_code)
end

function select_insecticide(df::DataFrame, insecticide_name::Union{AbstractString, Int})
    if typeof(insecticide_name)== AbstractString
        insecticide_name = string(uppercase(insecticide_name))
        @where(df, :chem_name .== insecticide_name)
    elseif typeof(insecticide_name)==Int
        df = dropmissing!(df)
        @where(df, :chem_code .== insecticide_name)
    else throw("Wrong insecticide name or code")
    end
end

# Function to append the data from different years
function append_years(county_name::AbstractString, first_year::Int, last_year::Union{Nothing,Int},
                      dropmissing::Bool = true)
    if last_year==nothing
        total = load_pur(county_name, first_year, false)
        return total
    end

    if dropmissing == false
        total = load_pur(county_name, first_year, false)
        for i in first_year+1:last_year
            pur = load_pur(county_name, i, false)
            total= append!(total, pur)
        end
        return total
    end

    total = load_pur(county_name, first_year)
    for i in first_year+1:last_year
        pur = load_pur(county_name, i)
        total= append!(total, pur)
    end
    return total
end

# Function to obtain the information in insecticide use from a given commodity and county
function insecticide_use(county_name::AbstractString,
                         insecticide_name::Union{AbstractString, Int},
                         site_code::Int,
                         first_year::Int, last_year::Union{Nothing,Int},
                         dropmissing::Bool=true)
    total = append_years(county_name, first_year, last_year, dropmissing)
    total = select_site_code(total, site_code)
    total = select_insecticide(total, insecticide_name)
    total = by(total, [:year]) do df
    DataFrame(lbs_chm_used = sum(df[:lbs_chm_used]),
              N = size(df,1))
    end
    return total
end

function plot_insecticide_use(county_name::AbstractString,
                         insecticide_name::AbstractString,
                         site_code::Int,
                         first_year::Int, last_year::Union{Nothing,Int},
                         dropmissing::Bool=true)
    result = insecticide_use(county_name, insecticide_name, site_code, first_year, last_year,
                             dropmissing)
    plot(result, x = :year, y = :lbs_chm_used, Geom.line)
end

plot(abamectin_monterey_154, x = :year, y = :lbs_chm_used, Geom.line, Theme(background_color="white"))


text:
using Random
using JuMP, BilevelJuMP

function bench_forecast(prods, samples, optimizer, mode, seed = 1234)

    rng = Random.MersenneTwister(seed)

    SampleSize = samples
    Products = prods

    Samples = 2:SampleSize

    phi1 = 0.7 * ones(Products)
    phi0 = 3 * ones(Products)

    mu = 3.0/(1.0 - 0.7)

    demand = ones(SampleSize, Products)

    for p in 1:Products
        for t in Samples
            demand[t, p] = max(0.0, phi0[p] + phi1[p] * demand[t-1, p] + 0.2 * randn(rng))
        end
    end

    p_s = 5 .+ 0.1*collect(1:Products)
    p_b = 4 .+ 0.1*collect(1:Products)
    p_r = 3 .+ 0.1*collect(1:Products)

    u = 8*Products

    ctrs = []
    vars = []

    # MOI.empty!(optimizer)
    model = BilevelModel(optimizer, mode = mode)
    try
        JuMP.set_time_limit_sec(model, MAX_TIME)
    catch e
        @show e
        @show "failed to set limit time"
    end
    # parameters
    v = @variable(Upper(model),  0 <= φ0[1:Products] <= 10)
    push!(vars, vec(v))
    v = @variable(Upper(model), -1 <= φ1[1:Products] <= 1)
    push!(vars, vec(v))

    #=
        lower model
    =#

    # buying quantity is the only
    v = @variable(Lower(model), 0 <= q_b[p=1:Products, t=Samples] <= u)
    push!(vars, vec(v))

    #
    v = @variable(Lower(model), 0 <= q_s[p=1:Products, t=Samples] <= u)
    push!(vars, vec(v))
    v = @variable(Lower(model), 0 <= q_r[p=1:Products, t=Samples] <= u)
    push!(vars, vec(v))

    c = @constraint(Lower(model),[p=1:Products,t=Samples],
        q_s[p,t] + q_r[p,t] <= q_b[p,t])
    push!(ctrs, vec(c))

    c = @constraint(Lower(model),[p=1:Products,t=Samples],
        q_s[p,t] <= φ0[p] + φ1[p] * demand[t-1, p])
    push!(ctrs, vec(c))

    c = @constraint(Lower(model),[t=Samples],
        sum(q_b[p,t] for p in 1:Products) <= u)
    push!(ctrs, vec(c))

    @objective(Lower(model), Min,
        sum(p_b[p] * q_b[p,t] - p_s[p] * q_s[p,t] - p_r[p] * q_r[p,t]
            for t in Samples, p in 1:Products))

    #=
        upper model
    =#

    #
    v = @variable(Upper(model), 0 <= q2_s[p=1:Products, t=Samples] <= u)
    push!(vars, vec(v))
    v = @variable(Upper(model), 0 <= q2_r[p=1:Products, t=Samples] <= u)
    push!(vars, vec(v))

    c = @constraint(Upper(model),[p=1:Products,t=Samples],
        q2_s[p,t] + q2_r[p,t] <= q_b[p,t])
    push!(ctrs, vec(c))

    c = @constraint(Upper(model),[p=1:Products,t=Samples],
        q2_s[p,t] <= demand[t, p])
    push!(ctrs, vec(c))

    @objective(Upper(model), Min,
        sum(p_b[p] * q_b[p,t] - p_s[p] * q2_s[p,t] - p_r[p] * q2_r[p,t]
            for t in Samples, p in 1:Products))

    #=
        Optimize
    =#

    optimize!(model)

    #=
    for v in vars
        val = value(v) # >= 0
    end

    for c in ctrs
        val = value(c) - normalized_rhs(c) # >= 0
    end
    =#

    @show primal_st = primal_status(model)
    @show term_st = termination_status(model) #in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED, MOI.ALMOST_LOCALLY_SOLVED]

    solve_t = solve_time(model)
    build_t = BilevelJuMP.build_time(model)

    obj_l = try
        objective_value(Lower(model))
    catch
        NaN
    end
    obj_u = try
        objective_value(Upper(model))
    catch
        NaN
    end
    gap = try
        bound = objective_bound(Upper(model))
        abs(obj_u - bound)/max(abs(bound), 1e-8)
    catch
        NaN
    end

    return primal_st, term_st, solve_t, build_t, obj_l, obj_u, gap
end


text:
module VortexHelperPufferBarrier
using ..Ahorn, Maple

@mapdef Entity "VortexHelper/PufferBarrier" PufferBarrier(x::Integer, y::Integer, width::Integer = 8, height::Integer = 8)

const placements = Ahorn.PlacementDict(
    "Puffer Barrier (Vortex Helper)" => Ahorn.EntityPlacement(
        PufferBarrier,
        "rectangle"
    ),
)

Ahorn.minimumSize(entity::PufferBarrier) = 8, 8
Ahorn.resizable(entity::PufferBarrier) = true, true

const color = (255, 189, 74, 180) ./ 255

function Ahorn.selection(entity::PufferBarrier)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return Ahorn.Rectangle(x, y, width, height)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::PufferBarrier, room::Maple.Room)
    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    Ahorn.drawRectangle(ctx, 0, 0, width, height, color, (0.0, 0.0, 0.0, 0.0))
end

end

max_stars_repo_name: catapillie/VortexHelper

max_issues_repo_name: catapillie/VortexHelper

max_forks_repo_name: catapillie/VortexHelper


text:
using Random, LinearAlgebra

import DifferentialEquations: ODEProblem, solve
import Distributions: MvNormal, rand

export generate_data

"""
    generate_data( ssm, u0; seed=42)

from StateSpace generate:
 - true state (xt)
 - partial/noisy observations (yo)
 - catalog
"""
function generate_data(ssm::StateSpaceModel, u0::Vector{Float64}, seed = 42)

    rng = MersenneTwister(seed)

    try
        @assert ssm.dt_states < ssm.dt_obs
    catch
        @error " ssm.dt_obs must be bigger than ssm.dt_states"
    end

    try
        @assert mod(ssm.dt_obs, ssm.dt_states) == 0.0
    catch
        @error " ssm.dt_obs must be a multiple of ssm.dt_states "
    end

    tspan = (0.0, ssm.nb_loop_test)

    prob = ODEProblem(ssm.model, u0, tspan, ssm.params)

    # generSate true state (xt)
    sol = solve(prob, saveat = ssm.dt_states * ssm.dt_integration)
    xt = TimeSeries(sol.t, sol.u)

    # generate  partial/noisy observations (yo)
    nt = xt.nt
    nv = xt.nv

    yo = TimeSeries(xt.t, xt.u .* NaN)
    step = ssm.dt_obs ÷ ssm.dt_states
    nt = length(xt.t)

    d = MvNormal(ssm.sigma2_obs .* Matrix(I, nv, nv))
    ε = rand(rng, d, nt)
    for j = 1:step:nt
        for i in ssm.var_obs
            yo.u[j][i] = xt.u[j][i] + ε[i, j]
        end
    end

    # generate catalog
    u0 = last(sol)
    tspan = (0.0, ssm.nb_loop_train)
    prob = ODEProblem(ssm.model, u0, tspan, ssm.params)
    sol = solve(prob, saveat = ssm.dt_integration)
    n = length(sol.t)

    catalog_tmp = sol.u

    xt, yo, Catalog(hcat(catalog_tmp...), ssm)

end


text:
mutable struct GLLayer{
    T<:Number,T1<:Number,Tm1<:AbstractMatrix{T1},
    T2<:Number,Tm2<:AbstractMatrix{T2},
    Talg<:OrdinaryDiffEq.OrdinaryDiffEqCompositeAlgorithm
    }
    input_size::Int

    Win::Tm1

    G::SimpleGraph{Int64}
    W::Tm2
    γ::T
    Γ::T
    g::T
    D::T

    alg::Talg
    tinput::Float64
    toutput::Float64
end

function GLLayer(tinput::Real, toutput::Real, input_size::Int,
    γ::Number, Γ::Number, g::Number, G::SimpleGraph{Int64},
    D::Number=0.0, W=missing, Win=missing;
    alg::OrdinaryDiffEq.OrdinaryDiffEqCompositeAlgorithm=AutoTsit5(Rosenbrock23())
    )
    N = nv(G)
    @assert input_size > 0
    T = any(map(x->typeof(x)<:Union{Float32, ComplexF32}, [γ, Γ, g, D])) ? Float32 : typeof(γ)
    T = any(map(x->typeof(x)<:Complex, [γ, Γ, g, D])) ? complex(T) : T;

    _W = if ismissing(W)
        rand_couplings(N, T)
    else
        @assert size(W,1) == size(W,2) == N
        W
    end

    _Win = if ismissing(Win)
        rand_input_matrix(N, input_size, T)
    else
        @assert size(Win, 1) == N
        @assert size(Win, 2) == input_size
        Win
    end

    return GLLayer{T,eltype(_Win),typeof(_Win),eltype(_W),typeof(_W), typeof(alg)}(
          input_size, _Win, G, _W, T(γ), T(Γ), T(g), T(D),
          alg, Float64(tinput), Float64(toutput)
          )
end

function rand_couplings(N::Int,T::DataType)
    W = randn(T, N, N) ./ 2T(sqrt(N))
    for i in 1:N
        W[i,i] = zero(T)
    end
    (W .+ W') ./ 2
end

function rand_input_matrix(nr, nc, T)
    Win = rand(T, nr, nc)
    x = rand(T, nc)
    x ./= norm(x)
    Win ./= norm(Win * x)
    Win
end

function (m::GLLayer)(u::InputSignal)
    tspan = (m.tinput, m.toutput)

    function dψ_det!(dψ, ψ, p, t)
        mul!(dψ, m.Win, u(t))
        @. @inbounds dψ += (m.γ - (m.Γ + im*m.g)*abs2(ψ)) * ψ

        for v in vertices(m.G)
            Nv = neighbors(m.G, v)
            @inbounds @views dψ[v] -= im .* 10dot(m.W[v,Nv],ψ[Nv])
        end

        return dψ
    end

    function dψ_stoch!(dψ, ψ, p, t)
        @. @inbounds dψ = m.D

        return dψ
    end

    T = first(typeof(m).parameters)
    ψ0 = zeros(complex(T), nv(m.G))

    pb = SDEProblem(dψ_det!, dψ_stoch!, ψ0, tspan)
    solve(pb, m.alg; dtmax=0.1u.Δt, save_start=false, save_on=false).u[]
end

function (m::GLLayer)(u::AbstractArray{I}, args...; kwargs...) where I<:InputSignal
    tspan = (m.tinput, m.toutput)

    function dψ_det!(dψ, ψ, i, t)
        mul!(dψ, m.Win, u[i][1](t))
        @. @inbounds dψ += (m.γ - (m.Γ + im*m.g)*abs2(ψ)) * ψ

        for v in vertices(m.G)
            Nv = neighbors(m.G, v)
            @inbounds @views dψ[v] -= im .* 10dot(m.W[v,Nv],ψ[Nv])
        end
    end

    function dψ_stoch!(dψ, ψ, i, t)
        @. @inbounds dψ = m.D

        return dψ
    end

    T = first(typeof(m).parameters)
    ψ0 = zeros(complex(T), nv(m.G))

    pb = SDEProblem(dψ_det!, dψ_stoch!, ψ0, tspan, [1])
    function prob_func(prob, i, repeat)
        prob.p[1] = i
        prob
    end
    Δt = minimum(map(x->x.Δt, u))
    sol = solve(EnsembleProblem(pb, prob_func=prob_func), m.alg, args...; kwargs..., trajectories=length(u), dtmax=0.1Δt, save_start=false, save_on=false)
    reshape([sol[i].u[:][] for i in 1:length(u)], size(u))
end


text:
using Plots     # importing the package
gr()            # faster backend

# defining the range
x = 0:0.005:1

m = 1
b = 0
f = m * x .+ b  # linear
g = x .^ 2      # quadratic
h = x .^ 3      # cubic

# plot out the first function
# place legend in top left corner
plot(x, f, label="linear", legend=:topleft)
# edit existing plot
plot!(x, g, label="quadratic")
# increase resolution with dpi parameter
plot!(x, h, label="cubic", dpi=300)

# adding axis and title labels
xaxis!("x [-]")
yaxis!("f(x) [-]")
title!("polynomials")

# save the figure with string name
savefig("myawesomeplot.png")

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
dynamics(n::Neuron) = n.dynamics
has_dynamics(n::Neuron, species) = haskey(dynamics(n), species)

# defined for reversal and conductance since they belong to a neuron. not conductance (belongs to Link)
has_dynamics(n::Neuron, ::Type{Current{I}}) where {I} = true
# has_dynamics(n::Neuron, ::Type{Reversal{I}}) where {I} = haskey(dynamics(n), I)

struct EmptyNeuron{F<:Number} <: Neuron
    somatic_parameters::Dict{DataType,F}
end

EmptyNeuron() = EmptyNeuron(
    Dict(Voltage => -60.0)
)
dynamics(::EmptyNeuron) = Dict{DataType,SpeciesDynamics}()

struct BasicNeuron{G<:Geometry,C<:FlowChannel,F<:Real,S<:SpeciesDynamics} <: Neuron
    geometry::G
    dynamics::Dict{DataType,S}
    somatic_parameters::Dict{DataType,F}
    channels::Vector{C}
    name::Symbol
end

"""
building neuron
- geometric and electrical (ie capacitance) parameters are unique and should be defined by user
- number of hooks has to be pre-specified
- reversals are defined by the ion channels to which it is connected, same for currents
- channels which have a PlasticityRule (need extra sensors) are treated inside b the same way as those that don't
- update equations: voltage equation gets added synapses...for now. TODO Make synapses like channels

# DO LATER: to make it easier for the user, add an extra input which is var (the var in question). so they dont have to find eg voltage by searching through vars and varnames
# add b as input to channels so they can sense whether their inputs are parmeters or not, using b.dynamics
"""

function (b::BasicNeuron)(; incoming_connections::Integer=0)
    #shared -> hooks
    # e.g. species = Voltage or species = Potassium
    has_dynamics(species) = haskey(b.dynamics, species)
    # track union of things sensed by the connected channels
    tracked_names = vcat(
        Voltage,
        b.channels .|> sensed |> Iterators.flatten |> unique,
        [keys(b.dynamics)...]
    ) |> unique!

    state_indices = findall(has_dynamics, tracked_names)
    param_indices = findall(!has_dynamics, tracked_names)

    # build state/param ModelingToolkit variables for each of these tracked species
    tracked = zeros(Num, length(tracked_names))

    tracked[state_indices] .= reduce(vcat,
        tracked_names[state_indices] .|> shorthand_name
        |> instantiate_variables)

    tracked[param_indices] .= reduce(vcat,
        tracked_names[param_indices] .|> shorthand_name
        |> instantiate_parameters)

    syns = [@variables $el(t) for el in [Symbol(:Isyn, i) for i = 1:incoming_connections]]
    my_sum(syns) = incoming_connections == 0 ? 0.0 : sum(reduce(vcat, syns))
    !(incoming_connections == 0) && (syns = reduce(vcat, syns))
    chs = [ch(b) for ch in b.channels]

    tracked_fluxes = map(tracked_names[state_indices]) do thing
        sum(zip(b.channels, chs)) do (channel, channel_sys)
            get_actuator(channel, channel_sys, thing)
        end
    end

    # UGLY, make better. maybe add generality for syns:. Like hooks = {Voltage, 3}
    tracked_fluxes[findfirst(x -> x == Voltage, tracked_names)] += my_sum(syns)

    # hook connections between channel sensors and tracked variables.
    outward_connections = reduce(vcat, map(tracked_names, tracked) do species, variable
        [(variable ~ get_sensor(ch, ch_sys, species))::Equation for (ch, ch_sys) in zip(b.channels, chs) if (get_sensor(ch, ch_sys, species) !== nothing)]::Vector{Equation}
    end)

    outward_species = reduce(vcat, map(tracked_names) do species
        [species for (ch, ch_sys) in zip(b.channels, chs) if (get_sensor(ch, ch_sys, species) !== nothing)]
    end
    )

    outward_connection_indices = findall(outward_connections) do el
        !ModelingToolkit.isparameter(el.rhs)
    end

    # add currents and inputs to the dynamics of tracked variables
    inward_connections = [
        b.dynamics[species](b, tracked, tracked_names, flux) for (flux, species) in zip(tracked_fluxes, tracked_names[state_indices])
    ]

    # default initial conditions for neural state variables
    state_defaults = Dict(tracked[el] => b.somatic_parameters[tracked_names[el]] for el in state_indices)

    # default values for hooked parameter values in channels (their internal states/parameters already have defaults)
    outward_param_indices = findall(outward_connections) do el
        ModelingToolkit.isparameter(el.rhs)
    end
    parameter_defaults = Dict(outward_connections[el].rhs => b.somatic_parameters[outward_species[el]] for el in unique(outward_param_indices))

    # if there are any internal state variables from the b.dynamics laws, store and default them.
    somatic_states = get_from(b.dynamics, get_states)
    somatic_state_defaults = mapreduce(merge, values(b.dynamics)) do s
        default_states(s, b, tracked, tracked_names)
    end

    # if there are any internal parameter variables from the b.dynamics laws, store and default them.
    somatic_params = get_from(b.dynamics, get_parameters)
    somatic_param_defaults = mapreduce(merge, values(b.dynamics)) do s
        default_params(s, b, tracked, tracked_names)
    end

    sys = ODESystem(
        vcat(inward_connections, outward_connections[outward_connection_indices]),
        t,
        vcat(tracked[state_indices], somatic_states, syns),
        somatic_params
        ;
        systems=chs,
        defaults=merge(state_defaults, parameter_defaults, somatic_state_defaults, somatic_param_defaults),
        name=b.name
    )
    return (incoming_connections == 0) ? structural_simplify(sys) : sys
end

export get_from, get_states, default_states, default_params, get_parameters


text:
# These tests use BenchmarkTools, and are unreliable in an
# uncontrolled environment (like Travis)

using LazyArrays, BenchmarkTools

@testset "concat" begin
    A = Vcat(Vector(1:10), Vector(1:20))
    b = Array{Int}(undef, 30)
    @test @belapsed(copyto!($b,$A)) < @belapsed(vcat($A.arrays...))

    A = Vcat(1:10, 1:20)
    b = Array{Int}(undef, 30)
    @test @belapsed(copyto!($b,$A)) < @belapsed(vcat($A.arrays...))

    A = Vcat(randn(2,10), randn(4,10))
    b = Array{Float64}(undef, 6,10)
    @test @belapsed(copyto!($b,$A)) < @belapsed(vcat($A.arrays...))

    A = Hcat(1:10, 2:11)
    b = Array{Int}(undef, 10, 2)
    @test @belapsed(copyto!($b,$A)) < @belapsed(hcat($A.arrays...))
end


text:
const BESSEL_ORDER_INFO = """
derivatives of Bessel functions with respect to the order are not implemented currently:
https://github.com/JuliaMath/SpecialFunctions.jl/issues/160
"""

const INCOMPLETE_GAMMA_INFO = """
derivatives of the incomplete Gamma functions with respect to parameter `a` are not
implemented currently:
https://github.com/JuliaMath/SpecialFunctions.jl/issues/317
"""

const INCOMPLETE_EXPINT_INFO = """
derivatives of the exponential integral with respect to parameter `ν` are not implemented
currently:
https://github.com/JuliaMath/SpecialFunctions.jl/issues/321
"""

ChainRulesCore.@scalar_rule(airyai(x), airyaiprime(x))
ChainRulesCore.@scalar_rule(airyaix(x), airyaiprimex(x) + sqrt(x) * Ω)
ChainRulesCore.@scalar_rule(airyaiprime(x), x * airyai(x))
ChainRulesCore.@scalar_rule(airyaiprimex(x), x * airyaix(x) + sqrt(x) * Ω)
ChainRulesCore.@scalar_rule(airybi(x), airybiprime(x))
ChainRulesCore.@scalar_rule(airybiprime(x), x * airybi(x))
ChainRulesCore.@scalar_rule(besselj0(x), -besselj1(x))
ChainRulesCore.@scalar_rule(
    besselj1(x),
    (besselj0(x) - besselj(2, x)) / 2,
)
ChainRulesCore.@scalar_rule(bessely0(x), -bessely1(x))
ChainRulesCore.@scalar_rule(
    bessely1(x),
    (bessely0(x) - bessely(2, x)) / 2,
)
ChainRulesCore.@scalar_rule(dawson(x), 1 - (2 * x * Ω))
ChainRulesCore.@scalar_rule(digamma(x), trigamma(x))

# TODO: use `invsqrtπ` if it is added to IrrationalConstants
ChainRulesCore.@scalar_rule(erf(x), (2 * exp(-x^2)) / sqrtπ)
ChainRulesCore.@scalar_rule(erf(x, y), (- (2 * exp(-x^2)) / sqrtπ, (2 * exp(-y^2)) / sqrtπ))
ChainRulesCore.@scalar_rule(erfc(x), - (2 * exp(-x^2)) / sqrtπ)
ChainRulesCore.@scalar_rule(logerfc(x), - (2 * exp(-x^2 - Ω)) / sqrtπ)
ChainRulesCore.@scalar_rule(erfcinv(x), - (sqrtπ * (exp(Ω^2) / 2)))
ChainRulesCore.@scalar_rule(erfcx(x), 2 * (x * Ω - inv(oftype(Ω, sqrtπ))))
ChainRulesCore.@scalar_rule(logerfcx(x), 2 * (x - exp(-Ω) / sqrtπ))
ChainRulesCore.@scalar_rule(erfi(x), (2 * exp(x^2)) / sqrtπ)
ChainRulesCore.@scalar_rule(erfinv(x), sqrtπ * (exp(Ω^2) / 2))

ChainRulesCore.@scalar_rule(gamma(x), Ω * digamma(x))
ChainRulesCore.@scalar_rule(
    gamma(a, x),
    (
        ChainRulesCore.@not_implemented(INCOMPLETE_GAMMA_INFO),
        - exp(-x) * x^(a - 1),
    ),
)
ChainRulesCore.@scalar_rule(
    gamma_inc(a, x, IND),
    @setup(z = exp(-x) * x^(a - 1) / gamma(a)),
    (
        ChainRulesCore.@not_implemented(INCOMPLETE_GAMMA_INFO),
        z,
        ChainRulesCore.NoTangent(),
    ),
    (
        ChainRulesCore.@not_implemented(INCOMPLETE_GAMMA_INFO),
        -z,
        ChainRulesCore.NoTangent(),
    ),
)
ChainRulesCore.@scalar_rule(
    invdigamma(x),
    inv(trigamma(invdigamma(x))),
)
ChainRulesCore.@scalar_rule(trigamma(x), polygamma(2, x))

# Bessel functions
ChainRulesCore.@scalar_rule(
    besselj(ν, x),
    (
        ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO),
        (besselj(ν - 1, x) - besselj(ν + 1, x)) / 2
    ),
)
ChainRulesCore.@scalar_rule(
    besseli(ν, x),
    (
        ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO),
        (besseli(ν - 1, x) + besseli(ν + 1, x)) / 2,
    ),
)
ChainRulesCore.@scalar_rule(
    bessely(ν, x),
    (
        ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO),
        (bessely(ν - 1, x) - bessely(ν + 1, x)) / 2,
    ),
)
ChainRulesCore.@scalar_rule(
    besselk(ν, x),
    (
        ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO),
        -(besselk(ν - 1, x) + besselk(ν + 1, x)) / 2,
    ),
)
ChainRulesCore.@scalar_rule(
    besselkx(ν, x),
    (
        ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO),
        -(besselkx(ν - 1, x) + besselkx(ν + 1, x)) / 2 + Ω,
    ),
)
ChainRulesCore.@scalar_rule(
    hankelh1(ν, x),
    (
        ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO),
        (hankelh1(ν - 1, x) - hankelh1(ν + 1, x)) / 2,
    ),
)
ChainRulesCore.@scalar_rule(
    hankelh1x(ν, x),
    (
        ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO),
        (hankelh1x(ν - 1, x) - hankelh1x(ν + 1, x)) / 2 - im * Ω,
    ),
)
ChainRulesCore.@scalar_rule(
    hankelh2(ν, x),
    (
        ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO),
        (hankelh2(ν - 1, x) - hankelh2(ν + 1, x)) / 2,
    ),
)
ChainRulesCore.@scalar_rule(
    hankelh2x(ν, x),
    (
        ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO),
        (hankelh2x(ν - 1, x) - hankelh2x(ν + 1, x)) / 2 + im * Ω,
    ),
)

ChainRulesCore.@scalar_rule(
    polygamma(m, x),
    (
        ChainRulesCore.NoTangent(),
        polygamma(m + 1, x),
    ),
)
# todo: setup for common expr
ChainRulesCore.@scalar_rule(
    beta(a, b),
    (Ω*(digamma(a) - digamma(a + b)),
     Ω*(digamma(b) - digamma(a + b)),)
)
ChainRulesCore.@scalar_rule(
    logbeta(a, b),
    (digamma(a) - digamma(a + b),
     digamma(b) - digamma(a + b),)
)

# actually is the absolute value of the logorithm of gamma paired with sign gamma
ChainRulesCore.@scalar_rule(
    logabsgamma(x),
    digamma(x),
    ChainRulesCore.ZeroTangent()
)

ChainRulesCore.@scalar_rule(loggamma(x), digamma(x))
ChainRulesCore.@scalar_rule(
    loggamma(a, x),
    (
        ChainRulesCore.@not_implemented(INCOMPLETE_GAMMA_INFO),
        -exp(- (x + Ω)) * x^(a - 1),
    )
)

# exponential integrals
ChainRulesCore.@scalar_rule(expint(z), - exp(-z) / z)
ChainRulesCore.@scalar_rule(
    expint(ν, z),
    (
        ChainRulesCore.@not_implemented(INCOMPLETE_EXPINT_INFO),
        - expint(ν - 1, z),
    )
)
ChainRulesCore.@scalar_rule(expintx(z), Ω - inv(z))
ChainRulesCore.@scalar_rule(
    expintx(ν, z),
    (
        ChainRulesCore.@not_implemented(INCOMPLETE_EXPINT_INFO),
        Ω - expintx(ν - 1, z),
    )
)
ChainRulesCore.@scalar_rule(expinti(x), exp(x) / x)
ChainRulesCore.@scalar_rule(sinint(x), sinc(invπ * x))
ChainRulesCore.@scalar_rule(cosint(x), cos(x) / x)

# non-holomorphic functions
function ChainRulesCore.frule((_, Δν, Δx), ::typeof(besselix), ν::Number, x::Number)
    # primal
    Ω = besselix(ν, x)

    # derivative
    ∂Ω_∂ν = ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO)
    a = (besselix(ν - 1, x) + besselix(ν + 1, x)) / 2
    ΔΩ = if Δx isa Real
        muladd(muladd(-sign(real(x)), Ω, a), Δx, ∂Ω_∂ν * Δν)
    else
        muladd(a, Δx, muladd(-sign(real(x)) * real(Δx), Ω, ∂Ω_∂ν * Δν))
    end

    return Ω, ΔΩ
end
function ChainRulesCore.rrule(::typeof(besselix), ν::Number, x::Number)
    Ω = besselix(ν, x)
    project_x = ChainRulesCore.ProjectTo(x)
    function besselix_pullback(ΔΩ)
        ν̄ = ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO)
        a = (besselix(ν - 1, x) + besselix(ν + 1, x)) / 2
        x̄ = project_x(muladd(conj(a), ΔΩ, - sign(real(x)) * real(conj(Ω) * ΔΩ)))
        return ChainRulesCore.NoTangent(), ν̄, x̄
    end
    return Ω, besselix_pullback
end

function ChainRulesCore.frule((_, Δν, Δx), ::typeof(besseljx), ν::Number, x::Number)
    # primal
    Ω = besseljx(ν, x)

    # derivative
    ∂Ω_∂ν = ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO)
    a = (besseljx(ν - 1, x) - besseljx(ν + 1, x)) / 2
    ΔΩ = if Δx isa Real
        muladd(a, Δx, ∂Ω_∂ν * Δν)
    else
        muladd(a, Δx, muladd(-sign(imag(x)) * imag(Δx), Ω, ∂Ω_∂ν * Δν))
    end

    return Ω, ΔΩ
end
function ChainRulesCore.rrule(::typeof(besseljx), ν::Number, x::Number)
    Ω = besseljx(ν, x)
    project_x = ChainRulesCore.ProjectTo(x)
    function besseljx_pullback(ΔΩ)
        ν̄ = ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO)
        a = (besseljx(ν - 1, x) - besseljx(ν + 1, x)) / 2
        x̄ = if x isa Real
            project_x(a * ΔΩ)
        else
            project_x(muladd(conj(a), ΔΩ, - sign(imag(x)) * real(conj(Ω) * ΔΩ) * im))
        end
        return ChainRulesCore.NoTangent(), ν̄, x̄
    end
    return Ω, besseljx_pullback
end

function ChainRulesCore.frule((_, Δν, Δx), ::typeof(besselyx), ν::Number, x::Number)
    # primal
    Ω = besselyx(ν, x)

    # derivative
    ∂Ω_∂ν = ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO)
    a = (besselyx(ν - 1, x) - besselyx(ν + 1, x)) / 2
    ΔΩ = if Δx isa Real
        muladd(a, Δx, ∂Ω_∂ν * Δν)
    else
        muladd(a, Δx, muladd(-sign(imag(x)) * imag(Δx), Ω, ∂Ω_∂ν * Δν))
    end

    return Ω, ΔΩ
end
function ChainRulesCore.rrule(::typeof(besselyx), ν::Number, x::Number)
    Ω = besselyx(ν, x)
    project_x = ChainRulesCore.ProjectTo(x)
    function besselyx_pullback(ΔΩ)
        ν̄ = ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO)
        a = (besselyx(ν - 1, x) - besselyx(ν + 1, x)) / 2
        x̄ = if x isa Real
            project_x(a * ΔΩ)
        else
            project_x(muladd(conj(a), ΔΩ, - sign(imag(x)) * real(conj(Ω) * ΔΩ) * im))
        end
        return ChainRulesCore.NoTangent(), ν̄, x̄
    end
    return Ω, besselyx_pullback
end