text:
@testset "goto symbols" begin
    @testset "goto local symbols" begin
        let str = """
            function localgotoitem(word, path, column, row, startRow, context) # L0
              position = row - startRow                                        # L1
              ls = locals(context, position, column)                           # L2
              filter!(ls) do l                                                 # L3
                l[:name] == word &&                                            # L4
                l[:line] < position                                            # L5
              end                                                              # L6
              map(ls) do l # there should be zero or one element in `ls`       # L7
                text = l[:name]                                                # L8
                line = startRow + l[:line] - 1                                 # L9
                gotoitem(text, path, line)                                     # L10
              end                                                              # L11
            end                                                                # L12
            """,
            localgotoitem(word, line) = Atom.localgotoitem(word, "path", typemax(Int), line + 1, 0, str)[1]

            let item = localgotoitem("row", 2)
                @test item.line === 0
                @test item.text == "row"
                @test item.file == "path"
            end
            @test localgotoitem("position", 2).line === 1
            @test localgotoitem("l", 4).line === 3
            @test localgotoitem("l", 8).line === 7
        end

        # ignore dot accessors
        let str = """
            function withdots(expr::CSTParser.EXPR)
                bind = CSTParser.bindingof(expr.args[1])
                val = bind.val
                return val
            end
            """,
            localgotoitem(word, line) = Atom.localgotoitem(word, "path", typemax(Int), line + 1, 0, str)[1]

            @test localgotoitem("expr.args", 1).line === 0
            @test localgotoitem("bind.val", 2).line === 1
        end

        # don't error on fallback case
        @test (Atom.localgotoitem("word", nothing, 1, 1, 0, ""); true)
    end

    @testset "strip preceding ! operator" begin
        using Atom: gotosymbol
        let
            ret = gotosymbol("!isnothing")
            @test isdefined(ret, :items)
            @test !isempty(ret.items)
        end

        let
            ret = gotosymbol("!")
            @test isdefined(ret, :items)
            @test !isempty(ret.items)
        end
    end

    @testset "goto global symbols" begin
        using Atom: globalgotoitems, toplevelgotoitems, SYMBOLSCACHE, updatesymbols,
                    clearsymbols, regeneratesymbols, methodgotoitems

        ## strip a dot-accessed modules
        let
            path = joinpath′(@__DIR__, "..", "src", "comm.jl")
            text = read(path, String)
            items = globalgotoitems("Atom.handlers", Atom, path, text)
            @test !isempty(items)
            @test items[1].file == path
            @test items[1].text == "handlers"
            items = globalgotoitems("Main.Atom.handlers", Atom, path, text)
            @test !isempty(items)
            @test items[1].file == path
            @test items[1].text == "handlers"

            # can access the non-exported (non-method) bindings in the other module
            path = joinpath′(@__DIR__, "..", "src", "goto.jl")
            text = read(@__FILE__, String)
            items = globalgotoitems("Atom.SYMBOLSCACHE", Main, @__FILE__, text)
            @test !isempty(items)
            @test items[1].file == path
            @test items[1].text == "SYMBOLSCACHE"
        end

        @testset "goto modules" begin
            let item = globalgotoitems("Atom", Main, nothing, "")[1]
                @test item.file == joinpath′(atomsrcdir, "Atom.jl")
                @test item.line == 3
            end
            let item = globalgotoitems("SubJunk", Junk, nothing, "")[1]
                @test item.file == subjunkspath
                @test item.line == 3
            end
        end

        @testset "goto toplevel symbols" begin
            ## where Revise approach works, i.e.: precompiled modules
            let path = joinpath′(atomsrcdir, "comm.jl")
                mod = Atom
                key = "Atom"
                word = "handlers"

                # basic
                let items = toplevelgotoitems(word, mod, path, "")
                    @test !isempty(items)
                    @test items[1].file == path
                    @test items[1].text == word
                end

                # check caching works
                @test haskey(SYMBOLSCACHE, key)

                # check the Revise-like approach finds all files in Atom module
                @test length(SYMBOLSCACHE[key]) == length(atommodfiles)

                # when `path` isn't given, i.e. via docpane / workspace
                let items = toplevelgotoitems(word, mod, nothing, "")
                    @test !isempty(items)
                    @test items[1].file == path
                    @test items[1].text == word
                end

                # same as above, but without any previous cache -- falls back to CSTPraser-based module-walk
                delete!(SYMBOLSCACHE, key)

                let items = toplevelgotoitems(word, mod, nothing, "")
                    @test !isempty(items)
                    @test items[1].file == path
                    @test items[1].text == word
                end

                # check CSTPraser-based module-walk finds all the included files
                # NOTE: webio.jl is excluded since `include("webio.jl")` is a toplevel call
                @test length(SYMBOLSCACHE[key]) == length(atommodfiles)
            end

            ## where the Revise-like approach doesn't work, e.g. non-precompiled modules
            let path = junkpath
                mod = Main.Junk
                key = "Main.Junk"
                word = "toplevelval"

                # basic -- no need to pass a buffer text
                let items = toplevelgotoitems(word, mod, path, "")
                    @test !isempty(items)
                    @test items[1].file == path
                    @test items[1].line == 14
                    @test items[1].text == word
                end

                # check caching works
                @test haskey(Atom.SYMBOLSCACHE, key)

                # when `path` isn't given, i.e.: via docpane / workspace
                let items = toplevelgotoitems(word, mod, nothing, "")
                    @test !isempty(items)
                    @test items[1].file == path
                    @test items[1].line == 14
                    @test items[1].text == word
                end
            end

            ## escape recursive `include` loop
            let
                @test (Atom._collecttoplevelitems_static(nothing, joinpath′(fixturedir, "self_recur.jl"); inmod = true); true)
                @test (Atom._collecttoplevelitems_static(nothing, joinpath′(fixturedir, "mutual_recur.jl"); inmod = true); true)
            end

            ## don't include bindings outside of a module
            let path = subjunkspath
                text = read(subjunkspath, String)
                mod = Main.Junk.SubJunk
                word = "imwithdoc"

                items = toplevelgotoitems(word, mod, path, text)
                @test length(items) === 1
                if length(items) === 1
                    @test items[1].file == path
                    @test items[1].line == 6
                    @test items[1].text == word
                end
            end

            ## `Main` module -- use a passed buffer text
            let path = joinpath′(@__DIR__, "runtests.jl")
                text = read(path, String)
                mod = Main
                word = "atomsrcdir"

                items = toplevelgotoitems(word, mod, path, text)
                @test !isempty(items)
                @test items[1].file == path
                @test items[1].text == word
            end

            # don't mix a function and a macro with the same name
            let items = toplevelgotoitems("samefoo", Main.Junk, nothing, "")
                @test length(items) == 1
                @test items[1].name == "samefoo"
                @test items[1].text == "samefoo(args)"
                @test items[1].file == junkpath
                @test items[1].line == 22
            end
            let items = toplevelgotoitems("@samefoo", Main.Junk, nothing, "")
                @test length(items) === 1
                @test items[1].name == "@samefoo"
                @test items[1].text == "samefoo(ex)"
                @test items[1].file == junkpath
                @test items[1].line == 23
            end
        end

        @testset "updating toplevel symbols" begin
            # check there is no cache before updating
            mod = Main.Junk
            key = "Main.Junk"
            path = junkpath
            text = read(path, String)
            @test filter(SYMBOLSCACHE[key][path]) do item
                item.name == "toplevelval2"
            end |> isempty

            # mock updatesymbol handler
            originallines = readlines(path)
            newtext = join(originallines[1:end - 1], '\n')
            word = "toplevelval2"
            newtext *= "\n$word = :youshoulderaseme\nend"
            updatesymbols(key, path, newtext)

            # check the cache is updated
            @test filter(SYMBOLSCACHE[key][path]) do item
                item.name == word
            end |> !isempty

            let items = toplevelgotoitems(word, mod, path, newtext)
                @test !isempty(items)
                @test items[1].file == path
                @test items[1].text == "toplevelval2"
            end

            # re-update the cache
            updatesymbols(key, path, text)
            @test filter(SYMBOLSCACHE[key][path]) do item
                item.name == word
            end |> isempty

            # don't differentiate between `Mod` and `Main.Mod`
            @test !haskey(SYMBOLSCACHE, "Main.Atom")
            prevstate = SYMBOLSCACHE["Atom"][atomjlfile]
            updatesymbols("Main.Atom", atomjlfile, "")
            @test !haskey(SYMBOLSCACHE, "Main.Atom")
            @test isempty(SYMBOLSCACHE["Atom"][atomjlfile])
            SYMBOLSCACHE["Atom"][atomjlfile] = prevstate

            # don't error on fallback case
            @test (updatesymbols(key, nothing, text); true)
        end

        @static VERSION ≥ v"1.5" && @testset "regenerating toplevel symbols" begin
            try
                Pkg.develop(path = normpath(@__DIR__, "fixtures", "Example"))

                regeneratesymbols()

                # cache symbols in loaded modules
                @test haskey(SYMBOLSCACHE, "Base")
                @test length(keys(SYMBOLSCACHE["Base"])) > 100
                @test haskey(SYMBOLSCACHE, "Atom")
                @test length(keys(SYMBOLSCACHE["Atom"])) === length(atommodfiles)

                # cache symbols even if not loaded
                @test haskey(SYMBOLSCACHE, "Example")
            finally
                Pkg.rm("Example")
            end
        end

        @testset "clear toplevel symbols" begin
            clearsymbols()

            @test length(keys(SYMBOLSCACHE)) === 0
        end

        # NOTE: test after the cache gets cleared
        @testset "goto for unloaded packages" begin
            exampledir = normpath(@__DIR__, "fixtures", "Example")
            rootfile = normpath(exampledir, "src", "Example.jl")
            @test !isempty(globalgotoitems("hello", "Example", rootfile, read(rootfile, String)))
        end

        @testset "goto methods" begin
            ## basic
            let ms = methods(Atom.handlemsg)
                @test length(methodgotoitems(ms)) === length(ms)
            end

            ## aggregate methods with default params
            let m = eval(:(module $(gensym(:atomjltestmodule)) end))
                @eval m function funcwithdefaultargs(args, defarg = "default") end
                let items = methodgotoitems(methods(m.funcwithdefaultargs))
                    # should be handled as an unique method
                    @test length(items) === 1
                    # show a method with full arguments
                    @test "funcwithdefaultargs(args, defarg)" in map(i -> i.text, items)
                end

                @eval m function funcwithdefaultargs(args::String, defarg = "default") end
                let items = methodgotoitems(methods(m.funcwithdefaultargs))
                    # should be handled as different methods
                    @test length(items) === 2
                    # show methods with full arguments
                    @test "funcwithdefaultargs(args, defarg)" in map(i -> i.text, items)
                    @test "funcwithdefaultargs(args::String, defarg)" in map(i -> i.text, items)
                end
            end
        end

        ## both the original methods and the toplevel bindings that are overloaded in a context module should be shown
        let items = globalgotoitems("isconst", Main.Junk, nothing, "")
            @test length(items) === 2
            @test "isconst(m::Module, s::Symbol)" in map(item -> item.text, items) # from Base
            @test "Base.isconst(::JunkType)" in map(item -> item.text, items) # from Junk
        end

        ## don't error on the fallback case
        @test (globalgotoitems("word", Main, nothing, ""); true)
    end
end


text:
using Plots
using Statistics
using DSP

include("prony.jl")

function computePronyErr(line,N,p=2,q=2)
    err_test=[]
    for k = 1:length(line)-N
        line_start = k
        line_end = line_start + N - 1
        b, a, err = prony(line[line_start:line_end], p, q)

        append!(err_test,err)

        # myPoly = PolynomialRatio(b, a)
        # append!(mean_test,mean(impz(myPoly, N)))
        # append!(std_test,std(impz(myPoly, N)))
        # append!(max_test,maximum(impz(myPoly, N)))
    end
    return err_test
end

function searchSynchro(line,N,threshold=0.1,p=2,q=2)
    err_test=computePronyErr(line,N,p,q)
    first_err=findfirst(err_test.>threshold)-1
    first_impact=first_err.+(N-1)

    return err_test, first_err, first_impact
end

function plotSearchSynchro(line,err_test,first_err,first_impact,experiment,save_filename)
    p1=plot(line,label="signal",xlims = (0,length(line)),linecolor=:blue,title="Temporal signal")
    p1=plot!(p1,[first_err;first_err],[0;maximum(line)],linecolor=:red,label="First err")
    p1=plot!(p1,[first_impact;first_impact],[0;maximum(line)],linecolor=:black,label="First impact",linewidth = 2)

    p3=plot(err_test,label="err",xlims = (0,length(line)),linecolor=:blue,title="Curve of prony's error")
    p3=plot!([first_err;first_err],[0;maximum(err_test)],linecolor=:red,label="First err")
    p3=plot!([first_impact;first_impact],[0;maximum(err_test)],linecolor=:black,label="First impact")

    pTotal=plot(p1, p3, layout = (2, 1),title=string("Experiment : ",experiment))
    savefig(pTotal,save_filename)
end

function searchAllSynchroMTiXsens(AccFree_MTiXsens,filenameMTiXsens,sensorName,N=20,threshold=0.1,p=2,q=2,f_coupure=15,f_sampling=100)
    index_test=[]
    # experiment=1
    for experiment=1:length(filenameMTiXsens[:,1]) # =1:length(filenameMTiXsens[:,1]))
        @info experiment
        xline = AccFree_MTiXsens[filenameMTiXsens[experiment,1]][sensorName][:,:FreeAcc_X]
        yline = AccFree_MTiXsens[filenameMTiXsens[experiment,1]][sensorName][:,:FreeAcc_Y]

        designmethod = Butterworth(5)
        ff = digitalfilter(Lowpass(f_coupure/f_sampling),designmethod)
        xline=filtfilt(ff,xline)
        yline=filtfilt(ff,yline)

        line = xline.^2 .+yline.^2
        # line=line./maximum(line)

        err_test, first_err, first_impact = searchSynchro(line,N,threshold,p,q)
        # print(err_test)
        append!(index_test,AccFree_MTiXsens[filenameMTiXsens[experiment,1]][sensorName][first_impact,:PacketCounter])

        save_filename=string("save_figure/figure_prony/MTiXsens_",experiment,".png")
        plotSearchSynchro(line,err_test,first_err,first_impact,filenameMTiXsens[experiment,1],save_filename)
    end
    println("search Synchro MTiXsens Finished")

    return index_test
