text:
import ThirdWay: ActivationPoint, Event
import Base.Collections.heappop!

facts("ActivationPoint") do
    context("jitter is randomized") do
        jitters = Set([ActivationPoint(0.0, 1).jitter for _ in 1:10])

        @fact length(jitters) > 1 --> true
    end

    context("(+) increments the time") do
        ap = ActivationPoint(0.0, 1)
        @fact (ap + 1.0).time == ap.time + 1.0 --> true
    end

    context("ordering") do
        @fact ActivationPoint(0.0, 0) < ActivationPoint(1.0, 0) --> true
        @fact ActivationPoint(0.0, 2) > ActivationPoint(0.0, 1) --> true
        @fact ActivationPoint(5.0, 5) > ActivationPoint(3.0, 1) --> true

        @fact ActivationPoint(1.0, 1, 5) < ActivationPoint(1.0, 1, 10) --> true
    end

    context("rejittering") do
        ap = ActivationPoint(10.0, 10)

        rejitters = Set([ActivationPoint(ap).jitter for _ in 1:10])

        @fact length(rejitters) > 1 --> true
    end
end

facts("Event") do
    context("uses ActivationPoint as an ordering") do
        a = Event(ActivationPoint(5.0, 5), println)
        b = Event(ActivationPoint(3.0, 3), println)

        @fact a > b --> true
    end
end

noop(env) = env

type MockCallable
    x::Int
end

Base.call(mc::MockCallable, env, sched::Schedule) = mc.x

type HelloAgent
    name::AbstractString
end

Base.call(agent::HelloAgent, env, sched) = push!(env, agent.name)

facts("Schedule") do
    context("push! onto heap ordered queue") do
        schedule = Schedule()

        times = [1.0, 5.0, 100.0, 33.0]
        orderings = [5, 3, 1, 2]

        all_pairs = Vector{Tuple{Float64, Int64}}()
        for t in times, o in orderings
            push!(schedule, ActivationPoint(t, o)) do env, _
                push!(env, (t, o))
            end
            push!(all_pairs, (t, o))
        end

        env = Vector{Tuple{Float64, Int64}}()
        while !isempty(schedule)
            heappop!(schedule.queue).action(env, schedule)
        end

        @fact env --> sort(all_pairs)
    end

    context("push! works on any callable") do
        schedule = Schedule()
        mc = MockCallable(10)
        @fact push!(schedule, mc, ActivationPoint(1.0, 0)) --> true

        @fact heappop!(schedule.queue).action("env", schedule) --> 10
    end

    context("will not schedule an item after ThirdWay.AFTER_SIMULATION") do
        schedule = Schedule()
        @fact push!(schedule, () -> true, ActivationPoint(1000.0, 0)) --> true

        ap = ActivationPoint(ThirdWay.AFTER_SIMULATION + 1, 0)
        @fact push!(schedule, () -> true, ap) --> false
    end

    context("will not accept an event prior to EPOCH") do
        schedule = Schedule()
        ap = ActivationPoint(ThirdWay.EPOCH - eps(ThirdWay.EPOCH), 0)
        @fact_throws ArgumentError push!(schedule, () -> true, ap)
    end

    context("will not accept an event with NaN as the time") do
        schedule = Schedule()
        @fact_throws ArgumentError push!(schedule, ActivationPoint(NaN, 0)) do
            # pass
        end
    end

    context("will not accept an event in the past") do
        schedule = Schedule()
        @fact push!(schedule, () -> true, ActivationPoint(5.0, 0)) --> true

        schedule.time = 10.0
        @fact_throws ArgumentError push!(schedule, ActivationPoint(5.0, 0)) do
            # pass
        end
    end

    context("can be clear!()ed") do
        schedule = Schedule()
        for i in 1:10
            push!(schedule, ActivationPoint(0.0, 0)) do env end
        end

        @fact length(schedule) --> 10
        @fact clear!(schedule) --> 10
        @fact length(schedule) --> 0
    end

    context("can be reset!()") do
        schedule = Schedule()
        schedule.time = 100.0
        schedule.steps = 100
        for i in 1:10
            push!(schedule, ActivationPoint(100.0, 0)) do env end
        end
        @fact length(schedule) --> 10

        reset!(schedule)

        @fact length(schedule) --> 0
        @fact schedule.time --> ThirdWay.BEFORE_SIMULATION
        @fact schedule.steps --> 0
    end

    context("can be copied") do
        schedule = Schedule()

        for name in ["Alice", "Bob", "Carol", "David"]
            push!(schedule, HelloAgent(name), ActivationPoint(0.0, 1))
        end

        cloned_sched = deepcopy(schedule)

        env_a = Vector{AbstractString}()
        run!(schedule, env_a)

        env_b = Vector{AbstractString}()
        run!(cloned_sched, env_b)

        @fact env_a --> env_b
    end

    context("can be reshuffle()ed") do
        schedule = Schedule()
        names = [randstring(5) for _ in 1:100]

        for name in names
            push!(schedule, HelloAgent(name), ActivationPoint(0.0, 1))
        end
        cloned_b = deepcopy(schedule)
        cloned_c = deepcopy(schedule)

        env_a = Vector{AbstractString}()
        run!(schedule, env_a)
        @fact is_complete(schedule) --> true

        env_b = Vector{AbstractString}()
        run!(cloned_b, env_b)
        @fact env_a == env_b --> true

        env_c = Vector{AbstractString}()
        reshuffle!(cloned_c)
        run!(cloned_c, env_c)

        @fact env_a != env_c --> true
    end

    context("schedule_once!") do
        schedule = Schedule()
        hit = false

        schedule_once!(schedule, 10.0, 3) do env, sch
            @fact sch.time --> 10.0
            @fact hit --> true
        end

        schedule_once!(schedule, 10.0, 1) do env, sch
            @fact sch.time --> 10.0
            @fact hit --> false
            hit = true
        end

        run!(schedule, "env")
    end

    context("schedule_once_in!") do
        schedule = Schedule()
        schedule.time = 100.0
        hit = false

        schedule_once_in!(schedule, 10.0, 3) do env, sch
            @fact sch.time --> 110.0
            @fact hit --> true
            @fact schedule.steps --> 0
        end

        schedule_once_in!(schedule, 10.0, 1) do env, sch
            @fact sch.time --> 110.0
            @fact hit --> false
            hit = true
            @fact schedule.steps --> 0
        end

        run!(schedule, "env")

        @fact schedule.steps --> 1
    end

    context("can be run! a specific number of steps") do
        schedule = Schedule()

        # 100 steps
        for delta in 1:100
            schedule_once_in!(schedule, Float64(delta)) do env, sch end
        end

        run!(schedule, Any[], Int64(7))

        @fact schedule.steps --> 7
    end

    context("can be run! with a callback") do
        schedule = Schedule()

        # 100 steps
        for delta in 1:100
            schedule_once_in!(schedule, Float64(delta)) do env, sch end
        end

        i = 1
        run!(schedule, Any[]) do env, schedule
            @fact schedule.steps --> i
            i += 1
            schedule.steps < 15
        end

        @fact schedule.steps --> 15
    end

    context("can be merge()d with another schedule") do
        schedule_a = Schedule()
        schedule_b = Schedule()
        times = shuffle(float(collect(1:10)))

        for t in times[1:5]
            schedule_once!(schedule_a, t, 0) do env, _
                push!(env, t)
            end
        end

        for t in times[6:10]
            schedule_once!(schedule_b, t, 0) do env, _
                push!(env, t)
            end
        end

        @fact length(schedule_a) --> 5
        @fact length(schedule_b) --> 5

        merge!(schedule_a, schedule_b)
        @fact length(schedule_a) --> 10

        env = Vector{Int64}()
        run!(schedule_a, env)
        @fact env --> sort(times)
    end

    context("cannot be merge()d with items in the past") do
        schedule_a = Schedule()
        schedule_a.time = 0.5
        schedule_b = Schedule()
        times = shuffle(float(collect(1:10)))

        for t in times[1:5]
            schedule_once!(schedule_a, t, 0) do env, _
                push!(env, t)
            end
        end

        schedule_once!(schedule_b, 0.0, 0) do env, _
            push!(env, t)
        end

        @fact_throws ErrorException merge!(schedule_a, schedule_b)
    end

    context("shows a nice human readable format") do
        schedule = Schedule()
        schedule.time = 100.0
        schedule.steps = 5
        schedule_once!(schedule, 200.0, 0) do env, _
            # pass
        end

        buf = IOBuffer()
        show(buf, schedule)
        expected = "Schedule at time=100.0 and steps=5 with 1 events."
        @fact UTF8String(buf.data) --> expected
    end

    context("schedule_repeating! repeats until stop()d") do
        schedule = Schedule()

        # Starting at t=0.0, repeat this action every 5.0 deltas.
        stoppable = schedule_repeating!(schedule, 0.0, 5.0, 0) do env, sched
            push!(env, schedule.time)
        end

        xs = Float64[]
        for step in 1:20
            step!(schedule, xs)
            @fact step --> schedule.steps

            # At this point, kill the action.
            if step == 10
                stop!(stoppable)
            end

            # One after the killed action, the schedule is exhausted.
            if is_complete(schedule)
                @fact step --> 11
                break
            end
        end

        @fact xs --> float(collect(0:5:45))
    end

    context("repeating action has a positive interval") do
        @fact_throws ArgumentError RepeatingAction(-1, 0) do env, sched
        end
    end
end

facts("SequenceOfActions") do
    context("executes each actionable in the specified order") do
        seq = SequenceOfActions([
            (env, sch) -> push!(env, "a"),
            (env, sch) -> push!(env, "b"),
            (env, sch) -> push!(env, "c")
        ])

        xs = AbstractString[]
        seq(xs, Schedule())
        @fact xs --> ["a", "b", "c"]
    end
end

facts("ShuffledActions") do
    context("shuffles then executes a set of actions") do
        seq = ShuffledActions([
            (env, sch) -> push!(env, "a"),
            (env, sch) -> push!(env, "b"),
            (env, sch) -> push!(env, "c")
        ])

        xs = AbstractString[]
        for _ in 1:10
            seq(xs, Schedule())
        end
        result = join(xs)
        repeated_str = repeat("abc", 10)
        @fact length(result) --> length(repeated_str)
        @fact repeated_str != result --> true
    end
end

facts("TentativeAction") do
    context("does not execute when stop!()d") do
        bad_f(env, schedule) = assert(false)
        @fact_throws AssertionError bad_f([], Schedule())

        tentative = TentativeAction(bad_f)
        stop!(tentative)
        tentative([], Schedule())
    end
end


text:
using RAFF
using Random
using Distributed
using Printf
using ArgParse

# Set Debug for Logging
using Logging
using Base.CoreLogging

global_logger(ConsoleLogger(stdout, Logging.Error))

s = ArgParseSettings()

@add_arg_table s begin

    "type"
    help = "Sequential(s) or parallel(p)"
    arg_type = String
    required = true

    "np"
    arg_type = Int
    help = "Number of points"
    required = true

    "nw"
    arg_type = Int
    help = "Number of workers"

end

parsed_args = parse_args(ARGS, s)

n = 2

np = parsed_args["np"]

p = Int(0.7 * np)

answer = [2.0, -0.5]

if parsed_args["type"] == "p"

    ###################
    # Distributed run #
    ###################

    nw = parsed_args["nw"]

    addprocs(nw)

    @everywhere using RAFF

    # Set Debug for Logging
    @everywhere using Logging
    @everywhere using Base.CoreLogging

    @everywhere global_logger(ConsoleLogger(stdout, Logging.Error))

    @everywhere gmodel!(x, t_, g) = begin
        g[1] = exp(t_ * x[2])
        g[2] = t_ * x[1] * exp(t_ * x[2]);
    end

    @everywhere model(x, t) = x[1] * exp(t * x[2])

    # First run
    praff(model, gmodel!, zeros(1, 3), n)

    # Set the seed for generating the same data
    Random.seed!(123456789)

    data, xSol = RAFF.generateNoisyData(model, n, np, p, answer)

    val, t, bytes, gctime, memallocs = try

        @timed praff(model, gmodel!, data, n)

    catch e

        ([1.0e+99, 1.0e+99], 1.0e+99, -1), -1.0, 0, -1.0, nothing

    end

    sol, fsol, psol = val

    @printf("%2d %10d %10.5s %15.8e %15.8e %15.8e %6d\n", nw, np,
            t, sol[1], sol[2], fsol, psol)

    rmprocs(workers())

else

    ##############
    # Serial run #
    ##############

    gmodel!(x, t_, g) = begin
        g[1] = exp(t_ * x[2])
        g[2] = t_ * x[1] * exp(t_ * x[2]);
    end

    model(x, t) = x[1] * exp(t * x[2])

    # First run
    raff(model, gmodel!, zeros(1, 3), n)

    # Set the seed for generating the same data
    Random.seed!(123456789)

    data, xSol = RAFF.generateNoisyData(model, n, np, p, answer)

    val, t, bytes, gctime, memallocs = try

        @timed raff(model, gmodel!, data, n)

    catch e

        ([1.0e+99, 1.0e+99], 1.0e+99, -1), -1.0, 0, -1.0, nothing

    end

    sol, fsol, psol = val

    @printf("%2d %10d %10.5s %15.8e %15.8e %15.8e %6d\n", 0, np,
            t, sol[1], sol[2], fsol, psol)

end

ext: jl


text:
"""
    showblockinterpretable(backend, encodings, block, data)

Decode `block` successively by applying `encodings` until a block is gotten
that can be shown by `backend`. Useful to visualize encoded data that
is not directly interpretable, for example an `Image{2}` representing an
encoded `Image`.

## Examples

```julia
using FastAI
encodings = (ImagePreprocessing(),)
block = ImageTensor{2}(3)
x = FastAI.mockblock(block)

showblockinterpretable(ShowText(), encodings, block, x)  # will decode to an `Image`
@test_throws showblock(ShowText(), encodings, block, x)  # will error
```

"""
function showblockinterpretable(backend::ShowBackend, encodings, block, data)
    res = decodewhile(
        block -> !isshowable(backend, block),
        encodings,
        Validation(),
        block,
        data)
    isnothing(res) && error("Could not decode to an interpretable block representation!")
    block_, data_ = res
    showblock(backend, block_, data_)
