text:
@testset "Summarizer" begin
d = StringDocument("""
    Discount retailer Poundworld has appointed administrators, putting 5,100 jobs at risk.
    The move came after talks with a potential buyer, R Capital, collapsed leaving Poundworld with no option other than administration.
    Poundworld, which serves two million customers a week from 355 stores, also trades under the Bargain Buys name.
    Administrators Deloitte stress the stores will continue to trade as normal with no redundancies at this time.
    It said in a statement: Like many high street retailers, Poundworld has suffered from high product cost inflation, decreasing footfall, weaker consumer confidence and an increasingly competitive discount retail market.
    Clare Boardman, joint administrator at Deloitte, said: The retail trading environment in the UK remains extremely challenging and Poundworld has been seeking to address this through a restructure of its business.
    Unfortunately, this has not been possible.
    She said Deloitte believed a buyer could be found for the business, or at least part of it.
    A spokesperson for Poundworlds owner, TPG said filing for administration had been a difficult decision.
    Despite investing resources to strengthen the business, the decline in UK retail and changing consumer behaviour affected Poundworld significantly, they added.
    """)

s = summarize(d)
@test length(s) == 5
@test s[1] == "Discount retailer Poundworld has appointed administrators, putting 5,100 jobs at risk."

s = summarize(d; ns=2)
@test length(s) == 2
end


text:
using Discord
include("base.jl")

function attack(client::Client,m::Message)
    global management_dict

    no_reserve = true

    channel_id_str = string(m.channel_id)

    if !haskey(management_dict,channel_id_str)
        return reply(client,m,"初期設定が完了しておりません．`.set`コマンドにて初期設定を完了させてください")
    end

    if !haskey(management_dict[channel_id_str],"remain_attack_channel")
        return reply(client,m,"残凸把握板の登録が完了しておりません．`.set`コマンドにて初期設定を完了させてください")
    end

    if  !haskey(management_dict[channel_id_str],"output_channel")
        return reply(client,m,"予約確認板の登録が完了しておりません．`.set`コマンドにて初期設定を完了させてください")
    end

    boss_number = string(management_dict[channel_id_str]["boss_number"])
    reserve_list =  management_dict[channel_id_str]["reserve"][boss_number]["reserve_list"]

    #予約中の人を探してフラグを付け替える
    for i = 1:length(reserve_list)
        if reserve_list[i]["id"] == m.author.id
            management_dict[channel_id_str]["reserve"][boss_number]["reserve_list"][i]["attack"] = true
            no_reserve = false
            break
        end
    end

    #予約していなかったら怒る
    if no_reserve
        return reply(client,m,"予約されていません.")
    end

    output_message_content = create_output_message(management_dict[channel_id_str])
    edit_message(client,management_dict[channel_id_str]["output_channel"],management_dict[channel_id_str]["output_message"],content=output_message_content)

    dt_now =now(tz"Asia/Tokyo") - Hour(5)
    dt_now_str = "$(month(dt_now))月$(day(dt_now))日"
    if dt_now_str != management_dict[channel_id_str]["today"]#日付が変わっていれば残凸状況を修正する
        management_dict[channel_id_str]["today"] = dt_now_str
        management_dict[channel_id_str]["remain_attack"] =initiate_remain_attack_dict(client,management_dict[channel_id_str])
        remain_attack_message_content = make_remain_attack_message_content(management_dict[channel_id_str]["remain_attack"],management_dict[channel_id_str]["today"])
        remain_attack_message = fetchval(create_message(client,management_dict[channel_id_str]["remain_attack_channel"],content=remain_attack_message_content))

        management_dict[channel_id_str]["remain_attack_message"] = remain_attack_message.id
    end

    save_management_dict(management_dict)
    create(client, Reaction, m, '👍')
end

ext: jl

max_stars_repo_name: hibibol/clan-battle

max_issues_repo_name: hibibol/clan-battle

max_forks_repo_name: hibibol/clan-battle


text:
#Tests des packages
include("../ALAS/include.jl")
include("AMtoMPCC.jl")

using MPCCmod
using OutputRelaxationmod
using Penalty
using DDirection
using LineSearch
using AlgoSetmod

using NLPModels
using AmplNLReader

using BenchmarkProfiles

name_pb=["bard1","bard1m","bard2","bard2m","bard3","bard3m","bilevel1",
         "bilevel1m","bilevel2","bilevel2m","bilevel3","bilin",
         "dempe","desilva","df1","ex9.1.1","ex9.1.2","ex9.1.3","ex9.1.4",
         "ex9.1.5","ex9.1.6","ex9.1.7","ex9.1.8","ex9.1.9","ex9.1.10",
         "ex9.2.1","ex9.2.2","ex9.2.3","ex9.2.4","ex9.2.5",
         "ex9.2.6","ex9.2.7","ex9.2.8","ex9.2.9","flp2","gauvin","hakonsen",
         "hs044-i","jr1","jr2","kth1","kth2","kth3","monteiro","monteiroB",
         "outrata31","outrata32","outrata33","outrata34","qpec1","qpec2",
         "ralph1","ralph2","scholtes1","scholtes2","scholtes3","scholtes4",
         "scholtes5","sl1","stackelberg1"]
fmpec=[]

#10 juillet 17 ++
#good:60,59,57,56,52,49,48,47,46,43,42,41,40,39,36,35,15,14,12, 2, 1
#feas:58,55,54,51,33,32,31,30,29,28,26,25,23,19,17,16,13,11, 8, 7, 5
#nbdd:53, 3
# bad:50,45,44,38,37,34,27,24,22,21,20,18,10, 9, 6, 4
#to check the feasibility 34,27

npb=length(name_pb);npb=10

#output variables:
output=-ones(npb,1)
T=-ones(npb,1)

for i=1:npb

 AMnlp = AmplModel(string("benchmark/macmpec/",name_pb[i],".nl"))

 nlp,G,H=AMtoMPCC(AMnlp)

 @printf("MacMPEC%i: %s \n",i,AMnlp.meta.name)

 #resolution avec ALAS
 exemple_nlp=MPCCmod.MPCC(nlp,G,H)

 @time xb,fb,orb,nb_eval = MPCCsolve.solve(exemple_nlp)
 @printf("%s ",orb.solve_message)
 println(nb_eval)

 T[i,1]=orb.solve_message=="Success"?sum(nb_eval):Inf

 #manque les unbounded
 if orb.solve_message=="Success"
  output[i,1]=0
 elseif orb.solve_message=="Success (with sub-pb failure)"
  output[i,1]=0
 elseif orb.solve_message=="Infeasible"
  output[i,1]=2
 elseif orb.solve_message=="Feasible, but not optimal"
  output[i,1]=1
 else
  output[i,1]=3
 end

 amplmodel_finalize(AMnlp)
end


text:
module CellValuesMocks

using CellwiseValues

export TestIterCellValue
export TestIndexCellValue
import Base: iterate
import Base: length
import Base: size
import Base: getindex
import CellwiseValues: evaluate

struct TestIndexCellValue{T} <: IndexCellValue{T,1}
  a::T
  l::Int
end

size(self::TestIndexCellValue) = (self.l,)

getindex(self::TestIndexCellValue,cell::Int) = self.a

struct TestIterCellValue{T} <: IterCellValue{T}
  a::T
  l::Int
end

length(v::TestIterCellValue) = v.l

@inline function iterate(v::TestIterCellValue)
  i = 0
  if v.l <= i; return nothing; end
  (v.a,i+1)
end

@inline function iterate(v::TestIterCellValue,state)
  i = state
  if v.l <= i; return nothing; end
  (v.a,i+1)
end

end # module


text:
using FLuxSymbolEquationsV2
using Test

@testset "FLuxSymbolEquationsV2.jl" begin
    # Write your tests here.
end


text:
using Weave

let start_dir = pwd()
    cd(@__DIR__)

    weave("../examples/FIR_design.jmd", doctype = "pandoc", out_path = "build/examples")
    weave("../examples/FIR_design.jmd", doctype = "md2html", out_path = "build/examples")
    weave("../examples/FIR_design_plots.jl", doctype = "md2html", out_path = "build/examples")

    # PDF outputs
    if haskey(ENV, "TRAVIS")
        # in Travis, just cp already generated PDFs
        cp("assets/FIR_design.pdf", "build/examples/FIR_design.pdf", force = true)
        cp("assets/FIR_design_plots.pdf", "build/examples/FIR_design_plots.pdf", force = true)
    else
        # otherwise try to generate them
        try
            weave("../examples/FIR_design.jmd", doctype = "md2pdf", out_path = "assets")
            weave("../examples/FIR_design_plots.jl", doctype = "md2pdf", out_path = "assets")
        catch err
            @error err
        end
    end

    cp("../examples/FIR_design.jmd", "build/examples/FIR_design.jmd", force = true)
    cp("build/examples/FIR_design.md", "build/examples/FIR_design.txt", force = true)
    cp("../examples/FIR_design_plots.jl", "build/examples/FIR_design_plots.jl", force = true)

    cd(start_dir)
end


text:
using BinaryBuilder

name = "MacOSLibc"
version = v"10.10"