end

function searchAllSynchro(data,m,filename,N=20,threshold=0.1,p=2,q=2)
    index_test=[]
    experiment=1
    for experiment=1:length(filename[:,1]) # =1:length(filenameMTiXsens[:,1]))
        @info experiment
        line = data[filename[experiment,1]][:,m]
        line=line.-mean(line[1:25])
        # line=line./maximum(line)

        f_coupure=25
        f_sampling=50

        designmethod = Butterworth(5)
        ff = digitalfilter(Lowpass(f_coupure/f_sampling),designmethod)
        line=filtfilt(ff,line)

        err_test, first_err, first_impact = searchSynchro(line,N,threshold,p,q)
        # print(err_test)
        append!(index_test,first_impact)

        save_filename=string("save_figure/figure_prony/",replace(filename[experiment,1]," "=> "_"),".png")
        plotSearchSynchro(line,err_test,first_err,first_impact,filename[experiment,1],save_filename)
    end
    println("search Synchro Finished")

    return index_test
end

ext: jl


text:
using SearchLight
using SearchLightPostgreSQL
SearchLight.Configuration.load() |> SearchLight.connect
SearchLight.Migrations.create_migrations_table()
SearchLight.Migrations.up()
using BitemporalPostgres, Dates, SearchLight, Test, TimeZones
#=workflow w1 blue rectangle
=#
w1 = Workflow()
w1.tsw_validfrom = ZonedDateTime(2014, 5, 30, 21, 0, 1, 1, tz"Africa/Porto-Novo")
SearchLight.Configuration.load() |> SearchLight.connect
t = TestDummyComponent()
tr = TestDummyComponentRevision(description = "blue")
ts = TestDummySubComponent(ref_super = t.id)
tsr = TestDummySubComponentRevision(description = "green")
create_entity!(w1)
create_component!(t, tr, w1)
println(tr)
create_subcomponent!(t, ts, tsr, w1)
w1.ref_history != Nothing
@testset "create entity" begin
    @test w1.is_committed == 0
    @test w1.ref_version == tr.ref_validfrom
    @test w1.ref_version == tsr.ref_validfrom
end
#=
committing w1 blue rectangle
=#
@testset "commit create entity" begin
    commit_workflow!(w1)
    @test w1.is_committed == 1
end
#=
workflow w2 yellow rectangle overlaps blue one
=#
w2 = Workflow(
    ref_history = w1.ref_history,
    tsw_validfrom = ZonedDateTime(2015, 5, 30, 21, 0, 1, 1, tz"Africa/Porto-Novo"),
)
tr2 = copy(tr)
tr2.description = "yellow"
update_entity!(w2)
update_component!(tr, tr2, w2)
println(tr2)
@testset "update entity" begin
    @test w2.ref_version == tr2.ref_validfrom
    @test w2.ref_version == tr.ref_invalidfrom
end
#=
committing w2 yellow rectangle
=#
commit_workflow!(w2)
@testset "commit update entity" begin
    @test w2.is_committed == 1
end
#=
workflow w3 red rectangle shadows red and overlaps blue one
=#
w3 = Workflow(
    ref_history = w1.ref_history,
    tsw_validfrom = ZonedDateTime(2014, 11, 30, 21, 0, 1, 1, tz"Africa/Porto-Novo"),
)
tr3 = copy(tr2)
tr3.description = "red"
update_entity!(w3)
update_component!(tr2, tr3, w3)
delete_component!(ts, w3)
tsr = find(TestDummySubComponentRevision, SQLWhereExpression("id=?", 1))[1]
println(tr3)
@testset "retrospective update entity " begin
    @test w3.ref_version == tr3.ref_validfrom
    @test w3.ref_version == tr2.ref_invalidfrom
    @test w3.ref_version == tsr.ref_invalidfrom
end
#= committing w3 red rectangle
=#

commit_workflow!(w3)
@testset "commit retrospective update entity " begin
    @test w3.is_committed == 1
end
#=
Testing reading
=#
@testset "reading tests" begin
    v1 = findversion(w1.ref_history, w1.tsdb_validfrom, w1.tsw_validfrom)

    @test findcomponentrevision(TestDummyComponentRevision, t.id, v1)[1].description ==
          "blue"
    @test findcomponentrevision(TestDummySubComponentRevision, ts.id, v1)[1].description ==
          "green"

    v2 = findversion(w2.ref_history, w2.tsdb_validfrom, w2.tsw_validfrom)

    @test findcomponentrevision(TestDummyComponentRevision, t.id, v2)[1].description ==
          "yellow"
    @test findcomponentrevision(TestDummySubComponentRevision, ts.id, v2)[1].description ==
          "green"

    v2a = findversion(w2.ref_history, w2.tsdb_validfrom, w2.tsw_validfrom - Dates.Second(1))

    @test findcomponentrevision(TestDummyComponentRevision, t.id, v2a)[1].description ==
          "blue"
    @test findcomponentrevision(TestDummySubComponentRevision, ts.id, v2a)[1].description ==
          "green"

    v3 = findversion(w3.ref_history, w3.tsdb_validfrom, w3.tsw_validfrom)

    @test findcomponentrevision(TestDummyComponentRevision, t.id, v3)[1].description ==
          "red"
    @test isempty(findcomponentrevision(TestDummySubComponentRevision, ts.id, v3))

    v3a = findversion(w3.ref_history, w3.tsdb_validfrom, w3.tsw_validfrom - Dates.Second(1))

    @test findcomponentrevision(TestDummyComponentRevision, t.id, v3a)[1].description ==
          "blue"
    @test findcomponentrevision(TestDummySubComponentRevision, ts.id, v3a)[1].description ==
          "green"

end

w4 = Workflow(
    ref_history = w1.ref_history,
    tsw_validfrom = ZonedDateTime(2017, 11, 30, 21, 0, 1, 1, tz"Africa/Porto-Novo"),
)
tr4 = copy(tr3)
tr4.description = "green"
t1 = TestDummyComponent()
tr5 = TestDummyComponentRevision(description = "pink")

# @test tr3.ref_invalidfrom==MaxVersion

update_entity!(w4)
update_component!(tr3, tr4, w4)
println(tr4)
create_component!(t1, tr5, w4)
println(tr5)

@testset "pending transaction tests" begin
    v4 = find(Version, SQLWhereExpression("id=?", w4.ref_version))[1].id
    @test findcomponentrevision(TestDummyComponentRevision, t.id, v4)[1].description ==
          "green"
    @test findcomponentrevision(TestDummyComponentRevision, t1.id, v4)[1].description ==
          "pink"
    @test w4.ref_version == tr3.ref_invalidfrom
    @test w4.ref_version == tr4.ref_validfrom
    @test MaxVersion == tr4.ref_invalidfrom
    @test w4.ref_version == tr5.ref_validfrom
end
println(tr4)
println(tr3)

currentVersion = find(Version, SQLWhereExpression("id=?", w4.ref_version))[1]
currentInterval =
    find(ValidityInterval, SQLWhereExpression("ref_version=?", w4.ref_version))[1]
tr3b = find(TestDummyComponentRevision, SQLWhereExpression("id=?", tr3.id))[1]
# delete(currentVersion)
rollback_workflow!(w4)
tr3b = find(TestDummyComponentRevision, SQLWhereExpression("id=?", tr3.id))[1]
println(tr3b)
@testset "rollback transaction tests" begin
    @test tr3b.ref_invalidfrom == MaxVersion
    @test findcomponentrevision(TestDummyComponentRevision, t.id, DbId(4))[1].description ==
          "red"
    @test isempty(find(TestDummyComponentRevision, SQLWhereExpression("id=?", tr4.id)))
    @test isempty(find(TestDummyComponent, SQLWhereExpression("id=?", t1.id)))
    @test isempty(find(Version, SQLWhereExpression("id=?", currentVersion.id)))
    @test isempty(find(ValidityInterval, SQLWhereExpression("id=?", currentInterval.id)))

end

@testset "mkforest" begin
    hforest =
        mkforest(DbId(1), MaxDate, ZonedDateTime(1900, 1, 1, 0, 0, 0, 0, tz"UTC"), MaxDate)
    print_tree(hforest)
    @test hforest[1].interval.ref_version.value == 3
    @test hforest[1].shadowed[1].interval.ref_version.value == 2
    @test hforest[2].interval.ref_version.value == 1
    @test isempty(hforest[2].shadowed)
end


text:
#=
The prime 41, can be written as the sum of six consecutive primes:

41 = 2 + 3 + 5 + 7 + 11 + 13
This is the longest sum of consecutive primes that adds to a prime below one-hundred.

The longest sum of consecutive primes below one-thousand that adds to a prime, contains 21 terms, and is equal to 953.

Which prime, below one-million, can be written as the sum of the most consecutive primes?
=#
using Primes

function calc()

  x = 10^6
  p = primes(x-1)
  largestrun = largestprime = 0

  for i in 1:length(p)
    currentsum = 0
    for j in i:length(p)
      currentsum += p[j]
      if currentsum >= x
        break
      end
      if isprime(currentsum) && j-i > largestrun
        largestrun = j-i
        largestprime = currentsum
      end
    end
  end
  largestprime
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
export estimate_rs

"""
This is a MLE for fitting the selection parameters to a time trace, given the neutral parameters
"""

##
"""
The log-likelihood difference is computed using multiple importance sampling
with a mixture of the neutral-beta and the binomial-reweighted neutral sample
1/2 * B(Dq, D(1-q)) + 1/2 * B(Dq+k, D(1-q) + N-k).
This give the benefit of using the same samples for both likelihoods, greatly reducing the noise.

"""

logbinomial(j, k) =  loggamma(1+j+k) - loggamma(1+j) - loggamma(1+k)
##
function randlogGamma(a,n;scale = 1)
    if a > .2
        return log.(rand(Gamma(a),n))
    else
        L = 1 / a - 1
        w = a / MathConstants.e / (1 - a)
        ww = 1 / (1 + w)
        η = z -> ( z >= 0 ? exp(-z) : w * L * exp(L * z))
        h = z -> exp(-z - exp(-z / a))
        function rh(a)
            while true
                U = rand()
                z = (U <= ww) ? -log(U / ww) : log(rand()) / L
                h(z) / η(z) > rand() && return(z)
            end
        end
        Z = Array{Float64}(undef,n)
        for ii in eachindex(Z)
            Z[ii] = log(scale) - rh(a) / a
        end
        return Z
    end
end

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

# Likelihood parameterized strictly in terms of exponents
function σ_likelihood(θ_wt, θ_mut, k_wt, k_mut) # mutant and wildtype counts
    n_samp = 1000 # hyperparameter
    loggam_a = vcat(
        randlogGamma(θ_wt, n_samp),
        randlogGamma(θ_wt+k_wt, n_samp)
    )
    loggam_b = vcat(
        randlogGamma(θ_mut, n_samp),
        randlogGamma(θ_mut+k_mut, n_samp)
    )
    xs = 1 ./(1 .+ exp.(loggam_b .- loggam_a))
    δ_z = logbeta(θ_wt, θ_mut) - logbeta(θ_wt+k_wt, θ_mut +k_mut)
    δ_ls = map((a,b) -> k_wt*a + k_mut*b - (k_mut + k_wt)*(max(a,b) + log1p(exp(-abs(a-b)))) + δ_z, loggam_a, loggam_b) # log(p(d+k)/p(d))
    σ -> log(sum( exp(x*σ) / (1 + exp(δ_l)) for (x,δ_l) in zip(xs,δ_ls))) -
        log(sum( exp(x*σ) / (1 + exp(-δ_l)) for (x,δ_l) in zip(xs,δ_ls))) -
        logbinomial(k_wt,k_mut)
end

# defines an importance sampler structure for defining a likelihood

struct LikelihoodSample
    w_θ::Array{Float64,1}
    w_k::Array{Float64,1}
    x::Array{Float64,1}
    θ_ts::Float64
end

function LikelihoodSample(μ_wt, μ_mut, θ_ts, k_wt, k_mut; n_samp = 1000)
    (θ_wt, θ_mut) = θ_ts .* (μ_wt, μ_mut)
    # quenched and variance reduced sampler
    loggam_a = vcat(
        randlogGamma(θ_wt, n_samp),
        randlogGamma(θ_wt+k_wt, n_samp)
    )
    loggam_b = vcat(
        randlogGamma(θ_mut, n_samp),
        randlogGamma(θ_mut+k_mut, n_samp)
    )
    δ_z = logbeta(θ_wt, θ_mut) - logbeta(θ_wt+k_wt, θ_mut +k_mut)
    δs = map((a,b) -> k_wt*a + k_mut*b - (k_mut + k_wt)*(max(a,b) + log1p(exp(-abs(a-b)))) + δ_z,
        loggam_a, loggam_b) # log(p(d+k)/p(d))
    LikelihoodSample(
        2 .* sigmoid.(-δs) ,
        2 .* sigmoid.(δs) ,
        sigmoid.(loggam_a .- loggam_b),θ_ts)
end

function like_rs(rs, ls::LikelihoodSample)
    θ_ts = ls.θ_ts
    z_θ = mean(w_θ*exp(θ_ts*rs*x) for (w_θ, x) in zip(ls.w_θ, ls.x))
    z_k = mean(w_k*exp(θ_ts*rs*x) for (w_k, x) in zip(ls.w_k, ls.x))
    return log(z_k/z_θ)
end

function grad_hess_rs(rs, ls::LikelihoodSample)
    θ_ts = ls.θ_ts
    z_θ = mean(w_θ*exp(θ_ts*rs*x) for (w_θ, x) in zip(ls.w_θ, ls.x))
    z_k = mean(w_k*exp(θ_ts*rs*x) for (w_k, x) in zip(ls.w_k, ls.x))
    x_θ = mean(w_θ*exp(θ_ts*rs*x)*x for (w_θ, x) in zip(ls.w_θ, ls.x)) / z_θ
    x_k = mean(w_k*exp(θ_ts*rs*x)*x for (w_k, x) in zip(ls.w_k, ls.x)) / z_k
    xx_θ = mean(w_θ*exp(θ_ts*rs*x)*x^2 for (w_θ, x) in zip(ls.w_θ, ls.x)) / z_θ
    xx_k = mean(w_k*exp(θ_ts*rs*x)*x^2 for (w_k, x) in zip(ls.w_k, ls.x)) / z_k
    g = (x_k - x_θ)*θ_ts
    h = ((xx_k - x_k^2) - (xx_θ - x_θ^2)) *θ_ts^2
    return (g,h)