end

"""
    showblocksinterpretable(backend, encodings, block, datas)

Multi-sample version [`showblockinterpretable`](#).
"""
function showblocksinterpretable(backend::ShowBackend, encodings, block, datas::AbstractVector)
    blockdatas = [decodewhile(
        block -> !isshowable(backend, block),
        encodings,
        Validation(),
        block,
        data) for data in datas]
    any(isnothing, blockdatas) && error("Could not decode to an interpretable block representation!")
    block_ = first(first(blockdatas))
    datas_ = last.(blockdatas)
    showblocks(backend, block_, datas_)
end

# Helpers

function isshowable(backend::S, block::B) where {S<:ShowBackend, B<:AbstractBlock}
    hasmethod(FastAI.showblock!, (Any, S, B, Any))
end

"""
    decodewhile(f, encodings, ctx, block, data) -> (block', data')

Decode `block` by successively applying `encodings` to decode in
reverse order until `f(block') == false`.
"""
function decodewhile(f, encodings, ctx, block::AbstractBlock, data)
    encodings === () && return nothing
    if f(block)
        return decodewhile(
            f,
            encodings[1:end-1],
            ctx,
            decodedblock(encodings[end], block, true),
            decode(encodings[end], ctx, block, data),
        )
    else
        return (block, data)
    end
end

function decodewhile(f, encodings, ctx, blocks::Tuple, datas::Tuple)
    encodings === () && return nothing
    results = Tuple(decodewhile(f, encodings, ctx, block, data)
                for (block, data) in zip(blocks, datas))
    any(isnothing, results) && return nothing
    blocks = first.(results)
    datas = last.(results)
    return blocks, datas
end

function decodewhile(f, encodings, ctx, (title, block)::Pair, data)
    encodings === () && return nothing
    res = decodewhile(f, encodings, ctx, block, data)
    isnothing(res) && return nothing
    block_, data_ = res
    (title => block_, data_)
end


text:
using Test
using Distributed
procs = addprocs(2)
@everywhere using ProgressMeter

@testset "map tests" begin
    println("Testing map functions")

    # basic
    vals = progress_map(1:10) do x
        sleep(0.1)
        return x^2
    end
    @test vals == map(x->x^2, 1:10)

    vals = progress_map(1:10, mapfun=pmap) do x
        sleep(0.1)
        return x^2
    end
    @test vals == map(x->x^2, 1:10)

    val = progress_map(1:10, mapfun=reduce) do x, y
        sleep(0.1)
        return x+y
    end
    @test val == reduce((x,y)->x+y, 1:10)

    # errors
    @test_throws ErrorException progress_map(1:10) do x
        if x > 3
            error("intentional error")
        end
        return x^2
    end
    println()

    @test_throws RemoteException progress_map(1:10, mapfun=pmap) do x
        if x > 3
            error("intentional error")
        end
        return x^2
    end
    println()

    @test_throws ErrorException progress_map(1:10, mapfun=reduce) do x, y
        if x > 3
            error("intentional error")
        end
        return x + y
    end
    println()

    # @showprogress
    vals = @showprogress map(1:10) do x
        return x^2
    end
    @test vals == map(x->x^2, 1:10)

    vals = @showprogress pmap(1:10) do x
        return x^2
    end
    @test vals == map(x->x^2, 1:10)

    val = @showprogress reduce(1:10) do x, y
        return x + y
    end
    @test val == reduce((x, y)->x+y, 1:10)

    # function passed by name
    function testfun(x)
        return x^2
    end
    vals = @showprogress map(testfun, 1:10)
    @test vals == map(testfun, 1:10)
    vals = @showprogress pmap(testfun, 1:10)
    @test vals == map(testfun, 1:10)
    val = @showprogress reduce(+, 1:10)
    @test val == reduce(+, 1:10)

    # #136: make sure mid progress shows up even without sleep
    println("Verify that intermediate progress is displayed:")
    @showprogress map(1:100) do i
        A = rand(10000,1000)
        sum(A)
    end

    # multiple args
    vals = @showprogress pmap((x,y)->x*y, 1:10, 2:11)
    @test vals == map((x,y)->x*y, 1:10, 2:11)

    # abstract worker pool arg
    wp = WorkerPool(procs)
    vals = @showprogress pmap(testfun, wp, 1:10)
    @test vals == map(testfun, 1:10)

    vals = @showprogress pmap(wp, 1:10) do x
        x^2
    end
    @test vals == map(testfun, 1:10)

    # Progress args
    vals = @showprogress 0.1 "Computing" pmap(testfun, 1:10)
    @test vals == map(testfun, 1:10)

    # named vector arg
    a = collect(1:10)
    vals = @showprogress pmap(x->x^2, a)
    @test vals == map(x->x^2, a)

    # global variable in do
    C = 10
    vals = @showprogress pmap(1:10) do x
        return C*x
    end
    @test vals == map(x->C*x, 1:10)

    # keyword arguments
    vals = @showprogress pmap(x->x^2, 1:100, batch_size=10)
    @test vals == map(x->x^2, 1:100)
    # with semicolon
    vals = @showprogress pmap(x->x^2, 1:100; batch_size=10)
    @test vals == map(x->x^2, 1:100)
end

rmprocs(procs)


text:
using DSP
include("chroma.jl")
include("dtw.jl")

"""
sndalign(
	signal1::Vector{Cdouble}, signal2::Vector{Cdouble},
	fs::Cdouble, wsize::Csize_t, overlap::Csize_t, γ::Cdouble, ϵ::Cdouble
)

Compute the set of indices `path1, path2` that align the two given signals.

The alignment is computed with the dtw algorithm over a set of chroma features
extracted from the signals' stft, therefore the returned vectors give the
alignment of the windows of the transform and not of the signals' samples
themselves.

**Parameters:**

`signal1`, `signal2`: The two signals to be aligned.

`fs`: The sampling rate of the signals.

`wsize`: The number of samples of the window to be used in the STFT.

`overlap`: The number of samples shared by two windows in the STFT.

`γ`: The compression factor to be used in the chroma features extraction.

`ϵ`: The smallest norm that a chroma vector needs to have to not be replaced by
a stub one.
"""
function sndalign(
	signal1::Vector{Cdouble},
	signal2::Vector{Cdouble},
	fs::Cdouble,
	wsize::Csize_t,
	overlap::Csize_t,
	γ::Cdouble,
	ϵ::Cdouble
)

   stft1 = stft(signal1, Int64(wsize), Int64(overlap);fs=fs,window=ones)
	chroma1 = chroma(stft1, fs, wsize, γ, ϵ)

	stft2 = stft(signal2, Int64(wsize), Int64(overlap);fs=fs,window=ones)
	chroma2 = chroma(stft2, fs, wsize, γ, ϵ)

	return dtw(chroma1, chroma2)

end

"""
c_sndalign(
	signal1::Vector{Cdouble}, signal2::Vector{Cdouble},
	fs::Cdouble, wsize::Csize_t, overlap::Csize_t, γ::Cdouble, ϵ::Cdouble
)

Compute the set of indices `path1, path2` that align the two given signals.

The alignment is computed with the dtw algorithm over a set of chroma features
extracted from the signals' stft, therefore the returned vectors give the
alignment of the windows of the transform and not of the signals' samples
themselves.

This is a wrapper that calls the C++ implementation of the alignment algorithm.

**Parameters:**

`signal1`, `signal2`: The two signals to be aligned.

`fs`: The sampling rate of the signals.

`wsize`: The number of samples of the window to be used in the STFT.

`overlap`: The number of samples shared by two windows in the STFT.

`γ`: The compression factor to be used in the chroma features extraction.

`ϵ`: The smallest norm that a chroma vector needs to have to not be replaced by
a stub one.
"""
function c_sndalign(
	signal1::Vector{Cdouble},
	signal2::Vector{Cdouble},
	fs::Cdouble,
	wsize::Csize_t,
	overlap::Csize_t,
	γ::Cdouble,
	ϵ::Cdouble
)

	alignmentLength = Array{Csize_t}(undef, 1)

	c_alignment = ccall(
		(:soundAlign, "./libsndalign.so"), Ptr{Csize_t},
		(Ptr{Cdouble}, Csize_t, Ptr{Cdouble}, Csize_t,
		Cdouble, Csize_t, Csize_t, Cdouble, Cdouble, Ptr{Csize_t}),
		signal1, length(signal1), signal2, length(signal2),
		fs, wsize, overlap, γ, ϵ, alignmentLength
	)

	alignmentLength = alignmentLength[] ÷ 2
	offset = alignmentLength * sizeof(Csize_t)
	alignment1 = Vector{Csize_t}(undef, alignmentLength)
	alignment2 = Vector{Csize_t}(undef, alignmentLength)
	unsafe_copyto!(pointer(alignment1), c_alignment, alignmentLength)
	unsafe_copyto!(pointer(alignment2), c_alignment + offset, alignmentLength)

	ccall(
		(:freeSndAligment, "./libsndalign.so"), Cvoid,
		(Ptr{Cvoid},),
		c_alignment
	)

	return (Int.(alignment1), Int.(alignment2))

end

max_stars_repo_name: gabriel-ss/interpretation-switcher

max_issues_repo_name: gabriel-ss/interpretation-switcher

max_forks_repo_name: gabriel-ss/interpretation-switcher


text:
using LoopVectorization
using LinearAlgebra
using Octavian

function RMP2_canonical_energy(ints::IntegralHelper{T,RIFIT,O}, Alg::RMP2b) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}
    Bvo = permutedims(ints["BOV"], (1,3,2))
    ϵo = ints["Fii"]
    ϵv = ints["Faa"]

    output(" Computing DF-MP2 Energy")
    output(" - Contraction engine: Octavian")
    v_size = length(ϵv)
    o_size = length(ϵo)

    # Pre-allocating arrays for threads
    BABs = [zeros(T, v_size, v_size) for _ = 1:Threads.nthreads()]

    # Vector containing the energy contribution computed by each thread
    ΔMP2s = zeros(T,Threads.nthreads())
    TWO = T(2.0)
    ONE = one(T)
    t = @elapsed begin
        @sync for i in 1:o_size
            Threads.@spawn begin
            id = Threads.threadid()

            Bab = BABs[id]
            @views Bi = Bvo[:,:,i]

            for j in i:o_size

                @views Bj = Bvo[:,:,j]
                matmul_serial!(Bab, transpose(Bi), Bj)

                eij = ϵo[i] + ϵo[j]
                E = zero(T)
                @turbo for a = eachindex(ϵv)
                    eija = eij - ϵv[a]
                    for b = eachindex(ϵv)
                        D = eija - ϵv[b]
                        E += Bab[a,b] * (TWO * Bab[a,b] - Bab[b,a]) / D
                    end
                end
                fac = i !== j ? TWO : ONE
                ΔMP2s[id] += fac * E
            end
        end # spawn
        end # sync
        Emp2 = sum(ΔMP2s)
    end # time
    output("Done in {:5.5f} seconds.", t)

    return Emp2
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
using Documenter, SEMFitting

makedocs()


text:
module DMM
using Libdl # for NI dlls
include("../INSTR_TYPES.jl")
include("../lib_locations.jl")
include("../visa/VISA.jl")
#
## include all DMM files
include("./Keithley2000/KE2000.jl")
# NI
include("NIDmm/NIDmm.jl")
include("NIDmm/NIDmm_Misc.jl")
end #endmodule


text:
@testset "2-area, 11-bus, 4-generator benchmark system Parsing " begin
    test_dir = mktempdir()
    ger4_raw_file = joinpath(PSSE_RAW_DIR, "Benchmark_4ger_33_2015.RAW")
    ger4_dyr_file = joinpath(PSSE_DYR_DIR, "Benchmark_4ger_33_2015.dyr")
    sys = PSB.build_system(PSB.PSSETestSystems, "psse_Benchmark_4ger_33_2015_sys")
    @test_throws PSY.DataFormatError System(ger4_raw_file, ger4_raw_file)
    @test_throws PSY.DataFormatError System(ger4_dyr_file, ger4_raw_file)
    dyn_injectors = get_components(DynamicInjection, sys)
    @test length(dyn_injectors) == 4
    for g in dyn_injectors
        m = PSY.get_machine(g)
        @test typeof(m) <: RoundRotorExponential
    end
    path = joinpath(test_dir, "test_dyn_system_serialization.json")
    to_json(sys, path)
    parsed_sys = System(path)
    dyn_injectors = get_components(DynamicInjection, parsed_sys)
    @test length(dyn_injectors) == 4
    for g in dyn_injectors
        m = PSY.get_machine(g)
        @test typeof(m) <: RoundRotorExponential
    end
end

@testset "2-area, 11-bus, 4-generator with renewables benchmark system Parsing " begin
    test_dir = mktempdir()
    sys = PSB.build_system(PSB.PSSETestSystems, "psse_renewable_parsing_1")
    dyn_injectors = get_components(DynamicInjection, sys)
    @test length(dyn_injectors) == 5
    for g in dyn_injectors
        if isa(g, DynamicGenerator)
            m = PSY.get_machine(g)
            @test typeof(m) <: RoundRotorExponential
        elseif isa(g, DynamicInverter)
            ic = PSY.get_inner_control(g)
            @test typeof(ic) <: RECurrentControlB
        else
            @error("Generator $g not supported")
        end
    end
    path = joinpath(test_dir, "test_dyn_system_serialization.json")
    to_json(sys, path)
    parsed_sys = System(path)
    dyn_injectors = get_components(DynamicInjection, parsed_sys)
    @test length(dyn_injectors) == 5
    for g in dyn_injectors
        if isa(g, DynamicGenerator)
            m = PSY.get_machine(g)
            @test typeof(m) <: RoundRotorExponential
        elseif isa(g, DynamicInverter)
            ic = PSY.get_inner_control(g)
            @test typeof(ic) <: RECurrentControlB
        else
            @error("Generator $g not supported")
        end
    end