# sources to build, such as mingw32, our patches, etc....
sources = [
    "https://github.com/phracker/MacOSX-SDKs/releases/download/10.13/MacOSX10.10.sdk.tar.xz" =>
    "4a08de46b8e96f6db7ad3202054e28d7b3d60a3d38cd56e61f08fb4863c488ce"
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/MacOSX*.sdk
sysroot="${prefix}/${target}/sys-root"

mkdir -p "${sysroot}"
mv * "${sysroot}"
"""

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

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(joinpath(prefix, "usr", "lib"), "libSystem", :libSystem),
]

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

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


text:
using BenchmarkTools
using LarSurf

# function grid_voxel_id_to_cart(grid_size::Array{Int, 1}, ind::Int64)
function grid_voxel_id_to_cart_calcul(grid_size, ind)
#     ind = (sz2 * sz3) * (i - 1)  + (j - 1) * sz3 + k
    sz1,sz2,sz3 = grid_size
    layer = sz2*sz3
    ir = div(ind - 1, layer) + 1
    rest1 = mod(ind - 1 , layer) + 1
#     println("rest1 ", rest1)
    row = sz3
    jr = div(rest1 -1, row) + 1
    rest2 = mod(rest1 - 1, row) + 1
    kr = rest2
#     println("rests: ", rest1, " ", rest2)
    return [kr, jr, ir]
end

function grid_voxel_id_to_cart_cartind(blocks_number_axis, block_i)
    a = Array{Int}(
        undef,
        blocks_number_axis[1],
        blocks_number_axis[2],
        blocks_number_axis[3]
    )
#     print("block_i: ", block_i)
    # TODO maybe the cartesian index is slow
    bsub = CartesianIndices(a)[block_i]
    return bsub
end

blocks_number_axis = [10,10,10]
_cart = Array{Int}(
    undef,
    blocks_number_axis[1],
    blocks_number_axis[2],
    blocks_number_axis[3]
)

function grid_voxel_id_to_cart_extern(blocks_number_axis, block_i)
#     print("block_i: ", block_i)
    # TODO maybe the cartesian index is slow
    bsub = CartesianIndices(_cart)[block_i]
    return bsub
end

bl1 = grid_voxel_id_to_cart_calcul([10,10,10], 125)
println(bl1)
bl2 = grid_voxel_id_to_cart_cartind([10,10,10], 125)
println(bl2)
# bl3 = LarSurf.grid_voxel_id_to_cart([10,10,10], 125)
bl3 = grid_voxel_id_to_cart_extern([10,10,10], 125)
println(bl3)

bm = @benchmark grid_voxel_id_to_cart_calcul([10,10,10], 472)
display(bm)

bm = @benchmark grid_voxel_id_to_cart_cartind([10,10,10], 472)
display(bm)

bm = @benchmark grid_voxel_id_to_cart_extern([10,10,10], 472)
display(bm)

# bm = @benchmark LarSurf.grid_voxel_id_to_cart([10,10,10], 472)
# display(bm)


text:
@testset "BSplineManifold" begin
    @testset "1dim" begin
        @testset "BSplineManifold-1dim" begin
            Random.seed!(42)

            P1 = BSplineSpace(1, Knots([0, 0, 1, 1]))
            n1 = dim(P1) # 2
            a = [Point(i, rand()) for i in 1:n1]  # n1 × n2 array of d̂-dim vector.
            M = BSplineManifold([P1], a)
            @test dim(M) == 1

            P1′ = BSplineSpace(2, Knots([-2, 0, 0, 1, 1, 2]))
            p₊ = [1]
            k₊ = [Knots()]

            @test P1 ⊑ P1′

            M′ = refinement(M, [P1′])
            M′′ = refinement(M, p₊=p₊, k₊=k₊)
            ts = [[rand()] for _ in 1:10]
            for t in ts
                @test M(t) ≈ M′(t)
                @test M(t) ≈ M′′(t)
            end
        end

        @testset "FastBSplineManifold-1dim" begin
            Random.seed!(42)

            P1 = FastBSplineSpace(1, Knots([0, 0, 1, 1]))
            n1 = dim(P1) # 2
            a = [Point(i, rand()) for i in 1:n1]  # n1 × n2 array of d̂-dim vector.
            M = FastBSplineManifold([P1], a)
            @test dim(M) == 1

            P1′ = FastBSplineSpace(2, Knots([-2, 0, 0, 1, 1, 2]))
            p₊ = [1]
            k₊ = [Knots()]

            @test P1 ⊑ P1′

            M′ = refinement(M, [P1′])
            M′′ = refinement(M, p₊=p₊, k₊=k₊)
            ts = [[rand()] for _ in 1:10]
            for t in ts
                @test M(t) ≈ M′(t)
                @test M(t) ≈ M′′(t)
            end
        end

        @testset "BSplineCurve" begin
            Random.seed!(42)

            P1 = FastBSplineSpace(1, Knots([0, 0, 1, 1]))
            n1 = dim(P1) # 2
            a = [Point(i, rand()) for i in 1:n1]  # n1 × n2 array of d̂-dim vector.
            M = BSplineCurve([P1], a)
            @test dim(M) == 1

            P1′ = FastBSplineSpace(2, Knots([-2, 0, 0, 1, 1, 2]))
            p₊ = [1]
            k₊ = [Knots()]

            @test P1 ⊑ P1′

            M′ = refinement(M, [P1′])
            M′′ = refinement(M, p₊=p₊, k₊=k₊)
            ts = [[rand()] for _ in 1:10]
            for t in ts
                @test M(t) ≈ M′(t)
                @test M(t) ≈ M′′(t)
            end
        end
    end

    @testset "2dim" begin
        @testset "BSplineManifold-2dim" begin
            Random.seed!(42)

            P1 = BSplineSpace(1, Knots([0, 0, 1, 1]))
            P2 = BSplineSpace(1, Knots([1, 1, 2, 3, 3]))
            n1 = dim(P1) # 2
            n2 = dim(P2) # 3
            a = [Point(i, j) for i in 1:n1, j in 1:n2]  # n1 × n2 array of d̂-dim vector.
            M = BSplineManifold([P1, P2], a)
            @test dim(M) == 2

            P1′ = BSplineSpace(2, Knots([-2, 0, 0, 1, 1, 2]))
            P2′ = BSplineSpace(1, Knots([-3, 1, 2, 1.45, 3, 4]))
            p₊ = [1, 0]
            k₊ = [Knots(), Knots(1.45)]

            @test P1 ⊑ P1′
            @test P2 ⊑ P2′

            M′ = refinement(M, [P1′, P2′])
            M′′ = refinement(M, p₊=p₊, k₊=k₊)
            ts = [[rand(), 1 + 2 * rand()] for _ in 1:10]
            for t in ts
                @test M(t) ≈ M′(t)
                @test M(t) ≈ M′′(t)
            end
        end

        @testset "FastBSplineManifold-2dim" begin
            Random.seed!(42)

            P1 = FastBSplineSpace(1, Knots([0, 0, 1, 1]))
            P2 = FastBSplineSpace(1, Knots([1, 1, 2, 3, 3]))
            n1 = dim(P1) # 2
            n2 = dim(P2) # 3
            a = [Point(i, j) for i in 1:n1, j in 1:n2]  # n1 × n2 array of d̂-dim vector.
            M = FastBSplineManifold([P1, P2], a)
            @test dim(M) == 2

            P1′ = FastBSplineSpace(2, Knots([0, 0, 0, 1, 1, 1]))
            P2′ = FastBSplineSpace(1, Knots([1, 1, 2, 1.45, 3, 3]))
            p₊ = [1, 0]
            k₊ = [Knots(), Knots(1.45)]

            @test P1 ⊆ P1′
            @test P2 ⊆ P2′

            M′ = refinement(M, [P1′, P2′])
            M′′ = refinement(M, p₊=p₊, k₊=k₊)
            ts = [[rand(), 1 + 2 * rand()] for _ in 1:10]
            for t in ts
                @test M(t) ≈ M′(t)
                @test M(t) ≈ M′′(t)
            end
        end

        @testset "BSplineSurface" begin
            Random.seed!(42)

            P1 = FastBSplineSpace(1, Knots([0, 0, 1, 1]))
            P2 = FastBSplineSpace(1, Knots([1, 1, 2, 3, 3]))
            n1 = dim(P1) # 2
            n2 = dim(P2) # 3
            a = [Point(i, j) for i in 1:n1, j in 1:n2]  # n1 × n2 array of d̂-dim vector.
            M = BSplineSurface([P1, P2], a)
            @test dim(M) == 2

            P1′ = FastBSplineSpace(2, Knots([0, 0, 0, 1, 1, 1]))
            P2′ = FastBSplineSpace(1, Knots([1, 1, 2, 1.45, 3, 3]))
            p₊ = [1, 0]
            k₊ = [Knots(), Knots(1.45)]

            @test P1 ⊆ P1′
            @test P2 ⊆ P2′

            M′ = refinement(M, [P1′, P2′])
            M′′ = refinement(M, p₊=p₊, k₊=k₊)
            ts = [[rand(), 1 + 2 * rand()] for _ in 1:10]
            for t in ts
                @test M(t) ≈ M′(t)
                @test M(t) ≈ M′′(t)
            end
        end
    end
end


text:
using Test, AxisKeys

@testset "offset" begin
    using OffsetArrays

    o = OffsetArray(rand(1:99, 5), -2:2)
    w = wrapdims(o, i='a':'e')
    @test axiskeys(w,1) isa OffsetArray
    @test w[i=-2] == w('a')
    @test_throws ArgumentError KeyedArray(o, i='a':'e')

    w′ = wrapdims(o)
    @test axiskeys(w′,1) == -2:2

end
@testset "unique" begin
    using UniqueVectors

    u = wrapdims(rand(Int8,5,1), UniqueVector, [:a, :b, :c, :d, :e], nothing)
    @test axiskeys(u,1) isa UniqueVector
    @test u(:b) == u[2,:]

    n = wrapdims(rand(2,100), UniqueVector, x=nothing, y=rand(Int,100))
    @test axiskeys(n,1) isa UniqueVector
    k = axiskeys(n, :y)[7]
    @test n(y=k) == n[:,7]

end
@testset "tables" begin
    using Tables

    @testset "source" begin
        R = wrapdims(rand(2,3), 11:12, 21:23)
        N = wrapdims(rand(2,3), a=[11, 12], b=[21, 22, 23.0])

        @test keys(first(Tables.rows(R))) == (:dim_1, :dim_2, :value)
        @test keys(first(Tables.rows(N))) == (:a, :b, :value)

        @test Tables.columns(N).a == [11, 12, 11, 12, 11, 12]
    end
    @testset "sink" begin
        A = KeyedArray(rand(24, 11, 3); time = 0:23, loc = -5:5, id = ["a", "b", "c"])
        table = Tables.columntable(A)

        # Test fully constructing from a table
        # Common when working with adhoc data
        B = wrapdims(table, :value, :time, :loc, :id)
        @test B == A

        # Test wrapping of key vectors, and wrong order:
        U = wrapdims(table, UniqueVector, :value, :id, :time, :loc)
        @test axiskeys(U, :time) isa UniqueVector
        @test U(time=3, id="b") == A(time=3, id="b")

        # Test populating an existing array (e.g., expected data based on calculated targets/offsets)
        C = KeyedArray(
            zeros(Float64, size(A));
            time = unique(table.time),
            loc = unique(table.loc),
            id = unique(table.id),
        )
        @test C != A
        AxisKeys.populate!(C, table, :value)
        @test C == A

        # Constructing a NamedDimsArray with different default value and table type
        # Partial populating
        r_table = Tables.rowtable(A)
        n = length(r_table)
        idx = rand(Bool, n)
        D = wrapdims(r_table[idx], :value, :time, :loc, :id; default=missing)
        # dimnames should still match, but we'll have missing values
        @test dimnames(D) == dimnames(A)
        @test any(ismissing, D)

        # BTW, this is why it's a method of wrapdims, not KeyedArray:
        # @code_warntype wrapdims(table, :value, :time, :loc, :id) # ::Any
        # @code_warntype wrapdims(r_table[idx], :value, :time, :loc, :id; default=missing)

        # Construction with invalid columns error as expected, but the specific error is
        # dependent on the table type.
        # ERROR: ArgumentError: wrong number of names, got (:q, :time, :loc, :id) with ndims(A) == 1
        @test_throws ArgumentError wrapdims(Tables.rowtable(A), :q, :time, :loc, :id)
        # ERROR: ArgumentError: wrong number of names, got (:value, :p, :loc, :id) with ndims(A) == 1
        @test_throws ArgumentError wrapdims(Tables.rowtable(A), :value, :p, :loc, :id)
        # ERROR: type NamedTuple has no field q
        @test_throws ErrorException wrapdims(Tables.columntable(A), :q, :time, :loc, :id)
        # ERROR: type NamedTuple has no field p
        @test_throws ErrorException wrapdims(Tables.columntable(A), :value, :p, :loc, :id)

        # Construction with duplicates
        # ERROR: ArgumentError: Key (Date("2019-01-01"), -5) is not unique
        @test_throws ArgumentError wrapdims(table, :value, :time, :loc)
        @test wrapdims(r_table, :value, :time, :loc; force=true) == C(:, :, Key("c"))
    end
end
@testset "stack" begin
    using LazyStack

    rin = [wrapdims(1:3, a='a':'c') for i=1:4]

    @test axiskeys(stack(rin), :a) == 'a':'c'
    @test axiskeys(stack(:b, rin...), :a) == 'a':'c' # tuple
    @test axiskeys(stack(z for z in rin), :a) == 'a':'c' # generator

    rout = wrapdims([[1,2], [3,4]], b=10:11)
    @test axiskeys(stack(rout), :b) == 10:11

    rboth = wrapdims(rin, b=10:13)
    @test axiskeys(stack(rboth), :a) == 'a':'c'
    @test axiskeys(stack(rboth), :b) == 10:13

    nts = [(i=i, j="j", k=33) for i=1:3]
    @test axiskeys(stack(nts), 1) == [:i, :j, :k]
    @test axiskeys(stack(:z, nts...), 1) == [:i, :j, :k]
    @test axiskeys(stack(n for n in nts), 1) == [:i, :j, :k]

end
@testset "dates" begin
    using Dates

    D = wrapdims(rand(2,53), row = [:one, :two], week = Date(2020):Week(1):Date(2021))
    w9 = axiskeys(D,:week)[9]
    @test w9 isa Date
    @test D(w9) == D[week=9]
    # But steps of Year(1) don't work, https://github.com/JuliaLang/julia/issues/35203

    @test D(==(Date(2020, 1, 8))) == D[:, 2:2]
    @test D(Near(Date(2020, 1, 10))) == D(Date(2020, 1, 8)) == D[:, 2]
    int = Interval(Date(2020, 1, 8), Date(2020, 1, 22))
    @test D(int) == D[:, 2:4]
    axiskeys(D[:,int],2) isa Array # ideally would be a range

end
@testset "inverted" begin
    using InvertedIndices

    K = wrapdims(rand(4,5))
    @test K[:, Not(4)] == K[:, vcat(1:3, 5)] == K(:, Base.Fix2(!=,4)) == K(:, Not(4))
    @test K[Not(1,4), :] == K[2:3, :] == K(r -> 2<=r<=3, :)

    N = wrapdims(rand(Int8, 2,3,4), a=[:one, :two], b='α':'γ', c=31:34)
    @test N[b=Not(2)] == N[:,[1,3],:] == N(b=Base.Fix2(!=,'β')) == N(:,['α','γ'],:)
    @test N[c=Not(2,4)] == N(c=Index[Not(2,4)]) == N(c=Not(2,4))

end
@testset "fourier" begin
    using FFTW

    times = 0.1:0.1:10
    data = rand(100) ./ 10; data[1:5:end] .= 1;
    A = KeyedArray(data, times)
    Atil = fft(A)
    @test axiskeys(Atil,1)[end] == -0.1

    Ashift = fftshift(Atil)
    @test axiskeys(Ashift,1)[end] == +4.9

    Ar = rfft(A) # different size, different freqencies
    @test axiskeys(Ar,1)[end] == 5.0

    A2 = ifft(ifftshift(Ashift))
    A ≈ A2
    @test all(mod.(axiskeys(A,1) .- axiskeys(A2,1),10) .≈ 0.1) # that's OK I think?

    data2 = randn(32,2);
    B = wrapdims(data2, time=100:10:410.0, col=[:a, :b])
    @test_skip Btil = fftshift(fft(B, :time), :time) # result does not have names
    Btil = fftshift(fft(B, :time), 1)
    @test parent(Btil) ≈ fftshift(fft(data2,1),1)

end
@testset "unitful fourier" begin
    using FFTW
    using Unitful: s

    times = 0.1s:0.1s:10s
    A = wrapdims(rand(100), t=times)

    ifft(fft(A))
    fft(ifft(A))
    irfft(rfft(A), 100)

    abs.(fft(A)) # keys remain Frequencies thanks to copy method

    fftshift(fft(A))
    ifftshift(A) # keys become a Vector here, unavoidable I think

    @test axiskeys(sortkeys(fft(A)),1) ≈ axiskeys(fftshift(fft(A)),1)
    @test_broken sortkeys(fft(A)) ≈ fftshift(fft(A)) # isapprox should be used for keys

end

@testset "statsbase" begin
    using CovarianceEstimation, StatsBase

    A = rand(4, 3)
    A_ka = KeyedArray(A, (0:3, [:a, :b, :c]))
    A_kanda = KeyedArray(A; time = 0:3, id = [:a, :b, :c])
    wv = aweights(rand(4))

    @testset "$f" for f in (mean, std, var)
        R = f == mean ? f(A, wv; dims=1) : f(A, wv, 1; corrected=true)

        R_ka = f(A_ka, wv; dims=1)
        R_kanda_int = f(A_kanda, wv; dims=1)
        R_kanda_sym = f(A_kanda, wv; dims=:time)
        expected_keys = (Base.OneTo(1), [:a, :b, :c])
        expected_names = (:time, :id)

        @test dimnames(R_kanda_int) == dimnames(R_kanda_sym) == expected_names
        @test axiskeys(R_ka) == axiskeys(R_kanda_int) == axiskeys(R_kanda_sym) == expected_keys
        @test parent(R_ka) ≈ parent(parent(R_kanda_int)) ≈ parent(parent(R_kanda_sym)) ≈ R
    end

    @testset "$f" for f in (cov, cor, scattermat)
        # Inconsistent statsbase behaviour
        kwargs = f === cov ? [:corrected => true] : []
        R = f === scattermat ? f(A, wv; dims=1, kwargs...) : f(A, wv, 1; kwargs...)

        R_ka = f(A_ka, wv; dims=1, kwargs...)
        R_kanda_int = f(A_kanda, wv; dims=1, kwargs...)
        R_kanda_sym = f(A_kanda, wv; dims=:time, kwargs...)
        expected_keys = ([:a, :b, :c], [:a, :b, :c])
        expected_names = (:id, :id)

        @test dimnames(R_kanda_int) == dimnames(R_kanda_sym) == expected_names
        @test axiskeys(R_ka) == axiskeys(R_kanda_int) == axiskeys(R_kanda_sym) == expected_keys
        @test parent(R_ka) ≈ parent(parent(R_kanda_int)) ≈ parent(parent(R_kanda_sym)) ≈ R
    end

    @testset "$f" for f in (mean_and_var, mean_and_std, mean_and_cov)
        R1, R2 = f(A, wv, 1; corrected=true)
        R1_ka, R2_ka = f(A_ka, wv; dims=1, corrected=true)
        R1_kanda_int, R2_kanda_int = f(A_kanda, wv; dims=1, corrected=true)
        R1_kanda_sym, R2_kanda_sym = f(A_kanda, wv; dims=:time, corrected=true)

        @test parent(R1_ka) ≈ parent(parent(R1_kanda_int)) ≈ parent(parent(R1_kanda_sym)) ≈ R1
        @test parent(R2_ka) ≈ parent(parent(R2_kanda_int)) ≈ parent(parent(R2_kanda_sym)) ≈ R2
    end

    @testset "conversions" begin
        @testset "cov2cor" begin
            @test cov2cor(cov(A_ka; dims=1), std(A_ka; dims=1)) ≈ cor(A_ka; dims=1)
            @test cov2cor(cov(A_ka; dims=2), std(A_ka; dims=2)) ≈ cor(A_ka; dims=2)
            @test cov2cor(cov(A_kanda; dims=:time), std(A_kanda; dims=:time)) ≈ cor(A_kanda; dims=:time)
            @test cov2cor(cov(A_kanda; dims=:id), std(A_kanda; dims=:id)) ≈ cor(A_kanda; dims=:id)
        end
        @testset "cor2cov" begin
            @test cor2cov(cor(A_ka; dims=1), std(A_ka; dims=1)) ≈ cov(A_ka; dims=1)
            @test cor2cov(cor(A_ka; dims=2), std(A_ka; dims=2)) ≈ cov(A_ka; dims=2)
            @test cor2cov(cor(A_kanda; dims=:time), std(A_kanda; dims=:time)) ≈ cov(A_kanda; dims=:time)
            @test cor2cov(cor(A_kanda; dims=:id), std(A_kanda, dims=:id)) ≈ cov(A_kanda; dims=:id)
        end
    end

    @testset "covariance estimation" begin
        ce = SimpleCovariance()

        @testset "unweighted" begin
            R = cov(ce, A; dims=1)
            R_ka = cov(ce, A_ka; dims=1)
            R_kanda_int = cov(ce, A_kanda; dims=1)
            R_kanda_sym = cov(ce, A_kanda; dims=:time)
            expected_keys = ([:a, :b, :c], [:a, :b, :c])
            expected_names = (:id, :id)

            @test dimnames(R_kanda_int) == dimnames(R_kanda_sym) == expected_names
            @test axiskeys(R_ka) == axiskeys(R_kanda_int) == axiskeys(R_kanda_sym) == expected_keys
            @test parent(R_ka) ≈ parent(parent(R_kanda_int)) ≈ parent(parent(R_kanda_sym)) ≈ R
        end

        @testset "weighted" begin
            R = cov(ce, A, wv; dims=1)
            R_ka = cov(ce, A_ka, wv; dims=1)
            R_kanda_int = cov(ce, A_kanda, wv; dims=1)
            R_kanda_sym = cov(ce, A_kanda, wv; dims=:time)
            expected_keys = ([:a, :b, :c], [:a, :b, :c])
            expected_names = (:id, :id)

            @test dimnames(R_kanda_int) == dimnames(R_kanda_sym) == expected_names
            @test axiskeys(R_ka) == axiskeys(R_kanda_int) == axiskeys(R_kanda_sym) == expected_keys
            @test parent(R_ka) ≈ parent(parent(R_kanda_int)) ≈ parent(parent(R_kanda_sym)) ≈ R
        end
    end
end


text:
module LadderDGA
#Base.Experimental.@optlevel 3

include("DepsInit.jl")

export sum_q, sum_freq, get_sum_helper
export ModelParameters, SimulationParameters, FreqGrid, EnvironmentVars
export LocalQuantities, NonLocalQuantities, ΓT, BubbleT, GνqT
export readConfig, setup_LDGA, calc_bubble, calc_χ_trilex, calc_Σ, full_run
export λ_correction!,calc_λsp_rhs_usable, calc_λsp_correction!
export calculate_Σ_ladder, writeFortranΣ

end


text:
module InheritanceTest
using ConcreteInterfaces.ComputedFieldTypes
using Base.Test

abstract type Bar end

@computed struct Parametric{T} <: Bar
    x::Base.promote_op(+, T, Float64)
end

@computed struct NonParametric <: Bar
    x::Base.promote_op(+, Int, Float64)
end

@computed struct NonParametricNoParent
    x::Base.promote_op(+, Int, Float64)
end

@testset "inheritance" begin
    @test isa(@inferred(Parametric{Int}(1)), Parametric{Int, Float64})
    @test isa(Parametric{Int}(1), Bar)
    @test isa(@inferred(Parametric{Float64}(1.0)), Parametric{Float64, Float64})

    @test isa(@inferred(NonParametric(1)), NonParametric{Float64})
    @test isa(NonParametric(1), Bar)

    @test isa(@inferred(NonParametricNoParent(1)), NonParametricNoParent{Float64})
    @test !isa(NonParametricNoParent(1), Bar)
end
end


text:
cd("/Users/adrienpajon/Documents/all.nosync/git_repo/pulse_analysis/2020-11-12 xp pulse all seuil/Julia_analysis")

## Open Libraries
using JLD

using CSV

using DataFrames

using Plots

using StatsPlots

using LinearAlgebra

# using FFTW #pwelch ? spectrum PSD ? (matlab)#
using DSP

# using Pyplot
pyplot()

## Constant
Ts_MTiXsens=0.01;
#Ts_Xsensor=0.02; #idéalement 50Hz
Ts_Xsensor=1/55.2; #correctif dans la pratique on est plus vers 55.2Hz
Ts_IBS = 0.02
convertlbpo2Nm = 0.11298482933333

## MTiXsens
#filenameMTiXsens : namefile of each experiment
#filenameMTiXsensID : name ID of each sensor with keys that refers to the position on the body
#sensorXsensNames : keys of each sensor equivalent to their position on the body
# include("imports/Import_filename_MTiXsens.jl")
filenameMTiXsens=load("imports/import_filename_MTiXsens_pneu.jld")["filenameMTiXsens"]
sensorXsensNames=load("imports/import_filename_MTiXsens_pneu.jld")["sensorXsensNames"]
filenameMTiXsensID=load("imports/import_filename_MTiXsens_pneu.jld")["filenameMTiXsensID"]

## Xsensor
#filenameXsensor : namefile of each experiment
# include("imports/Import_filename_Xsensor.jl")
filenameXsensor=load("imports/import_filename_Xsensor_pneu.jld")["filenameXsensor"]

## Instrumented_backseat
#filenameIBS : namefile of each experiment
# include("imports/Import_filename_IBS.jl")
filenameIBS=load("imports/import_filename_IBS_pneu.jld")["filenameIBS"]

# ## delete bad file's filename
# include("functions/deleterows.jl")
# badFileMTiXsens=[43]
# filenameMTiXsens=deleterows(filenameMTiXsens,badFileMTiXsens)
# # filenameXsensor=deleterows(filenameMTiXsens,badFile)
# badFileIBS=[41 63 67]
# filenameIBS=deleterows(filenameIBS,badFileIBS)

println("FileName import finished")

####################################################################
#script to read CSV and store in
# MTiXsens
# Xsensor
# IBS
encoderFolderPath="../measures/pneu/encoder_file/"
IBSFolderPath="../measures/pneu/IBS/"
XsensorFolderPath="../measures/pneu/Xsensor/"
MTiXsensFolderPath="../measures/pneu/MTiXsens/"
include("scripts/script_read_CSV.jl")

####################################################################
#script to extract data from MTiXsens and Xsensor in
# AccFree_MTiXsens
# Xsensor_reduce
include("scripts/script_data_extract.jl")
# include("scripts/script_reduce_IBS.jl")

####################################################################
#script example to find synchro values
#include("Main-search_synchro.jl")

# include("imports/Import_synchro.jl")
sync=load("imports/import_synchro_pneu.jld")["sync"]
println("Synchro import finished")

#script to synchronize datas in
#AccFree_MTiXsens_sync
#Xsensor_reduce_sync
#IBS_sync
include("scripts/script_synchro_datas.jl")

# include("imports/Import_synchro_perturb.jl")
sync_=load("imports/Import_synchro_perturb_pneu.jld")["sync_"]
step_max_perturb=load("imports/Import_synchro_perturb_pneu.jld")["step_max_perturb"]
println("Synchro perturb import finished")

include("scripts/script_synchro_perturb_datas.jl")

####################################################################
# #Plot datas
#
# #Plot MTiXsens acc in plot_Acc_MTiXsens
# include("scripts/figure_MTiXsensAcc.jl")
#
# #Plot Xsensor load and COM in plot_Xsensor
# include("scripts/figure_Xsensor.jl")
#
# #Plot Optotrak position in plot_Optotrak
# include("scripts/figure_Optotrak.jl")

####################################################################
#include("functions/lowpassf.jl")
#AccSignalFreeFilt=Dict()
#for n in sensorXsensNames
#    AccSignalFreeFilt[n]=[lowpassf(AccSignalFree[n][:,1],Ts,0.011) lowpassf(AccSignalFree[n][:,2],Ts,0.011) lowpassf(AccSignalFree[n][:,3],Ts,0.011)]
#end
#
#println("lowpass filter on signal finished")

####################################################################
#plot of AccSignalFreeFilt
#include("figure_AccFilt.jl")

####################################################################
#plot of AccSignalFreeFilt
#include("figure_AccFilt.jl")

####################################################################

println("Computation finished")


text:
module TIQI

include("colormaps.jl")
include("ionizer.jl")
include("plots.jl")
include("thresholding.jl")

end


text:
using .LinearSpace

struct LongDelays{T} end

import Base: transpose, +, -, *

abstract type AbstractOperator end

#@linearspace AbstractOperator{T} T

"""
*Atomic operator*: one that assemblechunk can deal with
"""
abstract type Operator <: AbstractOperator end

mutable struct TransposedOperator <: Operator
    op::Operator
end

scalartype(op::TransposedOperator) = scalartype(op.op)

mutable struct LinearCombinationOfOperators{T} <: AbstractOperator
    coeffs::Vector{T}
    ops::Vector
end

function scalartype(op::LinearCombinationOfOperators{T}) where T
    W = promote_type(T, scalartype(op.ops[1]))
    for i in 2:length(op.ops)
        W = promote_type(W, scalartype(op.ops[i]))
    end
    W
end

function +(a::LinearCombinationOfOperators, b::Operator)
    LinearCombinationOfOperators(
        [a.coeffs;[1.0]],
        [a.ops;[b]]
    )
end

function +(a::LinearCombinationOfOperators, b::LinearCombinationOfOperators)
    LinearCombinationOfOperators(
        [a.coeffs; b.coeffs],
        [a.ops; b.ops]
    )
end

+(a::Operator, b::LinearCombinationOfOperators) = b + a
+(a::Operator, b::Number) = a + (b * Identity())
+(a::Number, b::Operator) = b + a

*(a::Number, b::Operator) = LinearCombinationOfOperators([a], [b])
*(a::Number, b::LinearCombinationOfOperators) = LinearCombinationOfOperators(a * b.coeffs, b.ops)
-(a::AbstractOperator, b::AbstractOperator) = a + (-1.0) * b

function +(a::Operator, b::Operator)
    LinearCombinationOfOperators(
        [1.0, 1.0],
        [a, b]
    )
end

transpose(op::TransposedOperator) = op.op
transpose(op::Operator) = TransposedOperator(op)

function assemble(operator::AbstractOperator, test_functions, trial_functions,
    storage_policy = Val{:bandedstorage},
    long_delays_policy = LongDelays{:ignore})
    # This is a convenience function whose only job is to allocate
    # the storage for the interaction matrix. Further dispatch on
    # operator and space types is handled by the 4-argument version
    Z, store = allocatestorage(operator, test_functions, trial_functions,
        storage_policy, long_delays_policy)
    assemble!(operator, test_functions, trial_functions, store)
    sdata(Z)
end

function assemblerow(operator::AbstractOperator, test_functions, trial_functions,
    storage_policy = Val{:bandedstorage},
    long_delays_policy = LongDelays{:ignore})

    Z, store = allocatestorage(operator, test_functions, trial_functions,
        storage_policy, long_delays_policy)
    assemblerow!(operator, test_functions, trial_functions, store)
    sdata(Z)
end

function assemblecol(operator::AbstractOperator, test_functions, trial_functions,
    storage_policy = Val{:bandestorage},
    long_delays_policy = LongDelays{:ignore})

    Z, store = allocatestorage(operator, test_functions, trial_functions,
        storage_policy, long_delays_policy)
    assemblecol!(operator, test_functions, trial_functions, store)
    sdata(Z)
end

function allocatestorage(operator::AbstractOperator, test_functions, trial_functions,
    ::Type{Val{:bandedstorage}},
    ::Type{LongDelays{:ignore}})

    T = promote_type(
        scalartype(operator)       ,
        scalartype(test_functions) ,
        scalartype(trial_functions),
    )
    Z = SharedArray{T}(
        numfunctions(test_functions)  ,
        numfunctions(trial_functions),
    )
    fill!(Z, 0)
    store(v,m,n) = (Z[m,n] += v)
    return Z, store
end

function allocatestorage(operator::LinearCombinationOfOperators,
        test_functions::SpaceTimeBasis, trial_functions::SpaceTimeBasis,
        storage_policy::Type{Val{:bandedstorage}},
        long_delays_policy::Type{LongDelays{:ignore}})

    # TODO: remove this ugly, ugly patch
    Z, store = allocatestorage(operator.ops[end], test_functions, trial_functions,
        storage_policy, long_delays_policy)
end

function allocatestorage(operator::LinearCombinationOfOperators,
        test_functions::SpaceTimeBasis, trial_functions::SpaceTimeBasis,
        storage_policy::Type{Val{S}},
        long_delays_policy::Type{LongDelays{L}}) where {L,S}

    # TODO: remove this ugly, ugly patch
    Z, store = allocatestorage(operator.ops[end], test_functions, trial_functions,
        storage_policy, long_delays_policy)
end

function assemble!(operator::Operator, test_functions::Space, trial_functions::Space, store)

    P = Threads.nthreads()
    numchunks = P
    @assert numchunks >= 1
    splits = [round(Int,s) for s in range(0, stop=numfunctions(test_functions), length=numchunks+1)]

    Threads.@threads for i in 1:P
        lo, hi = splits[i]+1, splits[i+1]
        lo <= hi || continue
        test_functions_p = subset(test_functions, lo:hi)
        store1 = (v,m,n) -> store(v,lo+m-1,n)
        assemblechunk!(operator, test_functions_p, trial_functions, store1)
    end

end

# function assemble_st!(operator::Operator, test_functions::Space, trial_functions::Space, store)
#
#     # This method should only be called for `atomic` discrete operators, this means
#     # in particular that the spaces of test and trial functions are fully conforming
#     # to the Space concept and that the operator/space combinations conform the
#     # concept implicitly defined by the assemble! function in itegralop.jl
#     # (no more transposes or repositioning in a larger system for example)
#
#     @warn "Parallel assembly support disabled."
#
#     P = procs()
#
#     @assert length(P) == 1
#
#     if length(P) > 1; P = P[2:end]; end
#     numchunks = length(P)
#     @assert numchunks >= 1
#     splits = [round(Int,s) for s in linspace(0, numfunctions(test_functions), numchunks+1)]
#
#     T = typeof(test_functions)
#     S = eltype(test_functions.fns)
#
#     assemblechunk!(operator, test_functions, trial_functions, store)
# end

function assemble!(op::TransposedOperator, tfs::Space, bfs::Space, store)

    store1(v,m,n) = store(v,n,m)
    assemble!(op.op, bfs, tfs, store1)
end

function assemble!(op::LinearCombinationOfOperators, tfs::AbstractSpace, bfs::AbstractSpace, store)
    for (a,A) in zip(op.coeffs, op.ops)
        store1(v,m,n) = store(a*v,m,n)
        assemble!(A, tfs, bfs, store1)
    end
end

# Support for direct product spaces
function assemble!(op::Operator, tfs::DirectProductSpace, bfs::Space, store)
    I = Int[0]
    for s in tfs.factors push!(I, last(I) + numfunctions(s)) end
    for (i,s) in enumerate(tfs.factors)
        store1(v,m,n) = store(v,m + I[i], n)
        assemble!(op, s, bfs, store1)
    end
end

function assemble!(op::Operator, tfs::Space, bfs::DirectProductSpace, store)
    J = Int[0]
    for s in bfs.factors push!(J, last(J) + numfunctions(s)) end
    for (j,s) in enumerate(bfs.factors)
        store1(v,m,n) = store(v,m,n + J[j])
        assemble!(op, tfs, s, store1)
    end
end

function assemble!(op::Operator, tfs::DirectProductSpace, bfs::DirectProductSpace, store)
    I = Int[0]
    for s in tfs.factors push!(I, last(I) + numfunctions(s)) end
    for (i,s) in enumerate(tfs.factors)
        store1(v,m,n) = store(v,m + I[i],n)
        assemble!(op, s, bfs, store1)
    end
end


text:
using LoggingExtras
using LokiLogger
using Dates

# https://github.com/JuliaLang/julia/blob/1b93d53fc4bb59350ada898038ed4de2994cce33/base/logging.jl#L142-L151
function Base.parse(::Type{LogLevel}, s::String)
    if     s == string(Logging.BelowMinLevel) Logging.BelowMinLevel
    elseif s == string(Logging.Debug) Logging.Debug
    elseif s == string(Logging.Info) Logging.Info
    elseif s == string(Logging.Warn) Logging.Warn
    elseif s == string(Logging.Error) Logging.Error
    elseif s == string(Logging.AboveMaxLevel) Logging.AboveMaxLevel
    else
        m = match(r"LogLevel\((?<level>-?[1-9]\d*)\)", s)
        if isnothing(m)
            throw(ArgumentError("unknown log level"))
        else
            Logging.LogLevel(parse(Int, m[:level]))
        end
    end

end

function create_log_transformer(date_format)
    function transformer(log)
        merge(
            log,
            (
                datetime = Dates.format(now(), date_format),
                from=self(),
                myid=myid(),
            )
        )
    end
end

function create_default_fmt(with_color=false, is_expand_stack_trace=false)
    function default_fmt(iob, args)
        level, message, _module, group, id, file, line, kw = args
        color, prefix, suffix = Logging.default_metafmt(
            level, _module, group, id, file, line
        )
        ignore_fields = (:datetime, :path, :myid, :from)
        if with_color
            printstyled(iob, "$(kw.datetime) "; color=:light_black)
            printstyled(iob, prefix; bold=true, color=color)
            printstyled(iob, "[$(kw.from)@$(kw.myid)]"; color=:green)
            print(iob, message)
            for (k,v) in pairs(kw)
                if k ∉ ignore_fields
                    print(iob, " ")
                    printstyled(iob, k; color=:yellow)
                    printstyled(iob, "="; color=:light_black)
                    print(iob, v)
                end
            end
            !isempty(suffix) && printstyled(iob, " ($suffix)"; color=:light_black)
            println(iob)
        else
            print(iob, "$(kw.datetime) $prefix[$(kw.from)@$(kw.myid)]$message")
            for (k,v) in pairs(kw)
                if k ∉ ignore_fields
                    print(iob, " $k=$v")
                end
            end
            !isempty(suffix) && print(iob, " ($suffix)")
            println(iob)
        end
    end
end

function create_logger(config::Config)
    sinks = []

    if !isnothing(config.logging.loki_logger)
        push!(sinks, LokiLogger.Logger(config.logging.loki_logger.url))
    end

    if !isnothing(config.logging.driver_logger)
        driver_sinks = []
        console_logger_config = config.logging.driver_logger.console_logger
        if !isnothing(console_logger_config)
            push!(
                driver_sinks,
                FormatLogger(
                    create_default_fmt(
                        config.color,
                        console_logger_config.is_expand_stack_trace
                    )
                )
            )
        end
        rotating_logger_config = config.logging.driver_logger.rotating_logger
        if !isnothing(rotating_logger_config)
            mkpath(rotating_logger_config.path)
            push!(
                driver_sinks,
                DatetimeRotatingFileLogger(
                    create_default_fmt(),
                    rotating_logger_config.path,
                    rotating_logger_config.file_format,
                )
            )
        end
        if isempty(driver_sinks)
            push!(driver_sinks, current_logger())
        end
        push!(sinks, DriverLogger(TeeLogger(driver_sinks...)))
    end

    if isempty(sinks)
        push!(sinks, current_logger())
    end

    TeeLogger(
        (
            MinLevelLogger(
                TransformerLogger(
                    create_log_transformer(config.logging.date_format),
                    s
                ),
                parse(Logging.LogLevel, config.logging.log_level)
            )
            for s in sinks
        )...
    )
end

#####

Base.@kwdef struct DriverLogger <: AbstractLogger
    logger::TeeLogger
end

Logging.shouldlog(::DriverLogger, args...) = true
Logging.min_enabled_level(::DriverLogger) = Logging.BelowMinLevel
Logging.catch_exceptions(::DriverLogger) = true

struct LogMsg
    args
    kw
end

Logging.handle_message(logger::DriverLogger, args...; kw...) = LogMsg(args, kw) |> LOGGER

function (L::DriverLogger)(msg::LogMsg)
    handle_message(L.logger, msg.args...;msg.kw...)
end


text:
import HDF5

type SpikeWaveforms
    waveforms::Array{Float64,3}
    timestamps::Array{Float64,1}
end

abstract SpikeSorterResult
abstract FeatureModel

type PCAFeatureModel <: FeatureModel
    features::Array{Float64,2}
    model::PCA
end

transform(m::PCAFeatureModel, y) = transform(m.model, y)

type FeatureSorter{T1 <: FeatureModel, T2} <: SpikeSorterResult
    waveforms::SpikeWaveforms
    featuremodel::T1
    clusterid::Array{Int64,1}
    clustermodel::T2
end

function SpikeWaveforms(waveforms::Array{Float64,2}, timestamps::Array{Float64,1})
    npoints,nspikes = size(waveforms)
    _waveforms = zeros(npoints, 1, nspikes)
    _waveforms[:,1,:] = waveforms
    SpikeWaveforms(_waveforms, timestamps)
end

type TemplateFile
	templates::Array{Float64,3}
	cinv::Array{Float64,2}
	ntemplates::Int64
	nchannels::Int64
end

type SortedFile
	templates::Array{Float64,3}
	cinv::Array{Float64,2}
	ntemplates::Int64
	nchannels::Int64
	mlseq::Array{UInt8,2}
end

type Features
	spike_width::Float64
	spikes_in_bursts::Float64
	low_isi::Float64
	dv_ratio::Float64
end

@compat function TemplateFile(fname::AbstractString;verbose::Integer=0)
	if !HDF5.ishdf5(fname)
		verbose > 0 && println("$fname is not a valid hdf5 file")
		return nothing
	end
	ff = HDF5.h5open(fname)
	if !("spikeForms" in names(ff))
		verbose > 0 && println("File does not contain any spike forms")
		close(ff)
		return nothing
	end
	templates = HDF5.read(ff, "spikeForms")
	nspikes = size(templates,3)
	nchannels = size(templates,2)
	if nchannels > 1
		cinv = HDF5.read(ff, "cinv")
	else
		cinv = zeros(1,1)
		cinv[1,1] = HDF5.read(ff, "cinv")[:][1]
	end
	close(ff)
	return TemplateFile(templates, cinv,nspikes,nchannels)
end


text:
# This function reads in the necessary Earth orientation parameters for
# rotating gps ephemeris from the ITRF14 reference frame to the GCRF
# (Geocentric Celestial Reference Frame) often just refered to as
# the J2000 reference frame (i.e. in SPICE)
function readERPs(gpsWeekStart, weekDayStart, gpsWeekEnd, weekDayEnd)
    # If weekDayStart is 0, we need the final row of last weeks eop in addition to gpsWeekStart week
    needLastWeek = (weekDayStart == 0)

    # If weekDayEnd is 6, we need the first row in next weeks eop in addition to gpsWeekEnd week
    needNextWeek = (weekDayEnd == 6)

    # Compute number of weeks in requested data
    numReqWeeks = gpsWeekEnd - gpsWeekStart + 1

    # Size EOP Matrix
    eops = zeros(numReqWeeks*7 + Int(needLastWeek) + Int(needNextWeek), 18)

    # Read in data
    lineNum = 0

    # Read in last week data if necessary
    if needLastWeek
        # Open file
        #erpFile = datadir("igs") * "\\igs" * string(gpsWeekStart - 1) * "7.erp"
        erpFile = datadir("igs", "igs" * string(gpsWeekStart - 1) * "7.erp")

        if !isfile(erpFile)
            throw(ErrorException("The file " * erpFile * " does not exist and must be downloaded from CDDIS.\n" *
            "The file can be found at 'https://cddis.nasa.gov/archive/gnss/products/'"))
        end
        f = open(erpFile)

        # Get last line in file
        lastLine = ""
        for line in readlines(f)
            if length(line) < 2
                break
            end
            lastLine = line
        end

        # Put data in Matrix
        lineNum += 1
        splitLine = split(lastLine, " "; keepempty = false)
        for i in 1:18
            eops[lineNum, i] = parse(Float64, splitLine[i])
        end
        close(f)
    end

    # Read in span of weeks
    weeks = gpsWeekStart:gpsWeekEnd
    for week in weeks
        # Open file
        #erpFile = datadir("igs") * "\\igs" * string(week) * "7.erp"
        erpFile = datadir("igs", "igs" * string(week) * "7.erp")
        if !isfile(erpFile)
            throw(ErrorException("The file 'igs"*string(week)*"7.erp' does not exist in the folder ./data/igs and must be downloaded from CDDIS.\n" *
            "The file can be found at 'https://cddis.nasa.gov/archive/gnss/products/'"))
        end
        f = open(erpFile)

        # Read in data
        for line in readlines(f)
            if length(line) < 2
                break
            elseif isdigit(line[1])
                lineNum += 1
                splitLine = split(line, " "; keepempty = false)
                for i in 1:18
                   eops[lineNum, i] = parse(Float64, splitLine[i])
                end
            end
        end
        close(f)
    end

    # Read in last week data if necessary
    if needNextWeek
        # Open file
        #erpFile = datadir("igs") * "\\igs" * string(gpsWeekEnd + 1) * "7.erp"
        erpFile = datadir("igs", "igs" * string(gpsWeekEnd + 1) * "7.erp")

        if !isfile(erpFile)
            throw(ErrorException("The file 'igs"*string(gpsWeekEnd + 1)*"7.erp' does not exist in the folder ./data/igs and must be downloaded from CDDIS.\n" *
            "The file can be found at 'https://cddis.nasa.gov/archive/gnss/products/'"))
        end
        f = open(erpFile)

        # Get first line of data in file
        firstDataLine = ""
        for line in readlines(f)
            if isdigit(line[1])
                firstDataLine = line
                break
            end
        end

        # Put data in Matrix
        lineNum += 1
        splitLine = split(firstDataLine, " "; keepempty = false)
        for i in 1:18
            eops[lineNum, i] = parse(Float64, splitLine[i])
        end
        close(f)
    end

    return eops
end

max_stars_repo_name: GrantHecht/OptimalEstimationProject


text:
import FiniteDiff
import ForwardDiff
import ReverseDiff
import Zygote

using Test

F(x, θ, deg) = (θ[1] - x[1])^deg + θ[2] * (x[2] - x[1]^deg)^deg
F_idx_val(ca) = F(ca[Val(:x)], ca[Val(:θ)], ca[Val(:deg)])
F_idx_sym(ca) = F(ca[:x], ca[:θ], ca[:deg])
F_view_val(ca) = F(@view(ca[Val(:x)]), @view(ca[Val(:θ)]), ca[Val(:deg)])
F_view_sym(ca) = F(@view(ca[:x]), @view(ca[:θ]), ca[:deg])
F_prop(ca) = F(ca.x, ca.θ, ca.deg)

ca = ComponentArray(x = [1, 2], θ = [1.0, 100.0], deg = 2)
truth = [-400, 200]

@testset "$(nameof(F_))" for F_ in (F_idx_val, F_idx_sym, F_view_val, F_view_sym, F_prop)
    finite = FiniteDiff.finite_difference_gradient(ca -> F_(ca), ca).x
    @test finite ≈ truth

    forward = ForwardDiff.gradient(ca -> F_(ca), ca).x
    @test forward ≈ truth

    reverse = ReverseDiff.gradient(ca -> F_(ca), ca).x
    @test reverse ≈ truth

    zygote_full = Zygote.gradient(ca -> F_(ca), ca)[1]
    if F_ == F_prop && VERSION < v"1.3"
        @test_broken zygote_full.x ≈ truth
    else
        @test zygote_full.x ≈ truth
    end
end


text:
using Reexport

@reexport module InterpolationAndGrids
    using Interpolations

    include("findall_nan_chunks.jl")
    include("fill_nans.jl")

    include("grids/RectangularGrid.jl")
    include("methods/interpolation.jl")
    include("interpolate_and_bin.jl")
    include("binning.jl")
end


text:
function __init__()
    if !haskey(ENV, "SCALAR_FDFD_COLOR_THEME")
        SCALAR_FDFD_COLOR_THEME = :default
    else
        SCALAR_FDFD_COLOR_THEME = Symbol(ENV["SCALAR_FDFD_COLOR_THEME"])
    end

    global COLOR_SCHEME = SCALAR_FDFD_COLOR_THEME
    set_colors_by_scheme()
end

function set_colors_by_scheme()
    if COLOR_SCHEME == :default
        global BAND_COLOR = :lightgrey
        global BAND_WIDTH = 1.2
        global BAND_STYLE = :solid
        global GAP_COLOR = :lightgrey
        global GAP_WIDTH = 1
        global GAP_STYLE = :dash
        global DISPERSION_WIDTH = 3
        global DISPERSION_STYLE = :solid
    elseif COLOR_SCHEME == :dark
        global BAND_COLOR = :lightgrey
        global BAND_WIDTH = 1.2
        global BAND_STYLE = :solid
        global GAP_COLOR = :lightgrey
        global GAP_WIDTH = 1
        global GAP_STYLE = :dash
        global DISPERSION_WIDTH = 3
        global DISPERSION_STYLE = :solid
    elseif COLOR_SCHEME == :ggplot2
        global BAND_COLOR = :lightgrey
        global BAND_WIDTH = 1.2
        global BAND_STYLE = :solid
        global GAP_COLOR = :lightgrey
        global GAP_WIDTH = 1
        global GAP_STYLE = :dash
        global DISPERSION_WIDTH = 3
        global DISPERSION_STYLE = :solid
    elseif COLOR_SCHEME == :juno
        global BAND_COLOR = :lightgrey
        global BAND_WIDTH = 1.2
        global BAND_STYLE = :solid
        global GAP_COLOR = :lightgrey
        global GAP_WIDTH = 1
        global GAP_STYLE = :dash
        global DISPERSION_WIDTH = 3
        global DISPERSION_STYLE = :solid
    elseif COLOR_SCHEME == :lime
        global BAND_COLOR = :lightgrey
        global BAND_WIDTH = 1.2
        global BAND_STYLE = :solid
        global GAP_COLOR = :lightgrey
        global GAP_WIDTH = 1
        global GAP_STYLE = :dash
        global DISPERSION_WIDTH = 3
        global DISPERSION_STYLE = :solid
    elseif COLOR_SCHEME == :orange
        global BAND_COLOR = :lightgrey
        global BAND_WIDTH = 1.2
        global BAND_STYLE = :solid
        global GAP_COLOR = :lightgrey
        global GAP_WIDTH = 1
        global GAP_STYLE = :dash
        global DISPERSION_WIDTH = 3
        global DISPERSION_STYLE = :solid
    elseif COLOR_SCHEME == :sand
        global BAND_COLOR = :lightgrey
        global BAND_WIDTH = 1.2
        global BAND_STYLE = :solid
        global GAP_COLOR = :lightgrey
        global GAP_WIDTH = 1
        global GAP_STYLE = :dash
        global DISPERSION_WIDTH = 3
        global DISPERSION_STYLE = :solid
    elseif COLOR_SCHEME == :solarized
        global BAND_COLOR = :lightgrey
        global BAND_WIDTH = 1.2
        global BAND_STYLE = :solid
        global GAP_COLOR = :lightgrey
        global GAP_WIDTH = 1
        global GAP_STYLE = :dash
        global DISPERSION_WIDTH = 3
        global DISPERSION_STYLE = :solid
    elseif COLOR_SCHEME == :solarized_light
        global BAND_COLOR = :lightgrey
        global BAND_WIDTH = 1.2
        global BAND_STYLE = :solid
        global GAP_COLOR = :lightgrey
        global GAP_WIDTH = 1
        global GAP_STYLE = :dash
        global DISPERSION_WIDTH = 3
        global DISPERSION_STYLE = :solid
    elseif COLOR_SCHEME == :wong
        global BAND_COLOR = :lightgrey
        global BAND_WIDTH = 1.2
        global BAND_STYLE = :solid
        global GAP_COLOR = :lightgrey
        global GAP_WIDTH = 1
        global GAP_STYLE = :dash
        global DISPERSION_WIDTH = 3
        global DISPERSION_STYLE = :solid
    elseif COLOR_SCHEME == :wong2
        global BAND_COLOR = :lightgrey
        global BAND_WIDTH = 1.2
        global BAND_STYLE = :solid
        global GAP_COLOR = :lightgrey
        global GAP_WIDTH = 1
        global GAP_STYLE = :dash
        global DISPERSION_WIDTH = 3
        global DISPERSION_STYLE = :solid
    else
        global BAND_COLOR = :lightgrey
        global BAND_WIDTH = 1.2
        global BAND_STYLE = :solid
        global GAP_COLOR = :lightgrey
        global GAP_WIDTH = 1
        global GAP_STYLE = :dash
        global DISPERSION_WIDTH = 3
        global DISPERSION_STYLE = :solid
    end
    return nothing
end

"""
    cmapc, cmapk, cmapsim1, cmapsim2, n_mult, F_sign = fix_colormap(theme)
"""
function fix_colormap(theme)

    if theme ∈ [:dark, :juno]
        cmapc=:bkr
        cmapk=:ice
        cmapsim1=:dimgray
        cmapsim2=:bky
    elseif theme ∈ [:solarized]
        cmapc=:bky
        cmapk=:solar
        cmapsim1=:viridis
        cmapsim2=:bky
    elseif theme ∈ [:orange]
        cmapc=:bky
        cmapk=:haline
        cmapsim1=:inferno
        cmapsim2=:bky
    elseif theme ∈ [:lime]
        cmapc=:bky
        cmapk=:haline
        cmapsim1=:inferno
        cmapsim2=:bky
    elseif theme ∈ [:solarized_light]
        cmapc=:YlOrBr
        cmapk=:Greys
        cmapsim1=:YlOrBr
        cmapsim2=:RdYlBu
    else
        cmapc=:RdBu
        cmapk=:Greys
        cmapsim1=:Greys
        cmapsim2=:RdGy
    end

    return cmapc, cmapk, cmapsim1, cmapsim2
end


text:
function generate_bloobs(n_features = 2, n_centers = 3, n_samples = 100)
    centers = @parallel (hcat) for i = 1:n_features
        rand(Uniform(-10, 10), n_centers)
    end

    cluster_std = @parallel (hcat) for i = 1:n_features
      rand(Uniform(0.5, 4), n_centers)
    end

    n_samples_per_center = ones(Integer, n_centers) * int(n_samples / n_centers)

    for i = 1:(n_centers % n_samples)
        n_samples_per_center[i] += 1
    end

    X = @parallel (hcat) for i = 1:n_centers
      rand(MvNormal(ones(n_features) .* centers[i], eye(n_features) .* cluster_std[i]), n_samples_per_center[i])
    end

    X = X';

    Y = @parallel (vcat) for i = 1:n_centers
      y = ones(n_samples_per_center[i]) * i
    end

    ids = [1:size(X)[1]]
    shuffle!(ids)

    X = X[ids,:]
    Y = Y[ids];

    return (X, Y)
end

function generateBarsDataset(img_size::Int, noise_level::Float64, num_per_mixture::Array{Float64}, num_group::Int, num_data::Int)

  numdim   = img_size^2;
  numbars  = img_size*2

  mix_theta = zeros(numdim, numbars)

  # add horizontal bars
  for i in 1:img_size
    img = ones(img_size, img_size) .* noise_level / img_size^2;
    img[:,i] = img[:,i] + 1/img_size;
    mix_theta[:,i] = img[:] / sum(img[:]);
  end

  # add vertical bars
  for i = 1:img_size
    img = ones(img_size, img_size) .* noise_level / img_size^2;
    img[i,:] = img[i,:] + 1/img_size;
    mix_theta[:,img_size+i] = img[:]/sum(img[:]);
  end

  # generate samples
  # Note: strange behavior of cumsum!!!
  cumbar = cumsum(num_per_mixture, 2)

  samples = Vector{Vector{Int}}(num_group)

  for g in 1:num_group

      # random number of mixture components
      n_bars = 1 + sum(rand() .> cumbar)
      k = randperm(numbars)[1:n_bars]

      # get weigths
      theta = mean(mix_theta[:,k], 2)

      # get samples (e.g. words)
      samples[g] = vec( 1 + sum(repmat(rand(num_data)', numdim, 1) .> repmat(cumsum(theta, 1), 1, num_data), 1) )
  end

  # return samples and bars
  return (samples, mix_theta)

end


text:
function getGrid(minongrid, maxongrid, GridPoints, method)

#NB: This code borrows some code from Chris Carroll of John's Hopkins University
#His website is here http://www.econ2.jhu.edu/people/ccarroll/ and his
#notes on solving dynamic models are very good - they can be found
#here: http://www.econ2.jhu.edu/people/ccarroll/SolvingMicroDSOPs.pdf

#% ------------------------------------------------------------------------
# Explanation

# We need to a grid from a to b. A basic approach could involve spacing out
# the point equally. The following line of code would achieve this:
#     grid= linspace(a, b, GridPoints);

# If we want to space them out so that the growth rate of the distance
# between spaces is equal we achieve this by spacing out the logs of grid
# equally. We could achieve this by the following line of code:
#     grid= exp(linspace(log(a), log(b), GridPoints));

# This is problematic if a<=0. The approach taken here is to get a grid
# from log(1) to log(b - a + 1), exponentiate and then subrtact 1
# so that we have a grid from 0 to b-a. If we have add a to each point we
# get the desired result - a log-spaced grid from a to b

span = maxongrid - minongrid;     # b - a

if method=="equalsteps"
    grid= linspace(min, span, GridPoints);
elseif method== "logsteps"
  loggrid = linspace(log(1), log(1+span), GridPoints);
  grid = exp(loggrid)-1;
elseif method=="3logsteps"
  loggrid = linspace(log(1+log(1+log(1))), log(1+log(1+log(1+span))), GridPoints);
  grid = exp(exp(exp(loggrid)-1)-1)-1;
elseif method=="5logsteps"
  loggrid = linspace(log(1+log(1+log(1+log(1+log(1))))), log(1+log(1+log(1+log(1+log(1+span))))), GridPoints);
  grid = exp(exp(exp(exp(exp(loggrid)-1)-1)-1)-1)-1;
elseif method=="10logsteps"
  loggrid = linspace(log(1+log(1+log(1+log(1+log(1+log(1+log(1+log(1+log(1+log(1)))))))))), log(1+log(1+log(1+log(1+log(1+log(1+log(1+log(1+log(1+log(1+span)))))))))), GridPoints);
  grid = exp(exp(exp(exp(exp(exp(exp(exp(exp(exp(loggrid)-1)-1)-1)-1)-1)-1)-1)-1)-1)-1;
else
    error("Error in getgrid. You have entered an invalid method for choosing the distance between grid points. Method must be one of equalsteps, logsteps, 3logsteps, 5logsteps or 10logsteps.");
end

grid = grid + minongrid*ones(Float64,GridPoints,1);

return grid;

end

max_stars_repo_name: floswald/ucl-econ-julia

max_issues_repo_name: floswald/ucl-econ-julia

max_forks_repo_name: floswald/ucl-econ-julia


text:
const pubchem_features = [:acceptor, :anion, :cation, :donor, :hydrophobe, :rings, :aromaticrings]

"""
    feat = atoms_to_feature(atoms, dirs=nothing)

Combines Gaussian distributions specified by `atoms` to create a feature repressented by a
single Gaussian. The feature is centered at the center of mass of the distributions, and the
width `σ` of the feature is the average of the distance of individual features from the center
added to their widths.

Geometric constraints for the feature can be specified by `dirs`. The returned feature has no
direction by default.
"""
function atoms_to_feature(atoms::AbstractVector{<:AbstractIsotropicGaussian}, dirs=nothing)
    if length(atoms)==1
        μ = atoms[1].μ
        σ = atoms[1].σ
        ϕ = atoms[1].ϕ
    else
        μ = center_of_mass(atoms)
        # σ = mean([norm(μ-a.μ) + a.σ for (i,a) in enumerate(atoms)])
        σ = sum([norm(μ.-a.μ) for (i,a) in enumerate(atoms)])/length(atoms)
        ϕ = sum([a.ϕ for a in atoms])/length(atoms)
    end
    if isnothing(dirs)
        dirs = typeof(μ)[]
    end

    return IsotropicGaussian(μ, σ, ϕ, dirs)
end
atoms_to_feature(gmm::AbstractSingleGMM, dirs=nothing) = atoms_to_feature(gmm.gaussians, dirs)
atoms_to_feature(mol::UndirectedGraph, nodeset, dirs=nothing; kwargs...) = atoms_to_feature(MolGMM(nodesubgraph(mol,nodeset); kwargs...), dirs)
atoms_to_feature(submol::SubgraphView, dirs=nothing; kwargs...) = atoms_to_feature(submol.graph, nodeset(submol), dirs; kwargs...)

#TODO: fit single Gaussian to points sampled from sum of the features?
function combine_features(feats::AbstractVector{<:AbstractIsotropicGaussian}, weights=[f.ϕ for f in feats])
    ϕ = sum(weights)
    μ = center_of_mass(feats, weights)
    σ = sum([f.σ * weights[i]/ϕ for (i,f) in enumerate(feats)])
    return IsotropicGaussian(μ, σ, ϕ)
end


text:
export CompositeFactor, @composite

abstract type CompositeFactor <: FactorNode end

"""
The `@composite` macro allows for defining custom (composite) nodes. Composite
nodes allow for implementating of custom update rules that may be computationally
more efficient or convenient. A composite node can be defined with or without an
internal model. For detailed usage instructions we refer to the `composite_nodes` demo.
"""
macro composite(name::Symbol, exposed_vars::Expr, model::Expr)
    (exposed_vars.head == :tuple) || error("Exposed variables should be passed as Tuple")

    exposed_var_arguments = join(["$varname::Variable" for varname in exposed_vars.args], ", ")
    n_vars = length(exposed_vars.args)

    # Code for constructing interfaces
    interface_definitions = ""
    for idx = 1:n_vars
        varname = exposed_vars.args[idx]
        interface_definitions *= "self.i[:$varname] = self.interfaces[$idx] = ForneyLab.associate!(Interface(self), $varname)\n"
    end

    # Code for constructing exposed variables in inner graph
    exposed_var_definitions = ""
    terminal_definitions = ""
    for idx = 1:n_vars
        varname = exposed_vars.args[idx]
        exposed_var_definitions *= "local $varname = Variable(id=:$varname)\n"
        terminal_definitions *= "push!(self.terminals, Terminal($varname, self.interfaces[$idx], id=:$varname))\n"
        terminal_definitions *= "self.interface2terminal[self.interfaces[$idx]] = self.terminals[$idx]\n"
    end

    expr = parse("""
    mutable struct $name <: CompositeFactor
        id::Symbol
        interfaces::Vector{Interface}
        i::Dict{Symbol, Interface}
        inner_graph::FactorGraph
        terminals::Vector{Terminal}
        interface2terminal::Dict{Interface,Terminal}

        function $name($exposed_var_arguments; id=ForneyLab.generateId($name))
            self = new(id, Array{Interface}(undef, $n_vars), Dict{Symbol,Interface}())
            outer_graph = currentGraph()
            ForneyLab.addNode!(outer_graph, self)
            $interface_definitions

            # Build internal graph
            self.inner_graph = FactorGraph()
            self.terminals = Terminal[]
            self.interface2terminal = Dict{Interface,Terminal}()
            let
                $exposed_var_definitions
                $model
                $terminal_definitions
            end
            setCurrentGraph(outer_graph)

            return self
        end
    end
    """)

    return esc(expr)
end


text:
@sumProductRule(:node_type     => Nonlinear{Extended},
                :outbound_type => Message{GaussianMeanVariance},
                :inbound_types => (Nothing, Message{Gaussian}),
                :name          => SPNonlinearEOutNG)

@sumProductRule(:node_type     => Nonlinear{Extended},
                :outbound_type => Message{Gaussian},
                :inbound_types => (Message{Gaussian}, Nothing),
                :name          => SPNonlinearEIn1GG)

mutable struct SPNonlinearEOutNGX <: SumProductRule{Nonlinear{Extended}} end
outboundType(::Type{SPNonlinearEOutNGX}) = Message{GaussianMeanVariance}
function isApplicable(::Type{SPNonlinearEOutNGX}, input_types::Vector{<:Type})
    total_inputs = length(input_types)
    (total_inputs > 2) || return false
    (input_types[1] == Nothing) || return false

    for input_type in input_types[2:end]
        matches(input_type, Message{Gaussian}) || return false
    end

    return true
end

mutable struct SPNonlinearEInGX <: SumProductRule{Nonlinear{Extended}} end
outboundType(::Type{SPNonlinearEInGX}) = Message{Gaussian}
function isApplicable(::Type{SPNonlinearEInGX}, input_types::Vector{<:Type})
    total_inputs = length(input_types)
    (total_inputs > 2) || return false
    (input_types[1] != Nothing) || return false

    nothing_inputs = 0
    gaussian_inputs = 0
    for input_type in input_types
        if input_type == Nothing
            nothing_inputs += 1
        elseif input_type << Message{Gaussian}
            gaussian_inputs += 1
        end
    end

    return (nothing_inputs == 1) && (gaussian_inputs == total_inputs - 1)
end

mutable struct MNonlinearEInGX <: MarginalRule{Nonlinear{Extended}} end
function isApplicable(::Type{MNonlinearEInGX}, input_types::Vector{<:Type})
    total_inputs = length(input_types)
    (total_inputs > 2) || return false
    (input_types[1] == Nothing) || return false # Indicates marginalization over outbound variable

    for input_type in input_types[2:end]
        matches(input_type, Message{Gaussian}) || return false
    end

    return true
end


text:
#####
##### Dummy runs are used to ensure the absence of runtime errors in the code
##### before launching a training session
#####

using Setfield

dummy_run_mcts(p::MctsParams) = @set p.num_iters_per_turn = 4

dummy_run_player(p::Benchmark.Player) = p
dummy_run_player(p::Benchmark.Full) = @set p.params = dummy_run_mcts(p.params)
dummy_run_player(p::Benchmark.MctsRollouts) = @set p.params = dummy_run_mcts(p.params)

function dummy_run_sim(s::SimParams)
  s = @set s.num_games = ConstSchedule(1)
  s = @set s.num_workers = ConstSchedule(1)
  s = @set s.batch_size = ConstSchedule(1)
  s = @set s.use_gpu = false
  return s
end

# Returned modified parameters where all num_games fields are set to 1.
# The number of iterations is set to 2.
function dummy_run_params(params)
  params = @set params.self_play.sim = dummy_run_sim(params.self_play.sim)
  params = @set params.self_play.mcts = dummy_run_mcts(params.self_play.mcts)
  if !isnothing(params.arena)
    params = @set params.arena.sim = dummy_run_sim(params.arena.sim)
    params = @set params.arena.mcts = dummy_run_mcts(params.arena.mcts)
  end
  params = @set params.learning.max_batches_per_checkpoint = 2
  params = @set params.learning.num_checkpoints = min(params.learning.num_checkpoints, 2)
  params = @set params.num_iters = 2
  return params
end

function dummy_run_evaluation(d::Benchmark.Single)
  d = @set d.sim = dummy_run_sim(d.sim)
  d = @set d.player = dummy_run_player(d.player)
  return d
end

function dummy_run_evaluation(d::Benchmark.Duel)
  d = @set d.sim = dummy_run_sim(d.sim)
  d = @set d.player = dummy_run_player(d.player)
  d = @set d.baseline = dummy_run_player(d.baseline)
  return d
end

dummy_run_benchmark(benchmark) = map(dummy_run_evaluation, benchmark)

function dummy_run_experiment(e::Experiment)
  params = dummy_run_params(e.params)
  benchmark = dummy_run_benchmark(e.benchmark)
  return Experiment(e.name, e.gspec, params, e.mknet, e.netparams, benchmark)
end

"""

    dummy_run(experiment; [dir, nostdout])

Launch a training session where hyperparameters are altered so that training
finishes as quickly as possible.

This is useful to ensure the absence of runtime errors before
a real training session is started.
"""
function dummy_run(experiment::Experiment; dir=nothing, nostdout=false)
  experiment = dummy_run_experiment(experiment)
  isnothing(dir) && (dir = "sessions/dummy-$(experiment.name)")
  rm(dir, force=true, recursive=true)
  session = Session(experiment; dir, nostdout)
  resume!(session)
  return
end


text:
using POMDPs, Crux, Flux, POMDPGym

## Cartpole - V0
mdp = GymPOMDP(:CartPole, version = :v0)
as = actions(mdp)
S = state_space(mdp)

A() = DiscreteNetwork(Chain(Dense(Crux.dim(S)..., 64, relu), Dense(64, 64, relu), Dense(64, length(as))), as)
V() = ContinuousNetwork(Chain(Dense(Crux.dim(S)..., 64, relu), Dense(64, 64, relu), Dense(64, 1)))

# Solve with REINFORCE (~2 seconds)
𝒮_reinforce = REINFORCE(π=A(), S=S, N=10000, ΔN=500, a_opt=(epochs=5,))
@time π_reinforce = solve(𝒮_reinforce, mdp)

# Solve with A2C (~8 seconds)
𝒮_a2c = A2C(π=ActorCritic(A(), V()), S=S, N=10000, ΔN=500)
@time π_a2c = solve(𝒮_a2c, mdp)

# Solve with PPO (~15 seconds)
𝒮_ppo = PPO(π=ActorCritic(A(), V()), S=S, N=10000, ΔN=500)
@time π_ppo = solve(𝒮_ppo, mdp)

# Solve with DQN (~12 seconds)
𝒮_dqn = DQN(π=A(), S=S, N=10000)
𝒮_dqn.c_opt
@time π_dqn = solve(𝒮_dqn, mdp)

# Plot the learning curve
p = plot_learning([𝒮_reinforce, 𝒮_a2c, 𝒮_ppo, 𝒮_dqn], title = "CartPole-V0 Training Curves", labels = ["REINFORCE", "A2C", "PPO", "DQN"])

# Produce a gif with the final policy
gif(mdp, π_ppo, "cartpole_policy.gif", max_steps=100)

## Optional - Save data for imitation learning
# using BSON
# s = Sampler(mdp, 𝒮_dqn.π, max_steps=100, required_columns=[:t])
#
# data = steps!(s, Nsteps=10000)
# sum(data[:r])/100
# data[:expert_val] = ones(Float32, 1, 10000)
# data[:a]
#
# data = ExperienceBuffer(data)
# BSON.@save "examples/il/expert_data/cartpole.bson" data


text:
using Distributed

const Ncpu = 1
const totalStep = 1e7
const Repeat = 4

addprocs(Ncpu)

@everywhere using QuantumStatistics, LinearAlgebra, Random, Printf, StaticArrays, Statistics, BenchmarkTools, InteractiveUtils, Parameters

# claim all globals to be constant, otherwise, global variables could impact the efficiency
@everywhere const kF, m = 1.919, 0.5
@everywhere const β = 25.0 / kF^2

@everywhere const Nk,Nt = 8, 8
@everywhere const extQ = [@SVector [q, 0.0, 0.0] for q in range(0.0, stop=3.0 * kF, length=Nk)]
@everywhere const extT = range(0.0, stop=β, length=Nt)

@everywhere const obs1, obs2 = zeros(Float64, (Nk,Nk,Nt,Nt)), zeros(Float64, (Nk,Nk,Nt,Nt))

@everywhere function WRPA(τ,q)
    # temporarily using a toy model
    g=1.0

    factor=1.0
    q2=dot(q,q)
    if τ<-β
        τ=τ+2*β
    elseif τ<0
        τ=τ+β
        factor=-1.0
    end
    sq2g=sqrt(q2+g)
    return -factor*4*π/(q2+g)/sq2g*(exp(-sq2g*τ)+exp(-sq2g*(β-τ)))/(1-exp(-sq2g*β))
end

@everywhere function integrand(config)
    if config.curr.id == 1
        return 1.0
    elseif config.curr.id == 2
            return eval2(config)
        else
            return 0.0
    end
end

@everywhere function eval2(config)
    T, K, ExtQ, ExtT, Theta = config.var[1], config.var[2], config.var[3], config.var[4], config.var[5]
        # In case the compiler is too stupid, it is a good idea to explicitly specify the type here
    q = K[1]
    t = T[1]
    θ = Theta[1]

    k1 = extQ[ExtQ[1]] # external momentum
    k2 = extQ[ExtQ[2]] # external momentum
    t1 = extT[ExtT[3]] # external momentum
    t2 = extT[ExtT[4]] # external momentum

    k2 = norm(k2) * @SVector[cos(θ), sin(θ), 0]

    ω1 = (dot(q-k1, q-k1) - kF^2) * β
    τ1 = (t1-t)/β
    g1 = Spectral.kernelFermiT(τ1, ω1)

    ω2 = (dot(q-k2, q-k2) - kF^2) * β
    τ2 = (β-t1-t)/β
    g2 = Spectral.kernelFermiT(τ2, ω2)

    W1 = WRPA(t, q-k1-k2)
    W2 = WRPA(t-t1-t2,q)

    legendre=1.0
    factor = 1.0
    return g1 * g2 * W1 * W2 * factor * legendre
end

@everywhere function measure(config)
    diag = config.curr
    factor = 1.0 / diag.reWeightFactor
    extqidx = [config.var[3][1],config.var[3][2]]
    exttidx = [config.var[4][1],config.var[4][2]]
    extidx = CartesianIndex(Tuple(x for x in vcat(extqidx, exttidx)))
    if diag.id == 1
        obs1[extidx] += factor
    elseif diag.id == 2
        weight = integrand(config)
        obs2[extidx] += weight / abs(weight) * factor
    else
        return
    end
end

@everywhere function MC(totalStep, pid)
    rng = MersenneTwister(pid)

    K = MonteCarlo.FermiK(3, kF, 0.2 * kF, 10.0 * kF)
    T = MonteCarlo.Tau(β, β / 2.0)
    Theta = MonteCarlo.Angle(0.5)
    ExtQ = MonteCarlo.Discrete(1, length(extQ)) # external variable is specified
    ExtT = MonteCarlo.Discrete(1, length(extT)) # external variable is specified
    diag1 = MonteCarlo.Diagram(1, 0, [0, 0, 2, 2, 0])
    diag2 = MonteCarlo.Diagram(2, 1, [1, 1, 2, 2, 1])

    config = MonteCarlo.Configuration(totalStep, (diag1, diag2), (T, K, ExtQ, ExtT, Theta); pid=pid, rng=rng)

    # @code_warntype MonteCarlo.Configuration(totalStep, (diag1, diag2), (T, K, Ext); pid=pid, rng=rng)
    # @code_warntype MonteCarlo.Diagram(2, 1, [2, 1, 1], zeros(Float64, Ext.size))
    # @code_warntype MonteCarlo.increaseOrder(config, integrand)
    # @code_warntype integrand(config)
    # @code_warntype eval2(config)
    # @code_warntype measure(config)
    # @code_lowered MonteCarlo.changeVar(config, integrand)
    # exit()

    MonteCarlo.montecarlo(config, integrand, measure)

    return obs2 / sum(obs1) * ExtQ.size * ExtT.size
end

function run(repeat, totalStep)
    if Ncpu > 1
        observable = pmap((x) -> MC(totalStep, rand(1:10000)), 1:repeat)
    else
        observable = map((x) -> MC(totalStep, rand(1:10000)), 1:repeat)
    end

    obs = mean(observable)
    obserr = std(observable) / sqrt(length(observable))

    for (idx, t) in enumerate(extT)
        t = t[1]
        @printf("%10.6f  %10.6f ± %10.6f\n", t, obs[idx], obserr[idx])
    end

end

# @btime run(1, 10)
run(Repeat, totalStep)
# @time run(Repeat, totalStep)


text:
# ### Converter temperatura
# Leia um valor de temperatura em graus Celsius e converta-o em graus Fahrenheit e Kelvin
# * F = ((9C) / 5 + 32) ou F = (9C + 160) /5
# * K = C + 273

println("Informe a temperatura em Celsius")
temperatura = parse(Float16, readline())

temperatura += 273
println("Temperatura em Kelvin ",temperatura)

temperatura -= 273
temperatura = (( (temperatura * 9) / 5) + 32)
println("Temperatura em graus Fahrenheit ",temperatura)

max_stars_repo_name: codinginbrazil/functional-programming

max_issues_repo_name: codinginbrazil/functional-programming

max_forks_repo_name: codinginbrazil/functional-programming


text:
some_elm(::Type{String}) = ""

"""
Saving a String
"""
compress_then_write(b::Vector{String}, io) = begin
    # TODO compare whether StringArray is better
    # return compress_then_write(StringArray(b), io)

    # write the string one by one
    # do a Run-length encoding (RLE)
    previous_b = b[1]
    cnt = 1
    lens = Int[]
    str_lens = Int[]
    for i = 2:length(b)
        if b[i] != previous_b
            push!(str_lens, write(io, previous_b))
            push!(lens, cnt)
            cnt = 0
            previous_b = b[i]
        end
        cnt += 1
    end

    # reach the end: two situation
    # 1) it's a new element, so write it
    # 2) it's an existing element. Also write it
    push!(str_lens, write(io, previous_b))
    push!(lens, cnt)

    @assert sum(lens) == length(b)

    str_lens_compressed = Blosc.compress(Vector{UInt32}(str_lens))
    str_lens_bytes = write(io, str_lens_compressed)

    lens_compressed = Blosc.compress(Vector{UInt64}(lens))
    rle_bytes = write(io, lens_compressed)

    # return metadata
    return (
        string_compressed_bytes = sum(str_lens),
        string_len_bytes = str_lens_bytes,
        rle_bytes = rle_bytes,
        rle_len = length(str_lens),
        type = String,
        len = max(sum(str_lens), str_lens_bytes, rle_bytes),
    )
end

# load a string column
"""
    metadata should consists of length, compressed byte size of string-lengths,
    string content lengths
"""
column_loader!(_, ::Type{String}, io, metadata) = begin
    column_loader(String, io, metadata)
end

column_loader(::Type{String}, io, metadata) = begin
    buffer = Vector{UInt8}(undef, metadata.string_compressed_bytes)
    readbytes!(io, buffer, metadata.string_compressed_bytes)
    #return String(buffer)

    # read the string-lengths
    buffer2 = Vector{UInt8}(undef, metadata.string_len_bytes)
    readbytes!(io, buffer2, metadata.string_len_bytes)

    buffer3 = Vector{UInt8}(undef, metadata.rle_bytes)
    readbytes!(io, buffer3, metadata.rle_bytes)

    counts = Blosc.decompress(UInt64, buffer3)

    str_lens = Blosc.decompress(UInt32, buffer2)

    #return (String(buffer), str_lens, counts)

    lengths = inverse_rle(str_lens, counts)
    offsets = inverse_rle(vcat(0, cumsum(str_lens[1:end-1])), counts)

    #res = StringArray{String, 1}(buffer, vcat(1, cumsum(Blosc.decompress(UInt64, buffer3))[1:end-1]) .-1,  )
    res = StringArray{String,1}(buffer, offsets, lengths)
end


text:
using OdinSon

#TODO: add support for HSV so I can do the equivalent of HSV(rand())
render(Canvas([PolyLine(rand(2, 2), style=Style(stroke=RGB(rand(3)...))) for i in 1:100]))

max_stars_repo_name: gabrielgellner/OdinSon

max_issues_repo_name: gabrielgellner/OdinSon

max_forks_repo_name: gabrielgellner/OdinSon


text:
# Machines store data for training as *arguments*. Eg, in the
# supervised case, the first two arguments are always `X` and `y`. An
# argument could be raw data (such as a table) or, in the case of a
# learning network, a "promise" of data (aka "dynamic" data) -
# formally a `Node` object, which is *callable*. For uniformity of
# interface, even raw data is stored in a wrapper that can be called
# to return the object wrapped. The name of the wrapper type is
# `Source`, as these constitute the source nodes of learning networks.

# Here we define the `Source` wrapper, as well as some methods they
# will share with the `Node` type for use in learning networks.

## SOURCE TYPE

abstract type AbstractNode <: MLJType end
abstract type CallableReturning{K} end # scitype for sources and nodes

mutable struct Source <: AbstractNode
    data     # training data
    scitype::DataType
end

"""
    Xs = source(X=nothing)

Define, a learning network `Source` object, wrapping some input data
`X`, which can be `nothing` for purposes of exporting the network as
stand-alone model. For training and testing the unexported
network, appropriate vectors, tables, or other data containers are
expected.

The calling behaviour of a `Source` object is this:

    Xs() = X
    Xs(rows=r) = selectrows(X, r)  # eg, X[r,:] for a DataFrame
    Xs(Xnew) = Xnew

See also: [`@from_network`](@ref], [`sources`](@ref),
[`origins`](@ref), [`node`](@ref).

"""
source(X) = Source(X, scitype(X))
source() = source(nothing)

source(Xs::Source; args...) = Xs

ScientificTypes.scitype(X::Source) = CallableReturning{X.scitype}
ScientificTypes.elscitype(X::Source) = X.scitype
nodes(X::Source) = [X, ]
Base.isempty(X::Source) = X.data === nothing
nrows_at_source(X::Source) = nrows(X.data)

color(::Source) = :yellow

# make source nodes callable:
function (X::Source)(; rows=:)
    rows == (:) && return X.data
    return selectrows(X.data, rows)
end
(X::Source)(Xnew) = Xnew

# return a string of diagnostics for the call `X(input...; kwargs...)`
diagnostic_table_sources(X::AbstractNode) =
    "Learning network sources:\n"*
    "source\tscitype\n"*
    "-------------------------------------------\n"*
    reduce(*, ("$s\t$(scitype(s()))\n" for s in sources(X)))

function diagnostics(X::AbstractNode, input...; kwargs...)
    raw_args = map(X.args) do arg
        arg(input...; kwargs...)
    end
    _sources = sources(X)
    scitypes = scitype.(raw_args)
    mach = X.machine
    model = mach.model
    _input = input_scitype(model)
    _target = target_scitype(model)
    _output = output_scitype(model)

    table1 = "Incoming data:\n"*
    "arg of $(X.operation)\tscitype\n"*
    "-------------------------------------------\n"*
    reduce(*, ("$(X.args[j])\t$(scitypes[j])\n" for j in eachindex(X.args)))

    table2 =  diagnostic_table_sources(X)
    return """
    Model ($model):
    input_scitype = $_input
    target_scitype =$_target
    output_scitype =$_output

    $table1
    $table2"""
end

"""
    rebind!(s, X)

Attach new data `X` to an existing source node `s`. Not a public
method.

"""
function rebind!(s::Source, X)
    s.data = X
    s.scitype = scitype(X)
    return s
end

origins(s::Source) = [s,]

## DISPLAY FOR SOURCES AND OTHER ABSTRACT NODES

# show within other objects:
function Base.show(stream::IO, object::AbstractNode)
    str = simple_repr(typeof(object))
    show_handle(object) && (str *= " $(handle(object))")
    if false
        printstyled(IOContext(stream, :color=> SHOW_COLOR),
                    str, bold=false, color=:blue)
    else
        print(stream, str)
    end
    return nothing
end
show_handle(::Source) = true

# show when alone:
function Base.show(stream::IO, ::MIME"text/plain", source::Source)
    show(stream, source)
    print(stream, " \u23CE `$(elscitype(source))`")
    return nothing
end