end

function grad_hess_pop(rs,ls::Array{Array{LikelihoodSample,1},1}; average = true)
    g = 0
    h = 0
    for lsp in ls
        gp = 0
        hp = 0
        for lspt in lsp
            (gpt,hpt) =  grad_hess_rs(rs,lspt)
            gp += gpt
            hp += hpt
        end
        if average
            g += gp / length(lsp) #
            h += hp /length(lsp) # time averaging can be turned off here by commenting "/length(lsp)"
        else
            g += gp  #
            h += hp  # time averaging can be turned off here by commenting "/length(lsp)"
        end
    end
    return (g,h)
end

function rs_newton_search(ls; average = true)
    rs = 0
    while true
        (g,h) = grad_hess_pop(rs,ls ; average = average)
        if |(g^2/abs(h) < 10^(-8), abs(rs)>10^4)
            return (rs, g, h)
        end
        rs -= g/h
    end
end

function like_rs_pop(rs, ls; avg = 0)
    # ls is an array of arrays of likelihood sampers, [patients[timpoints]], where each time point
    # avg is a number between 0 and 1
    l = 0
    for lsp in ls
        lp = 0
        for lspt in lsp
            lp += like_rs(rs,lspt)
        end
        l += ( 1/length(lsp) * avg + 1 - avg)*lp
    end
    return l
end

function baysian_pop_rs(ls; samples = 10^4, burn_in = 10^3, avg = 0) # ls is an array of likelihood sampers
    ll = rs -> like_rs_pop(rs, ls; avg = avg)
    baysian_sample_rs(ll; samples = samples, burn_in = burn_in)
end

function baysian_sample_rs(ll; samples = 10^4, burn_in = 50)
    # takes a log-likelihood and returns the sampled distribuiton with
    run  = Float64[]
    rs = 100
    for ii in 1:(samples + burn_in)
        rs = mh_step(ll,rs, σ = 50)
        push!(run,rs)
    end
    return run[burn_in+1:end]
end

function mh_step(ll, θ; σ = 50) # 1-d mh_stepper
    oldll = ll(θ)
    new = θ + randn().* σ
    newll = ll(new)
    rand() < exp(newll - oldll) ? out = new : out = θ # ll difference, always accept if positive
    return out
end

# function get_countdata(df, locus, time_points)
#     ind = map( t->findfirst(x -> x>t, df.time), time_points)
#     ind = ind[ind .!= nothing]
#     Ds = df.D[ind] # population diversity
#     # subsampling element (this may be neccesary to limit the large number of genomes)
#     Ns = df.pop_size[ind] # total number of indivudals
#     Ks = round.(Int,map(x -> x[locus], df.freq[ind]) .* Ns) #number of mutants.
#     Js = Ns .- Ks
#     return (Ds, Js, Ks)
# end

function make_likelihood(df; locus = 1, timepoints = 100:100:1000, sub_n = 100)
    # Return the likelihood function for a single time trace
    ind = map( t->findfirst(x -> x>t, df.time),timepoints)
    ind = ind[ind .!= nothing]
    if length(ind) <= 1
        return (x->0)
    end
    θs = df.θ[ind] # population diversity
    # subsampling element (respresenting incomplete information in a typical genetic data point)
    Ns = df.pop_size[ind] # total number of indivudals
    Ks = round.(Int,map(x -> x[locus], df.freq[ind]) .* Ns) #number of mutants.
    Js = Ns .- Ks
    ks = map((k,j) -> rand(Hypergeometric(k, j, min(sub_n,k+j))) .* [1,-1] .+ [0,  min(sub_n,k+j)] , Ks, Js)
    # if pop.size is smaller than the subsample size (sub_n), than default to pop.size
    fn_list = map((θ_ref,k)->rs_likelihood(θ_ref, [1,1], k...), θs, ks)
    fn_normalize = rs_likelihood(10^-5, [1,1], 1,1)
    # this function creates a minimum at finite Δ even if there is no mutants observed
    # D sets rough scale at which fitness effect will be invisible.
    rs -> sum(f(rs) for f in fn_list) + 0.01*fn_normalize(rs)
end

function estimate_rs(dfs::Vector; locus = 1, timepoints = 100:100:1000, sub_n = 100)
    fnlist = map(x-> make_likelihood(x; locus = locus, timepoints = timepoints, sub_n = sub_n),dfs)
    like(rs) = mapreduce(f->f(rs...), +, fnlist)
    s0 = [0.0]
    y = Optim.minimizer(optimize(like, s0, LBFGS(), autodiff = :forward))
    return -y[1]
end

function get_Δ(dfs::Vector; locus = 1, timepoints = 100:100:1000, sub_n = 100)
    fnlist = map(x-> make_likelihood(x; locus = locus, timepoints = timepoints, sub_n = sub_n),dfs)
    like(Δ) = mapreduce(f->f(Δ...), +, fnlist)
    return like
end

max_stars_repo_name: StatPhysBio/Tomoko

max_issues_repo_name: StatPhysBio/Tomoko

max_forks_repo_name: StatPhysBio/Tomoko


text:
"Structure representing a real number that is computable to arbitrary precision."
mutable struct ExactReal <: Real
    f::Function
    largest_n::Int  # largest n for which computed
    largest_value::BigInt  # value at that n
    msd::Int  # minimum significant digit
end

ExactReal(f) = ExactReal(f, -1, -1, -1)

"""
Memoized evaluation of the nth value in the representation
"""
function (x::ExactReal)(n::Integer)

    if n == x.largest_n
        return x.largest_value

    elseif n < x.largest_n
        return x.largest_value >> (x.largest_n - n)

    end

    value = (x.f)(n)

    x.largest_n = n
    x.largest_value = value

    return value
end

const B = big(2)  # base of the expansion