end

@testset "240 Bus WECC system Parsing " begin
    test_dir = mktempdir()
    sys = PSB.build_system(PSB.PSSETestSystems, "psse_240_case_renewable_sys")
    dyn_injectors = get_components(DynamicInjection, sys)
    @test length(dyn_injectors) == 146
    for g in dyn_injectors
        if isa(g, DynamicGenerator)
            m = PSY.get_machine(g)
            @test typeof(m) <: RoundRotorQuadratic
        elseif isa(g, DynamicInverter)
            ic = PSY.get_inner_control(g)
            @test typeof(ic) <: RECurrentControlB
        else
            @error("Generator $g not supported")
        end
    end
    path = joinpath(test_dir, "test_dyn_system_serialization.json")
    to_json(sys, path)
    parsed_sys = System(path)
    dyn_injectors = get_components(DynamicInjection, parsed_sys)
    @test length(dyn_injectors) == 146
    for g in dyn_injectors
        if isa(g, DynamicGenerator)
            m = PSY.get_machine(g)
            @test typeof(m) <: RoundRotorQuadratic
        elseif isa(g, DynamicInverter)
            ic = PSY.get_inner_control(g)
            @test typeof(ic) <: RECurrentControlB
        else
            @error("Generator $g not supported")
        end
    end
end

@testset "GENCLS dyr parsing" begin
    threebus_raw_file = joinpath(PSSE_TEST_DIR, "ThreeBusNetwork.raw")
    gencls_dyr_file = joinpath(PSSE_TEST_DIR, "TestGENCLS.dyr")
    nogencls_dyr_file = joinpath(PSSE_TEST_DIR, "Test-NoCLS.dyr")
    sexs_dyr_file = joinpath(PSSE_TEST_DIR, "Test_SEXS.dyr")
    sys_gencls = PSB.build_system(PSB.PSSETestSystems, "psse_3bus_gen_cls_sys")
    sys_nogencls = PSB.build_system(PSB.PSSETestSystems, "psse_3bus_no_cls_sys")
    sys_sexs = PSB.build_system(PSB.PSSETestSystems, "psse_3bus_SEXS_sys")

    #Check that generator at bus 102 (H = 0) is a Source, and not ThermalStandard.
    @test isnothing(get_component(ThermalStandard, sys_gencls, "generator-102-1"))
    @test typeof(get_component(StaticInjection, sys_gencls, "generator-102-1")) <: Source
    #Check that generator at bus 102 (not provided in dyr data) is a Source.
    @test isnothing(get_component(ThermalStandard, sys_nogencls, "generator-102-1"))
    @test typeof(get_component(StaticInjection, sys_nogencls, "generator-102-1")) <: Source

    #Check Parameters of GENROE raw file
    genroe_gen = get_component(DynamicInjection, sys_gencls, "generator-101-1")
    m = PSY.get_machine(genroe_gen)
    s = PSY.get_shaft(genroe_gen)
    @test typeof(m) <: RoundRotorExponential

    @test get_R(m) == 0.0
    @test get_Td0_p(m) == 8.0
    @test get_Td0_pp(m) == 0.03
    @test get_Tq0_p(m) == 0.4
    @test get_Tq0_pp(m) == 0.05
    @test get_Xd(m) == 1.8
    @test get_Xq(m) == 1.7
    @test get_Xd_p(m) == 0.3
    @test get_Xd_pp(m) == 0.25
    @test get_Xl(m) == 0.2
    @test get_H(s) == 6.5
    @test get_D(s) == 0.0

    #Check Parameters of other GENCLS file
    gencls_gen = get_component(DynamicInjection, sys_gencls, "generator-103-1")
    m = PSY.get_machine(gencls_gen)
    s = PSY.get_shaft(gencls_gen)
    @test typeof(m) <: BaseMachine

    @test get_R(m) == 0.0
    @test get_Xd_p(m) == 0.2
    @test get_H(s) == 3.1
    @test get_D(s) == 2.0

    #Check Parameters of other GENCLS file
    gencls_gen = get_component(DynamicInjection, sys_nogencls, "generator-103-1")
    m = PSY.get_machine(gencls_gen)
    s = PSY.get_shaft(gencls_gen)
    @test typeof(m) <: BaseMachine

    @test get_R(m) == 0.0
    @test get_Xd_p(m) == 0.2
    @test get_H(s) == 3.1
    @test get_D(s) == 2.0

    #Check Parameters of SEXS file
    genroe_sexs = get_component(DynamicInjection, sys_sexs, "generator-101-1")
    a = PSY.get_avr(genroe_sexs)
    @test typeof(a) <: AVR
    @test get_Ta_Tb(a) == 0.4
    @test get_Tb(a) == 5.0
    @test get_K(a) == 20.0
    @test get_Te(a) == 1.0
    E_min, E_max = get_V_lim(a)
    @test E_min == -50.0
    @test E_max == 50.0
end

@testset "2000-Bus Parsing" begin
    test_dir = mktempdir()
    texas2000_raw_file = joinpath(TAMU_DIR, "ACTIVSg2000.RAW")
    texas2000_dyr_file = joinpath(PSSE_DYR_DIR, "ACTIVSg2000_dynamics.dyr")
    sys = PSB.build_system(PSB.PSSETestSystems, "psse_ACTIVSg2000_sys")
    for g in get_components(ThermalStandard, sys)
        @test !isnothing(get_dynamic_injector(g))
    end
    path = joinpath(test_dir, "test_dyn_system_serialization_2000.json")
    to_json(sys, path)
    parsed_sys = System(path)
    dyn_injectors = get_components(DynamicInjection, parsed_sys)
    @test length(dyn_injectors) == 435
    for g in get_components(ThermalStandard, parsed_sys)
        @test !isnothing(get_dynamic_injector(g))
    end
end


text:
function champs_2(t, fft_ubar, fft_vbar, fft_ug, fft_vg,
                  A1, A2, matr, conjmatr, sigma, llambda,
                  ktau, epsilon, Ntaumm)

    champu, champv = ftau(t, fft_ubar, fft_vbar,
                     A1, A2, matr, conjmatr, sigma, llambda)

    champu_fft = fft(champu, 1)
    champv_fft = fft(champv, 1)

    champu_fft[1, :] .= 0.0
    champv_fft[1, :] .= 0.0

    dtauh1u = ifft(champu_fft, 1)
    dtauh1v = ifft(champv_fft, 1)

    h1u = epsilon * ifft(champu_fft ./ (1im * ktau), 1)
    h1v = epsilon * ifft(champv_fft ./ (1im * ktau), 1)

    C1u = fft_ubar .+ h1u
    C1v = fft_vbar .+ h1v

    ffu, ffv = ftau(t,
                    C1u .+ fft_ug,
                    C1v .+ fft_vg,
                    A1, A2, matr, conjmatr,
                    sigma, llambda)

    champu, champv = ftau(t, C1u, C1v, A1, A2, matr, conjmatr, sigma, llambda)

    champu_fft = fft(champu, 1)
    champv_fft = fft(champv, 1)

    champubaru = champu_fft[1, :] / Ntaumm
    champubarv = champv_fft[1, :] / Ntaumm

    champu1, champv1 = duftau(t, fft_ubar, fft_vbar,
                              champubaru, champubarv,
                              A1, A2, matr, conjmatr, sigma, llambda,
                              Ntaumm)

    champu2, champv2 = dtftau(t, fft_ubar, fft_vbar,
                              A1, A2, matr, conjmatr,
                              sigma, llambda, Ntaumm)

    champu_fft = fft(champu1 .+ champu2, 1)
    champv_fft = fft(champv1 .+ champv2, 1)

    champu_fft[1, :] .= 0.0
    champv_fft[1, :] .= 0.0

    dtC1u = champubaru .+ epsilon * ifft(champu_fft / (1im * ktau), 1)
    dtC1v = champubarv .+ epsilon * ifft(champv_fft / (1im * ktau), 1)

    champgu = ffu .- dtauh1u .- dtC1u
    champgv = ffv .- dtauh1v .- dtC1v

    champgu_fft = fft(champgu, 1)
    champgv_fft = fft(champgv, 1)

    champmoyu = champgu_fft[1, :] / Ntaumm
    champmoyv = champgv_fft[1, :] / Ntaumm

    champgu_fft[1, :] .= 0.0
    champgv_fft[1, :] .= 0.0

    ichampgu = ifft(champgu_fft ./ (1im * ktau), 1)
    ichampgv = ifft(champgv_fft ./ (1im * ktau), 1)

    champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv

end


text:
@testset "Entropy" begin
    @test_throws MethodError Entropy()
    @test_throws MethodError Entropy(0.0)
    @test_throws MethodError Entropy(0, 1.0)
    for b in -2:0
        @test_throws ArgumentError Entropy(b, 2)
        @test_throws ArgumentError Entropy(2, b)
    end

    let dist = Entropy(2)
        @test dist.data == zeros(Int, 2)
        @test dist.bs == (2,)
        @test dist.N == 0
    end

    let dist = Entropy(2, 3)
        @test dist.data == zeros(Int, 6)
        @test dist.bs == (2, 3)
        @test dist.N == 0
    end

    let dist = Entropy(3)
        observe!(dist, [1,2,3,2,3])
        @test dist.data == [1,2,2]
        @test dist.bs == (3,)
        @test dist.N == 5
    end

    let dist = Entropy(3)
        observe!(dist, [1 2 3 2 3])
        @test dist.data == [1,2,2]
        @test dist.bs == (3,)
        @test dist.N == 5
    end

    let dist = Entropy(2, 2)
        observe!(dist, [2 1 2 1; 1 1 2 2])
        @test dist.data == [1,1,1,1]
        @test dist.bs == (2, 2)
        @test dist.N == 4
    end

    let dist = Entropy(2, 3)
        data = zeros(Int, 2, 4, 3)
        data[:,:,1] = [1 1 1 1; 2 2 2 2]
        data[:,:,2] = [1 2 2 2; 1 3 3 2]
        data[:,:,3] = [1 1 1 1; 2 2 2 2]

        observe!(dist, data)

        @test dist.data == [1,0,8,1,0,2]
        @test dist.bs == (2, 3)
        @test dist.N == 12
    end

    let dist = Entropy(2, 3)
        observe!(dist, [1 1 1 1; 2 2 2 2])
        observe!(dist, [1 2 2 2; 1 3 3 2])
        @test dist.data == [1,0,4,1,0,2]
        @test dist.bs == (2, 3)
        @test dist.N == 8
    end

    @test_throws ArgumentError Entropy(Int[])
    @test_throws ArgumentError Entropy(zeros(Int, 0, 0))
    @test_throws ArgumentError Entropy([1,2,0])
    @test_throws ArgumentError Entropy([0 1 1; 1 1 2])
    @test_throws ArgumentError Entropy([2 1 1; 0 1 2])

    let dist = Entropy([1 1 1 1 1 2 2 2; 2 2 2 2 1 1 2 2])
        @test dist.data == [1,1,4,2]
        @test dist.bs == (2, 2)
        @test dist.N == 8
    end

    let dist = Entropy([1 1 1 1 1 1 1 1; 2 2 2 2 1 1 2 2])
        @test dist.data == [2,6]
        @test dist.bs == (1,2)
        @test dist.N == 8
    end

    let data = zeros(Int, 2, 5, 3)
        data[:,:,1] = [1 1 1 1 1; 2 2 2 2 2]
        data[:,:,2] = [1 2 2 1 1; 3 3 2 2 1]
        data[:,:,3] = [1 1 1 3 3; 1 1 1 2 2]

        dist = Entropy(data)

        @test dist.data == [4,0,0,6,1,2,1,1,0]
        @test dist.bs == (3, 3)
        @test dist.N == 15
    end

    @test_throws BoundsError observe!(Entropy(2), [1 2 3])
    @test_throws BoundsError observe!(Entropy(2), [1,2,3])
    @test_throws BoundsError observe!(Entropy(2, 2), [1 2 3; 1 2 3])

    let dist = observe!(Entropy(2), [1,2,2,1,1])
        clear!(dist)
        @test dist.data == zeros(Int, 2)
        @test dist.bs == (2,)
        @test dist.N == 0
    end

    let dist = observe!(Entropy(2, 2), [1 1 2 2; 1 2 1 2])
        clear!(dist)
        @test dist.data == zeros(Int, 4)
        @test dist.bs == (2,2)
        @test dist.N == 0
    end

    let data = zeros(Int, 2, 5, 3)
        data[:,:,1] = [1 1 1 1 1; 2 2 2 2 2]
        data[:,:,2] = [1 2 2 1 1; 3 3 2 2 1]
        data[:,:,3] = [1 1 1 3 3; 1 1 1 2 2]

        dist = Entropy(data)
        clear!(dist)

        @test dist.data == zeros(Int, 9)
        @test dist.bs == (3, 3)
        @test dist.N == 0
    end

    let xs = [1 1 1 1 2 2 2 2;
              2 2 2 2 1 1 1 1]
        @test entropy(xs) ≈ 1.0 atol=1e-6
    end

    let xs = [1 1 2 2 2 2 1 1 1;
              2 2 1 1 1 1 2 2 2]
        @test entropy(xs) ≈ 0.991076 atol=1e-6
    end

    let xs = [2 2 1 2 1 2 2 2 1;
              2 2 1 1 1 2 1 2 2]
        @test entropy(xs) ≈ 1.836592 atol=1e-6
    end

    let xs = [1 1 1 1 1 1 1 1 1;
              2 2 2 1 1 1 2 2 2]
        @test entropy(xs) ≈ 0.918296 atol=1e-6
    end

    let xs = [2 2 2 2 1 1 1 1 2;
              2 2 2 1 1 1 2 2 2]
        @test entropy(xs) ≈ 1.836592 atol=1e-6
    end

    let xs = [2 2 1 1 2 2 1 1 2;
              2 2 2 1 1 1 2 2 2]
        @test entropy(xs) ≈ 1.891061 atol=1e-6
    end

    let xs = [1 2 1 2 1 2 2 2 2]
        @test entropy(xs) ≈ 0.918296 atol=1e-6
    end

    let xs = [1,2,1,2,1,2,2,2,2]
        @test entropy(xs) ≈ 0.918296 atol=1e-6
    end

    let xs = [1,2,3,1,2,3,1,1,3]
        @test entropy(xs) ≈ 1.530493 atol=1e-6
    end

    let xs = zeros(Int, 1, 5, 2)
        xs[:,:,1] = [1,2,1,2,1]
        xs[:,:,2] = [1,2,3,2,1]

        @test entropy(xs) ≈ 1.360964 atol=1e-6
    end

    let xs = zeros(Int, 2, 9, 2)
        xs[:,:,1] = [2 2 1 1 2 2 1 1 2;
                     2 2 2 1 1 1 2 2 2]
        xs[:,:,2] = [1 1 1 1 1 1 1 1 1;
                     2 2 2 1 1 1 2 2 2]

        @test entropy(xs) ≈ 1.765246 atol=1e-6
    end
end


text:
using Test
using GmshSDK
import WavePropBase as WPB

@testset "Basic tests" begin
    @gmsh begin
        WPB.clear_entities!()
        # try some basic things just to make sure it does not error
        gmsh.model.geo.addPoint(0,0,0)    # test native CAO
        gmsh.model.occ.addSphere(0,0,0,1) # test occ is available
        gmsh.model.occ.synchronize()
        gmsh.model.mesh.generate(3)
        gmsh.model.mesh.partition(5) # test metis is available
        @test true == true #
    end
    # read domain from a model
    @gmsh begin
        WPB.clear_entities!()
        gmsh.model.occ.addSphere(0,0,0,1)
        gmsh.model.occ.synchronize()
        Ω = GmshSDK.domain()
        @test WPB.geometric_dimension(Ω) == 3
    end
    # read domain and mesh from a model
    Ω,M = @gmsh begin
        WPB.clear_entities!()
        gmsh.model.occ.addDisk(0,0,0,1,1)
        gmsh.model.occ.synchronize()
        gmsh.model.mesh.generate(2)
        Ω = GmshSDK.domain(dim=2)
        M = GmshSDK.meshgen(Ω,dim=2)
        return Ω,M
    end
    Ω,M = @gmsh begin
        WPB.clear_entities!()
        gmsh.model.occ.addRectangle(0,0,0,1,1)
        gmsh.model.occ.synchronize()
        gmsh.model.mesh.generate(2)
        Ω = GmshSDK.domain(dim=2)
        M = GmshSDK.meshgen(Ω,dim=2)
        return Ω,M
    end
    # read domain from a file
    @gmsh begin
        fname = joinpath(GmshSDK.PROJECT_ROOT,"test","circle.geo")
        Ω =  GmshSDK.read_geo(fname;dim=2)
    end
    # read domain and mesh from a file
    Ω, msh = @gmsh begin
        WPB.clear_entities!()
        dir = @__DIR__
        fname = joinpath(dir,"circle.msh")
        GmshSDK.read_msh(fname;dim=2)
    end
end

max_stars_repo_name: maltezfaria/GmshSDK

max_issues_repo_name: maltezfaria/GmshSDK

max_forks_repo_name: maltezfaria/GmshSDK


text:
# This file is a part of TextSearch.jl

export neardup

function neardup(X::AbstractVector{T}, α=0.1) where T
    idx = InvertedFile()
    res = KnnResult(1)
    L = zeros(Int, length(X))
    D = zeros(Float64, length(X))
    L[1] = 1
    D[1] = 0.0
    push!(idx, 1 => X[1])
    @inbounds for i in 2:length(X)
        empty!(res)
        x = X[i]
        search(idx, x, res)
        if length(res) == 0 || minimum(res) > α
            push!(idx, i => x)
            L[i] = i
            D[i] = 0.0
        else
            L[i], D[i] = first(res)
        end
    end

    (idx=idx, nn=L, dist=D)
end


text:
using Base: @deprecate

@deprecate intersection(pt1::Point, pt2::Point, pt3::Point, pt4::Point) intersectionlines(pt1::Point, pt2::Point, pt3::Point, pt4::Point)

@deprecate seednoise(a::Array{Int64}) seedpatentednoise()


text:
using SafeTestsets

@time begin
    @time @safetestset "Quantum Harmonic oscilator tests" begin include("harmonic_oscilator.jl") end
    @time @safetestset "Schwarzschild quasinormal modes tests" begin include("schwarzschild_qnm.jl") end
end


text:
using Pkg
Pkg.add("DataStructures")
Pkg.add("Revise")

max_stars_repo_name: juhamyllari/minipascal-julia

max_issues_repo_name: juhamyllari/minipascal-julia

max_forks_repo_name: juhamyllari/minipascal-julia


text:
# functions for building various standard types of graphs

export Complete, Path, Cycle, RandomGraph, RandomRegular, RandomSBM
export RandomTree, code_to_tree
export Grid, Wheel, Cube, BuckyBall, Johnson
export Petersen, Kneser, Paley, Knight, Frucht, Hoffman, HoffmanSingleton

"""
`Complete(n)` returns a complete graph with `n` vertices `1:n`.

`Complete(n,m)` returns a complete bipartite graph with `n` vertices
in one part and `m` vertices in the other.

`Complete([n1,n2,...,nt])` returns a complete multipartite graph with
parts of size `n1`, `n2`, ..., `nt`.
"""
function Complete(n::Int)
    G = IntGraph(n)

    for k=1:n-1
        for j=k+1:n
            add!(G,j,k)
        end
    end
    name(G,"Complete graph K($n)")
    return G
end

# Create a complete bipartite graph
function Complete(n::Int, m::Int)
    G = IntGraph(n+m)
    for u=1:n
        for v=n+1:n+m
            add!(G,u,v)
        end
    end
    name(G,"Complete bipartite graph K($n,$m)")
    return G
end

# Create the complete multipartite graph with given part sizes
function Complete(parts::Array{Int,1})
    # check all part sizes are positive
    for p in parts
        if p < 1
            error("All part sizes must be positive")
        end
    end

    n = sum(parts)
    G = IntGraph(n)

    np = length(parts)
    if np < 2
        return G
    end

    # create table of part ranges
    ranges = Array{Int}(undef,np,2)   # old ranges = Array(Int,np,2)
    ranges[1,1] = 1
    ranges[1,2] = parts[1]
    for k = 2:np
        ranges[k,1] = ranges[k-1,2] + 1
        ranges[k,2] = ranges[k,1] + parts[k] - 1
    end

    # Add all edges between all parts
    for i=1:np-1
        for j=i+1:np
            for u=ranges[i,1]:ranges[i,2]
                for v=ranges[j,1]:ranges[j,2]
                    add!(G,u,v)
                end
            end
        end
    end
    name(G,"Complete multipartite graph K$parts")
    return G
end

# Create a path graph on n vertices
"""
`Path(n)` creates a path graph with `n` vertices named `1:n`.

`Path(array)` creates a path graph with vertices `array[1]`,
`array[2]`, etc.
"""
function Path(n::Int)
    G = IntGraph(n)
    for v = 1:n-1
        add!(G,v,v+1)
    end
    name(G,"Path graph P($n)")
    return G
end

# Create a path graph from a list of vertices
function Path(verts::Array{T}) where {T}
    G = SimpleGraph{T}()
    n = length(verts)

    if n==1
        add!(G,verts[1])
    end
    for k = 1:n-1
        add!(G,verts[k],verts[k+1])
    end
    name(G,"Path graph with $n vertices")
    return G
end

# Create a cycle graph on n vertices
"""
`Cycle(n)` creates a cycle with vertex set `1:n`.
"""
function Cycle(n::Int)
    if n<3
        error("Cycle requires 3 or more vertices")
    end
    G = Path(n)
    add!(G,1,n)
    name(G,"Cycle graph C($n)")
    return G
end

# Create the wheel graph on n vertices: a cycle on n-1 vertices plus
# an additional vertex adjacent to all the vertices on the wheel.
"""
`Wheel(n)` creates a wheel graph with `n` vertices. That is, a cycle
with `n-1` vertices `1:(n-1)` all adjacent to a common single vertex,
`n`.
"""
function Wheel(n::Int)
    if n < 4
        error("Wheel graphs must have at least 4 vertices")
    end
    G = Cycle(n-1)
    for k=1:n-1
        add!(G,k,n)
    end
    name(G,"Wheel graph with $n vertices")
    return G
end

# Create a grid graph
"""
`Grid(n,m)` creates an `n`-by-`m` grid graph. For other grids, we
suggest `Path(n1)*Path(n2)*Path(n3)` optionally wrapped in
`relabel`. See also: `Cube`.
"""
function Grid(n::Int, m::Int)
    G = SimpleGraph{Tuple{Int,Int}}()

    # add the vertices
    for u=1:n
        for v=1:m
            add!(G,(u,v))
        end
    end

    #horizontal edges
    for u=1:n
        for v=1:m-1
            add!(G,(u,v),(u,v+1))
        end
    end

    # vertical edges
    for v=1:m
        for u=1:n-1
            add!(G,(u,v),(u+1,v))
        end
    end

    # embedding
    d = Dict{Any, Array{Float64,1}}()
    for v in G.V
        x,y = v
        d[v] = [Float64(x),Float64(y)]
    end
    embed(G,d)
    recenter(G)

    name(G,"Grid graph of size $n-by-$m")
    return G
end

# Create an Erdos-Renyi random graph
"""
`RandomGraph(n,p=0.5)` creates an Erdos-Renyi random graph with `n`
vertices and edge probability `p`.
"""
function RandomGraph(n::Int, p::Real=0.5)
    G = IntGraph(n)

    # guess the size of the edge set to preallocate storage
    m = round(Int,n*n*p)+1

    # generate the edges
    for v=1:n-1
        for w=v+1:n
            if (rand() < p)
                add!(G,v,w)
            end
        end
    end
    name(G,"Erdos-Renyi random graph, p=$p")
    return G
end

# Generate a random tree on vertex set 1:n. All n^(n-2) trees are
# equally likely.

"""
`RandomTree(n)` creates a random tree on `n` vertices each with
probability `1/n^(n-2)`.
"""
function RandomTree(n::Int)
    if n<0   # but we allow n==0 to give empty graph
        error("Number of vertices cannot be negative")
    end

    if n<2
        return IntGraph(n)
    end

    code = [ mod(rand(Int),n)+1 for _ in 1:n-2 ]
    G =  code_to_tree(code)
    name(G,"Random tree")
    return G
end

# This is a helper function for RandomTree that converts a Prufer code
# to a tree. No checks are done on the array fed into this function.
function code_to_tree(code::Array{Int,1})
    n = length(code)+2
    G = IntGraph(n)
    degree = ones(Int,n)  # initially all 1s

    #every time a vertex appears in code[], up its degree by 1
    for c in code
        degree[c]+=1
    end

    for u in code
        for v in 1:n
            if degree[v]==1
                add!(G,u,v)
                degree[u] -= 1
                degree[v] -= 1
                break
            end
        end
    end

    last = findall(x->x!=0,degree)
    add!(G,last[1],last[2])

    return G
end

# Create the Cube graph with 2^n vertices
"""
`Cube(n)` creates the `n`-dimensional cube graph. This graph has `2^n`
vertices named by all possible length-`n` strings of 0s and 1s. Two
vertices are adjacent iff they differ in exactly one position.
"""
function Cube(n::Integer=3)
    G = StringGraph()
    for u=0:2^n-1
        for shift=0:n-1
            v = xor( (1<<shift), u)
            add!(G,string(u,base=2,pad=n), string(v,base=2,pad=n))
        end
    end
    name(G,"Cube graph Q($n)")
    return G
end

# Create the BuckyBall graph

"""
`BuckyBall()` returns the Bucky ball graph.
"""
function BuckyBall()
    G = IntGraph()
    edges = [(1,3), (1,49), (1,60), (2,4), (2,10), (2,59),
	     (3,4), (3,37), (4,18), (5,7), (5,9), (5,13),
	     (6,8), (6,10), (6,17), (7,8), (7,21), (8,22),
	     (9,10), (9,57), (11,12), (11,13), (11,21), (12,28),
	     (12,48), (13,14), (14,47), (14,55), (15,16), (15,17),
	     (15,22), (16,26), (16,42), (17,18), (18,41), (19,20),
	     (19,21), (19,27), (20,22), (20,25), (23,24), (23,32),
	     (23,35), (24,26), (24,39), (25,26), (25,31), (27,28),
	     (27,31), (28,30), (29,30), (29,32), (29,36), (30,45),
	     (31,32), (33,35), (33,40), (33,51), (34,36), (34,46),
	     (34,52), (35,36), (37,38), (37,41), (38,40), (38,53),
	     (39,40), (39,42), (41,42), (43,44), (43,47), (43,56),
	     (44,46), (44,54), (45,46), (45,48), (47,48), (49,50),
	     (49,53), (50,54), (50,58), (51,52), (51,53), (52,54),
	     (55,56), (55,57), (56,58), (57,59), (58,60), (59,60),
	     ]
    for e in edges
        add!(G,e[1],e[2])
    end
    name(G,"Buckyball graph")
    return G
end

# The Kneser graph Kneser(n,k) has C(n,k) vertices that are the
# k-element subsets of 1:n in which two vertices are adjacent if (as
# sets) they are disjoint. The Petersen graph is Kneser(5,2).

# import IterTools.subsets

"""
`Kneser(n,m)` creates the Kneser graph whose vertices are all the
`m`-element subsets of `1:n` in which two vertices are adjacent iff
they are disjoint.
"""
function Kneser(n::Int,k::Int)
    A = collect(1:n)
    vtcs = [Set(v) for v in IterTools.subsets(A,k)]
    G = SimpleGraph{Set{Int}}()

    for v in vtcs
        add!(G,v)
    end

    n = length(vtcs)
    for i=1:n-1
        u = vtcs[i]
        for j=i+1:n
            v = vtcs[j]
            if length(intersect(u,v))==0
                add!(G,u,v)
            end
        end
    end
    name(G,"Kneser($n,$k) graph")
    return G
end

"""
`Johnson(n,k)` creates the Johnson graph whose vertices
and the `k`-element subsets of `{1,2,...,n}`. Vertices `v`
and `w` are adjacent if their intersection has size `k-1`.
"""
function Johnson(n::Int, k::Int)
    if k<0 || n<0
        error("n,k must be nonnegative")
    end

    if k>n
        error("k must not be greater than n")
    end

    A = collect(1:n)
    vtcs = [Set{Int}(v) for v in IterTools.subsets(A,k)]
    G = SimpleGraph{Set{Int}}()

    for v in vtcs
        add!(G,v)

        for j=1:n
            if !in(j,v)
                for x in v
                    w = deepcopy(v)
                    push!(w,j)
                    pop!(w,x)
                    add!(G,v,w)
                end
            end
        end
    end

    name(G,"Johnson($n,$k) graph")
    return G
end

# Create the Petersen graph.
"""
`Petersen()` returns the Petersen graph. The vertices are labeled as
the 2-element subsets of `1:5`. Wrap in `relabel` to have vertices
named `1:10`. See also: `Kneser`.
"""
function Petersen()
  G = Kneser(5,2)
  name(G,"Petersen graph")
  return G
end

# Create Paley graphs

"""
`Paley(p)` creates the Paley graph with `p` vertices named
`0:(p-1)`. Here `p` must be a prime with `p%4==1`. Vertices `u` and
`v` are adjacent iff `u-v` is a quadratic residue (perfect square)
modulo `p`.
"""
function Paley(p::Int)
    if mod(p,4) != 1 || ~isprime(p)
        error("p must be a prime congruent to 1 mod 4")
    end

    # Quadratic residues mod p
    qrlist = unique( [ mod(k*k,p) for k=1:p ] )

    G = IntGraph()
    for u = 0:p-1
        for k in qrlist
            v = mod(u+k,p)
            add!(G,u,v)
        end
    end
    name(G,"Paley($p) graph")
    return G
end

"""
`Frucht()` returns the Frucht graph: A 12-vertex, 3-regular
graph with no non-nontrivial automorphisms.
"""
function Frucht()
  G = Cycle(12)
  more_edges = [ (1,6), (2,4), (3,11), (5,7), (8,10), (9,12) ]
  add_edges!(G,more_edges)
  name(G,"Frucht graph")
  return G
end

# Called by RandomRegular ... one step
function RandomRegularBuilder(n::Int, d::Int)
    # caller has already checked the values of n,d are legit
    vlist = randperm(n*d)
    G = IntGraph(n*d)
    for v=1:2:n*d
        add!(G,vlist[v], vlist[v+1])
    end

    for v = n:-1:1
        mushlist = collect( d*(v-1)+1 : v*d )
        for k=d-1:-1:1
            contract!(G,mushlist[k],mushlist[k+1])
        end
    end
    return relabel(G)
end

"""
`RandomRegular(n,d)` creates a random `d`-regular graph on `n`
vertices. This can take a while especially if the arguments are
large. Call with an optional third argument to activate verbose
progress reports: `RandomRegular(n,p,true)`.
"""
function RandomRegular(n::Int, d::Int, verbose::Bool=false)
    # sanity checks
    if n<1 || d<1 || (n*d)%2==1
        error("n,d must be positive integers and n*d even")
    end
    if verbose
        println("Trying to build ", d, "-regular graph on ",
                n, " vertices")
        count::Int = 0
    end

    while true
        if verbose
            count += 1
            println("Attempt ", count)
            tic()
        end
        g = RandomRegularBuilder(n,d)
        if verbose
            toc();
        end
        dlist = deg(g)
        if dlist[1] == dlist[n]
            if verbose
                println("Success")
            end
            name(g,"Random regular graph, d=$d")
            return g
        end
        if verbose
            println("Failed; trying again")
        end
    end
end

"""
`RandomSBM(bmap,pmat)` creates a random stochastic block model random graph.
The vector `bmap` is a list of `n` positive integers giving the block number
of vertices `1:n`. The `i,j`-entry of the matrix `pmat` gives the probability
of an edge from a vertex in block `i` to a vertex in block `j`.

`RandomSBM(n,pvec,pmat)` creates such a graph with `n` vertices. The vector
`pvec` gives the probabilities that vertices fall into a given block.
"""
function RandomSBM(bmap::Vector{Int}, pmat::Array{S,2}) where {S<:Real}
  n = length(bmap)    # no of vertices
  b = maximum(bmap)   # no of blocks

  @assert minimum(bmap)>0 "Block numbers must be positive"
  @assert pmat == pmat' "Probability matrix must be square and symmetrical"
  r = size(pmat,1)
  @assert b<=r "Insufficient rows/cols in probability matrix"
  @assert minimum(pmat)>=0 && maximum(pmat)<=1 "Entries in probability matrix out of range"

  G = IntGraph(n)
  for u=1:n-1
    bu = bmap[u]
    for v=u+1:n
      bv = bmap[v]
      p = pmat[bu,bv]
      if rand() <= p
        add!(G,u,v)
      end
    end
  end
  name(G,"Random stochastic block model graph")
  return G
end

function RandomSBM(n::Int, pvec::Vector{S}, pmat::Array{T,2}) where {S<:Real,T<:Real}
  @assert minimum(pvec)>=0 "Entries in pvec must be nonnegative"
  @assert sum(pvec)==1 "Entries in pvec must sum to 1"
  bmap = [ random_choice(pvec) for _=1:n]
  return RandomSBM(bmap,pmat)
end

"""
`Knight(r::Int=8,c::Int=8)` creates a Knight's Moves graph on a
`r`-by-`c` grid. That is, the vertices of this graph are the squares
of an `r`-by-`c` chess board. Two vertices are adjacent if a Knight
can go from one of these squares to the other in a single move.
"""
function Knight(r::Int=8,c::Int=8)
    vtcs = collect(Base.Iterators.product(1:r,1:c))
    G = SimpleGraph{Tuple{Int64,Int64}}()
    for v in vtcs
        add!(G,v)
    end

    for v in vtcs
        for w in vtcs
            xv = collect(v)
            xw = collect(w)
            z = sort(map(abs,xv - xw))
            if z==[1,2]
                add!(G,v,w)
            end
        end
    end
    name(G,"Knight's tour graph on $r-by-$c chess board")
    return G
end

"""
`HoffmanSingleton()` creates the Hoffman-Singleton graph. This is
a 7-regular graph whose diameter is 2 and whose girth is 5.
"""
function HoffmanSingleton()
  G = StringGraph()
  # P-pentagons
  for i=0:4
    for j=0:4
      jj = mod(j+1,5)
      v = "P$i$j"
      w = "P$i$jj"
      add!(G,v,w)
    end
  end

  # Q-pentagrams
  for i=0:4
    for j=0:4
      jj = mod(j+2,5)
      v = "Q$i$j"
      w = "Q$i$jj"
      add!(G,v,w)
    end
  end

  # Connections
  for i=0:4
    for j=0:4
      for k=0:4
        v = "P$i$j"
        x = mod(i*k+j,5)
        w = "Q$k$x"
        add!(G,v,w)
      end
    end
  end
  name(G,"Hoffman-Singleton graph")
  return G
end

"""
`Hoffman()` creates the Hoffman graph which is cospectral with,
but not isomorphic to, `Cube(4)`.
"""
function Hoffman()
  D = [
  1 1 1 1 0 0 0 0
  1 1 1 0 1 0 0 0
  1 0 0 1 0 1 1 0
  0 1 0 1 0 1 0 1
  0 0 1 1 0 0 1 1
  1 0 0 0 1 1 1 0
  0 1 0 0 1 1 0 1
  0 0 1 0 1 0 1 1]

  A = [0*D D; D' 0*D]

  G = SimpleGraph(A)
  name(G,"Hoffman graph")
  return G
end


text:
import Base: Meta
const Self = @__MODULE__
const Template_mtb = IdDict(
    GlobalRef(Base, :materialize) => GlobalRef(Self, :mtb_materialize),
    GlobalRef(Base, :materialize!) => GlobalRef(Self, :mtb_materialize!),
    GlobalRef(Self, :tab_materialize) => GlobalRef(Self, :mtab_materialize),
)
const Template_tab = IdDict(
    GlobalRef(Base, :materialize) => GlobalRef(Self, :tab_materialize),
    GlobalRef(Self, :mtb_materialize) => GlobalRef(Self, :mtab_materialize),
)
const Template_mtab = IdDict(
    GlobalRef(Base, :materialize!) => GlobalRef(Self, :mtb_materialize!),
    GlobalRef(Base, :materialize) => GlobalRef(Self, :mtab_materialize),
)

const Template_lzb = IdDict(
    GlobalRef(Base, :materialize) => GlobalRef(Base.Broadcast, :instantiate),
)
## general function
using MacroTools
# replace return
keepreturn(ex) = Meta.isexpr(ex, :return) ? Expr(:call, :𝐫𝐞𝐭𝐮𝐫𝐧, ex.args[1]) : ex
function expandbroadcast(ex)
    Meta.isexpr(ex, :call) || return ex
    if ex.args[1] === :broadcast
        ex.args[1] = GlobalRef(Base, :broadcasted)
        return :($(GlobalRef(Base, :materialize))($ex))
    elseif ex.args[1] === :broadcast!
        ex′ = Expr(:call, GlobalRef(Base, :broadcasted), ex.args[2], ex.args[4:end]...)
        return :($(GlobalRef(Base, :materialize!))($(ex.args[3]), $ex′))
    end
    ex
end
replacetemp(ex, Temp::Symbol) = ex
replacetemp(ex::Core.Slot, Temp::Symbol) = Symbol(Temp, "_slot_", ex.id)
replacetemp(ex::Core.SSAValue, Temp::Symbol) = Symbol(Temp, "_ssa_", ex.id)
replacetemp(ex::Core.GotoIfNot, Temp::Symbol) = Core.GotoIfNot(replacetemp(ex.cond, Temp), ex.dest)
replacetemp(ex::Core.ReturnNode, Temp::Symbol) = Core.ReturnNode(replacetemp(ex.val, Temp))
replacetemp(ex::Expr, Temp::Symbol) = (ex.args .= replacetemp.(ex.args, Temp); ex)

function lowerhack(mod::Module, ex::Expr, temps, addarg = nothing)
    ex = MacroTools.postwalk(keepreturn ∘ expandbroadcast, ex)
    ci = Meta.lower(mod, ex)
    ci.head === :error && error(ci.args[1])
    code = only(ci.args).code
    Temp, tag = gensym(:𝐭𝐞𝐦𝐩), gensym(:𝐠𝐨𝐭𝐨)
    labels = Int[]
    for I in eachindex(code)
        ex = replacetemp(code[I], Temp)
        if ex isa Core.ReturnNode
            code[I] = :($(ex.val))
        elseif ex isa Core.GotoNode
            label = Symbol(tag, ex.label)
            code[I] = :(@goto $label)
            push!(labels, ex.label)
        elseif ex isa Core.GotoIfNot
            dest = Symbol(tag, ex.dest)
            code[I] = :($(ex.cond) || @goto $dest)
            push!(labels, ex.dest)
        elseif Meta.isexpr(ex, :call)
            func = ex.args[1]
            if func === :𝐫𝐞𝐭𝐮𝐫𝐧
                code[I] = :(return $(ex.args[2]))
            else
                ex.args[1] = get(temps, func, func)
                !isnothing(addarg) && ex.args[1] !== func && push!(ex.args, addarg)
                Temp′ = replacetemp(Core.SSAValue(I), Temp)
                code[I] = :(local $Temp′ = $ex)
            end
        elseif !Meta.isexpr(ex, :(=))
            Temp′ = replacetemp(Core.SSAValue(I), Temp)
            code[I] = :(local $Temp′ = $ex)
        end
    end
    for I in unique!(labels)
        code[I] = :(@label $(Symbol(tag, I)); $(code[I]))
    end
    exblock = Expr(:block)
    exblock.args = code
    exblock
end

"""
**Example:**
```julia
@mtb @. a = sin(c)
@mtb a = sin.(c)
@mtb broadcast!(sin,a,c)
@mtb a = broadcast!(sin,c)
@mtb 2 a = broadcast!(sin,c) # use 2 threads
```
"""
macro mtb(args...)
    global Template_mtb
    na = length(args)
    if na == 1
        nthread, ex = num_threads(), args[1]
    elseif na == 2
        nthread, ex = args
    else
        error("Invalid input")
    end
    (!isa(nthread, Integer) || nthread <= 1) && return esc(ex)
    nthread = min(nthread, Threads.nthreads())
    return esc(lowerhack(Self, ex, Template_mtb, nthread))
end

"""
**Example:**
```julia
@tab @. (a,b) = sincos(c)
@tab @. a,b = sincos(c)
@tab (a,b) = @. sincos(c)
```
**Note:**
@tab is needed only for non in-place situation.
"""
macro tab(ex)
    global Template_tab
    esc(lowerhack(Self, ex, Template_tab))
end

"""
    @mtab [n] ex = @mtb [n] @tab ex
"""
macro mtab(args...)
    global Template_mtab
    global Template_tab
    na = length(args)
    if na == 1
        nthread, ex = num_threads(), args[1]
    elseif na == 2
        nthread, ex = args
    else
        error("Invalid input")
    end
    (!isa(nthread, Integer) || nthread <= 1) &&
        return esc(lowerhack(Self, ex, Template_tab))
    nthread = min(nthread, Threads.nthreads())
    return esc(lowerhack(Self, ex, Template_mtab, nthread))
end

"""
    @lzb  return lazy Broadcasted object
"""
macro lzb(ex)
    global Template_lzb
    esc(lowerhack(Self, ex, Template_lzb))
end


text:
# ----------------------------------------------------
## --- Channel definition
# ----------------------------------------------------
# Object with all channel attribute (but no realisations)
"""
---
Channel model object
# --- Syntax
- profile		  : Name of profile [String]
- carrierFreq	  : Carrier frequency [Union{Int,Float64}]
- samplingFreq	  : Sampling frequency [Union{Int,Float64}]
- speed			  : Desired speed (km/h) [Union{Int,Float64}]
- powerProfile	  : Distribution of power values in dB [Array{Float64}]
- delayProfile	  : Distribution of delay values in s [Array{Float64}] (same size as powerProfile)
- dopplerFreq	  : Doppler frequency (inherited from samplingFreq  and speed)
- delaySpread	  : Max support of CIR in samples [Int]
# ---
"""
struct ChannelModel
	profile::String;					  # Name of profile
	carrierFreq::Union{Int,Float64};	  # Carrier frequency (Hz)
	samplingFreq::Union{Int,Float64};	  # Sampling frequency (Hz)
	speed::Union{Int,Float64};			  # Velocity (km/h)
	powerProfile::Union{Array{Int},Array{Float64}};		  # Vector of power value (dB)
	delayProfile::Union{Array{Int},Array{Float64}};		  # Vector of delay value (s or samples)
	dopplerFreq::Union{Int,Float64};	  # Doppler frequency in Hz
	#randSeed::Int;					  	  # Seed for random variables %
	delaySpread::Int;					  # Delay spread value (samples)
end

#
"""
---
Object with channel realisations
# --- Syntax
- timeVarying	: Flag for constant vs time varying channel [Int]
- cir			: CIR matrix (nbTap x nbChannel) [Union{Array{Int},Array{Float64},Array{Complex{Float64}}}]
- channelModel  : Model use for generation [DigitalComm.ChannelModel]
- powerLin		: Power distribution with the rayleigh distribution [Array{Complex{Float64}}]
- randSeed		: Seed use for rayleigh generation
# ---
# v 1.0 - Robin Gerzaguet.
"""
struct ChannelImpl
	timeVarying::Int;					  # Time varing or constant channel
	cir::Union{Array{Int},Array{Float64},Array{Complex{Float64}}}; # CIR matrix (nbTap x nbReal)
	channelModel::ChannelModel;			  # Channel model baseline
	powerLin::Array{Complex{Float64}};	  # Power distribution (with random rayleigh application)
	randSeed::Int;						  # Used seed
end

# ----------------------------------------------------
## ---  Channel initialisation
# ----------------------------------------------------
"""
---
 Create a channel object associated to physical parameters and propagation profile.
# --- Syntax
channObj  = initChannel(profile,carrierFreq,samplingFreq,dopplerFreq;powerProfile=[],delayProfile=[],randSeed=-1)
# --- Input parameters
- profile : Multipath profile (see below) [String]
- carrierFreq : Carrier frequency in Hz [Union{Int,Float64}]
- samplingFreq: baseband sampling frequency [Union{Int,Float64}]
- speed : Velocity (km/h) [Float64]
- powerProfile	: User defined power profile (default []). In case when it is defined, delayProfile should also be defined and profile is not used. [Array{Float}]
- delayProfile	: User defined delay profile in second [Array{Float64}]
- randSeed : Seed if necessary (default: -1) [Float64]
# --- Output parameters
- channel : Julia Channel Object [Channel]
# --- Channel models
-	ETU	  : Extended Typical Urban
-	EVA	  : Extended Vehicular area
-	TDLC  :
-	Rayleigh : Single tap rayleigh channel model
# ---
"""
function initChannel(profile::String,carrierFreq::Union{Int,Float64},samplingFreq::Union{Int,Float64},speed::Union{Int,Float64};powerProfile=[],delayProfile=[],tau=1)
	# ----------------------------------------------------
	# --- Init parameters
	# ----------------------------------------------------
	# --- Creating profile without case sensitive
	profile = lowercase(profile);
	# ----------------------------------------------------
	# --- Channel model switch
	# ----------------------------------------------------
	if profile == "multipath"
		# ----------------------------------------------------
		# --- Custom multipath channel profile
		# ----------------------------------------------------
		# In this case, the profile is direclty inherited from powerProfile and delayProfile array
		# --- Calculating doppler frequency
		dopplerFreq	  = speed * 1e3 / 3600 * carrierFreq / 3e8;
        # --- Estimating delay profile
        delaySpread = Int( ceil( maximum(delayProfile) * samplingFreq))
		# --- Creating object
		return ChannelModel("multipath",carrierFreq,samplingFreq,speed,powerProfile,delayProfile,dopplerFreq,delaySpread);
	elseif profile == "none" || profile == "awgn"
		# ----------------------------------------------------
		# --- AWGN channel model
		# ----------------------------------------------------
		# AWGN channel model, constant with no doppler frequency
		return ChannelModel("none",carrierFreq,samplingFreq,0,[1.0],[1.0],0,1);
	elseif profile == "exponential"
		# ----------------------------------------------------
		# --- Exponential rayleigh profile
		# ----------------------------------------------------
		# --- Setting doppler frequency
		dopplerFreq	  = speed * 1e3 / 3600 * carrierFreq / 3e8;
		# --- Creating random profile
		# We use tau as a specific parameter for exponential decay profile
		# This is a constant time --> Channel spread on 5tau
		# tau can be given in s or in samples. We assume that if tau > 1, it is in samples.
		if tau < 1
			# --- second to samples conversion
			tau	  = tau * samplingFreq;
		end
		# --- Channel delay spread
		delaySpread	  = 5 * tau;
		# --- Enveloppe
		enveloppe	  = 10*log10.(exp.(-(0:1:delaySpread-1) ./ tau));
		# --- Creating object
		return ChannelModel("exponential",carrierFreq,samplingFreq,speed,enveloppe,collect(1:delaySpread),dopplerFreq,delaySpread);
	elseif profile == "flatrayleigh"
		# ----------------------------------------------------
		# --- Flat fading channel
		# ----------------------------------------------------
		# --- Calculating doppler frequency
		dopplerFreq	  = speed * 1e3 / 3600 * carrierFreq / 3e8;
		# --- Creating object
		return ChannelModel("flatRayleigh",carrierFreq,samplingFreq,speed,[1],[1],dopplerFreq,1);
	elseif profile == "eva"
		# ----------------------------------------------------
		# ---EVA channel
		# ----------------------------------------------------
		# --- Calculating doppler frequency
		dopplerFreq	    = speed * 1e3 / 3600 * carrierFreq / 3e8;
		# --- EVA power and delay profile
		# second
		powerProfile			= [0.0; -1.5; -1.4; -3.6; -0.6; -9.1; -7.0; -12.0; -16.9];
		tapDelaySecond		= [0; 30; 150; 310; 370; 710; 1090; 1730; 2510] * 1e-9;
		delaySpread		    = ceil(maximum(tapDelaySecond) * samplingFreq);
		return ChannelModel("EVA",carrierFreq,samplingFreq,speed,powerProfile,tapDelaySecond,dopplerFreq,delaySpread);
	elseif profile == "consteva"
		# ----------------------------------------------------
		# ---EVA channel
		# ----------------------------------------------------
		# --- EVA power and delay profile
		# second
		powerProfile			= [0.0; -1.5; -1.4; -3.6; -0.6; -9.1; -7.0; -12.0; -16.9];
		tapDelaySecond		= [0; 30; 150; 310; 370; 710; 1090; 1730; 2510] * 1e-9;
		delaySpread		    = ceil(maximum(tapDelaySecond) * samplingFreq);
		return ChannelModel("EVA",carrierFreq,samplingFreq,0,powerProfile,tapDelaySecond,0,delaySpread);

	elseif profile == "etu"
		# ----------------------------------------------------
		# --- ETU fading channel
		# ----------------------------------------------------
		# --- Calculating doppler frequency
		dopplerFreq	    = speed * 1e3 / 3600 * carrierFreq / 3e8;
		# --- EVA power and delay profile
		# second
		powerProfile		= [-1.0; -1.0; -1.0; +0.0; +0.0; +0.0; -3.0; -5.0; -7.0];
		tapDelaySecond		= [0;   50;  120;  200;  230;  500; 1600;  2300; 5000] * 1e-9;
		delaySpread		    = ceil(maximum(tapDelaySecond) * samplingFreq);
		return ChannelModel("ETU",carrierFreq,samplingFreq,speed,powerProfile,tapDelaySecond,dopplerFreq,delaySpread);
	elseif profile == "constetu"
		# ----------------------------------------------------
		# --- ETU fading channel: w/o modification even with speed
		# ----------------------------------------------------
		# --- EVA power and delay profile
		# second
		powerProfile			= [-1.0;-1.0; -1.0; +0.0; +0.0; +0.0; -3.0; -5.0; -7.0];
		tapDelaySecond		= [0;   50;  120;  200 ; 230;  500; 1600;  2300; 5000] * 1e-9;
		delaySpread		    = ceil(maximum(tapDelaySecond) * samplingFreq);
		return ChannelModel("ETU",carrierFreq,samplingFreq,0,powerProfile,tapDelaySecond,0,delaySpread);

	end
end

"""
---
Returns a complete FIR response based on the power profile (in linear scale), the delay profile (in s) and the sampling frequency.
# --- Syntax
cir = getFIRResponse(delayProfile:,powerProfile,freq,sincSupport=5,interpSystem=0)
# --- Input parameters
- delayProfile	: Vector of delay in second [Union{Array{Int},Array{Float64}}]
- powerProfile	: Vector of attenuation with same size of delayProfile in linear scale with rayleigh distribution
- freq			: Sampling frequency
- sincSupport	: Size of interpolator (default 5)
- interpSystem	: Forcing extra delay for all samples to have proper interpolation of the FIR beginning (default 0)
# --- Ouput parameters
- cir			: Finite impulse response [Array{Complex{Float64}},max(delayProfile)+sincSupport]
# ---
"""
function getFIRResponse(delayProfile::Union{Array{Int},Array{Float64}},powerProfile,freq,sincSupport=5,interpSystem=0)
	# --- Instantiate grid
	# FIXME: A signal processing question here.
	# % We should do fractional interpolation as the delay provided in s
	# % will not match the time grid (and depends on the sampling frequency).
	# % I propose Shannon interpolation (and not Lagragian here) in order to have something
	# % slightly different with the farrow interpolation that can be done in practical Rx
	# % However, it leads to the definition of a sinc support as the spread of the interpolator is infinite.
	# % 5 additional taps for sinc should be sufficient
	# % Two ways to to this:
	# %	  (0) - Forcing 5 taps before and after: proper interpolator but all delays will be shifted by 5Te
	# %	  (1) - Keep all the same delays values: it means that the interpolation will not be perfect @ the first samples
	# % interpSystem = 0 or 1
	if Bool(interpSystem)
		# No additional delay
		sS		= sincSupport;			# Half support (only append @ end)
		fdGap	= 0;					# We keep the delay grid
	else
		# Grid is shifted to enable FD (fractional delay) interpolation
		sS		= 2*sincSupport;		# Support is doubled (before and after delay grid)
		fdGap	= sincSupport / freq;	# We shift all the delay value to be sure first samples are FD.
	end
	# Final CIR with size of nbTaps + overhead due to sinc. interpolation
	nbTaps	  = length(delayProfile);
	delaySpread = Int(ceil(delayProfile[end]*freq));
	cir		  = zeros(Complex{Float64}, delaySpread + sS);
	# --- Ensuring additional delay for proper interpolation
	delay	  = delayProfile .+ fdGap;
	# sinc. interpolation due to delay in s (not in grid)
	for k = 1 : 1 : delaySpread + sS
		# --- Each index is superimposition of all timed compondent
		for ∂ = 1 : 1 : nbTaps
			# --- Shannon interpolation with sinc. on sincSupport
			cir[k] += powerProfile[∂] * sinc.(pi*(k/freq - delay[∂]) * freq);
		end
	end
	return cir;
end

# ----------------------------------------------------
# --- Create and instanciate channel model
# ----------------------------------------------------
"""
---
Create a channelImpl based on desired channel model and number of realisations
# --- Syntax
channelImpl	  = getChannel(nbSamples,channelModel,randSeed=-1)
# --- Input parameters
- nbSamples		: Number of samples on which channel will be applied [Int]
- channelModel	: Channel object [ChannelModel]
- randSeed		: Desired seed (default -1)
# --- Output parameters
- channelImpl	: Channel implementation
# ---
"""
function getChannel(nbSamples::Int,channelModel::ChannelModel,randSeed=-1)
	# ----------------------------------------------------
	# --- Constant or varying channel
	# ----------------------------------------------------
	if channelModel.speed == 0 || channelModel.dopplerFreq == 0
		timeVarying = 0;
	else timeVarying = 1;
	end
	# ----------------------------------------------------
	# --- Other parameters
	# ----------------------------------------------------
	sincSupport	  = 5;	  # Extra support for sinc pulse shape (in samples)
	# --- Setting the random seed if given
	if randSeed != -1
		# --- Seed as the second parameter
		Random.seed!(randSeed);
	end
	if channelModel.profile == "none"
		# ----------------------------------------------------
		# --- Special case of none channel
		# ----------------------------------------------------
		# Channel is unitary, no power loss (no rayleigh distrib)
		cir = [1.0];
		return ChannelImpl(0,cir,channelModel,[1.0],randSeed);
	else
		if Bool(timeVarying)
			# ----------------------------------------------------
			# --- Time varying channel
			# ----------------------------------------------------
			# FIXME: A signal processing question here.
			# % We should do fractional interpolation as the delay provided in s
			# % will not match the time grid (and depends on the sampling frequency).
			# % I propose Shannon interpolation (and not Lagragian here) in order to have something
			# % slightly different with the farrow interpolation that can be done in practical Rx
			# % However, it leads to the definition of a sinc support as the spread of the interpolator is infinite.
			# % 5 additional taps for sinc should be sufficient
			# % Two ways to to this:
			# %	  (0) - Forcing 5 taps before and after: proper interpolator but all delays will be shifted by 5Te
			# %	  (1) - Keep all the same delays values: it means that the interpolation will not be perfect in the first samples
			# % interpSystem = 0 or 1
			interpSystem  = 1;
			if Bool(interpSystem)
				# No additional delay
				sS		= sincSupport;			# Half support (only append @ end)
				fdGap	= 0;					# We keep the delay grid
			else
				# Grid is shifted to enable FD (fractional delay) interpolation
				sS		= 2*sincSupport;		# Support is doubled (before and after delay grid)
				fdGap	= sincSupport / freq;	# We shift all the delay value to be sure first samples are FD.
			end
			# --- Generate realisation of rayleigh distribution
			nbTaps	  = length(channelModel.powerProfile);
			delaySpread = Int(ceil(channelModel.delayProfile[end]*channelModel.samplingFreq));
			# --- Creating amplitudes
			powerLin  = 10 .^(channelModel.powerProfile / 10);
			# --- Correlated rayleigh amplitude
			raylAmpl  = rayleighChan(nbTaps,nbSamples,channelModel.samplingFreq,channelModel.dopplerFreq,randSeed);
			# --- Getting sampling frequency
			freq	  = channelModel.samplingFreq;
			# --- Ensuring additional delay for proper interpolation
			delay	  = channelModel.delayProfile .+ fdGap;
            # --- Shannon interpolation with sinc. on sincSupport
            shannonInterp = sinc.(pi.*(transpose(1:delaySpread+sS)./freq .- delay) .* freq);
            # --- Each index is superimposition of all timed compondent
            cir = transpose(powerLin) .* raylAmpl[1:nbSamples] * shannonInterp;
			return ChannelImpl(1,cir,channelModel,powerLin,randSeed);
		else
			# ----------------------------------------------------
			# --- Constant FIR channel
			# ----------------------------------------------------
			# --- Generate realisation of rayleigh distribution
			nbTaps	  = length(channelModel.powerProfile);
			raylComp  = 1/sqrt(2)*(randn(nbTaps) + 1im .* randn(nbTaps));
			# --- Creating amplitudes
			powerLin  = 10 .^(channelModel.powerProfile / 10) .* raylComp;
			# --- Create FIR
			cir	  = getFIRResponse(channelModel.delayProfile,powerLin,channelModel.samplingFreq);
			# --- Create object
			return ChannelImpl(0,cir,channelModel,powerLin,randSeed);
		end
	end
end

# --- Channel model
"""
---
Apply a channel implementation  to an input signal
# --- Syntax
sigChan	  : applyChannel(sigId,channelImpl)
# --- Input parameters
- sigId	  : Input signal [Array{Any}]
- channelImpl : Channel implementation [ChannelImpl]
# --- Output parameters
- sigChan	: Output signal [Array{Any}]
# ---
# v 1.0 - Robin Gerzaguet.
"""

function applyChannel(sigId,channelImpl)

	if !Bool(channelImpl.timeVarying)
		# --- Classic convolution
		return sigChan = conv(sigId,channelImpl.cir);
	else
		# --- Time vraying convolution
		# Check size
		if length(sigId) != size(channelImpl.cir,1)
			error("Time varying channel should have same size as incoming signal");
		end
		# Init outut vector: Output is a convolution
		nbTaps	=size(channelImpl.cir,2);
		nbOut	=length(sigId)+nbTaps
		sigChan = zeros(Complex{Float64},nbOut)
        timeVaryingConv!(sigChan, channelImpl.cir, sigId, nbTaps)

        # At the end we have zero
		buffEnd = zeros(Complex{Float64},nbTaps *2);
		buffEnd[1:nbTaps] = sigChan[1:nbTaps];
        timeVaryingConvTail!(sigChan, channelImpl.cir, buffEnd, length(sigId), nbTaps)

		return sigChan;
	end
end

function timeVaryingConv!(sigChan, cir, sigId, nbTaps)
    # n is time index of output, δ is lag index
    for n = 1:length(sigId)
        # --- At beginning, not possible to have all index
        nbMaxTap = min(n,nbTaps);
        # --- Classic time domain convolution
        for δ = 1:nbMaxTap
            @inbounds sigChan[n] += cir[n,δ] * sigId[n-δ+1];
        end
    end
end

function timeVaryingConvTail!(sigChan, cir, buffEnd, nbSigSamples, nbTaps)
    for n = 1:nbTaps
        # --- Classic time domain convolution
        for δ = 1 : 1 : nbTaps
            @inbounds sigChan[nbSigSamples+n] += cir[end,δ] * buffEnd[n+nbTaps-δ+1];
        end
    end
end


text:
function test()::Int64
    a::Array{Int64} = [1, 2, 3]
    return a[1+1]
end

function main()
    b = test()
    @assert(b == 2)
    println(join(["OK"], " "))
end

main()


text:
facts("Goldencross trading logic") do
  context("Market state") do
    fm = TradingLogic.goldencrossmktstate
    @fact fm(120.0, 50.0) --> :trendup
    @fact fm(20.0, 50.0) --> :trenddown
    @fact fm(20.0, 20.0) --> :undefined
    @fact fm(NaN, 20.0) --> :undefined
    @fact fm(20.0, NaN) --> :undefined
    @fact fm(NaN, NaN) --> :undefined
  end
  context("Target position") do
    tq = 100
    ft(mkt, posnow) = TradingLogic.goldencrossposlogic(
      mkt, tq, [posnow])
    stlim = Array(Float64, 0)

    # up, zero position
    @fact ft(:trendup, 0) --> (tq, stlim)
    # up, position less than target (e.g. partial fill)
    p = round(Int64, tq/2)
    @fact ft(:trendup, p) --> (tq - p, stlim)

    # down: sell
    @fact ft(:trenddown, tq) --> (-tq, stlim)
    @fact ft(:trenddown, p) --> (-p, stlim)
    # down: nothing left to sell
    @fact ft(:trenddown, 0) --> (0, stlim)

    # hold position in line with the market state
    @fact ft(:trendup, tq) --> (0, stlim)

    # undefined: wait
    @fact ft(:undefined, tq) --> (0, stlim)

    # negative position: should not happen, close it
    @fact ft(:trendup, -5) --> (5, stlim)
    @fact ft(:trenddown, -5) --> (5, stlim)
    @fact ft(:undefined, -5) --> (5, stlim)
  end
end

facts("Goldencross strategy backtesting") do
  context("Boeing stock over 50 years vs. quantstrat") do
    # quantstrat input: OHLC data
    ohlc_BA = TimeSeries.readtimearray(
      rel("quantstrat/goldencross/data/OHLC_BA_2.csv"))
    # parameters should match
    #  test/quantstrat/goldencross/quantstrat_goldencross.R
    mafast = 50
    maslow = 200
    targetqty = 100
    # restrict between start and end dates
    #println(ohlc_BA)
    ##ohlc_test = ohlc_BA[Date(1961,12,31):Date(2010,1,1)]
    # in reality: quantstrat transactions are present
    # past the end date
    date_final = Date(2012,8,31)
    ohlc_test = ohlc_BA[Date(1961,12,31):date_final]
    #println(ohlc_test)
    # NOTE: excluding period around 2012-09-06 where
    # quantstrat behavior differs from TradingLogic by design

    # quantstrat output: transactions
    txnsdf = DataFrames.readtable(
      rel("quantstrat/goldencross/transactions.csv"),
      header = true,
      names = [:datestr, :qty, :prc, :fees, :val, :avgcost, :pl],
      eltypes = [String, Int64, Float64, Float64,
                 Float64, Float64, Float64])[2:end,:]
    # vectors to verify
    vqty = convert(Array, txnsdf[:qty])
    vprc = convert(Array, txnsdf[:prc])
    vpnlcum = cumsum(convert(Array, txnsdf[:pl]))
    # NOTE: quantstrat records transaction times when
    #  signal is fired not when open fill-price is taken
    # adjusting for that
    vdate = Date(DateTime(convert(Array, txnsdf[:datestr]),
                          "yyyy-mm-dd HH:MM:SS"))
    oneday = Day(1)
    for i = 1:length(vdate)
      vdate[i] = vdate[i] + oneday
    end

    s_ohlc = Reactive.Signal((Dates.DateTime(ohlc_test.timestamp[1]),
                              vec(ohlc_test.values[1,:])))
    ohlc_inds = Dict{Symbol,Int64}()
    ohlc_inds[:open] = 1
    ohlc_inds[:close] = 4

    # backtest at next-open price
    # quantstrat fills tracsactions at next open on enter-signal
    s_pnow = Reactive.map(s -> s[2][ohlc_inds[:open]], s_ohlc, typ=Float64)
    blotter = TradingLogic.emptyblotter()

    s_status = TradingLogic.runtrading!(
      blotter, true, s_ohlc, ohlc_inds, s_pnow, 0,
      TradingLogic.goldencrosstarget, targetqty, mafast, maslow)

    s_perf = TradingLogic.tradeperfcurr(s_status)

    for i = 2:length(ohlc_test)
      push!(s_ohlc, (Dates.DateTime(ohlc_test.timestamp[i]),
                     vec(ohlc_test.values[i,:])))
      #println(s_perf.value)
    end
    # no errors
    @fact s_status.value[1] --> true

    #TradingLogic.printblotter(STDOUT, blotter)
    metr = [:DDown]
    vt, perfm = TradingLogic.tradeperf(blotter, metr)
    #println(vt)
    #println(perfm[:PnL])

    # transaction matching
    @fact length(perfm[:Qty]) --> length(txnsdf[:datestr])
    @fact Date(vt) --> vdate
    @fact perfm[:Qty] --> vqty
    @fact perfm[:FillPrice] --> roughly(vprc)

    # profit loss over time
    @fact perfm[:PnL] --> roughly(vpnlcum)
    #println(perfm)

    # quantstrat/goldencross/results_summary.txt
    pnlnet = 2211.0 # Net.Trading.PL
    ddownmax = 17374.0 # Max.Drawdown

    # not final PnL without exit yet
    @fact TradingLogic.tradepnlfinal(blotter) - pnlnet > 10.0 --> true
    # but final PnL if exit price is given: last timestep close price
    pfinal = s_ohlc.value[2][ohlc_inds[:close]]
    @fact TradingLogic.tradepnlfinal(blotter, pfinal) --> roughly(pnlnet)

    # add exit to blotter for cumulative statistics at the end date
    blotter[DateTime(date_final)] = (-sum(perfm[:Qty]), pfinal)

    @fact TradingLogic.tradepnlfinal(blotter) --> roughly(pnlnet)
    # perfm[:DDown] does not have the true max. drawdown
    # as it is only blotter timesteps based;
    # use performance metrics signal for that
    @fact s_perf.value[2] --> roughly(ddownmax)
  end
end


text:
"""
Delegates loading functions & traits for the HARPS spectrograph

Author: Alex Wise and collaborators
Created: April 2021
"""

"""
Module providing types and traits and customized functions for the HARPS Spectrograph.
- https://www.eso.org/sci/facilities/lasilla/instruments/harps.html
"""
module HARPS
using RvSpectMLBase
import RvSpectMLBase: AbstractInstrument, AbstractInstrument1D, AbstractInstrument2D
import ..EchelleInstruments
import ..EchelleInstruments: read_manifest, read_metadata_from_fits, check_metadata_fields_expected_present
#import ..EchelleInstruments: read_header
import ..EchelleInstruments: default_Δv_to_avoid_tellurics

using CSV, DataFrames, Query
using FITSIO
using Interpolations
using NaNMath, Missings
using Polynomials

""" Trait for 1D Extracted spectra from HARPS """
struct HARPS1D <: AbstractInstrument1D end

""" Trait for 2D Extracted spectra from HARPS """
struct HARPS2D <: AbstractInstrument2D end

""" Trait to specify any 1D or 2D Extracted spectra from HARPS """
const AnyHARPS = Union{HARPS1D,HARPS2D}

export HARPS, HARPS1D, HARPS2D, AnyHARPS

# traits.jl imports from RvSpectMBase on its own
include("traits.jl")
export min_order, max_order, min_pixel_in_order, max_pixel_in_order
export orders_to_use_default, min_col_default, max_col_default
export metadata_symbols_default, metadata_strings_default
export default_ccf_mask_v_width
export get_inst_module
export bad_col_ranges

include("io.jl")
# not exported, since don't have parameter that specializes to instrument

#import RvSpectMLBase: continuum_normalize_spectrum, continuum_normalize_spectra  # Not yet implemented
#import RvSpectMLBase: filter_line_list
#import RvSpectMLBase: find_worst_telluric_in_each_chunk # Not yet implemented
import RvSpectMLBase: make_clean_line_list_from_tellurics, choose_obs_idx_for_init_guess
include("util.jl") # TODO
#export continuum_normalize_spectrum!, continuum_normalize_spectra! # Not yet implemented
export filter_line_list
export make_clean_line_list_from_tellurics
export make_λ_list_for_bad_columns
#export choose_obs_idx_for_init_guess

end

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
__precompile__()
module SeisRemoveEQ

include("utils.jl")
include("map_removeEQ.jl")
using .Utils
using .Map_removeEQ

using Distributed
@everywhere using SeisIO, Dates, Printf, JLD2, FileIO

export seisremoveEQ

function seisremoveEQ(InputDict::Dict)

	Utils.initlogo()

	jldopen(InputDict["finame"], "r") do file
		InputDict["DLtimestamplist"] = file["info/DLtimestamplist"];
		InputDict["stationlist"] = file["info/stationlist"];
	end

	fodir = InputDict["fodir"]
	foname = InputDict["foname"]

	if foname == InputDict["finame"]
		error("Input name and output name is identical; which may cause overwrite. Please change the output name. Abort.")
	end

	mkpath(fodir)
    fopath=joinpath(fodir, foname*".jld2")
	InputDict["fopath"] = fopath

	tmppath = joinpath(fodir, "./seisremoveeq_tmp")
	InputDict["tmppath"] = tmppath
	mkpath(tmppath)

	if InputDict["IsOutputRemovalFrac"]
		removal_fractionpath = joinpath(fodir, "./removal_fraction")
		InputDict["removal_fractionpath"] = removal_fractionpath
		mkpath(removal_fractionpath)
	end

	Utils.defaultinputdict!(InputDict)

	#print parameters
	printparams(InputDict)

	printstyled("---Start removing EQ---\n"; color=:cyan, bold=true)

	# choose converting time window
	mapidlist = []
	for i = 1:length(InputDict["DLtimestamplist"])
		y, jd = parse.(Int64, split(InputDict["DLtimestamplist"][i], ".")[1:2])
		m, d = j2md(y,jd)
		curdate=DateTime(y, m, d)
		if !InputDict["IsStartendtime"]
			# convert all time stamp
			push!(mapidlist, i)
		else
			if curdate >= InputDict["starttime"] && curdate <= InputDict["endtime"]
				push!(mapidlist, i)
			else
				continue
			end
		end
 	end

	if isempty(mapidlist)
		error("no data is within start and endtime. abort.")
	end

	InputDict["NumofTimestamp"] = length(mapidlist)

	t_removeeq = @elapsed pmap(x -> map_removeEQ(x, InputDict), mapidlist)

	# convert intermediate file to prescibed file format (JLD2, ASDF, ...)
	InputDict["DLtimestamplist_selected"] = InputDict["DLtimestamplist"][mapidlist]

	t_convert = @elapsed convert_tmpfile(InputDict)

	printstyled("---Summary---\n"; color=:cyan, bold=true)
	println("time to remove EQ    =$(t_removeeq)[s]")
	println("time to convert      =$(t_convert)[s]")
	print("\nprocess done at:")
	println(now())
	str = "EQ is successfully removed from raw data:\nOutput = $fopath\n"
	printstyled(str; color=:green, bold=true)

end

end # module

max_stars_repo_name: kura-okubo/SeisRemoveEQ

max_issues_repo_name: kura-okubo/SeisRemoveEQ

max_forks_repo_name: kura-okubo/SeisRemoveEQ


text:
module Lists

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element

export list, item, item_section, itemsection, item_label, itemlabel

function __init__()
  Genie.Renderer.Html.register_normal_element("q__list", context = Genie.Renderer.Html)
  Genie.Renderer.Html.register_normal_element("q__item", context = Genie.Renderer.Html)
  Genie.Renderer.Html.register_normal_element("q__item__section", context = Genie.Renderer.Html)
  Genie.Renderer.Html.register_normal_element("q__item__label", context = Genie.Renderer.Html)
end

function list(args...;
              wrap::Function = StippleUI.DEFAULT_WRAPPER,
              kwargs...)
  wrap() do
    Genie.Renderer.Html.q__list(args...; kwargs...)
  end
end

function item(args...; kwargs...)
  Genie.Renderer.Html.q__item(args...; attributes([kwargs...], StippleUI.API.ATTRIBUTES_MAPPINGS)...)
end

function item_section(args...; kwargs...)
  Genie.Renderer.Html.q__item__section(args...; attributes([kwargs...], StippleUI.API.ATTRIBUTES_MAPPINGS)...)
end
const itemsection = item_section

function item_label(args...; kwargs...)
  Genie.Renderer.Html.q__item__label(args...; kwargs...)
end
const itemlabel = item_label

end


text:
## Initial parameter estimators ——————————————————————
"
$SIGNATURES

Convert autoregressive prediction parameters `a_ar` estimated from regularly sampled
data at interval `Tr` to SDE parameters `a_sde`. If `rm_negreal` is true,
remove negative real roots; otherwise, throw an error.

Returns SDE parameters `a_sde`, and `negreal_count`, the number of negative real roots.
If `remove_negreal` is `false`, only `a_sde` is returned.
"
function ar2sdepar(a_ar,Tr,remove_negreal)
  p = [one(eltype(a_ar));-a_ar]
  pol = Polynomial(p[end:-1:1])
  rd = roots(pol)
  idx_negreal = map(x->isreal(x) && real(x)<0,rd)
  negreal_count = count(idx_negreal)
  if !remove_negreal && negreal_count>0 error("Negative real roots") end
  rd = rd[.!idx_negreal]
  a_sde = if isempty(rd)
    Float64[]
  else
    r = log.(rd)/Tr
    roots2sdepar(r,def_tol)
  end
  return remove_negreal ? (a_sde,negreal_count) : a_sde
end

is_expanding(problem,proc_prevs′) = problem.p>maximum(order.(proc_prevs′))

abstract type InitialEstimate end

"
Add random real root from an interval

$FIELDS
"
struct AddRoot<:InitialEstimate
  "root magnitude interval"
  i::Interval{Float64}
  "Random number generator"
  rng::MersenneTwister
  function AddRoot(i,rng)
    if !(i≫0) error("Root magnitude interval must be ≫ 0") end
    return new(i,rng)
  end
end

Base.show(io::IO,x::AddRoot) = print(io,"root ~ $(x.i)")

const default_addroot_int = 10.0..100.0

const default_addroot_rng = MersenneTwister(463468)

AddRoot() = AddRoot(default_addroot_int,copy(default_addroot_rng))

randroots(m::AddRoot,n::Int) = -rand(m.rng,Uniform(m.i.first,m.i.last),n)

function initial_estimate(m::AddRoot,problem::ProblemLike,proc_prevs)
  if !is_expanding(problem,proc_prevs) return nothing end
  proc = last(proc_prevs)
  n = problem.p-order(proc)
  r = [SdeRoots(proc).r;
       randroots(m,n)]
  eq = SdeRoots(sortroots(r);σy=problem.σy)
  return SdeProcess(eq)
end

"Mode for adding complex roots"
@enum ComplexRootMode KEEP_PREVIOUS MAX_COMPLEX ALL_RANDOM

"
Add random complex root from an interval

$FIELDS
"
struct AddComplexRoot<:InitialEstimate
  "Real part"
  i_re::Interval{Float64}
  "Imaginary part"
  i_im::Interval{Float64}
  "Root mode"
  rootmode::ComplexRootMode
  "Random Number generator"
  rng::MersenneTwister
  function AddComplexRoot(i_re,i_im,rootmode,rng)
    if !(i_im≫0) error("Imaginary part interval must be ≫ 0") end
    return new(i_re,i_im,rootmode,rng)
  end
end

function Base.show(io::IO,x::AddComplexRoot)
  print(io,x.rootmode," root ~ $(x.i_re)+$(x.i_im)im")
end

const default_addimagroot_int = 1e-2..10

function AddComplexRoot(i_re=default_addroot_int,i_im=default_addimagroot_int)
  rng = copy(default_addroot_rng)
  return AddComplexRoot(i_re,i_im,MAX_COMPLEX,rng)
end

function randroots(m::AddComplexRoot,n::Integer)
  nc = n÷2
  rc = Vector{ComplexF64}(undef,2nc)
  for i in 1:nc
    a = -rand(m.rng,Uniform(m.i_re.first,m.i_re.last))
    b = rand(m.rng,Uniform(m.i_im.first,m.i_im.last))
    rc[2(i-1)+1] = a-b*im
    rc[2(i-1)+2] = a+b*im
  end
  nr = n-2nc
  rr = -rand(m.rng,Uniform(m.i_re.first,m.i_re.last),nr)
  return [rc;rr]
end

function initial_estimate(m::AddComplexRoot,problem::ProblemLike,proc_prevs)
  if !is_expanding(problem,proc_prevs) return nothing end
  eq = if m.rootmode==KEEP_PREVIOUS
    # Add a complex root(s) to a previous model, keeping both complex and real roots
    proc = nothing
    for pp in reverse(proc_prevs)
      p = order(pp)
      if iseven(problem.p-p)
        proc = pp; break
      end
    end
    if isnothing(proc) return nothing end
    n = problem.p-order(proc)
    r = [SdeRoots(proc).r;randroots(m,n)]
    σy = std(getfield.(problem.data,:y);mean=0) # std(proc)
    SdeRoots(sortroots(r);σy=σy)
  elseif m.rootmode==MAX_COMPLEX
    # Add to complex roots of previous model, thereby maximizing the number of complex roots
    proc = proc_prevs[end]
    r = copy(SdeRoots(proc).r)
    r = filter(x->!isreal(x),r)
    n_add = problem.p-length(r)
    append!(r,EstimateSDE.randroots(m,n_add))
    SdeRoots(sortroots(r);σy=problem.σy)
  elseif m.rootmode==ALL_RANDOM
    # Purely random complex roots
    r = randroots(m,problem.p)
    y = getfield.(problem.data,:y)
    SdeRoots(sortroots(r))
  else error("Unknown root mode $(m.rootmode)")
  end
  return SdeProcess(sde(problem,eq))
end

"
Initial model returned only when a model of `order(proc)` is requested,
`nothing` otherwise.

$FIELDS
"
struct InitialModel{T<:SdeProcess} <:InitialEstimate
  "process"
  proc::T
end

function initial_estimate(m::InitialModel,problem::Problem,proc_prevs′)
  proc = if is_expanding(problem,proc_prevs′) && order(m.proc)==problem.p
    m.proc else nothing end
  return proc
end

"
Autoregressive models from resampled signals. Return the model with the largest
log likelihood. Does not use `proc_prev`.
Negative real roots are replaced by SDE roots using `Addroot`.

$FIELDS
"
struct ResampleAr<:InitialEstimate
  "Model sampling intervals"
  Ts::Vector{Float64}
  "Interpolation method"
  imethod::Symbol
  "Additional roots"
  addroot::AddComplexRoot
end

function initial_estimate(m::ResampleAr,problem::ProblemLike,proc_prevs::Any)
  if !is_expanding(problem,proc_prevs) return nothing end
  p = problem.p
  t = [x.t for x in problem.data]; y = [x.y for x in problem.data]
  Tmin = minimum(m.Ts)
  _,yr0 = interpolate(t,y,Tmin;method=m.imethod)
  proc_prev = proc_prevs[end]

  # Autoregressive models from resampled data
  procs = map(m.Ts) do T
    yr = T==Tmin ? yr0 : resample(yr0,Tmin/T)
    arhat,P,_ = estimate(ArModelStruct(p),reshape(yr,length(yr),1))
    a = predpar(arhat)[1][1,1,:,1]
    a_sde,n = ar2sdepar(a,T,true)
    r = isempty(a_sde) ? ComplexF64[] : SdeRoots(Sde(a_sde)).r
    rn = randroots(m.addroot,n)
    eq = SdeRoots(sortroots([r;rn]),problem.σy)
    SdeProcess(eq)
  end

  # Select model with maximum likelihood
  L = map(x->logpdf(problem,x),procs)
  eq = procs[argmax(L)].eq

  return SdeProcess(sde(problem,eq))
end

struct RemoveRoots<:InitialEstimate
  addroot::AddRoot
end

function initial_estimate(m::RemoveRoots,problem::ProblemLike,proc_prevs)
  if is_expanding(problem,proc_prevs) return nothing end
  proc_prevs = filter(x->order(x)>problem.p,proc_prevs)
  Lmax = -Inf; eq = nothing;
  for proc in proc_prevs
    r = roots(proc.eq)
    rs = r[end-problem.p+1:end]
    if imag(rs[1])>0  rs[1] = randroots(m.addroot,1)[1] end
    eq_c = SdeRoots(sortroots(rs);σy=problem.σy)
    L = logpdf(problem,SdeProcess(eq_c))
    if L>Lmax  eq = eq_c; Lmax = L  end
  end
  if isnothing(eq) return nothing end
  proc = SdeProcess(sde(problem,eq))
  return proc
end

function default_init_tors(rng::MersenneTwister,Ts::Vector{<:Number},n_random_init=5)
  addcomplex = AddComplexRoot(default_addroot_int,default_addimagroot_int,MAX_COMPLEX,rng)
  addroot = AddRoot(default_addroot_int,rng)
  init_tors = [
    addroot,
    addcomplex,
    ResampleAr(Ts,:linear,addcomplex),
    ResampleAr(Ts,:nearest,addcomplex),
    RemoveRoots(addroot)
  ]
  append!(init_tors,default_random_init_tors(rng,n_random_init))
  return init_tors
end

"
$SIGNATURES

Default initial estimators for `problem`. `T` is the time scale of interest. If omitted,
it is derived from the problem data using `minimal_time`.
"
function default_init_tors(problem::Problem,T::Number=minimal_time(problem))
  Ts = [0.5,1,2]*T
  @warn "Using system randomness for random number generator"
  return default_init_tors(MersenneTwister(),Ts)
end

"
$SIGNATURES

Default random initial estimators
"
function default_random_init_tors(rng::AbstractRNG,N::Int)
  init_tors = map(1:N) do _
    rng_c = MersenneTwister(rand(rng,UInt32,4))
    AddComplexRoot(default_addroot_int,default_addimagroot_int,ALL_RANDOM,rng_c)
  end
  return init_tors
end