ExactReal(q::Rational) = ExactReal(n -> ((B^n) * numerator(q)) ÷ denominator(q))
ExactReal(q::Integer) = ExactReal(q//1)

const exact_real_precision = big(100)

"""
Returns an interval (as a tuple) containing the exact real number;
this interval is of approximate width `2^(-n)`.
"""
function value_range(x::ExactReal, n=exact_real_precision)
    numerator = x(n)
    denominator = big(2)^n

    return (Float64( (numerator - 1) / denominator, RoundDown),
            Float64( (numerator + 1) / denominator, RoundUp))
end

function Base.show(io::IO, r::ExactReal)
    lo, hi = value_range(r)

    print(io, "ExactReal in ",
        lo, "..", hi)

end


text:
module debug

# Home to all of the routines that are useful for debugging purposes

using ..constants
using ..model
using ..image

import PyPlot
import PyPlot.plt
using LaTeXStrings

export format_grid, synthesize_image, plot_chmaps, plot_slices

"""Format grid values to a standard string"""
function format_grid(grid::Grid, npix::Int)
    #nr, ntheta, r_in, r_out, npix
    r_in = @sprintf("%.1f", grid.rs[1]/AU)
    r_out = @sprintf("%.1f", grid.rs[end]/AU)
    return "$(grid.nr)_$(grid.ntheta)_$(r_in)_$(r_out)_$npix"
end

"""Synthesize images for these parameters, grid, and npix values. dlams provides the frequencies at which to synthesize images. Assumes that you have run DJInitialize.jl inside of homedir. When done, copies the image to outfile."""
function synthesize_image(dlams::Array{Float64, 1}, pars::Parameters, grid::Grid, npix::Int, outfile::AbstractString, species::AbstractString)

    # Assumes that all of the necessary files are already in the current working directory.
    write_grid("", grid)

    # Compute parameter file using model.jl, write to disk in current directory
    write_model(pars, "", grid, species)

    # Doppler shift the dataset wavelengths to rest-frame wavelength
    nchan = length(dlams)
    beta = pars.vel/c_kms # relativistic Doppler formula
    lams = Array(Float64, nchan)
    for i=1:nchan
        lams[i] =  dlams[i] * sqrt((1. - beta) / (1. + beta)) # [microns]
    end

    write_lambda(lams, "") # write into current directory

    # Run RADMC-3D, redirect output to /dev/null
    run(pipeline(`radmc3d image incl $(pars.incl) posang $(pars.PA) npix $npix loadlambda`, DevNull))

    # Copy the output image to the new filename
    run(`mv image.out image_$(outfile).out`)
end

# Plot the channel maps using sky convention
function plot_chmaps(img::image.SkyImage, lam0, fname="channel_maps_sky.png")

    vels = c_kms * (img.lams .- lam0)/lam0

    (im_ny, im_nx) = size(img.data)[1:2] # y and x dimensions of the image

    # Image needs to be flipped along RA dimension
    ext = (img.ra[end], img.ra[1], img.dec[1], img.dec[end])

    nlam = length(img.lams)

    # Figure out how many plots we'll have.
    ncols = 8
    nrows = ceil(Int, nlam/ncols)

    fig, ax = plt[:subplots](nrows=nrows, ncols=ncols, figsize=(12, 1.5 * nrows))

    # max = maximum(log10(img.data))

    vvmax = maxabs(img.data)
    norm = PyPlot.matplotlib[:colors][:Normalize](0, vvmax)

    for row=1:nrows
        for col=1:ncols
            iframe = col + (row - 1) * ncols

            if col != 1 || row != nrows
                ax[row, col][:xaxis][:set_ticklabels]([])
                ax[row, col][:yaxis][:set_ticklabels]([])
            else
                ax[row, col][:set_xlabel](L"$\Delta \alpha$ ('')")
                ax[row, col][:set_ylabel](L"$\Delta \delta$ ('')")
            end

            if iframe > nlam
                # Plot a blank square if we run out of channels
                ax[row, col][:imshow](zeros((im_ny, im_nx)), cmap=plt[:get_cmap]("PuBu"), vmin=0, vmax=20, extent=ext, origin="lower")

            else
                #Flip the frame for Sky convention
                frame = flipdim(img.data[:,:,iframe], 2)

                # frame += 1e-99 #Add a tiny bit so that we don't have log10(0)
                # lframe = log10(frame)
                #
                # ax[row, col][:imshow](lframe, extent=ext, interpolation="none", vmin=max - 6, vmax=max, origin="lower", cmap=plt[:get_cmap]("PuBu"))

                ax[row, col][:imshow](frame, extent=ext, interpolation="none", origin="lower", cmap=plt[:get_cmap]("PuBu"), norm=norm)

                ax[row, col][:annotate](@sprintf("%.1f", vels[iframe]), (0.1, 0.8), xycoords="axes fraction", size=8)
            end
        end
    end
    fig[:subplots_adjust](hspace=0.06, wspace=0.01, top=0.9, bottom=0.1, left=0.05, right=0.95)
    plt[:savefig](fname)
end

function plot_slices(img::image.SkyImage, index, fname="")
end

end #module

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
function Newton_Spectral(
  nlp::AbstractNLPModel,
  x₀::AbstractVector;
  τ₀::Float64 = 0.0005,
  ϵ::Float64 = 1e-6,
  maxiter::Int = 200,
)
  x = copy(x₀)
  iter = 0
  f, g = obj(nlp, x), grad(nlp, x)

  while (norm(g, Inf) > ϵ) && (iter <= maxiter)
    H = Matrix(Symmetric(hess(nlp, x), :L))
    Δ, O = eigen(H)

    # Boost negative values of Δ to 1e-8
    D = Δ .+ max.((1e-8 .- Δ), 0.0)

    d = -O * diagm(1.0 ./ D) * O' * g

    # Simple Armijo backtracking
    hp0 = g' * d
    t = 1.0
    ft = obj(nlp, x + t * d)
    while ft > (f + τ₀ * t * hp0)
      t /= 2.0
      ft = obj(nlp, x + t * d)
    end
    x += t * d
    f, g = ft, grad(nlp, x)
    iter += 1
  end
  if iter > maxiter
    @warn "Iteration limit"
  end

  return x, f, norm(g), iter
end


text:
module CalibrationErrors

using DataStructures
using Distances
using Parameters
using StatsBase

using LinearAlgebra
using Statistics

# estimation
export calibrationerror

# estimators
export ECE, BiasedSKCE, LinearUnbiasedSKCE, QuadraticUnbiasedSKCE

# binning algorithms
export UniformBinning, MedianVarianceBinning

# scalar-valued kernels
export SquaredExponentialKernel, ExponentialKernel, GaussianKernel, LaplacianKernel

# matrix-valued kernels
export UniformScalingKernel, DiagonalKernel

include("generic.jl")
include("binning/generic.jl")
include("binning/uniform.jl")
include("binning/medianvariance.jl")
include("ece.jl")
include("kernels/generic.jl")
include("kernels/scalar.jl")
include("kernels/matrix.jl")
include("skce/generic.jl")
include("skce/biased.jl")
include("skce/unbiased.jl")

end # module


text:
using BenchmarkTools
using ForwardDiff
using Test
using Random
using StaticArrays
using LinearAlgebra
import TrajectoryOptimization: kinematics, ∇rotate, ∇composition1, ∇composition2, vector
const TO = TrajectoryOptimization

q1 = rand(UnitQuaternion)
q2 = rand(UnitQuaternion)
r = @SVector rand(3)
ω = @SVector rand(3)

ϕ = 0.1*@SVector [1,0,0]
@test logm(expm(ϕ)) ≈ ϕ
@test expm(logm(q1)) ≈ q1
@test angle(expm(ϕ)) ≈ 0.1
# @btime expm($ϕ)
# @btime logm($u1)
@test norm(q1 ⊕ ϕ) ≈ 1
@test q1 ⊖ q2 isa SVector{3}
@test (q1 ⊕ ϕ) ⊖ q1 ≈ ϕ

# Test inverses
q3 = q2*q1
@test q2\q3 ≈ q1
@test q3/q1 ≈ q2
@test inv(q1)*r ≈ q1\r
@test r ≈ q3\(q2*q1*r)

q = q1
rhat = UnitQuaternion(r)
@test q*r ≈ Vmat()*Lmult(q)*Rmult(q)'*Vmat()'r
@test q*r ≈ Vmat()*Lmult(q)*Lmult(rhat)*Tmat()*SVector(q)
@test q*r ≈ Vmat()*Rmult(q)'*Rmult(rhat)*SVector(q)

@test ForwardDiff.jacobian(q->UnitQuaternion{VectorPart}(q)*r,SVector(q)) ≈ ∇rotate(q,r)
# @btime ForwardDiff.jacobian(q->UnitQuaternion(q)*$r,SVector($q))
# @btime ∇rotate($q,$r)

@test ForwardDiff.jacobian(q->SVector(q2*UnitQuaternion{VectorPart}(q)),SVector(q1)) ≈ ∇composition1(q2,q1)
@test ForwardDiff.jacobian(q->SVector(UnitQuaternion{VectorPart}(q)*q1),SVector(q2)) ≈ ∇composition2(q2,q1)
# @btime ForwardDiff.jacobian(q->SVector($u2*UnitQuaternion(q)),SVector($u1))
# @btime ∇composition1($u2,$u1)
# @btime ForwardDiff.jacobian(q->SVector(UnitQuaternion(q)*$u1),SVector($u2))
# @btime ∇composition2($u2,$u1)

@test Lmult(q) ≈ ∇composition1(q,q2)

ϕ = @SVector zeros(3)
@test TO.∇differential(q) ≈ Lmult(q)*jacobian(VectorPart,ϕ)
@test TO.∇differential(q) ≈ Lmult(q)*jacobian(ExponentialMap,ϕ)
@test TO.∇differential(q) ≈ Lmult(q)*jacobian(CayleyMap,ϕ)
@test TO.∇differential(q) ≈ Lmult(q)*jacobian(MRPMap,ϕ)

############################################################################################
#                          ROLL, PITCH, YAW EULER ANGLES
############################################################################################

e0 = @SVector [deg2rad(45), deg2rad(60), deg2rad(20)]
e1 = RPY(e0...)
# e2 = RotXYZ(e0...)

# @test rotmat(e1) == RotX(e0[1])*RotY(e0[2])*RotZ(e0[3])
# @test rotmat(e1) ≈ e2

# @test e1*r ≈ e2*r
# @btime $e2*$r
# @btime $e1*$r

R = rotmat(e1)
e1_ = TO.from_rotmat(rotmat(e1))
@test e1_.θ ≈ e1.θ
@test e1_.ψ ≈ e1.ψ
@test e1_.ϕ ≈ e1.ϕ

# @test rotmat(e1*e1) ≈ RotXYZ(e2*e2)

e1 = RPY(rand(3)...)
e2 = RPY(rand(3)...)
R = rotmat(e2*e1)

# Test inverses
e1,e2 = rand(RPY), rand(RPY)
e3 = e2*e1
@test e2\e3 ≈ e1
@test e3/e1 ≈ e2
@test r ≈ e3\(e2*e1*r)

@test ∇rotate(e1,r) isa SMatrix{3,3}
@test ∇composition1(e2,e1) isa SMatrix{3,3}
@test ∇composition2(e2,e1) isa SMatrix{3,3}

# sϕ1,cϕ1 = sincos(e1.ϕ)
# sθ1,cθ1 = sincos(e1.θ)
# sψ1,cψ1 = sincos(e1.ψ)
#
# sϕ2,cϕ2 = sincos(e2.ϕ)
# sθ2,cθ2 = sincos(e2.θ)
# sψ2,cψ2 = sincos(e2.ψ)
#
# R11 = cθ2*cψ2*(cθ1*cψ1) - cθ2*sψ2*(sϕ1*sθ1*cψ1+cϕ1*sψ1) + sθ2*(-cϕ1*sθ1*cψ1+sϕ1*sψ1)
# R11 ≈ R[1,1]
# R12 = cθ2*cψ2*(-cθ1*sψ1) - cθ2*sψ2*(-sϕ1*sθ1*sψ1+cϕ1*cψ1) + sθ2*(cϕ1*sθ1*sψ1+sϕ1*cψ1)
# R12 ≈ R[1,2]
#
# dxϕ1 = -cθ2*sψ2*(cϕ1*sθ1*cψ1-sϕ1*sψ1) + sθ2*( sϕ1*sθ1*cψ1+cϕ1*sψ1)
# dxθ1 = -cθ2*cψ2*(sθ1*cψ1) - cθ2*sψ2*(sϕ1*cθ1*cψ1) + sθ2*(-cϕ1*cθ1*cψ1)
# dxψ1 = -cθ2*cψ2*(cθ1*sψ1) - cθ2*sψ2*(-sϕ1*sθ1*sψ1+cϕ1*cψ1) + sθ2*( cϕ1*sθ1*sψ1+sϕ1*cψ1)
#
# dyϕ1 = -cθ2*sψ2*(-sϕ1*sθ1*sψ1+cϕ1*cψ1) + sθ2*(cϕ1*sθ1*sψ1+sϕ1*cψ1)
#
# dψdϕ1 = -R11/(R11^2 + R12^2)*dϕ1 + (R12)/(R11^2 + R12^2)*dϕ1
#
# vals1 = @SVector [e1.ϕ, e1.θ, e1.ψ]
# R2 = rotmat(e2)
# f(e) = rotmat_to_rpy(R2*rotmat(e...))
# f(vals1)
# ForwardDiff.jacobian(f, vals1)

# @btime ∇rotation1($e2, $e1)
# @btime ∇rotation2($e2, $e1)

############################################################################################
#                              MODIFIED RODRIGUES PARAMETERS
############################################################################################
p1 = MRP(q1)
p2 = MRP(q2)
# @btime $p2*$p1
# @btime $p2*$r

# Test rotations and composition
@test p2*r ≈ q2*r
@test p2*p1*r ≈ q2*q1*r
@test rotmat(p2)*r ≈ p2*r
p3 = p2*p1
@test p2\p3 ≈ p1
@test p3/p1 ≈ p2

@test r ≈ p3\(p2*p1*r)

r = @SVector rand(3)
R = rotmat(q1)
@test R*r ≈ q1*r

p = MRP(p1)
R1 = rotmat(p)
@test R1*r ≈ R*r
@test p*r ≈ R*r

@test UnitQuaternion(p) ≈ q1

# Test composition jacobians
@test ForwardDiff.jacobian(x->SVector(p2*MRP(x)),SVector(p1)) ≈ ∇composition1(p2,p1)
@test ForwardDiff.jacobian(x->SVector(MRP(x)*p1),SVector(p2)) ≈ ∇composition2(p2,p1)
p0 = MRP(0,0,0)
@test TO.∇differential(p2) ≈ ∇composition1(p2,p0)

############################################################################################
#                              RODRIGUES PARAMETERS
############################################################################################

g1 = RodriguesParam(q1)
g2 = RodriguesParam(q2)
@test q2 ≈ -UnitQuaternion(g2)
@test g2 ≈ RodriguesParam(UnitQuaternion(g2))

# Test compostion and rotation
@test g1*r ≈ q1*r
@test g2*r ≈ q2*r
@test g2*g1*r ≈ q2*q1*r
@test rotmat(g2)*r ≈ g2*r

g3 = g2*g1
@test g2\g3 ≈ g1
@test g3/g1 ≈ g2
@test r ≈ g3\(g2*g1*r)
@test r ≈ (g3\g2*g1)*r

# Test Jacobians
@test ForwardDiff.jacobian(g->RodriguesParam(g)*r, SVector(g1)) ≈ ∇rotate(g1, r)

function compose(g2,g1)
    N = (g2+g1 + g2 × g1)
    D = 1/(1-g2'g1)
    return D*N
end
@test ForwardDiff.jacobian(g->compose(SVector(g2),g), SVector(g1)) ≈ ∇composition1(g2,g1)
@test ForwardDiff.jacobian(g->compose(g,SVector(g1)), SVector(g2)) ≈ ∇composition2(g2,g1)

g0 = RodriguesParam{Float64}(0,0,0)
@test ∇composition1(g2, g0) ≈ TO.∇differential(g2)

# Conversions
import TrajectoryOptimization: rotmat_to_quat
Random.seed!(1) # i = 3
q = rand(UnitQuaternion)
@test rotmat_to_quat(rotmat(q)) ≈ q

Random.seed!(2) # i = 4
q = rand(UnitQuaternion)
@test rotmat_to_quat(rotmat(q)) ≈ q

Random.seed!(3) # i = 2
q = rand(UnitQuaternion)
@test rotmat_to_quat(rotmat(q)) ≈ q

Random.seed!(5) # i = 1
q = rand(UnitQuaternion)
@test rotmat_to_quat(rotmat(q)) ≈ q

#
# # Quadrotors
# quad1 = Dynamics.Quadrotor()
# quad2 = Dynamics.Quadrotor2{Quat{VectorPart}}()
# x,u = rand(quad2)
# dynamics(quad1,x,u) ≈ dynamics(quad2,x,u)
# # @btime dynamics($quad1,$x,$u)
# # @btime dynamics($quad2,$x,$u)
# #
# # @btime state_diff($quad1,$x,$(1.1x))
# # @btime state_diff($quad2,$x,$(1.1x))
# # @btime state_diff_jacobian($quad1, $x)
# # @btime state_diff_jacobian($quad2, $x)
#
# a = @SVector rand(3)
# b = @SVector rand(3)
# @test skew(a)*b == cross(a,b)
# skew(b)*a == -skew(a)*b
#
# p = a
# r = b
# f(x) = MRP_to_DCM(x)*r
# @test ForwardDiff.jacobian(f,p) ≈ MRP_rotate_jacobian(p,r)
#
# # Test MRP compositio
# p1 = @SVector rand(3)
# p2 = @SVector rand(3)
# r = @SVector rand(3)
# R1 = MRP_to_DCM(p1)
# R2 = MRP_to_DCM(p2)
# R3 = R2*R1
# p3 = MRP_composition(p2,p1)
# @test MRP_rotate_vec(p1,r) ≈ R1*r
# @test MRP_rotate_vec(p2,r) ≈ R2*r
# @test MRP_rotate_vec(p3,r) ≈ R3*r
#
# ForwardDiff.jacobian(x->MRP_composition(x,p1), p2) ≈ MRP_composition_jacobian_p2(p2,p1)
# ForwardDiff.jacobian(x->MRP_composition(p2,x), p1) ≈ MRP_composition_jacobian_p1(p2,p1)


text:
abstract type AbstractAnalysis end

abstract type AbstractPlotSpecification <: AbstractAnalysis end
abstract type AbstractSpaceTimePlotSpecification <: AbstractPlotSpecification end

# Default that should work for most
function output_name(aps::AbstractPlotSpecification)
	aps.output_name
end


text:
"""
`MaternKernel([ρ=1.0,[ν=1.0]])`
The matern kernel is an isotropic Mercer kernel given by the formula:
```
    κ(x,y) = 2^{1-ν}/Γ(ν)*(√(2ν)‖x-y‖)^ν K_ν(√(2ν)‖x-y‖)
```
For `ν=n+1/2, n=0,1,2,...` it can be simplified and you should instead use [`ExponentialKernel`](@ref) for `n=0`, [`Matern32Kernel`](@ref), for `n=1`, [`Matern52Kernel`](@ref) for `n=2` and [`SqExponentialKernel`](@ref) for `n=∞`.
"""
struct MaternKernel{Tr<:Transform, Tν<:Real} <: Kernel{Tr}
    transform::Tr
    ν::Tν
    function MaternKernel{Tr, Tν}(t::Tr, ν::Tν) where {Tr, Tν}
        @check_args(MaternKernel, ν, ν > zero(Tν), "ν > 0")
        return new{Tr, Tν}(t, ν)
    end
end

MaternKernel(ρ::Real=1.0, ν::Real=1.5) = MaternKernel(ScaleTransform(ρ), ν)

MaternKernel(ρ::AbstractVector{<:Real},ν::Real=1.5) = MaternKernel(ARDTransform(ρ), ν)

MaternKernel(t::Tr, ν::T=1.5) where {Tr<:Transform, T<:Real} = MaternKernel{Tr, T}(t, ν)

params(k::MaternKernel) = (params(transform(k)),k.ν)
opt_params(k::MaternKernel) = (opt_params(transform(k)),k.ν)

@inline kappa(κ::MaternKernel, d::Real) = iszero(d) ? one(d) : exp((1.0-κ.ν)*logtwo-logabsgamma(κ.ν)[1] + κ.ν*log(sqrt(2κ.ν)*d)+log(besselk(κ.ν,sqrt(2κ.ν)*d)))

metric(::MaternKernel) = Euclidean()

"""
`Matern32Kernel([ρ=1.0])`
The matern 3/2 kernel is an isotropic Mercer kernel given by the formula:
```
    κ(x,y) = (1+√(3)ρ‖x-y‖)exp(-√(3)ρ‖x-y‖)
```
"""
struct Matern32Kernel{Tr} <: Kernel{Tr}
    transform::Tr
end

@inline kappa(κ::Matern32Kernel, d::Real) = (1+sqrt(3)*d)*exp(-sqrt(3)*d)

metric(::Matern32Kernel) = Euclidean()

"""
`Matern52Kernel([ρ=1.0])`
The matern 5/2 kernel is an isotropic Mercer kernel given by the formula:
```
    κ(x,y) = (1+√(5)ρ‖x-y‖ + 5ρ²‖x-y‖^2/3)exp(-√(5)ρ‖x-y‖)
```
"""
struct Matern52Kernel{Tr} <: Kernel{Tr}
    transform::Tr
end

@inline kappa(κ::Matern52Kernel, d::Real) = (1+sqrt(5)*d+5*d^2/3)*exp(-sqrt(5)*d)

metric(::Matern52Kernel) = Euclidean()


text:
# some views are in colorchannels.jl
using Colors, FixedPointNumbers, ImageCore, OffsetArrays, Test
using OffsetArrays: IdentityUnitRange

@testset "rawview" begin
    a = map(N0f8, rand(3,5))
    a[2,2] = N0f8(0.5)
    v = rawview(a)
    @test v[2,2] === a[2,2].i
    v[1,3] = 0xff
    @test a[1,3] === N0f8(1)
    v[1,3] = 0x01
    @test a[1,3] === N0f8(1/255)
    s = view(a, 1:2, 1:2)
    v = rawview(s)
    @test v[2,2] === a[2,2].i
    v[2,2] = 0x0f
    @test a[2,2].i == 0x0f
    @test rawview(v) === v
end

@testset "normedview" begin
    a = rand(UInt8, 3, 5)
    a[2,2] = 0x80
    v = normedview(a)
    @test v[2,2] === N0f8(0.5)
    v[1,3] = 1
    @test a[1,3] === 0xff
    v[1,3] = 1/255
    @test a[1,3] === 0x01
    s = view(a, 1:2, 1:2)
    v = normedview(s)
    @test v[2,2] === N0f8(0.5)
    v[2,2] = 15/255
    @test a[2,2] == 0x0f
    @test normedview(v) === v
    @test normedview(N0f8, v) === v
end

@testset "permuteddimsview" begin
    a = [1 3; 2 4]
    v = permuteddimsview(a, (1,2))
    @test v == a
    v = permuteddimsview(a, (2,1))
    @test v == a'
    a = rand(3,7,5)
    v = permuteddimsview(a, (2,3,1))
    @test v == permutedims(a, (2,3,1))
end

@testset "StackedView" begin
    for (A, B, T) = (([1 3;2 4], [-1 -5; -2 -3], Int),
                     ([1 3;2 4], [-1.0 -5.0; -2.0 -3.0], Float64))
        V = @inferred(StackedView(A, B))
        @test eltype(V) == T
        @test size(V) == (2, 2, 2)
        @test axes(V) === (Base.OneTo(2), Base.OneTo(2), Base.OneTo(2))
        @test @inferred(V[1,1,1]) === T(1)
        @test @inferred(V[2,1,1]) === T(-1)
        @test V[1,:,:] == A
        @test V[2,:,:] == B
        @test_throws BoundsError V[0,1,1]
        @test_throws BoundsError V[3,1,1]
        @test_throws BoundsError V[2,0,1]
        @test_throws BoundsError V[2,3,1]
        @test_throws BoundsError V[1,1,0]
        @test_throws BoundsError V[1,1,3]
        V32 = @inferred(StackedView{Float32}(A, B))
        @test eltype(V32) == Float32
        @test V32[1,1,2] == Float32(3)
        V[1,2,2] = 0
        @test A[2,2] == 0
        V[2,1,2] = 11
        @test B[1,2] == 11

        V = @inferred(StackedView(A, zeroarray, B))
        @test eltype(V) == T
        @test size(V) == (3, 2, 2)
        @test axes(V) === (Base.OneTo(3), Base.OneTo(2), Base.OneTo(2))
        @test V[1,:,:] == A
        @test all(V[2,:,:] .== 0)
        @test V[3,:,:] == B
        @test_throws ErrorException V[2,1,1] = 7
        V32 = @inferred(StackedView{Float32}(A, zeroarray, B))
        @test eltype(V32) == Float32
        @test V32[1,1,2] == Float32(3)
    end

    # With mixed grayscale/real arrays
    a, b = Gray{N0f8}[0.1 0.2; 0.3 0.4], [0.5 0.6; 0.7 0.8]
    V = @inferred(StackedView(a, b))
    @test eltype(V) == Float64
    @test V[1,1,1] === Float64(N0f8(0.1))
    @test V[2,1,1] === 0.5
    V = @inferred(StackedView{N0f8}(a, b))
    @test eltype(V) == N0f8
    @test V[1,1,1] === N0f8(0.1)
    @test V[2,1,1] === N0f8(0.5)
    @test b[1,1] === 0.5
    V = @inferred(StackedView(a, zeroarray, b))
    @test eltype(V) == Float64
    V = @inferred(StackedView{N0f8}(a, zeroarray, b))
    @test eltype(V) == N0f8

    # With colorview
    a = [0.1 0.2; 0.3 0.4]
    b = [0.5 0.6; 0.7 0.8]
    z = zeros(2,2)  # because setindex! won't work with zeroarray
    # GrayA
    v = @inferred(colorview(GrayA, a, b))
    @test eltype(v) == GrayA{Float64}
    @test @inferred(v[2,1]) === GrayA(0.3,0.7)
    v = @inferred(colorview(GrayA{N0f8}, a, b))
    @test @inferred(v[2,1]) === GrayA{N0f8}(0.3,0.7)
    v[1,2] = GrayA(0.25, 0.25)
    @test @inferred(v[1,2]) === GrayA{N0f8}(0.25, 0.25)
    v = @inferred(colorview(GrayA{N0f8}, a, zeroarray))
    @test @inferred(v[2,1]) === GrayA{N0f8}(0.3,0)
    @test_throws ErrorException (v[1,2] = GrayA(0.25, 0.25))
    # RGB
    v = @inferred(colorview(RGB{N0f8}, a, zeroarray, b))
    @test @inferred(v[2,1]) === RGB{N0f8}(0.3,0,0.7)
    v = colorview(RGB{N0f8}, a, z, b)
    @test @inferred(v[2,1]) === RGB{N0f8}(0.3,0,0.7)
    v[2,1] = RGB(0,0.9,0)
    @test @inferred(v[2,1]) === RGB{N0f8}(0,0.9,0)
    @test a[2,1] == b[2,1] == 0
    @test z[2,1] == N0f8(0.9)
    # RGBA
    v = @inferred(colorview(RGBA{N0f8}, a, zeroarray, b, b))
    @test @inferred(v[2,2]) === RGBA{N0f8}(0.4,0,0.8,0.8)
    v = colorview(RGBA{N0f8}, a, copy(z), b, copy(z))
    v[2,2] = RGBA(0.75,0.8,0.75,0.8)
    @test @inferred(v[2,2]) === RGBA{N0f8}(0.75,0.8,0.75,0.8)
    @test a[2,2] == b[2,2] == N0f8(0.75)

    @test eltype(zeroarray) == Union{}
    @test_throws ErrorException colorview(RGB{N0f8}, zeroarray, zeroarray, zeroarray)
    @test_throws DimensionMismatch StackedView(rand(2,3), rand(2,5))

    # With padding
    r = reshape([0.1,0.2], 2, 1)
    g = [false true]
    cv = colorview(RGB, paddedviews(0, r, g, zeroarray)...)
    @test cv[1,1] === RGB(0.1,0,0)
    @test cv[2,1] === RGB(0.2,0,0)
    @test cv[1,2] === RGB(0,1.0,0)
    @test cv[2,2] === RGB(0,0.0,0)
    @test axes(cv) === (Base.OneTo(2), Base.OneTo(2))

    a = [0.1 0.2; 0.3 0.4]
    b = OffsetArray([0.1 0.2; 0.3 0.4], 0:1, 2:3)
    @test_throws DimensionMismatch colorview(RGB, a, b, zeroarray)
    cv = colorview(RGB, paddedviews(0, a, b, zeroarray)...)
    @test axes(cv) == (0:2, 1:3)
    @test red.(cv[axes(a)...]) == a
    @test green.(cv[axes(b)...]) == b
    @test parent(copy(cv)) == [RGB(0,0,0)   RGB(0,0.1,0)   RGB(0,0.2,0);
                               RGB(0.1,0,0) RGB(0.2,0.3,0) RGB(0,0.4,0);
                               RGB(0.3,0,0) RGB(0.4,0,0)   RGB(0,0,0)]
    chanv = channelview(cv)
    @test @inferred(axes(chanv)) === (IdentityUnitRange(1:3), IdentityUnitRange(0:2), IdentityUnitRange(1:3))
    @test chanv[1,1,1] == 0.1
    @test chanv[2,1,2] == 0.3

    @test_throws ErrorException paddedviews(0, zeroarray, zeroarray)

    # Just to boost coverage. These methods are necessary to
    # make inference happy.
    @test_throws ErrorException ImageCore._unsafe_getindex(1, (1,1))
    @test_throws ErrorException ImageCore._unsafe_setindex!(1, (1,1), 0)
end

@testset "Multi-component colorview" begin
    r, g, b = rand(Gray{N0f16},5,5), rand(Gray{N0f16},5,5), rand(Gray{N0f16},5,5)
    A = colorview(RGB, r, g, b)
    @test A[1,1] == RGB(r[1,1], g[1,1], b[1,1])
    o = ones(2, 2)
    A = colorview(RGB, o, o, zeroarray)
    @test A[1,1] == RGB(1,1,0)
    @test size(A) == (2,2)
    @test axes(A) == (Base.OneTo(2),Base.OneTo(2))
    A = colorview(RGB, o, zeroarray, o)
    @test A[1,1] == RGB(1,0,1)
    @test size(A) == (2,2)
    @test axes(A) == (Base.OneTo(2),Base.OneTo(2))
    A = colorview(RGB, zeroarray, o, o)
    @test A[1,1] == RGB(0,1,1)
    @test size(A) == (2,2)
    @test axes(A) == (Base.OneTo(2),Base.OneTo(2))
end

nothing


text:
include("./slices_observer.jl")

my_file = 	"./data/2021-07-27/binary_bh_M_1.0_alpha_10.0_Rorbit_6.0_Rcavity_40.0_N_128_ti_0.0_tf_102.0.h5"
my_file = "./data/2021-07-27/testfile8.h5.h5"

slices_observer(my_file)


text:
using FileIO, ImageCore

function multsizes(dimspecs...)
    szs = Dims{length(dimspecs)}[]
    sz = map(first, dimspecs)
    mult = map(s->s[2], dimspecs)
    stop = map(last, dimspecs)
    while true
        push!(szs, sz)
        sz = sz .* mult
        any(sz .> stop) && break
    end
    return szs
end

function generate_random(workdir::AbstractString, ext=".tif"; # ext=Dict(RGB=>".ppm", Gray=>".pgm");
                         sizes = multsizes((32, 2, 2048), (32, 2, 2048)),
                         eltypes = (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, RGB{N0f8}, RGB{N0f16}, RGB{Float32},))
    if !isdir(workdir)
        mkpath(workdir)
    end
    for T in eltypes, sz in sizes
        Tstr = (T <: Gray ? "gray_" : "rgb_")*string(eltype(T))
        if isa(ext, AbstractDict)
            thisext = ext[T <: Gray ? Gray : RGB]
        else
            thisext = ext
        end
        if eltype(T) <: AbstractFloat && thisext ∈ (".pbm", ".pgm", ".ppm")
            @warn "extension $thisext does not support floating point, skipping"
            continue
        end
        szstr = join(sz, '_')
        img = rand(T, sz)
        save(joinpath(workdir, "random_$(Tstr)_$szstr$thisext"), img)
    end
end

generate_bool(sz) = rand(Bool, sz)

function generate_dblcone(sz)
    img = trues(sz)
    middle = map(axes(img)) do ax
        (first(ax) + last(ax)) ÷ 2
    end
    jmiddle = CartesianIndex(Base.tail(middle))
    jsz = Base.tail(size(img))
    for j in CartesianIndices(Base.tail(axes(img)))
        for i in axes(img, 1)
            img[i, j] = abs2((i - middle[1])/size(img, 1)) < sum(abs2, Tuple(j - jmiddle)./jsz)
        end
    end
    return img
end

function generate_spiral(sz)
    img = zeros(Gray{N0f8}, sz)
    middle = map(axes(img)) do ax
        (first(ax) + last(ax)) ÷ 2
    end
    corner = Tuple(first(CartesianIndices(img)))
    orthog = (last(axes(img, 1)), Base.tail(corner)...)
    I1 = CartesianIndex(ntuple(_->1, ndims(img)))
    Iend = last(CartesianIndices(img))
    # Archimedian spiral
    for θ in LinRange(0, 6π, 6*maximum(size(img)))
        r = 0.95*θ/(6π)
        yn, xn = r*cos(θ), r*sin(θ)
        c = CartesianIndex(round.(Int, middle .+ yn .* (corner .- middle) .+ xn .* (orthog .- middle)))
        I = max(CartesianIndex(corner), c-I1):min(Iend, c+I1)
        isempty(I) && break
        imgtmp = 0.5 .+ 0.5 .* rand(size(I)...)
        img[I] = imgtmp
    end
    return img
end

function generate_special(workdir::AbstractString, ext=".tif")
    if !isdir(workdir)
        mkpath(workdir)
    end
    save(joinpath(workdir, "components2d"*ext), generate_bool((2048, 2048)))
    save(joinpath(workdir, "components3d"*ext), generate_bool((512, 512, 50)))
    save(joinpath(workdir, "dblcone2d"*ext), generate_dblcone((2048, 2048)))
    save(joinpath(workdir, "dblcone3d"*ext), generate_dblcone((512, 512, 50)))
    save(joinpath(workdir, "spiral2d"*ext), generate_spiral((2048, 2048)))
    save(joinpath(workdir, "spiral3d"*ext), generate_spiral((2048, 2048, 50)))
end


text:
#struct for saving the state of the system in some half
struct HydroStatus
	rho::Float64 #density
	u::Float64 #velocity
	p::Float64 #pressure
	gamma::Float64 #adiabatic index
	c::Float64 #speed of cound
	A::Float64 #Toro constant (4.8)
	B::Float64 #Toro constant (4.8)

	function HydroStatus(rho, u, p, gamma)
		c = sqrt(gamma * p / rho)
		A = 2. / (gamma + 1.) / rho
		B = (gamma - 1.) / (gamma + 1.) * p
		new(rho, u, p, gamma, c, A, B)
	end
end

#like the HydroStatus but saves only the important parameters - density, pressure and velocity
struct HydroStatusSimple
	rho::Float64 #density
	u::Float64 #velocity
	p::Float64 #pressure
end
HydroStatusSimple(status::HydroStatus) = HydroStatusSimple(status.rho, status.u, status.p)

#different types to dispatch the first guess of p_star calculation
#see "PressureGuess.jl"
abstract type FirstGuessScheme end
struct TwoRarefaction <: FirstGuessScheme end
struct PrimitiveValue <: FirstGuessScheme end
struct TwoShock <: FirstGuessScheme end
struct MeanPressure <: FirstGuessScheme end

#the type of the wave that will be created on some side of the contact discontinuity
abstract type WaveType end
struct Rarefaction <: WaveType end
struct Shock <: WaveType end
struct Vaccum <: WaveType end #for vaccum generated in the middle

#the side in which the calculation is happening
abstract type Side end
struct LeftSide <: Side end
struct RightSide <: Side end

#in some equations there is factor difference between the left and right side calculations
side_factor(::LeftSide) = -1
side_factor(::RightSide) = 1

#choose the some variable based on the side given
choose_side(left, right, ::LeftSide) = left
choose_side(left, right, ::RightSide) = right

max_stars_repo_name: ryarazi/ExactRiemannProblemSolver

max_issues_repo_name: ryarazi/ExactRiemannProblemSolver

max_forks_repo_name: ryarazi/ExactRiemannProblemSolver


text:
include("utility.jl")
# these functions are for updating local counts in sparseSCVB0
function  burnInUpdate(sparseVBparams,documentTopicCounts,documents,burnInPerDoc,docInd,docInd_mini,docLength,tau2,kappa2,sparseCounts,topicMem,sampleSize,sparsity);
    Beta = sparseVBparams.Beta;
    Alpha = sparseVBparams.Alpha;
    BetaSum = sparseVBparams.BetaSum;
    numTopics = sparseVBparams.numTopics;
    wordTopicCounts=sparseVBparams.wordTopicCounts;
    topicCounts=sparseVBparams.topicCounts;
    numCachedSamples=sparseVBparams.numCachedSamples;
    cached_qw=sparseVBparams.cached_qw;
    cachedSamples=sparseVBparams.cachedSamples;
    sampleIndx=sparseVBparams.sampleIndx;
    cached_Qw=sparseVBparams.cached_Qw;

    numDistinctTokens = length(documents[docInd].tokenIndex);
    stepSizeCounter = 1;
    stepSize2 = (stepSizeCounter + tau2)^-kappa2;
    right_step = 0.0; # initializing right step size for local expected counts to make it soft local scope

    probs = Array{Float64}(numTopics);
    for burn = 1:burnInPerDoc
        for tokenInd = 1:numDistinctTokens
            word = documents[docInd].tokenIndex[tokenInd];
            sum_probs = 0.0; # total of sparse bucket distribution
            for topic = 1:numTopics
                probs[topic] = (wordTopicCounts[word,topic] + Beta[word]) .*
                           (documentTopicCounts[docInd_mini,topic] + Alpha[topic]) ./
                           (topicCounts[topic] + BetaSum);
                sum_probs += probs[topic];
            end

            # effective stepsize for local expected counts
            tokenCountInDoc = documents[docInd].tokenCount[tokenInd];
            left_step = (1 - stepSize2)^tokenCountInDoc;
            right_step = 1 - (1-stepSize2)^tokenCountInDoc;
            # update the local expected counts

            for topic = 1:numTopics
                probs[topic]/=sum_probs;
                documentTopicCounts[docInd_mini,topic] = left_step .* documentTopicCounts[docInd_mini,topic] + right_step .* docLength .* probs[topic];
            end

            stepSizeCounter +=tokenCountInDoc;
            stepSize2 = (stepSizeCounter + (burn-1).*docLength + tau2)^-kappa2;
        end
        #documentTopicCounts,sparseCounts = sparsificationHeuristic(right_step,docLength,docInd,sparseCounts,documentTopicCounts);
    end
    documentTopicCounts,sparseCounts = sparsificationHeuristic(right_step,docLength,docInd_mini,sparseCounts,documentTopicCounts,sparsity);
    return sparseVBparams,documentTopicCounts,sparseCounts,topicMem;
end

function  mainLoopUpdate(sparseVBparams,documentTopicCounts,documents,burnInPerDoc,docInd,docInd_mini,docLength,tau2,kappa2,sparseCounts,
                    topicMem,wordTopicCounts_hat,topicCounts_hat,sampleSize);

    Beta = sparseVBparams.Beta;
    Alpha = sparseVBparams.Alpha;
    BetaSum = sparseVBparams.BetaSum;
    numTopics = sparseVBparams.numTopics;
    wordTopicCounts=sparseVBparams.wordTopicCounts;
    topicCounts=sparseVBparams.topicCounts;
    numCachedSamples=sparseVBparams.numCachedSamples;
    cached_qw=sparseVBparams.cached_qw;
    cachedSamples=sparseVBparams.cachedSamples;
    sampleIndx=sparseVBparams.sampleIndx;
    cached_Qw=sparseVBparams.cached_Qw;

    numDistinctTokens = length(documents[docInd].tokenIndex);
    stepSizeCounter = 1;
    stepSize2 = (stepSizeCounter + burnInPerDoc.*docLength + tau2)^-kappa2;
    right_step = 0.0; # initializing right step size for local expected counts to make it soft local scope

    for tokenInd = 1:numDistinctTokens
        word = documents[docInd].tokenIndex[tokenInd];

        pseudoProb,sparseCounts,topicMem,cachedSamples,cached_qw,sampleIndx,cached_Qw = drawTopic(sparseCounts,docInd_mini,documentTopicCounts,wordTopicCounts,Beta,
                              topicCounts,BetaSum,Alpha,cached_qw,topicMem,tokenInd,word,sampleIndx,numCachedSamples,cachedSamples,numTopics,sampleSize,cached_Qw)
                              # effective stepsize for local expected counts

        # effective stepsize for local expected counts
        tokenCountInDoc = documents[docInd].tokenCount[tokenInd];
        left_step = (1 - stepSize2)^tokenCountInDoc;
        right_step = 1 - (1-stepSize2)^tokenCountInDoc;
        # update the local expected counts
        for topic = 1:length(sparseCounts[docInd_mini].nonzeroDocIndex)
            documentTopicCounts[docInd_mini,sparseCounts[docInd_mini].nonzeroDocIndex[topic]] = left_step .* documentTopicCounts[docInd_mini,sparseCounts[docInd_mini].nonzeroDocIndex[topic]];
        end
        # calculating the constant update part of local and estimate expected counts
        localUpdate = (right_step.*docLength)/sampleSize;
        estimateUpdate = tokenCountInDoc/sampleSize;
        for topic = 1:sampleSize
            documentTopicCounts[docInd_mini,pseudoProb[topic]]+=localUpdate;
            # update the estimate of global expected counts
            wordTopicCounts_hat[word,pseudoProb[topic]]+=estimateUpdate;
            topicCounts_hat[pseudoProb[topic]]+=estimateUpdate;
        end
        stepSizeCounter+=tokenCountInDoc;
        stepSize2 = (stepSizeCounter + burnInPerDoc.*docLength + tau2)^-kappa2;
    end
    # sparsification heuristic by a clever threshold
    #documentTopicCounts,sparseCounts = sparsificationHeuristic(right_step,docLength,docInd,sparseCounts,documentTopicCounts);

    return sparseVBparams,documentTopicCounts,wordTopicCounts_hat,topicCounts_hat,sparseCounts,topicMem;
end

function sparsificationHeuristic(right_step,docLength,docInd_mini,sparseCounts,documentTopicCounts,sparsity)
    SparseTh = sparsity*right_step.*docLength; #sparsification threshold
    trackZeroInd = Array{Int64}(0); # to keep track of zero index for sparse counts
    for zeroInd = 1:length(sparseCounts[docInd_mini].nonzeroDocIndex)
        if documentTopicCounts[docInd_mini,sparseCounts[docInd_mini].nonzeroDocIndex[zeroInd]] < SparseTh
            documentTopicCounts[docInd_mini,sparseCounts[docInd_mini].nonzeroDocIndex[zeroInd]]=0.0;
            push!(trackZeroInd,zeroInd);
        end
    end
    # deleteing zero valued index for nonzeroDocIndex field of sparse counts
    deleteat!(sparseCounts[docInd_mini].nonzeroDocIndex,trackZeroInd);

    return documentTopicCounts,sparseCounts;
end

function drawTopic(sparseCounts,docInd_mini,documentTopicCounts,wordTopicCounts,Beta,
                      topicCounts,BetaSum,Alpha,cached_qw,topicMem,tokenInd,word,
                      sampleIndx,numCachedSamples,cachedSamples,numTopics,sampleSize,cached_Qw)
    nonzeroLength = length(sparseCounts[docInd_mini].nonzeroDocIndex);
    # compute sparse bucket distribution
    Pw = 0.0; # total of sparse bucket distribution
    pw = Array{Float64}(nonzeroLength);
    for topic = 1:nonzeroLength
        pw[topic] = documentTopicCounts[docInd_mini,sparseCounts[docInd_mini].nonzeroDocIndex[topic]] .*
                (wordTopicCounts[word,sparseCounts[docInd_mini].nonzeroDocIndex[topic]] + Beta[word]) ./
                (topicCounts[sparseCounts[docInd_mini].nonzeroDocIndex[topic]] + BetaSum);
        Pw = Pw + pw[topic];
    end
    Qw = cached_Qw[word]; # total of dense bucket distribution
    # generating pseudo-variational distribution
    pseudoProb = Array{Int64}(0);
    for s = 1:sampleSize
        oldTopic = topicMem[docInd_mini].Assigned[tokenInd];
    # choose one of the buckets (sparse & dense) to draw a sample
        if rand() < Pw./(Pw+Qw)
            indxDiscrete = sampleFromDiscrete(pw,Pw,nonzeroLength);
            newTopic = sparseCounts[docInd_mini].nonzeroDocIndex[indxDiscrete];
        else
            if sampleIndx[word]>numCachedSamples
                cached_Qw[word] = 0.0;
                for topic = 1:numTopics
                    cached_qw[word,topic] = Alpha[topic] .* (wordTopicCounts[word,topic] + Beta[word]) ./ (topicCounts[topic] + BetaSum);
                    cached_Qw[word] += cached_qw[word,topic];
                end
                aliasTable = generateAlias(cached_qw[word,:],numTopics);
                for sample = 1:numCachedSamples
                    cachedSamples[word,sample] = sampleAlias(aliasTable,numTopics);
                end
                sampleIndx[word] = 1;
            end
            newTopic = cachedSamples[word,sampleIndx[word]];
            sampleIndx[word]+=1;
        end
        if oldTopic != newTopic
    # accept the new sample with a probability
            acceptanceRatio = ((documentTopicCounts[docInd_mini,newTopic]+Alpha[newTopic])./(documentTopicCounts[docInd_mini,oldTopic]+Alpha[oldTopic])).*
              ((wordTopicCounts[word,newTopic]+Beta[word])./(wordTopicCounts[word,oldTopic]+Beta[word])).*
              ((topicCounts[oldTopic]+BetaSum)./(topicCounts[newTopic]+BetaSum)).*
              (((documentTopicCounts[docInd_mini,oldTopic].*(wordTopicCounts[word,oldTopic]+Beta[word])./(topicCounts[oldTopic]+BetaSum))*Pw + Qw*cached_qw[word,oldTopic])./
              ((documentTopicCounts[docInd_mini,newTopic].*(wordTopicCounts[word,newTopic]+Beta[word])./(topicCounts[newTopic]+BetaSum))*Pw + Qw*cached_qw[word,newTopic]));
            if rand()<acceptanceRatio
                topicMem[docInd_mini].Assigned[tokenInd] = newTopic; # re-store updated topic
            else
                newTopic = oldTopic;
            end
        end
        #if (newTopic in sparseCounts[docInd].nonzeroDocIndex) == false
        if documentTopicCounts[docInd_mini,newTopic] == 0.0
            if (newTopic in sparseCounts[docInd_mini].nonzeroDocIndex) == false
                push!(sparseCounts[docInd_mini].nonzeroDocIndex,newTopic);
            end
        end
        push!(pseudoProb,newTopic);
    end
    return pseudoProb,sparseCounts,topicMem,cachedSamples,cached_qw,sampleIndx,cached_Qw;
end


text:
cd(@__DIR__)
using Pkg
Pkg.activate(".")
##

include("sp_model.jl")
include("plot_utils.jl")
#include("results_handling.jl")
include("evaluation_utils.jl")
##

using DataFrames
using CSV

offset = 6531
timesteps = 1:168

pv = CSV.read("timeseries/basic_example.csv", DataFrame)[timesteps .+ offset, 3]
wind = CSV.read("timeseries/basic_example.csv", DataFrame)[timesteps .+ offset, 4]
demand = CSV.read("timeseries/basic_example.csv", DataFrame)[timesteps .+ offset, 2]

##

pars = copy(default_es_pars)

##

using Statistics
average_hourly_demand = mean(demand)

pars[:recovery_time] = 24
pars[:c_storage] = 100.
pars[:c_pv] = 300.
pars[:c_wind] = 800.
pars[:c_sto_op] = 0.00001

heatdemand = copy(demand)./300.
heatdemand0 = zeros(length(demand))
es = define_energy_system(pv, wind, demand, heatdemand; p = pars, strict_flex = true)

##
n = 100
F_max = average_hourly_demand * 0.1 # Have ancillary services equal to 10% of our typical demand
t_max = length(pv) - es.parameters[2].defaults[:recovery_time]
scens = simple_flex_sampler(n, F_max, t_max)

##

using Cbc

sp = instantiate(es, scens, optimizer = Cbc.Optimizer)

##

optimize!(sp)

##

println("Objective value: $(objective_value(sp))")

println("PV: $(round(value.(sp[1, :u_pv]); digits = 2)) -- Wind: $(round(value.(sp[1, :u_wind]); digits = 2)) -- Battery: $(round(value.(sp[1, :u_storage]); digits = 2))")
println("Heat storage: $(round(value.(sp[1, :u_heat_storage]); digits = 2))")

# for s in 1:minimum((n, 10))
#     println("Second stage objective value in scenario $s: $(objective_value(sp, s))")
# end

##

plot_results(sp, pv, wind, demand, hd = heatdemand, s = 8, stage_1 = [:gci, :gco], stage_2 = [:gci2, :gco2])

plot_difference(sp, s=8)

##
cost_pos, pot_pos, cost_neg, pot_neg = test_decision(sp, 1:t_max)

plot_flexibility(1:t_max, cost_pos, pot_pos, cost_neg, pot_neg, objective_value(sp))

##
sp0 = instantiate(es, no_flex_pseudo_sampler(), optimizer = Cbc.Optimizer)

optimize!(sp0)

plot_results(sp0, pv, wind, demand, hd = heatdemand, s = 1, stage_1 = [:gci, :gco], stage_2 = [:gci2, :gco2])

cost_pos0, pot_pos0, cost_neg0, pot_neg0 = test_decision(sp0, 1:t_max)

plot_flexibility(1:t_max, cost_pos0, pot_pos0, cost_neg0, pot_neg0, objective_value(sp0))

##
cost_pos, pot_pos, cost_neg, pot_neg = test_decision(sp, 1:t_max)

plot_flexibility(1:t_max, cost_pos, pot_pos, cost_neg, pot_neg, objective_value(sp))

#label_av = ["+ flexibility in scenario-aware system", "+ flexibility in scenario-unaware system"]

plt_av = plot();
flexibility_availability(plt_av, pot_pos)
flexibility_availability(plt_av, pot_pos0)
flexibility_availability(plt_av, pot_neg)
flexibility_availability(plt_av, pot_neg0)
##


text:
type PsPotInfo
  nbeta::Int
  NLmax::Int
end


text:
# OpenCL.Device

immutable Device <: CLObject
    id :: CL_device_id
end

Base.pointer(d::Device) = d.id

function Base.show(io::IO, d::Device)
    strip_extra_whitespace = r"\s+"
    device_name = replace(d[:name], strip_extra_whitespace, " ")
    platform_name = replace(d[:platform][:name], strip_extra_whitespace, " ")
    ptr_val = convert(UInt, Base.pointer(d))
    ptr_address = "0x$(hex(ptr_val, Sys.WORD_SIZE>>2))"
    print(io, "OpenCL.Device($device_name on $platform_name @$ptr_address)")
end

Base.getindex(d::Device, dinfo::Symbol) = info(d, dinfo)

macro int_info(func, cl_device_info, return_type)
    quote
        function $(esc(func))(d::Device)
            result = Ref{$return_type}()
            @check api.clGetDeviceInfo(d.id, $cl_device_info,
                                       sizeof($return_type), result, C_NULL)
            return result[]
        end
    end
end

let profile(d::Device) = begin
        size = Ref{Csize_t}()
        @check api.clGetDeviceInfo(d.id, CL_DEVICE_PROFILE, 0, C_NULL, size)
        result = Vector{CL_char}(size[])
        @check api.clGetDeviceInfo(d.id, CL_DEVICE_PROFILE, size[], result, C_NULL)
        bs = CLString(result)
        return bs
    end

    version(d::Device) = begin
        size = Ref{Csize_t}()
        @check api.clGetDeviceInfo(d.id, CL_DEVICE_VERSION, 0, C_NULL, size)
        result = Vector{CL_char}(size[])
        @check api.clGetDeviceInfo(d.id, CL_DEVICE_VERSION, size[], result, C_NULL)
        bs = CLString(result)
        return bs
    end

    driver_version(d::Device) = begin
        size = Ref{Csize_t}()
        @check api.clGetDeviceInfo(d.id, CL_DRIVER_VERSION, 0, C_NULL, size)
        result = Vector{CL_char}(size[])
        @check api.clGetDeviceInfo(d.id, CL_DRIVER_VERSION, size[], result, C_NULL)
        bs = CLString(result)
        return string(replace(bs, r"\s+", " "))
    end

    extensions(d::Device) = begin
        size = Ref{Csize_t}()
        @check api.clGetDeviceInfo(d.id, CL_DEVICE_EXTENSIONS, 0, C_NULL, size)
        result = Vector{CL_char}(size[])
        @check api.clGetDeviceInfo(d.id, CL_DEVICE_EXTENSIONS, size[], result, C_NULL)
        bs = CLString(result)
        return String[string(s) for s in split(bs)]
    end

    platform(d::Device) = begin
        result = Ref{CL_platform_id}()
        @check api.clGetDeviceInfo(d.id, CL_DEVICE_PLATFORM,
                                   sizeof(CL_platform_id), result, C_NULL)
        return Platform(result[])
    end

    name(d::Device) = begin
        size = Ref{Csize_t}()
        @check api.clGetDeviceInfo(d.id, CL_DEVICE_NAME, 0, C_NULL, size)
        result = Vector{CL_char}(size[])
        @check api.clGetDeviceInfo(d.id, CL_DEVICE_NAME,
                                   size[] * sizeof(CL_char), result, C_NULL)
        n = CLString(result)
        return string(replace(n, r"\s+", " "))
    end

    device_type(d::Device) = begin
        result = Ref{CL_device_type}()
        @check api.clGetDeviceInfo(d.id, CL_DEVICE_TYPE,
                                   sizeof(CL_device_type), result, C_NULL)
        result = result[]
        if result == CL_DEVICE_TYPE_GPU
            return :gpu
        elseif result == CL_DEVICE_TYPE_CPU
            return :cpu
        elseif result == CL_DEVICE_TYPE_ACCELERATOR
            return :accelerator
        elseif result == CL_DEVICE_TYPE_CUSTOM
            return :custom
        else
            return :unknown
        end
    end

    has_image_support(d::Device) = begin
        has_support = Ref{CL_bool}(CL_FALSE)
        @check api.clGetDeviceInfo(d.id, CL_DEVICE_IMAGE_SUPPORT,
                                   sizeof(CL_bool), has_support, C_NULL)
        return has_support[] == CL_TRUE
    end

    @int_info(queue_properties, CL_DEVICE_QUEUE_PROPERTIES, CL_command_queue_properties)

    has_queue_out_of_order_exec(d::Device) =
        (queue_properties(d) & CL_QUEUE_OUT_OF_ORDER_EXEC_MODE_ENABLE) != 0

    has_queue_profiling(d::Device) =
        (queue_properties(d) & CL_QUEUE_PROFILING_ENABLE) != 0

    has_native_kernel(d::Device) = begin
        result = Ref{CL_device_exec_capabilities}()
        @check api.clGetDeviceInfo(d.id, CL_DEVICE_EXECUTION_CAPABILITIES,
                                   sizeof(CL_device_exec_capabilities), result, C_NULL)
        return (result[] & CL_EXEC_NATIVE_KERNEL) != 0
    end

    @int_info(vendor_id,             CL_DEVICE_VENDOR_ID,                CL_uint)
    @int_info(max_compute_units,     CL_DEVICE_MAX_COMPUTE_UNITS,        CL_uint)
    @int_info(max_work_item_dims,    CL_DEVICE_MAX_WORK_ITEM_DIMENSIONS, CL_uint)
    @int_info(max_clock_frequency,   CL_DEVICE_MAX_CLOCK_FREQUENCY,      CL_uint)
    @int_info(address_bits,          CL_DEVICE_ADDRESS_BITS,             CL_uint)
    @int_info(max_read_image_args,   CL_DEVICE_MAX_READ_IMAGE_ARGS,      CL_uint)
    @int_info(max_write_image_args,  CL_DEVICE_MAX_WRITE_IMAGE_ARGS,     CL_uint)
    @int_info(global_mem_size,       CL_DEVICE_GLOBAL_MEM_SIZE,          CL_ulong)
    @int_info(max_mem_alloc_size,    CL_DEVICE_MAX_MEM_ALLOC_SIZE,       CL_ulong)
    @int_info(max_const_buffer_size, CL_DEVICE_MAX_CONSTANT_BUFFER_SIZE, CL_ulong)
    @int_info(local_mem_size,        CL_DEVICE_LOCAL_MEM_SIZE,           CL_ulong)

    has_local_mem(d::Device) = begin
        result = Ref{CL_device_local_mem_type}()
        @check api.clGetDeviceInfo(d.id, CL_DEVICE_LOCAL_MEM_TYPE,
                                   sizeof(CL_device_local_mem_type), result, C_NULL)
        return result[] == CL_LOCAL
    end

    host_unified_memory(d::Device) = begin
        result = Ref{CL_bool}(CL_FALSE)
        @check api.clGetDeviceInfo(d.id, CL_DEVICE_HOST_UNIFIED_MEMORY,
                                   sizeof(CL_bool), result, C_NULL)
        return result[] == CL_TRUE
    end

    available(d::Device) = begin
        result = Ref{CL_bool}(CL_FALSE)
        @check api.clGetDeviceInfo(d.id, CL_DEVICE_AVAILABLE,
                                   sizeof(CL_bool), result, C_NULL)
        return result[] == CL_TRUE
    end

    compiler_available(d::Device) = begin
        result = Ref{CL_bool}(CL_FALSE)
        @check api.clGetDeviceInfo(d.id, CL_DEVICE_COMPILER_AVAILABLE,
                                   sizeof(CL_bool), result, C_NULL)
        return result[] == CL_TRUE
    end

    max_work_item_size(d::Device) = begin
        dims = Ref{CL_uint}()
        @check api.clGetDeviceInfo(d.id, CL_DEVICE_MAX_WORK_ITEM_DIMENSIONS,
                                   sizeof(CL_uint), dims, C_NULL)
        result = Vector{Csize_t}(dims[])
        @check api.clGetDeviceInfo(d.id, CL_DEVICE_MAX_WORK_ITEM_SIZES,
                                   sizeof(Csize_t) * dims[], result, C_NULL)
        return tuple([Int(r) for r in result]...)
    end

    @int_info(max_work_group_size, CL_DEVICE_MAX_WORK_GROUP_SIZE, Csize_t)
    @int_info(max_parameter_size, CL_DEVICE_MAX_PARAMETER_SIZE,  Csize_t)
    @int_info(profiling_timer_resolution, CL_DEVICE_PROFILING_TIMER_RESOLUTION, Csize_t)

    max_image2d_shape(d::Device) = begin
        width  = Ref{Csize_t}()
        height = Ref{Csize_t}()
        @check api.clGetDeviceInfo(d.id, CL_DEVICE_IMAGE2D_MAX_WIDTH,
                                   sizeof(Csize_t), width,  C_NULL)
        @check api.clGetDeviceInfo(d.id, CL_DEVICE_IMAGE2D_MAX_HEIGHT,
                                   sizeof(Csize_t), height, C_NULL)
        return (width[], height[])
    end

    max_image3d_shape(d::Device) = begin
        width  = Ref{Csize_t}()
        height = Ref{Csize_t}()
        depth =  Ref{Csize_t}()
        @check api.clGetDeviceInfo(d.id, CL_DEVICE_IMAGE3D_MAX_WIDTH,
                                   sizeof(Csize_t), width, C_NULL)
        @check api.clGetDeviceInfo(d.id, CL_DEVICE_IMAGE3D_MAX_HEIGHT,
                                   sizeof(Csize_t), height, C_NULL)
        @check api.clGetDeviceInfo(d.id, CL_DEVICE_IMAGE3D_MAX_DEPTH,
                                   sizeof(Csize_t), depth, C_NULL)
        return (width[], height[], depth[])
    end

    info_map = Dict{Symbol, Function}(
        :driver_version => driver_version,
        :version => version,
        :profile => profile,
        :extensions => extensions,
        :platform => platform,
        :name => name,
        :device_type => device_type,
        :has_image_support => has_image_support,
        :queue_properties => queue_properties,
        :has_queue_out_of_order_exec => has_queue_out_of_order_exec,
        :has_queue_profiling => has_queue_profiling,
        :has_native_kernel => has_native_kernel,
        :vendor_id => vendor_id,
        :max_compute_units => max_compute_units,
        :max_work_item_size => max_work_item_size,
        :max_clock_frequency => max_clock_frequency,
        :address_bits => address_bits,
        :max_read_image_args => max_read_image_args,
        :max_write_image_args => max_write_image_args,
        :global_mem_size => global_mem_size,
        :max_mem_alloc_size => max_mem_alloc_size,
        :max_const_buffer_size => max_const_buffer_size,
        :local_mem_size => local_mem_size,
        :has_local_mem => has_local_mem,
        :host_unified_memory => host_unified_memory,
        :available => available,
        :compiler_available => compiler_available,
        :max_work_group_size => max_work_group_size,
        :max_work_item_dims => max_work_item_dims,
        :max_parameter_size => max_parameter_size,
        :profiling_timer_resolution => profiling_timer_resolution,
        :max_image2d_shape => max_image2d_shape,
        :max_image3d_shape => max_image3d_shape
    )

    function info(d::Device, s::Symbol)
        try
            func = info_map[s]
            func(d)
        catch err
            if isa(err, KeyError)
                throw(ArgumentError("OpenCL.Device has no info for: $s"))
            else
                throw(err)
            end
        end
    end
end

function cl_device_type(dtype::Symbol)
    if dtype == :all
        cl_dtype = CL_DEVICE_TYPE_ALL
    elseif dtype == :cpu
        cl_dtype = CL_DEVICE_TYPE_CPU
    elseif dtype == :gpu
        cl_dtype = CL_DEVICE_TYPE_GPU
    elseif dtype == :accelerator
        cl_dtype = CL_DEVICE_TYPE_ACCELERATOR
    elseif dtype == :custom
        cl_dtype = CL_DEVICE_TYPE_CUSTOM
    elseif dtype == :default
        cl_dtype = CL_DEVICE_TYPE_DEFAULT
    else
        throw(ArgumentError("Unknown device type: $dtype"))
    end
    return cl_dtype
end


text:
using Test
using ADBCVUMPS
using ADBCVUMPS:indexperm_symmetrize
using LinearAlgebra

@testset "$(Ni)x$(Nj) bcipeps" for Nj = 1:3, Ni = 1:3
    bcipeps = Array{Array,2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        bcipeps[i,j] = rand(3,3,3,3,2)
    end
    @test SquareBCIPEPS(bcipeps) isa BCIPEPS{}
    bcipeps[1,1] = rand(3,3,4,3,2)
    @test_throws DimensionMismatch SquareBCIPEPS(bcipeps)

    bcipeps[1,1] = rand(3,3,3,3,2)
    bcipeps = indexperm_symmetrize(SquareBCIPEPS(bcipeps))
    for j = 1:Nj, i = 1:Ni
        @test bcipeps.bulk[i,j] ==
            permutedims(bcipeps.bulk[i,j], (1,4,3,2,5)) == # up-down
            permutedims(bcipeps.bulk[i,j], (3,2,1,4,5)) == # left-right
            permutedims(bcipeps.bulk[i,j], (2,1,4,3,5)) == # diagonal
            permutedims(bcipeps.bulk[i,j], (4,3,2,1,5))    # rotation
    end
end


text:
#=
    Port of SPPL model
    https://github.com/probcomp/sppl/blob/master/examples/fairness-hiring-model-1.ipynb

    The SPPL model gives the following exact probabilities:
        p_hire_given_minority:0.007131626828051439
        p_hire_given_majority:0.01945024229170891
        min/maj: 0.3666600508668959

    This Turing.jl model:
    Summary Statistics
            parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
                Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

           ethnicity    1.8531    0.3540     0.0035    0.0046   6492.4337    1.0000     1198.5294
    years_experience    7.5026    1.9193     0.0192    0.0230   6416.1713    0.9999     1184.4510
        college_rank   20.6405    7.0939     0.0709    0.0827   6775.2100    1.0000     1250.7310
                hire    0.0199    0.1397     0.0014    0.0018   7133.8377    1.0002     1316.9351

    Quantiles
            parameters      2.5%     25.0%     50.0%     75.0%     97.5%
                Symbol   Float64   Float64   Float64   Float64   Float64

           ethnicity    1.0000    2.0000    2.0000    2.0000    2.0000
    years_experience    4.0000    6.0000    8.0000    9.0000   11.0000
        college_rank    5.9163   16.9625   20.4955   24.4271   35.7348
                hire    0.0000    0.0000    0.0000    0.0000    0.0000

    Distributions of variable ethnicity
    majority   =>    8531  (0.853100)
    minority   =>    1469  (0.146900)
    ethnicity: 1 minority, 2: majority:
    Distributions of variable (num:0)
    [2, 0]  =>      8345 (0.8345)
    [1, 0]  =>      1456 (0.1456)
    [2, 1]  =>      186 (0.0186)
    [1, 1]  =>      13 (0.0013)
    hire | minority: 0.008849557522123894
    hire | majority: 0.021802836713163756
    min / maj: 0.4058901893614997

=#

using Turing
include("jl_utils.jl")

#=
SPPL model:
# Prior on ethnicity, experience, and college rank.
ethnicity ~= choice({'minority': .15, 'majority': .85})
years_experience ~= binom(n=15, p=.5)
if (ethnicity == 'minority'):
    college_rank ~= dlaplace(loc=25, a=1/5)
else:
    college_rank ~= dlaplace(loc=20, a=1/5)

# Top 50 colleges and at most 20 years of experience.
condition((0 <= college_rank) <= 50)
condition((0 <= years_experience) <= 20)

# Hiring decision.
if college_rank <= 5:
    hire ~= atomic(loc=1)
else:
    switch (years_experience) cases (years in range(0, 20)):
        if ((years - 5) > college_rank):
            hire ~= atomic(loc=1)
        else:
            hire ~= atomic(loc=0)
=#
@model function fairness_model1()
    minority = 1
    majority = 2

    ethnicity ~ Categorical([0.15,0.85])
    years_experience ~ Binomial(15,0.5)
    college_rank ~ ethnicity == minority ? Laplace(25,5) : Laplace(20,5)

    # Top 50 colleges and at most 20 years of experience.
    true ~ Dirac(college_rank <= 50)
    true ~ Dirac(years_experience <= 20)

    # Hiring decision.
    if college_rank <= 5
        hire ~ Dirac(1)
    else
        # hakank: I'm not sure that I understand this: Why years_experience related to college rank?
        # Some seconds later: Ah, that's from the underlying decision tree!
        if ((years_experience - 5) > college_rank)
            hire ~ Dirac(1)
        else
            hire ~ Dirac(0)
        end
    end

    return [ethnicity, hire]

end

model = fairness_model1()

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

display(chns)

show_var_dist_pct(chns,:ethnicity,["minority","majority"])
eth = group(chns,:ethnicity).value.data |> make_hash

println("ethnicity: 1 minority, 2: majority:")
chains_params = Turing.MCMCChains.get_sections(chns, :parameters)
genq = generated_quantities(model, chains_params)
show_var_dist_pct(genq)

hhh = genq[:,1] |> make_hash
hire_min = get(hhh,[1,1],0) / eth[1]
hire_maj = get(hhh,[2,1],0) / eth[2]
println("hire | minority: $hire_min")
println("hire | majority: $hire_maj")
println("min / maj: $(hire_maj > 0 ? hire_min / hire_maj : 0)")

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
const scales = Observable{Any}(Dict())

abstract type AbstractScale end

struct DiscreteScale
    options::Union{AbstractVector, Nothing}
    function DiscreteScale(x)
        options = x isa AbstractVector ? x : nothing
        return new(options)
    end
end
DiscreteScale(v::AbstractObservable) = DiscreteScale(v[])

rank(c) = levelcode(c)
rank(n::Integer) = n
rank(v::AbstractArray) = rank(only(v))

function get_attr(d::DiscreteScale, value, unique)
    v = strip_name(value)
    # Could be tricky in case of many datasets
    n = sum(t -> !isless(rank(v), rank(t)), unique)
    scale = d.options
    return scale === nothing ? n : scale[mod1(n, length(scale))]
end

struct ContinuousScale
    extrema::Union{NTuple{2, Float64}, Nothing}
    function ContinuousScale(x)
        extrema = x isa NTuple{2, Number} ? map(Float64, x) : nothing
        return new(extrema)
    end
end

function get_attr(c::ContinuousScale, value, extrema)
    scale = c.extrema
    scale === nothing && return value
    min, max = extrema
    smin, smax = scale
    @. smin + value * (smax - smin) / (max - min)
end

function get_extrema(specs::AbstractVector{<:AbstractElement})
    d = Dict{Symbol, NTuple{2, Float64}}()
    for spec in specs
        style = spec.style
        for (k, val) in pairs(style.value)
            a, b = get(d, k, (Inf, - Inf))
            a′, b′ = val isa AbstractVector{<:Number} ? extrema(val) : (Inf, -Inf)
            d[k] = (min(a, a′), max(b, b′))
        end
    end
    return d
end

function get_unique(specs::AbstractVector{<:AbstractElement})
    d = Dict{Symbol, Set{Any}}()
    for spec in specs
        pkeys = spec.pkeys
        for (k, val) in pairs(pkeys)
            grp = strip_name(val)
            v = get!(d, k, Set())
            push!(v, grp)
        end
    end
    return d
end

function computescales(specs)
    unique_dict = get_unique(specs)
    extrema_dict = get_extrema(specs)
    [computescales(spec, unique_dict, extrema_dict) for spec in specs]
end

function computescales(spec::Spec{T}, unique_dict, extrema_dict) where {T}
    # temporary! Should have a sensible default scales set,
    # both discrete and continuous
    scales[] = (; AbstractPlotting.current_default_theme()[:palette]...)
    l = (layout_x = nothing, layout_y = nothing)
    discrete_scales = map(DiscreteScale, merge(scales[], spec.options, l))
    continuous_scales = map(ContinuousScale, spec.options)
    disc_options = applytheme(discrete_scales, spec.pkeys, unique_dict)
    cont_options = applytheme(continuous_scales, spec.style, extrema_dict)
    options = foldl(merge, (spec.options, disc_options, cont_options))
    return Spec{T}(pkeys=spec.pkeys, style=spec.style, options=options)
end

applytheme(scales, style::Style, metadata) = applytheme(scales, style.value, metadata)
function applytheme(scales, grp::NamedTuple{names}, metadata) where names
    res = map(names) do key
        haskey(scales, key) ? get_attr(scales[key], grp[key], metadata[key]) : grp[key]
    end
    return NamedTuple{names}(res)
end


text:
println("Running WCSUtils tests.")

import WCS: world_to_pix, pix_to_world

const dat_dir = joinpath(Pkg.dir("SloanDigitalSkySurvey"), "dat")

"Test that the identity WCSTransform works as expected."
function test_id_wcs()
    rand_coord = rand(2, 10)
    @test pix_to_world(WCSUtils.wcs_id, rand_coord) == rand_coord
    @test world_to_pix(WCSUtils.wcs_id, rand_coord) == rand_coord
end

function test_pixel_deriv_to_world_deriv()
    field_dir = joinpath(dat_dir, "sample_field")
    run_num = "003900"
    camcol_num = "6"
    field_num = "0269"

    # The gain is wrong but it doesn't matter.
    wcs = SDSS.load_raw_field(field_dir, run_num, camcol_num, field_num, 1, 1.0)[7];

    function test_fun(pix_loc::Array{Float64, 1})
        pix_loc[1]^2 + 0.5 * pix_loc[2]
    end

    function test_fun_grad(pix_loc::Array{Float64, 1})
        Float64[2 * pix_loc[1], 0.5 ]
    end

    function test_fun_world(world_loc::Array{Float64, 1}, wcs)
        pix_loc = WCSUtils.world_to_pix(wcs, world_loc)
        test_fun(pix_loc)
    end

    pix_del = 1e-3
    world_del = 1e-9
    pix_loc = Float64[5, 5]
    pix_loc_1 = pix_loc + pix_del * [1, 0]
    pix_loc_2 = pix_loc + pix_del * [0, 1]
    world_loc = pix_to_world(wcs, pix_loc)
    world_loc_1 = world_loc + world_del * [1, 0]
    world_loc_2 = world_loc + world_del * [0, 1]

    @test_approx_eq_eps test_fun(pix_loc) test_fun_world(world_loc, wcs) 1e-8

    pix_deriv = test_fun_grad(pix_loc)
    world_deriv = Float64[ (test_fun_world(world_loc_1, wcs) -
                            test_fun_world(world_loc, wcs)) / world_del
                           (test_fun_world(world_loc_2, wcs) -
                            test_fun_world(world_loc, wcs)) / world_del ]

    relative_err = (WCSUtils.pixel_deriv_to_world_deriv(wcs, pix_deriv, pix_loc) -
                    world_deriv) ./ abs(world_deriv)
    @test_approx_eq_eps relative_err [ 0 0 ] 1e-3
end

function test_world_to_pix()
    field_dir = joinpath(dat_dir, "sample_field")
    run_num = "003900"
    camcol_num = "6"
    field_num = "0269"

    # The gain will not be used.
    nelec, calib_col, sky_grid, sky_x, sky_y, sky_image, wcs =
        SDSS.load_raw_field(field_dir, run_num, camcol_num, field_num, 3, 1.0);

    pix_center = Float64[0.5 * size(nelec, 1), 0.5 * size(nelec, 1)]
    pix_loc = pix_center + [5., 3.]
    world_center = pix_to_world(wcs, pix_center)
    world_loc = pix_to_world(wcs, pix_loc)

    function test_jacobian(wcs, pix_center, world_center)
      wcs_jacobian = WCSUtils.pixel_world_jacobian(wcs, pix_center);

      pix_loc_test1 = world_to_pix(wcs, world_loc)
      pix_loc_test2 = world_to_pix(wcs_jacobian, world_center, pix_center,
                                   world_loc)

      # Note that the accuracy of the linear approximation isn't great.
      @test_approx_eq(pix_loc_test1, pix_loc)
      @test_approx_eq_eps(pix_loc_test2, pix_loc, 1e-2)
    end

    @test WCSUtils.pixel_world_jacobian(WCSUtils.wcs_id, pix_center) == [1.0 0.0; 0.0 1.0];

    test_jacobian(wcs, pix_center, world_center)
end

test_id_wcs()
test_pixel_deriv_to_world_deriv()
test_world_to_pix()


text:
# ---
# title: 824. Goat Latin
# id: problem824
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/goat-latin/description/>
# hidden: true
# ---
#
# A sentence `S` is given, composed of words separated by spaces. Each word
# consists of lowercase and uppercase letters only.
#
# We would like to convert the sentence to " _Goat Latin "_ (a made-up language
# similar to Pig Latin.)
#
# The rules of Goat Latin are as follows:
#
#   * If a word begins with a vowel (a, e, i, o, or u), append `"ma"` to the end of the word.
# For example, the word 'apple' becomes 'applema'.
#
#
#   * If a word begins with a consonant (i.e. not a vowel), remove the first letter and append it to the end, then add `"ma"`.
# For example, the word `"goat"` becomes `"oatgma"`.
#
#
#   * Add one letter `'a'` to the end of each word per its word index in the sentence, starting with 1.
# For example, the first word gets `"a"` added to the end, the second word gets
# `"aa"` added to the end and so on.
#
# Return the final sentence representing the conversion from `S` to Goat Latin.
#
#
#
# **Example 1:**
#
#
#
#     Input: "I speak Goat Latin"
#     Output: "Imaa peaksmaaa oatGmaaaa atinLmaaaaa"
#
#
# **Example 2:**
#
#
#
#     Input: "The quick brown fox jumped over the lazy dog"
#     Output: "heTmaa uickqmaaa rownbmaaaa oxfmaaaaa umpedjmaaaaaa overmaaaaaaa hetmaaaaaaaa azylmaaaaaaaaa ogdmaaaaaaaaaa"
#
#
#
#
# Notes:
#
#   * `S` contains only uppercase, lowercase and spaces. Exactly one space between each word.
#   * `1 <= S.length <= 150`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Roots, DataFrames, DataFramesMeta, Query

function compute_marginal_CI(α, Λ_of_t, wms_of_t; mass = 0.95)
    qinfCI, qupCI = get_quantiles_from_mass(mass)
    K = length(α)
    times = Λ_of_t |> keys |> collect |> sort
    Iterators.product(times, 1:K) |>
        collect |>
        x -> map(pr -> DataFrame(time = pr[1],
            margin = pr[2],
            infCI = compute_quantile_beta_mixture(α, Λ_of_t[pr[1]], wms_of_t[pr[1]], pr[2], qinfCI),
            supCI = compute_quantile_beta_mixture(α, Λ_of_t[pr[1]], wms_of_t[pr[1]], pr[2], qupCI)), x) |>
        dflist -> vcat(dflist...)

end

function compute_quantile_beta_mixture(α, Λ, wms, margin, q::Float64)
    function Beta_mixture_cdf(x::Real)
        return sum(wms.*Float64[cdf(Beta(m[margin] + α[margin], sum(m .+ α) - m[margin] - α[margin]),x) for m in Λ])
    end
    return fzero(x -> Beta_mixture_cdf(x)-q, 0, 10^9)
end

function marginal_densities(α, Λ_of_t, wms_of_t)
    K = α |> length
    map(margin -> for_one_marginal(α, margin, Λ_of_t, wms_of_t), 1:K) |>
    dfs -> vcat(dfs...)
end

function for_one_marginal(α, margin, Λ_of_t, wms_of_t)
    times = Λ_of_t |> keys |> collect |> sort

    marginal_psi_t = [create_marginal_beta_mixture(α, Λ_of_t[t], wms_of_t[t], margin) for t in times]
    [DataFrame(x = range(0, stop = 1, length = 100)) |> df -> @transform(df, dens = f.(:x)) for f in marginal_psi_t] |>
dfs -> map((dfi, ti) -> dfi |> @mutate(time = ti, variable = "x$margin") |> DataFrame, dfs, times) |>
    dfs -> vcat(dfs...)
end

function create_marginal_beta_mixture(α, Λ, wms, margin)
    function Beta_mixture(x::Real)
        return sum(wms.*Float64[pdf(Beta(m[margin] + α[margin], sum(m .+ α) - m[margin] - α[margin]),x) for m in Λ])
    end
    return Beta_mixture
end

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering