text:
using BioAlignments
using MIToS.MSA
using MIToS.PDB
using StatsPlots
using Statistics
using DataFrames
using Plots
using CSV

function fasta2pir(aln)
    seq_pos = 0
    ref_pos = 0
    last_seq_pos = 0
    seq2ref = Dict{Int,Int}()
    for (seq_res, ref_res) in alignment(aln)
        if seq_res != '-'
            seq_pos += 1
        end
        if ref_res != '-'
            ref_pos += 1
        end
        if seq_pos != last_seq_pos
            seq2ref[seq_pos] = ref_pos
            last_seq_pos = seq_pos
        end
    end
    seq2ref
end

function get_exons(query_name)
    df = DataFrame(symbol = Char[], length = Int[])
    lines = readlines("$(query_name)_$(iteration).exons_lengths.txt")
    exons = [parse(Int, ss) for ss in split(lines[2])]
    # MYO1B example:
    # lines[1] e.g. ">P1;ENSG00000128641 ENST00000304164 %()+0123456789:<=?" *
    # "ABCDEFGHIJKLMNOPQRSTUVWXYZ[]^`abcdefghijklmnopqrstuvwxyz{|}~¼½¾¿ÀÁÂÃ" *
    # "ÄÅÆÇÈÉÊËÌÍÎÏÐÑÒÓÔÕÖ×ØÙÚÛÜÝÞ"
    exons_symbol = split(lines[1])[3]
    # exons_symbol e.g. "%()+0123456789:<=?ABCDEFGHIJKLMNOPQRSTUVWXYZ[]^`ab" *
    # "cdefghijklmnopqrstuvwxyz{|}~¼½¾¿ÀÁÂÃÄÅÆÇÈÉÊËÌÍÎÏÐÑÒÓÔÕÖ×ØÙÚÛÜÝÞ"
    for (i, exon) in enumerate(exons_symbol)
        if i == length(exons)
            break
        end
        push!(df, (exon, (exons[i+1] - exons[i])))
    end
    return df
end

function get_infos(exons, msa_total, alignment)
    new_f2p = Dict(v => k for (k, v) in f2p)
    fragment_templates = Int64[]
    fragment_identity_mean = Float64[]
    fragment_identity_max = Float64[]
    fragment_templates_names = []
    alltemplates = []
    first_pos_pir = Int64[]
    last_pos_pir = Int64[]
    index = 1
    loop_end = 0
    i = 1
    loop_start = i
    goon = true
    debut = 1
    while goon && i <= size(exons, 1)
#         println(i)
        identityliste = Float64[]
        templateWellAligned = 0
        templatesIndex = []
#         println(debut)
#         println(debut+exons.length[i]-1)

        fin = debut + exons.length[i] - 1
#         println(fin)

        # if returns a 0 element array, push 0 in everything
        fragment_index = [new_f2p[i] for i in debut:fin if i in keys(new_f2p)]
#                     println("fragment index = $(fragment_index)")
        if isempty(fragment_index)
            push!(first_pos_pir, -1)
            push!(last_pos_pir, -1)
            push!(templatesIndex, "None")
            push!(alltemplates, "None")
            push!(fragment_templates_names, 0)
            push!(fragment_templates, 0)
            push!(fragment_identity_mean, 0)
            push!(fragment_identity_max, 0)
        else
#             println(fragment_index)
            index_debut = fragment_index[1]
            index_end = fragment_index[end]
#             println("index_debut = $(index_debut)")
#             println(" end = $(index_end)")
#         if index_end == size(msa_total, 2)
#             goon= false
#             break
#         end
#         if f2p[debut] == f2p[debut+exons.length[i]-1]
#             if f2p[debut] != size(msa_total, 2)
#                 println("go to the next iteration")
#                 loop_start = i+1
#             end
#         else

            #fragment = msa_total[2:nsequences(msa_total), index_debut:index_end]
            fragment = msa_total[:, index_debut:index_end]
#             println(fragment)
            pir_positions = getcolumnmapping(fragment)
#             println("pir pos")
#             println(pir_positions)
            pir_first = pir_positions[1]
            pir_last = pir_positions[end]
            push!(first_pos_pir, pir_first)
            push!(last_pos_pir, pir_last)
            cov = 100 .* coverage(fragment) * size(fragment, 2) /
                  exons.length[i]
#             println(cov)
            for k = 2:length(cov)
                if cov[k] >= 50
                    templateWellAligned += 1
                    push!(templatesIndex, names(cov)[1][k])
                    push!(alltemplates, names(cov)[1][k])
                end
            end
            push!(fragment_templates_names, templatesIndex)
            push!(fragment_templates, length(templatesIndex))
            for l in templatesIndex
                # @show l
                # @show fragment[[1, l], :]
                # @show percentidentity(fragment[1, :], fragment[l, :])
                push!(
                    identityliste,
                    percentidentity(fragment[1, :], fragment[l, :]),
                )
            end
            if (isempty(identityliste)) # to avoid LOADERROR when trying to do mean() or max() of an empty vector
                push!(fragment_identity_mean, 0)
                push!(fragment_identity_max, 0)
            else
                max_pid = maximum(identityliste)
                # @show max_pid
                push!(fragment_identity_mean, mean(identityliste))
                push!(fragment_identity_max, max_pid)
            end
        end
        loop_end = i
        debut = debut + exons.length[i]
        i = i + 1
    end
    return fragment_templates,
        fragment_identity_mean,
        fragment_identity_max,
        fragment_templates_names,
        alltemplates,
        loop_end,
        loop_start,
        first_pos_pir,
        last_pos_pir
end

# query_name = "./test_258682/ENSG00000058404_ENST00000258682/ENSG00000058404_ENST00000258682"#ARGS[1]
# iteration = 1#ARGS[2]
# query_name = "./FMR1/ENSG00000102081_ENST00000370475/ENSG00000102081_ENST00000370475"#ARGS[1]
# iteration = 4#ARGS[2]
query_name = ARGS[1]
iteration = ARGS[2]

gap = Residue('-')
open("$(query_name)_$(iteration).fa", "r") do f
    global seq_lines = readlines(f)
end
sequence = seq_lines[2]
exons = DataFrame(symbol = Char[], length = Int[])
exons_filtered = DataFrame(symbol = Char[], length = Int[])
exons = get_exons(query_name)
# println(exons)
exons_filtered = filter(row -> row[:length] > 7, exons)
msa_total = read("$(query_name)_$(iteration).pir", PIR, deletefullgaps = false)
setreference!(msa_total, 1)
adjustreference!(msa_total)
costmodel = AffineGapScoreModel(BLOSUM62, gap_open = -10, gap_extend = -1)
aln = pairalign(
    SemiGlobalAlignment(),
    join(msa_total[1, :]),
    seq_lines[2],
    costmodel,
)
f2p = fasta2pir(aln)
#println(sequence)
# println(index_debut.start)
df = DataFrame(
    exonNumber = Int[],
    symbol = Char[],
    nbTemp = Int[],
    meanIdentity = Float64[],
    maxIdentity = Float64[],
    color = String[],
    legend = String[],
    first_pos_pir = Int64[],
    last_pos_pir = Int64[],
    templatesnames = [],
)
fragments_identity_mean = Float64[]
fragments_identity_max = Float64[]
mean_rmsd = Float64[]
fragments_templates = Int64[]
fragments_templates_names = []
alltemplates = []
first_pos_pir = Int64[]
last_pos_pir = Int64[]
fragments_templates,
fragments_identity_mean,
fragments_identity_max,
fragments_templates_names,
alltemplates,
loop_end,
loop_start,
first_pos_pir,
last_pos_pir = get_infos(exons, msa_total, f2p)
n = exons_filtered.symbol
global it = 0
for i = 1:loop_end-loop_start+1
#     println(i)
#      if exons.length[i] > 7
    if (fragments_identity_max[i] < 25)
        push!(
            df,
            (
             loop_start + it,
             Char(exons.symbol[loop_start+it]),
             fragments_templates[i],
             fragments_identity_mean[i],
             fragments_identity_max[i],
             "red",
             "0-24.99%",
             first_pos_pir[i],
             last_pos_pir[i],
             fragments_templates_names[i],
            ),
        )#, mean_rmsd[i]))
    elseif ((fragments_identity_max[i] >= 25) &&
            (fragments_identity_max[i] < 50))
        push!(
            df,
            (
             loop_start + it,
             Char(exons.symbol[loop_start+it]),
             fragments_templates[i],
             fragments_identity_mean[i],
             fragments_identity_max[i],
             "orange",
             "25-49.99%",
             first_pos_pir[i],
             last_pos_pir[i],
             fragments_templates_names[i],
            ),
        )#, mean_rmsd[i]))
    elseif ((fragments_identity_max[i] >= 50) &&
            (fragments_identity_max[i] < 75))
        push!(
            df,
            (
             loop_start + it,
             Char(exons.symbol[loop_start+it]),
             fragments_templates[i],
             fragments_identity_mean[i],
             fragments_identity_max[i],
             "blue",
             "50-74.99%",
             first_pos_pir[i],
             last_pos_pir[i],
             fragments_templates_names[i],
            ),
        )#, mean_rmsd[i]))
    else
        push!(
            df,
            (
             loop_start + it,
             Char(exons.symbol[loop_start+it]),
             fragments_templates[i],
             fragments_identity_mean[i],
             fragments_identity_max[i],
             "green",
             "75-100%",
             first_pos_pir[i],
             last_pos_pir[i],
             fragments_templates_names[i],
            ),
        )#, mean_rmsd[i]))
    end
    global it += 1
#      end
end
total_length = 0
iterative_number = DataFrame(symbol = Char[], number = Char[])
for i = 1:size(df, 1)
    if (exons.length[i] < 5)
        if (df.maxIdentity[i]) > 25
            push!(iterative_number, (exons.symbol[loop_start+i-1], '1'))
        elseif (df.first_pos_pir[i] == -1)
            push!(iterative_number, (exons.symbol[loop_start+i-1], '0'))
        else
            push!(iterative_number, (exons.symbol[loop_start+i-1], '3'))
        end
    elseif (df.maxIdentity[i] > 25)
        push!(iterative_number, (exons.symbol[loop_start+i-1], '1'))
    elseif (df.first_pos_pir[i] == -1)
        push!(iterative_number, (exons.symbol[loop_start+i-1], '0'))
    else
        push!(iterative_number, (exons.symbol[loop_start+i-1], '2'))
    end
end
df2 = DataFrame(name = [], nb = [])
for i in sequencenames(msa_total)
    counted = count(c -> c == i, collect(alltemplates))
    push!(df2, (i, counted))
end
a = groupby(df2, 2, sort = true)
df3 = DataFrame(numberOfExons = [], nbOfTemplates = [])
for i = 2:length(a)
    push!(df3, (a[i].nb[1], length(a[i].nb)))
end

lines = readlines("$(query_name)_annotated.pir")

collected_exons_symbol = collect(lines[2])
first = findfirst(isequal(iterative_number.symbol[1]), collected_exons_symbol)
last = findlast(
    isequal(iterative_number.symbol[end]),
    collected_exons_symbol,
)

global exons_numbers = ""
for i = 1:length(iterative_number.symbol)
    global exons_numbers = exons_numbers *
                           (iterative_number.number[i]^((findlast(
        isequal(iterative_number.symbol[i]),
        lines[2],
    ) + 1) - findfirst(
        isequal(iterative_number.symbol[i]),
        lines[2],
    )))
end
l = collect(lines[4])
global number = 1
for i = first:last
    # if already runned once but still not found, do not rerun
    if (l[i] == '2' && exons_numbers[number] == '2')
        l[i] = '3'
    # if runned once but not in the window search of hhblits, do not change the nuber
    elseif (l[i] == '2' && exons_numbers[number] == '0')
        l[i] = '2'
    else
        l[i] = exons_numbers[number]
    end
    global number = number + 1
end
# println(join(l))
open("$(query_name)_annotated.pir", "w") do f
    write(f, lines[1] * "\n")
    write(f, lines[2] * "\n")
    write(f, lines[3] * "\n")
    write(f, join(l))
end

# creation of every plots and saving it in <name query>.pdf
if !isempty(df)
    p1 = @df df bar(
        :nbTemp,
        color = :color,
        group = :legend,
        bar_width = 0.9,
        title = "Number of templates by exons",
        xlabel = "Exons",
        xlim = (0, size(df, 1) + 1),
        xticks = (1:1:size(df, 1)+2, :symbol),
        ylabel = "Number of templates",
        ylim = (0, maximum(:nbTemp) + 1),
        yticks = 0:5:maximum(:nbTemp),
        legend = :topright,
    )

    p4 = @df df3 bar(
        :numberOfExons,
        :nbOfTemplates,
        title = "number of exons covered by the number of templates",
        xlabel = "Number of exons covered",
        ylabel = "Number of templates",
        xlim = (0, loop_end + 1),
        legend = false,
    )
    p = plot(p1, p4, size = (1280, 720))
    savefig(p, "$(query_name)_$(iteration)_plots.pdf")

end
CSV.write("$(query_name)_$(iteration)_plot1.csv", df)
CSV.write("$(query_name)_$(iteration)_plot_2.csv", df2)

# p = plot(p1,p4,size = (1280, 720))

max_stars_repo_name: elolaine/PhyloSofS

max_issues_repo_name: elolaine/Phylosofs

max_forks_repo_name: PhyloSofS-Team/PhyloSofS


text:
"""
    AbstractSimData

Supertype for simulation data objects. Thes hold grids, settings other objects required
to run the simulation and potentially requireing access from rules.

An `AbstractSimData` object is accessable in [`applyrule`](@ref) as the first parameter.

Multiple grids can be indexed into using their key if you need to read
from arbitrary locations:

```julia
funciton applyrule(data::AbstractSimData, rule::SomeRule{Tuple{A,B}},W}, (a, b), I) where {A,B,W}
    grid_a = data[A]
    grid_b = data[B]
    ...
```

In single-grid simulations `AbstractSimData` objects can be indexed directly as
if they are a `Matrix`.

## Methods

- `currentframe(data)`: get the current frame number, an `Int`
- `currenttime(data)`: the current frame time, which `isa eltype(tspan)`
- `aux(data, args...)`: get the `aux` data `NamedTuple`, or `Nothing`.
    adding a `Symbol` or `Val{:symbol}` argument will get a field of aux.
- `tspan(data)`: get the simulation time span, an `AbstractRange`.
- `timestep(data)`: get the simulaiton time step.
- `boundary(data)` : returns the [`BoundaryCondition`](@ref) - `Remove` or `Wrap`.
- `padval(data)` : returns the value to use as grid border padding.

These are available, but you probably shouldn't use them and their behaviour
is not guaranteed in furture versions. They will mean rule is useful only
in specific contexts.

- `settings(data)`: get the simulaitons [`SimSettings`](@ref) object.
- `extent(data)` : get the simulation [`AbstractExtent`](@ref) object.
- `init(data)` : get the simulation init `AbstractArray`/`NamedTuple`
- `mask(data)` : get the simulation mask `AbstractArray`
- `source(data)` : get the `source` grid that is being read from.
- `dest(data)` : get the `dest` grid that is being written to.
- `radius(data)` : returns the `Int` radius used on the grid,
    which is also the amount of border padding.
"""
abstract type AbstractSimData{S} end

@propagate_inbounds Base.setindex!(d::AbstractSimData, x, I...) = setindex!(first(grids(d)), x, I...)
@propagate_inbounds Base.getindex(d::AbstractSimData, I...) = getindex(first(grids(d)), I...)
Base.keys(d::AbstractSimData) = keys(grids(d))
Base.values(d::AbstractSimData) = values(grids(d))
Base.first(d::AbstractSimData) = first(grids(d))
Base.last(d::AbstractSimData) = last(grids(d))

gridsize(d::AbstractSimData) = gridsize(first(d))
padval(d::AbstractSimData) = padval(extent(d))

extent(d::AbstractSimData) = d.extent
frames(d::AbstractSimData) = d.frames
grids(d::AbstractSimData) = d.grids
init(d::AbstractSimData) = init(extent(d))
mask(d::AbstractSimData) = mask(first(d))
aux(d::AbstractSimData, args...) = aux(extent(d), args...)
auxframe(d::AbstractSimData, key) = auxframe(d)[_unwrap(key)]
auxframe(d::AbstractSimData) = d.auxframe
tspan(d::AbstractSimData) = tspan(extent(d))
timestep(d::AbstractSimData) = step(tspan(d))
currentframe(d::AbstractSimData) = d.currentframe
currenttime(d::AbstractSimData) = tspan(d)[currentframe(d)]

# Getters forwarded to data
Base.getindex(d::AbstractSimData, key::Symbol) = getindex(grids(d), key)

# Get the actual current timestep, e.g. seconds instead of variable periods like Month
currenttimestep(d::AbstractSimData) = currenttime(d) + timestep(d) - currenttime(d)

# Uptate timestamp
function _updatetime(simdata::AbstractSimData, f::Integer)
    @set! simdata.currentframe = f
    @set simdata.auxframe = _calc_auxframe(simdata)
end

# Convert regular index to block index
@inline _indtoblock(x::Int, blocksize::Int) = (x - 1) ÷ blocksize + 1

# Convert block index to regular index
@inline _blocktoind(x, blocksize) = (x - 1) * blocksize + 1

"""
    SimData <: AbstractSimData

    SimData(extent::AbstractExtent, ruleset::AbstractRuleset)

Simulation dataset to hold all intermediate arrays, timesteps
and frame numbers for the current frame of the simulation.

Additional methods not found in `AbstractSimData`:

- `rules(d::SimData)` : get the simulation rules.
- `ruleset(d::SimData)` : get the simulation [`AbstractRuleset`](@ref).
"""
struct SimData{S<:Tuple,G<:NamedTuple,E,RS,F,CF,AF} <: AbstractSimData{S}
    grids::G
    extent::E
    ruleset::RS
    frames::F
    currentframe::CF
    auxframe::AF
end
function SimData{S}(
    grids::G, extent::E, ruleset::RS, frames::F, currentframe::CF, auxframe::AF
) where {S,G,E,RS,F,CF,AF}
    SimData{S,G,E,RS,F,CF,AF}(grids, extent, ruleset, frames, currentframe, auxframe)
end
SimData(o, ruleset::AbstractRuleset) = SimData(o, extent(o), ruleset)
function SimData(o, extent::AbstractExtent, ruleset::AbstractRuleset)
    frames_ = if hasdelay(rules(ruleset))
        isstored(o) || _notstorederror()
        frames(o)
    else
        nothing
    end
    SimData(extent, ruleset, frames_)
end
# Convert grids in extent to NamedTuple
SimData(extent::AbstractExtent, ruleset::AbstractRuleset, frames=nothing) =
    SimData(_asnamedtuple(extent), ruleset)
function SimData(extent::AbstractExtent{<:NamedTuple{Keys}}, ruleset::AbstractRuleset, frames=nothing) where Keys
    # Calculate the neighborhood radus (and grid padding) for each grid

    S = Val{Tuple{gridsize(extent)...}}()
    radii = map(k-> Val{get(radius(ruleset), k, 0)}(), Keys)
    radii = NamedTuple{Keys}(radii)
    grids = _buildgrids(extent, ruleset, S, radii)
    # Construct the SimData for each grid
    SimData(grids, extent, ruleset, frames)
end
function SimData(
    grids::G, extent::AbstractExtent, ruleset::AbstractRuleset, frames
) where {G<:Union{NamedTuple{<:Any,<:Tuple{GridData,Vararg}},GridData}}
    currentframe = 1; auxframe = nothing
    S = Tuple{gridsize(extent)...}
    # SimData is isbits-only
    s_extent = StaticExtent(extent)
    s_ruleset = StaticRuleset(ruleset)
    SimData{S}(grids, s_extent, s_ruleset, frames, currentframe, auxframe)
end

_buildgrids(extent, ruleset, S, radii::NamedTuple) =
    map((r, in, pv) -> _buildgrids(extent, ruleset, S, r, in, pv), radii, init(extent), padval(extent))
function _buildgrids(extent, ruleset, ::Val{S}, ::Val{R}, init, padval) where {S,R}
    ReadableGridData{S,R}(
        init, mask(extent), proc(ruleset), opt(ruleset), boundary(ruleset), padval
    )
end

ConstructionBase.constructorof(::Type{<:SimData{S}}) where S = SimData{S}

ruleset(d::SimData) = d.ruleset
rules(d::SimData) = rules(ruleset(d))
boundary(d::SimData) = boundary(ruleset(d))
proc(d::SimData) = proc(ruleset(d))
opt(d::SimData) = opt(ruleset(d))
settings(d::SimData) = settings(ruleset(d))

# When no simdata is passed in, create new AbstractSimData
function _initdata!(::Nothing, output, extent::AbstractExtent, ruleset::AbstractRuleset)
    SimData(output, extent, ruleset)
end
# Initialise a AbstractSimData object with a new `Extent` and `Ruleset`.
function _initdata!(
    simdata::SimData, output, extent::AbstractExtent, ruleset::AbstractRuleset
)
    # TODO: make sure this works with delays and new outputs?
    map(copy!, values(simdata), values(init(extent)))
    @set! simdata.extent = StaticExtent(extent)
    @set! simdata.ruleset = StaticRuleset(ruleset)
    if hasdelay(rules(ruleset))
        isstored(o) || _not_stored_delay_error()
        @set! simdata.frames = frames(o)
    end
    simdata
end

"""
    SimData <: AbstractSimData

    RuleData(extent::AbstractExtent, settings::SimSettings)

`AbstractSimData` object that is passed to rules. Basically
a trimmed-down version of [`SimData`](@ref).
"""
struct RuleData{S<:Tuple,G<:NamedTuple,E,Se,F,CF,AF} <: AbstractSimData{S}
    grids::G
    extent::E
    settings::Se
    frames::F
    currentframe::CF
    auxframe::AF
end
function RuleData{S}(
    grids::G, extent::E, settings::Se, frames::F, currentframe::CF, auxframe::AF
) where {S,G,E,Se,F,CF,AF}
    RuleData{S,G,E,Se,F,CF,AF}(grids, extent, settings, frames, currentframe, auxframe)
end
function RuleData(d::AbstractSimData{S}) where S
    RuleData{S}(grids(d), extent(d), settings(d), frames(d), currentframe(d), auxframe(d))
end

ConstructionBase.constructorof(::Type{<:RuleData{S}}) where {S} = RuleData{S}

settings(d::RuleData) = d.settings
boundary(d::RuleData) = boundary(settings(d))
proc(d::RuleData) = proc(settings(d))
opt(d::RuleData) = opt(settings(d))

ext: jl


text:
using AprilTags
using ImageCore
using FileIO
using ImageMagick
using ImageDraw
using ColorTypes
using FixedPointNumbers
using Test

@testset "AprilTags" begin
    # To test without extra dependancies
    # image = rand(UInt8, 480,640)
    # image[50:79,50:79]     = kron(reinterpret(UInt8,getAprilTagImage(0)), ones(UInt8, 3,3))
    # image[150:179,150:179] = kron(reinterpret(UInt8,getAprilTagImage(1)), ones(UInt8, 3,3))
    # image[250:279,250:279] = kron(reinterpret(UInt8,getAprilTagImage(2)), ones(UInt8, 3,3))
    # image = Gray.(reinterpret(N0f8, image))
    #
    #
    # imageCol = RGB.(image)
    # refpoints = [[63.9, 63.9],
    #             [163.9, 163.9],
    #             [263.9, 263.9]]

    image = load(dirname(Base.source_path()) *"/../data/tagtest.jpg")
    imageCol = load(dirname(Base.source_path()) *"/../data/colortag.jpg")
    refpoints = [[404.5, 176.1],
                 [134.0, 216.1],
                 [412.0, 130.1]]

    @testset "Low-level API" begin
        # test wrappers
        #create april tag detector
        td = apriltag_detector_create()

        #create tag family
        tf = tag36h11_create()

        #add family to detector
        apriltag_detector_add_family(td, tf)

        #create image8 object for april tags
        image8 = convert(AprilTags.image_u8_t, image)

        # run detector on image
        detections = apriltag_detector_detect(td, image8)

        # copy detections
        tags = getTagDetections(detections)

        #extract tag centres
        cpoints = map(tag->[tag.c[2],tag.c[1]],tags)
        @test cpoints ≈ refpoints atol=0.5

        # test convertions
        image8_from_u8 = convert(AprilTags.image_u8_t, reinterpret(UInt8, image)[:,:])
        @test image8_from_u8.width == image8.width
        @test image8_from_u8.height == image8.height
        @test image8_from_u8.stride == image8.stride
        #TODO: maybe add test for content of pointer

        apriltag_detections_destroy(detections)
        # Cleanup: free the detector and tag family when done.
        apriltag_detector_destroy(td)
        tag36h11_destroy(tf)
    end

    @testset "High-level API" begin
        # test the high level convenience functions
        detector = AprilTagDetector()
        tags2 = detector(image)

        @test length(detector(gray.(image))) == 3
        @test length(detector(reinterpret(UInt8,image)[:,:])) == 3

        tagsth = AprilTags.threadcalldetect(detector, image)
        @test length(tagsth) == 3

        #test on random image, should detect zero tags
        @test length(detector(rand(Gray{N0f8},100,100))) == 0

        #getters -- compare with default
        @test detector.nThreads == 1
        @test detector.quad_decimate == 2.0 # TODO: default changed to 2 - consider changing back
        @test detector.quad_sigma == 0.0
        @test detector.refine_edges == 1
        @test detector.decode_sharpening == 0.25

        #setters -- set new values
        detector.nThreads = 4
        detector.quad_decimate = 2.0
        detector.quad_sigma = 0.1
        detector.refine_edges = 0
        detector.decode_sharpening = 0.2
        # detector.refine_decode = 1
        # detector.refine_pose = 1

        #getters -- compare with set values just now
        @test detector.nThreads == 4
        @test detector.quad_decimate == 2.0
        @test detector.quad_sigma == 0.1f0
        @test detector.refine_edges == 0
        @test detector.decode_sharpening == 0.2
        # @test detector.refine_decode == 1
        # @test detector.refine_pose == 1

        #test @show overload
        # @test sprint((t,s)->show(t,"text/plain",s), detector) == "AprilTagDetector\nnThreads: 4\nquad_decimate: 1.0\nquad_sigma: 0.0\nrefine_edges: 1\nrefine_decode: 1\nrefine_pose: 1\n"
        @test sprint((t,s)->show(t,"text/plain",s), detector) == "AprilTagDetector\nnThreads: 4\nquad_decimate: 2.0\nquad_sigma: 0.1\nrefine_edges: 0\ndecode_sharpening: 0.2\n"

        cpoints = map(tag->[tag.c[2],tag.c[1]],tags2)
        freeDetector!(detector)
        @test cpoints ≈ refpoints atol=0.5

        pose = homography_to_pose(tags2[1].H, -520., 520., 320., 240.)
        # TODO create better ref pose
        refpose = [0.65  0.12 -0.75  -5.43;
                   0.39 -0.90  0.19  -6.20;
                  -0.65 -0.41 -0.63 -19.62;
                   0.0   0.0   0.0    1.0]
        @test pose[1:3,1:3] ≈ refpose[1:3,1:3] atol = 0.05
        @test pose[1:3,4] ≈ refpose[1:3,4] atol = 0.1

        #test drawing functions
        fx = 524.040
        fy = 524.040
        cx = 319.254
        cy = 251.227
        K = [fx 0  cx;
              0 fy cy]
        imCol = RGB.(image)
        foreach(tag->drawTagBox!(imCol, tag), tags2)
        # test one pixel to be correctly drawn
        t1xy = round.(Int,tags2[1].p[1])
        @test imCol[t1xy[2],t1xy[1]] == RGB{N0f8}(0.0, 1.0, 0.0)
        #thicker lines also
        foreach(tag->drawTagBox!(imCol,tag, width = 2, drawReticle = false), tags2)
        foreach(tag->drawTagBox!(imCol,tag, width = 3, drawReticle = true), tags2)
        #TODO: verify that drawing is correct

        # test for bounds error in drawing functions with thicker lines
        foreach(tag->drawTagBox!(imCol,tag, width = 1000, drawReticle = true), tags2)

        foreach(tag->drawTagAxes!(imCol,tag, K), tags2)
        #TODO: verify that drawing tag axis is correct

        # test constructors for other families
        detector2 = AprilTagDetector(AprilTags.tag25h9)
        freeDetector!(detector2)

        detector2 = AprilTagDetector(AprilTags.tag16h5)
        freeDetector!(detector2)

        detector2 = AprilTagDetector(AprilTags.tag36h11)
        freeDetector!(detector2)

        # detector2 = AprilTagDetector(AprilTags.tag36h10)
        # freeDetector!(detector2)

        reftag36h11_0 = Gray{N0f8}[ 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0;
                                    1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0;
                                    1.0 0.0 1.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0;
                                    1.0 0.0 0.0 1.0 1.0 1.0 0.0 1.0 0.0 1.0;
                                    1.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0 0.0 1.0;
                                    1.0 0.0 1.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0;
                                    1.0 0.0 0.0 1.0 0.0 1.0 1.0 0.0 0.0 1.0;
                                    1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0;
                                    1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0;
                                    1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0]
        @test reftag36h11_0 == getAprilTagImage(0)

        # reftag36h10_0 = Gray{N0f8}[ 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0;
        #                             1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0;
        #                             1.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 0.0 1.0;
        #                             1.0 0.0 0.0 0.0 1.0 0.0 1.0 0.0 0.0 1.0;
        #                             1.0 0.0 1.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0;
        #                             1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 0.0 1.0;
        #                             1.0 0.0 0.0 1.0 1.0 0.0 1.0 0.0 0.0 1.0;
        #                             1.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 0.0 1.0;
        #                             1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0;
        #                             1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0]
        # @test reftag36h10_0 == getAprilTagImage(0, AprilTags.tag36h10)

        reftag16h5_0 = Gray{N0f8}[  1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0;
                                    1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0;
                                    1.0 0.0 0.0 0.0 1.0 0.0 0.0 1.0;
                                    1.0 0.0 0.0 0.0 1.0 1.0 0.0 1.0;
                                    1.0 0.0 0.0 0.0 0.0 1.0 0.0 1.0;
                                    1.0 0.0 1.0 0.0 1.0 1.0 0.0 1.0;
                                    1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0;
                                    1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0]
        @test reftag16h5_0 == getAprilTagImage(0, AprilTags.tag16h5)

        reftag25h9_0 = Gray{N0f8}[  1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0;
                                    1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0;
                                    1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0;
                                    1.0 0.0 0.0 1.0 0.0 1.0 1.0 0.0 1.0;
                                    1.0 0.0 1.0 0.0 0.0 1.0 0.0 0.0 1.0;
                                    1.0 0.0 1.0 1.0 1.0 1.0 1.0 0.0 1.0;
                                    1.0 0.0 1.0 0.0 0.0 0.0 1.0 0.0 1.0;
                                    1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0;
                                    1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0]
        @test reftag25h9_0 == getAprilTagImage(0, AprilTags.tag25h9)

        # TODO test, just placeholder for now
        detector = AprilTagDetector()
        detector.quad_decimate = 1.0 #NOTE see line 84
        fx = 524.040
        fy = 524.040
        cx = 319.254
        cy = 251.227
        taglength = 0.172
        (tags, poses) = detectAndPose(detector, image, fx, fy, cx, cy, taglength)
        # TODO test here, this is just result used as test.
        @test all(isapprox.(poses[1], [ 0.630532  -0.391638  0.670111  -0.461887;
                                        0.118082   0.901717  0.415889   0.494238;
                                       -0.767128  -0.183103  0.614807   1.69053],
                                      atol = 0.01))
        freeDetector!(detector)

    end

    include("homography.jl")

    @testset "Errors" begin
        #testing freed detectors errors
        detector = AprilTagDetector()
        freeDetector!(detector)
        @test_throws ErrorException tags = detector(image)
        @test_throws ErrorException tags = AprilTags.threadcalldetect(detector, image)
        @test_throws ErrorException AprilTags.setnThreads(detector, 4)
        @test_throws ErrorException AprilTags.setquad_decimate(detector, 1.0)
        @test_throws ErrorException AprilTags.setquad_sigma(detector,0.0)
        @test_throws ErrorException AprilTags.setrefine_edges(detector,1)
        @test_throws ErrorException AprilTags.setdecode_sharpening(detector,0.2)
        # @test_throws ErrorException AprilTags.setrefine_decode(detector,1)
        # @test_throws ErrorException AprilTags.setrefine_pose(detector,1)

        @test_throws ErrorException AprilTags.getnThreads(detector)
        @test_throws ErrorException AprilTags.getquad_decimate(detector)
        @test_throws ErrorException AprilTags.getquad_sigma(detector)
        @test_throws ErrorException AprilTags.getrefine_edges(detector)
        @test_throws ErrorException AprilTags.getdecode_sharpening(detector)
        # @test_throws ErrorException AprilTags.getrefine_decode(detector)
        # @test_throws ErrorException AprilTags.getrefine_pose(detector)

        @test freeDetector!(detector) == nothing
        #testing NULL tag families errors
        detector = AprilTagDetector()
        detector.tf = C_NULL
        @test_throws ErrorException tags = detector(image)
        @test_throws ErrorException tags = AprilTags.threadcalldetect(detector, image)
        @test freeDetector!(detector) == nothing
    end

    @testset "Color Image Conversion" begin
        detector = AprilTagDetector()
        detector.quad_decimate = 1.0 #NOTE see line 84
        tags = detector(imageCol)
        @test length(tags) == 1
        freeDetector!(detector)
    end
end


text:
# Variables
Z3Real(name::String) = z3[:Real](name)
Z3Int(name::String) = z3[:Int](name)
Z3Bool(name::String) = z3[:Bool](name)

function variable(model::Z3Model,vartype::Symbol=:Cont,lb::Real=-Inf, ub::Real=Inf;name=nothing)
  model.numvars += 1
  numvars = model.numvars
  push!(model.varvalue,NaN)
  push!(model.varlower,lb)
  push!(model.varupper,ub)
  push!(model.vartype,vartype)
  if vartype == :Cont
    name = name == nothing ? "x$numvars" : name
    var = Z3Real(name)
  elseif vartype == :Int
    name = name == nothing ? "z$numvars" : name
    var = Z3Int(name)
  elseif vartype == :Bool || vartype == :Bin
    name = name == nothing ? "y$numvars" : name
      var = Z3Bool(name)
  else
      error("Invalid variable type. Valid types ar :Cont, :Int, :Bool")
  end
    push!(model.varname,name)
    push!(model.varptr,var)
    if vartype != :Bool
        lb != -Inf && untrackedconstraint(model, var >= lb)
        ub != Inf && untrackedconstraint(model, var <= ub)
    end
    return var
end

add_cvar!(model::Z3Model,lb::Real=-Inf,ub::Real=Inf) = variable(model,:Cont,lb,ub)
add_ivar!(model::Z3Model,lb::Real=-Inf,ub::Real=Inf) = variable(model,:Int,lb,ub)
add_bvar!(model::Z3Model) = variable(model,:Bool,0,1)


text:
__precompile__()

module Poltergeist

using Base, ApproxFun, ApproxFunBase, BandedMatrices, Compat, DualNumbers, DomainSets, StaticArrays, LightGraphs#, PyPlot#, FastTransforms
using LinearAlgebra

import Base: values,getindex,setindex!,*,+,-,==,<,<=,>,|,
  >=,/,^,\,∪,∘,transpose, size, length,  eltype, inv, mod, convert#issymmetric,
@compat import LinearAlgebra: eigvals, eigvecs
import ApproxFun: domainspace, rangespace, domain, RaggedMatrix,
  resizedata!, colstop, Infinity,
  fromcanonicalD, tocanonicalD, fromcanonical, tocanonical, space, eigs, cfstype, prectype
import ApproxFunBase: israggedbelow, CachedOperator
import DomainSets: AbstractInterval, UnionDomain

export (..), Interval, PeriodicSegment, Fun, eigs

include("general.jl")
include("maps/maps.jl")
include("HGraphs.jl")
include("hofbauerextension.jl")
include("InducedMap.jl")
include("BackSum.jl")
include("spectral/spectral.jl")
include("timeseries.jl")
include("poetry.jl")

end # module


text:
module RadialBasisFunctionModels #src

using Base: NamedTuple, promote_eltype, inner_mapslices!
export RBFModel, RBFInterpolationModel #src
export Multiquadric, InverseMultiquadric, Gaussian, Cubic, ThinPlateSpline #src

export RBFInterpolator
export RBFMachine, RBFMachineWithKernel, fit!, add_data!

export auto_grad, auto_jac, grad, jac, eval_and_auto_grad
export eval_and_auto_jac, eval_and_grad, eval_and_jac

# Dependencies of this module:
using StaticPolynomials
using ThreadSafeDicts
using Memoization: @memoize
using StaticArrays
using LinearAlgebra: norm
using Lazy: @forward
using Parameters: @with_kw

# We define evaluation functions, so that the right type is returned by `polys`. #src
for V in [:SizedVector, :MVector]
    @eval Base.@propagate_inbounds StaticPolynomials.evaluate( F :: PolynomialSystem, x :: $V ) = $V( StaticPolynomials._evaluate( F, x ) ) #src
    @eval Base.@propagate_inbounds StaticPolynomials.evaluate( F :: PolynomialSystem, x :: $V , p) = $V( StaticPolynomials._evaluate( F, x , p) ) #src
end

import Zygote as Zyg
using Zygote: Buffer

# TODO also set Flux.trainable to make inner parameters trainable #src

# # Radial Basis Function Models

# The module `RadialBasisFunctionModels` provides utilities to work with radial
# basis function [RBF] models.
# Given ``N`` data sites ``X = \{ x^1, …, x^N \} ⊂ ℝ^n`` and values
# ``Y = \{ y^1, …, y^N \} ⊂ ℝ``, an interpolating RBF model ``r\colon ℝ^n → ℝ``
# has the form
# ```math
# r(x) = \sum_{i=1}^N w_i φ( \| x - x^i \|_2 ) + p(x),
# ```
# where `p` is a multivariate polynomial.
# The radial function ``φ\colon [0, ∞) \to ℝ`` defines the RBF and we can solve for
# the coefficients ``w`` by solving the interpolation system
# ```math
# \begin{equation}
# r( x^i ) \stackrel{!}= y^i \quad \text{for all }i=1,…,N
# \label{eqn:coeff_basic}
# \end{equation}
# ```

# !!! note
#     See the section about **[Getting the Coefficients](@ref)** for how we actually solve the equation system.

# ## Radial Basis Function Sum.

# The function ``k(•) = φ(\|•\|_2)`` is radially symmetric around the origin.
# ``k`` is called the kernel of an RBF.
#
# We define an abstract super type for radial functions:
abstract type RadialFunction <: Function end

# Each Type that inherits from `RadialFunction` should implement
# an evaluation method.
# It takes the radius/distance ``ρ = ρ(x) = \| x - x^i \|`` from
# ``x`` to a specific center ``x^i``.
(φ :: RadialFunction )( ρ :: Real ) :: Real = Nothing;
# We also need the so called order of conditional positive definiteness:
cpd_order( φ :: RadialFunction) :: Int = nothing;
# The derivative can also be specified. It defaults to
df( φ :: RadialFunction, ρ ) = Zyg.gradient( φ, ρ )[1]

# The file `radial_funcs.jl` contains various radial function implementations.
include("radial_funcs.jl")

# From an `RadialFunction` and a vector we can define a shifted kernel function.
const NumberOrVector = Union{<:Real, AbstractVector{<:Real}}

struct ShiftedKernel{RT <: RadialFunction, CT <: AbstractVector{<:Real}} <: Function
    φ :: RT
    c :: CT
end

norm2( vec ) = norm(vec, 2)

"Evaluate kernel `k` at `x - k.c`."
function (k::ShiftedKernel)( x :: AbstractVector{<:Real} )
    return k.φ( norm2( x .- k.c ) )
end

# A vector of ``N`` kernels is a mapping ``ℝ^n → ℝ^N, \ x ↦ [ k₁(x), …, k_N(x)] ``.

_eval_vec_of_kernels( K, x ) = [k(x) for k ∈ K]

"Evaluate ``x ↦ [ k₁(x), …, k_{N_c}(x)]`` at `x`."
( K::AbstractVector{<:ShiftedKernel})( x ) = _eval_vec_of_kernels( K, x )

# Suppose, we have calculated the distances ``\|x - x^i\|`` beforehand.
# We can save redundant effort by passing them to the radial functions of the kernels.

"Evaluate `k.φ` for distance `ρ` where `ρ` should equal `x - k.c` for the argument `x`."
eval_at_dist( k :: ShiftedKernel , ρ :: Real ) = k.φ(ρ)

"Evaluate ``x ↦ [ k₁(x), …, k_{N_c}(x)]``, provided the distances ``[ ρ_1(x), …, ρ_{N_c}(x) ]``."
function eval_at_dist( K::AbstractVector{<:ShiftedKernel}, dists :: AbstractVector{<:Real})
    [ eval_at_dist(k,ρ) for (k,ρ) ∈ zip(K,dists) ]
end

# Provided we have solved the interpolation system, the weights for the radial basis function
# part of ``r`` are ``w``, where ``w`` is a vector of length ``N_c`` or a matrix in ``ℝ^{N_c \times k}``
# where k is the number of outputs.
# We treat the general case ``k\ge 1`` and always assume ``w`` to be a matrix.

struct RBFSum{
    KT <: AbstractVector{<:ShiftedKernel},
    WT <: AbstractMatrix{<:Real}
}
    kernels :: KT
    weights :: WT # can be a normal matrix or a SMatrix

    num_outputs :: Int
end

# Make it display nicely:
function Base.show( io :: IO, rbf :: RBFSum{KT,WT} ) where {KT, WT}
    compact = get(io, :compact, false)
    if compact
        print(io, "RBFSum{$(KT), $(WT)}")
    else
        n_out, n_kernels = size(rbf.weights)
        print(io, "RBFSum\n")
        print(io, "* with $(n_kernels) kernels in an array of type $(KT)\n")
        print(io, "* and a $(n_kernels)×$(n_out) weight matrix of type $(WT).")
    end
end

# We can easily evaluate the `ℓ`-th output of the `RBFPart`:
@doc "Evaluate outut `ℓ` of RBF sum `rbf::RBFSum`"
function (rbf :: RBFSum)(x :: AbstractVector, ℓ :: Int)
    return (rbf.weights[ℓ,:]'rbf.kernels(x))[1]
end

# The overall output is a vector, and we also get it via matrix multiplication.
_eval_rbfsum(rbf::RBFSum, x ) = rbf.weights*rbf.kernels(x)
"Evaluate `rbf::RBFSum` at `x`."
(rbf :: RBFSum)( x :: AbstractVector ) = _eval_rbfsum(rbf, x)

# We want to return the right type and use `_type_guard`:
_type_guard( x , :: Type{<:Vector}, :: Int ) = convert( Vector, x)
for V in [:SVector, :MVector, :SizedVector ]
    @eval _type_guard( x, ::Type{ <: $V }, n_out :: Int ) = convert($V{ n_out }, x)
end

(rbf :: RBFSum)( x :: Vector ) = _type_guard( _eval_rbfsum(rbf, x), Vector, rbf.num_outputs )
function (rbf :: RBFSum)( x :: T ) where T<:Union{SVector,MVector,SizedVector}
    return _type_guard( _eval_rbfsum(rbf, x), T, rbf.num_outputs )
end

# As before, we allow to pass precalculated distance vectors:
function eval_at_dist( rbf::RBFSum, dists :: AbstractVector{<:Real}, ℓ :: Int )
   rbf.weights[ℓ,:]'eval_at_dist( rbf.kernels, dists )
end

function eval_at_dist( rbf :: RBFSum, dists :: AbstractVector{<:Real})
   vec(rbf.weights*eval_at_dist(rbf.kernels, dists ))
end

# For the PolynomialTail do something similar and
# use a `StaticPolynomials.PolynomialSystem` with a weight matrix.
include("empty_poly_sys.jl")

# This allows for the `PolySum`. `polys` evaluates the polynomial basis and
# `weights` are determined during training/fitting.
struct PolySum{
        PS <: Union{EmptyPolySystem, PolynomialSystem},
        WT <: AbstractMatrix
    }
    polys :: PS
    weights :: WT       # n_out × n_polys matrix
    num_outputs :: Int

    function PolySum( polys :: PS, weights :: WT) where{PS, WT}
        n_out, n_polys = size(weights)
        @assert npolynomials(polys) == n_polys "Number of polynomials does not match."
        new{PS,WT}(polys, weights, n_out)
    end
end

eval_psum( p :: PolySum, x ) = p.weights * p.polys(x)
(p :: PolySum)(x :: AbstractVector ) = eval_psum( p, x )
(p :: PolySum)(x :: Vector ) = _type_guard(eval_psum(p,x), Vector, p.num_outputs )
(p :: PolySum)(x :: T) where T<:Union{SVector,MVector,SizedVector} = _type_guard( eval_psum(p,x), T, p.num_outputs)

(p :: PolySum)(x,ℓ::Int) = (p.weights[ℓ,:]'p.polys(x))[end]

# We now have all ingredients to define the model type.

"""
    RBFModel{V, RS, PS, M}

* `V` is `true` by default. It can be set to `false` only if the number
  of outputs is 1. Then scalars are returned.

"""
struct RBFModel{
        V,
        RS <: RBFSum,
        PS <: PolySum,
        M }
    rbf :: RS
    psum :: PS

    ## Information fields
    num_vars :: Int
    num_outputs :: Int
    num_centers :: Int

    meta :: M

    function RBFModel( rbf :: RS, psum :: PS, num_vars, num_outputs, num_centers, meta :: M = nothing; vec_output :: Bool = true) where{RS,PS,M}
        return new{vec_output, RS,PS, M}(rbf, psum, num_vars, num_outputs, num_centers, meta)
    end
end

# We want a model to be displayed in a sensible way:
function Base.show( io :: IO, mod :: RBFModel{V,RS,PS,M} ) where {V,RS,PS,M}
    compact = get(io, :compact, false)
    if compact
        print(io, "$(mod.num_vars)D$(mod.num_outputs)D-RBFModel{$(V)}")
    else
        print(io, "RBFModel{$(V),$(RS),$(PS)}\n")
        if V
            print(io, "\twith vector output ")
        else
            print(io, " with scalar output ")
        end
        print(io, "and $(mod.num_centers) centers, ")
        print(io, "mapping from ℝ^$(mod.num_vars) to ℝ^$(mod.num_outputs).")
    end
end

# Evaluation is easy. We accept an additional `::Nothing` argument that does nothing
# for now, but saves some typing later.
function vec_eval(mod :: RBFModel, x :: AbstractVector{<:Real}, :: Nothing)
    return mod.rbf(x) .+ mod.psum( x )
end

function scalar_eval(mod :: RBFModel, x :: AbstractVector{<:Real}, :: Nothing )
    return (mod.rbf(x) + mod.psum( x ))[1]
end

## @doc "Evaluate model `mod :: RBFModel` at vector `x`."
( mod :: RBFModel{true, RS, PS, M} where {RS,PS,M} )(x :: AbstractVector{<:Real}, ℓ :: Nothing = nothing ) = vec_eval(mod,x,ℓ)
( mod :: RBFModel{false, RS, PS, M} where {RS,PS,M} )(x :: AbstractVector{<:Real}, ℓ :: Nothing = nothing ) = scalar_eval(mod,x,ℓ)

"Evaluate scalar output `ℓ` of model `mod` at vector `x`."
function (mod :: RBFModel)( x :: AbstractVector{<:Real}, ℓ :: Int)
    return mod.rbf(x, ℓ) + mod.psum( x, ℓ )
end

## scalar input
const NothInt = Union{Nothing,Int}

function (mod :: RBFModel)(x :: Real, ℓ :: NothInt = nothing )
    @assert mod.num_vars == 1 "The model has more than 1 inputs. Provide a vector `x`, not a number."
    mod( [x,], ℓ)
end

include("derivatives.jl")
include("constructors.jl")
#-
include("mlj_interface.jl")
# [^wild_diss]: “Derivative-Free Optimization Algorithms For Computationally Expensive Functions”, Wild, 2009.
# [^wendland]: “Scattered Data Approximation”, Wendland
# [^adv_eco]: “Advanced Econometrics“, Takeshi Amemiya
end #src


text:
#compare sequence
seq = Sequence("AAAATTTTCCCCGGGGAAAATTTTCCCCGGGG")
same_seq = Sequence("AAAATTTTCCCCGGGGAAAATTTTCCCCGGGG")
reverse_seq = Sequence("GGGGCCCCTTTTAAAAGGGGCCCCTTTTAAAA")
compelment_seq = Sequence("TTTTAAAAGGGGCCCCTTTTAAAAGGGGCCCC")
reverse_compelment_seq = Sequence("CCCCGGGGAAAATTTTCCCCGGGGAAAATTTT")
qual = Quality("FFFF####AAAAFFFFFFFFFFFFFFFFFFFF")
seq_name = "this_is_a_sequence_name"
fq = FastqRead(seq_name, seq, "+", qual)
reverse_compelment_fq = FastqRead(seq_name, reverse_compelment_seq, "+", reverse(qual))
fa = FastaRead(seq_name, seq)
reverse_compelment_fa = FastaRead(seq_name, reverse_compelment_seq)

@test seq == same_seq
@test seq == -reverse_seq
@test seq == !compelment_seq
@test seq == ~reverse_compelment_seq
@test seq == -(-seq)
@test seq == !(!seq)
@test seq == ~(~seq)
@test complement(dna("ATCG"))== dna("TAGC")
@test reverse_compelment_seq == !(-seq)
@test reverse_compelment_seq == -(!seq)
@test reverse_compelment_fq.sequence == (~fq).sequence
@test reverse_compelment_fq.quality == (~fq).quality
@test reverse_compelment_fa.sequence == (~fa).sequence
@test dna("AAAATTTTCCCCGGGGAAAATTTTCCCCGGGG") == Sequence("AAAATTTTCCCCGGGGAAAATTTTCCCCGGGG", DNA_SEQ)
@test rna("AAAATTTTCCCCGGGGAAAATTTTCCCCGGGG") == Sequence("AAAATTTTCCCCGGGGAAAATTTTCCCCGGGG", RNA_SEQ)
@test aa("AAAATTTTCCCCGGGGAAAATTTTCCCCGGGG") == Sequence("AAAATTTTCCCCGGGGAAAATTTTCCCCGGGG", AA_SEQ)


text:
Base.@kwdef struct FittingParameter{S, T} <: AbstractFittingParameter{T}
    xviews::VectorRange{T}                              # references to a partial range of the global parameter vector
    pviews::Vector{RowRange{T}}                         # references to the internal array of model parameters (for the entire period of fixing)
    yviews::CartesianView{T}                            # references to the internal array of model parameters (at the beginning of the period of fixing)
    gviews::VectorRange{T}                              # references to a partial range of the global gradient vector
    bviews::VectorView{T}                               # references to auxiliary buffer 1
    bviews2::VectorView{T}                              # references to auxiliary buffer 1
    V::BitMatrix                                        # array of parameter fixing periods
    vviews::Vector{ColumnSliceBool}                     # references to columns of the array of the period of validity of the parameter
end

Base.@kwdef struct FittingCache{T<:AbstractFloat}
    xbuf::Vector{T}
    ybuf::Vector{T}
    gbuf::Vector{T}
    tbuf::Vector{T}
    tbuf2::Vector{T}
    cbuf::Vector{T}
    cbuf2::Vector{T}
end

Base.@kwdef struct FittingSet{T<:AbstractFloat, PS<:AbstractParametersScaling{T}, TP<:Tuple{Vararg{<:AbstractFittingParameter{T}}}}
    params::TP
    scale::PS
    cache::FittingCache{T}
end

function build_fitting_views(df::AbstractDataFrame, prob::NonlinearProblem, name::Symbol)
    P = getfield(prob.params, name)
    pviews, yviews = @with df begin
        view.(Ref(P), :Istart, UnitRange.(:Jstart, :Jstop)),
        view(P, CartesianIndex.(:Istart, :Jstart))
    end
    return pviews, yviews
end

function build_validity_matrix(df::AbstractDataFrame, N, M)
    V = falses(N, M)
    for (v, dfr) ∈ zip(eachrow.((V, df))...)
        v[UnitRange(dfr.Jstart, dfr.Jstop)] .= true
    end
    return V, collect(eachcol(V))
end

function FittingCache{T}(Nt, Nc, Nx) where {T}
    kwargs = (
        xbuf = Array{T}(undef, Nx),
        ybuf = Array{T}(undef, Nx),
        gbuf = Array{T}(undef, Nx),
        tbuf = Array{T}(undef, Nt),
        tbuf2 = Array{T}(undef, Nt),
        cbuf = Array{T}(undef, Nc),
        cbuf2 = Array{T}(undef, Nc),
    )
    FittingCache{T}(kwargs...)
end

function FittingParameter{S, T}(df::AbstractDataFrame, prob::NonlinearProblem{T}, cache::FittingCache{T}, rng) where {S, T}

    # Number of time steps
    Nd = length(prob.pviews)

    # Auxiliary buffers
    idx = df[:, :Istart]
    if S === :Tconn
        bviews = view(cache.cbuf, idx)
        bviews2 = view(cache.cbuf2, idx)
    else
        bviews = view(cache.tbuf, idx)
        bviews2 = view(cache.tbuf2, idx)
    end

    # References to global vectors
    xviews = view(cache.ybuf, rng)
    gviews = view(cache.gbuf, rng)

    # References to internal parameters of model
    pviews, yviews = build_fitting_views(df, prob, S)
    # Parameter fixing periods
    V, vviews = build_validity_matrix(df, nrow(df), Nd)

    FittingParameter{S, T}(; xviews, pviews, yviews, gviews, bviews, bviews2, V, vviews)
end

function FittingSet{T}(df::AbstractDataFrame, prob::NonlinearProblem{T}, scale::AbstractParametersScaling{T}) where {T}

    # Number of connections and tanks
    Nc, Nt = size(prob.C)

    # Allocate memory for buffers
    df_view = getparams(df, Val(:tanks), Val(:var))
    cache = FittingCache{T}(Nt, Nc, nrow(df_view))

    stop = 0
    # FIXED: By itself 'map' for 'GroupedDataFrame'
    # is not possible (reserved), but possible by 'pairs(..)'
    gd = groupby(df_view, :Parameter)
    params = map(pairs(gd)) do (key, df)
        sym = key.Parameter
        start = stop + 1
        stop += nrow(df)
        FittingParameter{sym, T}(df, prob, cache, start:stop)
    end |> Tuple

    FittingSet{T, typeof(scale), typeof(params)}(; params, scale, cache)
end

function getparams!(fset::FittingSet)
    @unpack params, scale = fset
    @unpack xbuf, ybuf = fset.cache

    # FIXED: Using 'map' instead of 'for' preserves 'type-stability'
    map(params) do param
        @unpack xviews, yviews = param
        copyto!(xviews, yviews)
    end
    scalex!(xbuf, ybuf, scale)

    return xbuf
end

function getparams!!(x, fset::FittingSet)
    copyto!(x, getparams!(fset))
end

function setparams!(fset::FittingSet, xnew)
    @unpack params, scale = fset
    @unpack xbuf, ybuf = fset.cache

    copyto!(xbuf, xnew)
    unscalex!(ybuf, xbuf, scale)

    # FIXED: Using 'map' instead of 'for' preserves 'type-stability'
    map(params) do param
        @unpack xviews, pviews = param
        # FIXED: Faster than 'fill!.(pviews, xviews)'
        @inbounds @simd for i = 1:length(xviews)
            fill!(pviews[i], xviews[i])
        end
    end

    return ybuf
end


text:
module _XSampEn
export XSampEn
using Statistics: std
    """
        XSamp, A, B = XSampEn(Sig)

    Returns the cross-sample entropy estimates (`XSamp`) and the number of
    matched vectors (m:B, m+1:A) for m = [0,1,2] estimated for the two
    univariate data sequences contained in 'Sig' using the default parameters:
    embedding dimension = 2, time delay = 1,
    radius distance threshold= 0.2*SD(`Sig`),  logarithm = natural

        XSamp, A, B = XSampEn(Sig::AbstractArray{T,2} where T<:Real; m::Int=2, tau::Int=1, r::Real=0.2*std(Sig,corrected=false), Logx::Real=exp(1))

    Returns the cross-sample entropy estimates (`XSamp`) for dimensions [0,1,...,m]
    estimated between the data sequences in `Sig` using the specified
    'keyword'  arguments:

    # Arguments:
    `m`     - Embedding Dimension, a positive integer  [default: 2]   \n
    `tau`   - Time Delay, a positive integer         [default: 1]   \n
    `r`     - Radius Distance Threshold, a positive scalar [default: 0.2*SD(`Sig`)] \n
    `Logx`  - Logarithm base, a positive scalar      [default: natural] \n

    See also `XFuzzEn`, `XApEn`, `SampEn`, `SampEn2D`, `XMSEn`, `ApEn`

    # References:
    	[1] Joshua S Richman and J. Randall Moorman.
            "Physiological time-series analysis using approximate entropy
            and sample entropy."
            American Journal of Physiology-Heart and Circulatory Physiology
            (2000)

    """
    function XSampEn(Sig::AbstractArray{T,2} where T<:Real; m::Int=2,
        tau::Int=1, r::Real=0.2*std(Sig,corrected=false), Logx::Real=exp(1))

    (size(Sig,2) > size(Sig,1)) ? Sig = transpose(Sig) : nothing

    N = size(Sig,1)
    (N>10 && size(Sig,2)==2) ? nothing :  error("Sig:   must be a 2-columns matrix")
    (m > 0) ? nothing : error("m:     must be an integer > 0")
    (tau>0) ? nothing : error("tau:   must be an integer > 0")
    (r >=0) ? nothing : error("r:     must be a positive value")
    (Logx>0) ? nothing : error("Logx:     must be a positive number > 0")

    S1 = Sig[:,1]; S2 = Sig[:,2]
    Counter = 1*(abs.(S1 .- transpose(S2)) .<= r)
    M = vcat(m*ones(Int,N-(m*tau)), repeat((m-1):-1:1,inner=tau))
    A = zeros(Int,m+1)
    B = zeros(Int,m+1)
    A[1] = sum(Counter);  B[1] = N*N;

    for n = 1:N - tau
        ix = findall(Counter[n, :] .== 1)
        for k = 1:M[n]
            ix = ix[ix .+ (k*tau) .<= N]
            isempty(ix) ? break : nothing
            p1 = repeat(transpose(S1[n:tau:n+(tau*k)]), size(ix,1))
            p2 = S2[ix .+ transpose(collect(0:tau:(k*tau)))]
            ix = ix[findall(maximum(abs.(p1 - p2),dims=2) .<= r)]
            Counter[n, ix] .+= 1
        end
    end

    for k = 1:m
        A[k+1] = sum(Counter.>k)
        B[k+1] = sum(Counter.>=k)
    end

    XSamp = -log.(Logx, A./B)
    return XSamp, A, B
    end

end

"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub
"""


text:
function r3mv!(a, x, y)
    y[1] = a[1,1]*x[1] + a[1,2]*x[2] + a[1,3]*x[3]
    y[2] = a[2,1]*x[1] + a[2,2]*x[2] + a[2,3]*x[3]
    y[3] = a[3,1]*x[1] + a[3,2]*x[2] + a[3,3]*x[3]
    return
end


text:
@everywhere function pushwithlock!(circdq, what)
    withlock(circdq.lck) do
        push!(circdq, what)
    end
end

@everywhere function shiftwithlock!(circdq)
    withlock(circdq.lck) do
        shift!(circdq)
    end
end

@everywhere function popwithlock!(circdq)
    withlock(circdq.lck) do
        pop!(circdq)
    end
end

function splicewithlock(circdq)
    withlock(circdq.lck) do
        push!(circdq, 1)
        push!(circdq, 2)
        push!(circdq, 3)
        push!(circdq, 4)
        push!(circdq, 5)

        @test typeof(circdq[1:4]) == Vector{Int}
        @test length(circdq[1:end]) == length(circdq)

        @test splice!(circdq, 3) == 3
        @test splice!(circdq, 3) == 4
        @test splice!(circdq, 1) == 1
        @test splice!(circdq, 1) == 2
        @test splice!(circdq, 1) == 5
        @test length(circdq) == 0
    end
end

function test_shared_circdq()
    path = pwd()
    println("testing circular deque")
    master_circdq = SharedCircularDeque{Int}(path, 10; create=true)
    println("    closing and...")
    close(master_circdq)
    println("    deleting and...")
    delete!(master_circdq)
    println("    ...creating new")
    master_circdq = SharedCircularDeque{Int}(path, 10; create=true)

    try
        W = workers()
        println("    workers: ", W)
        for w in W
            remotecall_wait(()->(global circdq=SharedCircularDeque{Int}(path, 10); nothing), w)
        end
        println("    created on workers")

        @test length(master_circdq) == 0
        splicewithlock(master_circdq)

        remotecall_wait(()->(pushwithlock!(circdq, 2); nothing), W[1])
        @test remotecall_fetch(()->length(circdq), W[1]) == 1
        @test length(master_circdq) == 1
        @test 2 in master_circdq
        @test !(3 in master_circdq)
        remotecall_wait(()->(pushwithlock!(circdq, 3); nothing), W[2])
        @test remotecall_fetch(()->length(circdq), W[2]) == 2
        @test length(master_circdq) == 2
        @test 3 in master_circdq
        @test !(4 in master_circdq)
        remotecall_wait(()->(shiftwithlock!(circdq); nothing), W[1])
        @test remotecall_fetch(()->length(circdq), W[1]) == 1
        @test length(master_circdq) == 1
        @test !(2 in master_circdq)
        remotecall_wait(()->(popwithlock!(circdq); nothing), W[2])
        @test remotecall_fetch(()->length(circdq), W[2]) == 0
        @test length(master_circdq) == 0
        @test !(3 in master_circdq)

        splicewithlock(master_circdq)

        println("    test done")
        println("    closing...")
        close(master_circdq)
        remotecall_wait(()->close(circdq), W[1])
        remotecall_wait(()->close(circdq), W[2])
    finally
        println("    deleting...")
        delete!(master_circdq)
        println("    done.")
    end
end


text:
base_dir = dirname(dirname(pathof(PowerSystems)))

sys = PowerSystems.parse_standard_files(joinpath(MATPOWER_DIR, "case5_re.m"))

@testset "Check bus index" begin
    @test sort([b.number for b in collect(get_components(Bus, sys))]) == [1, 2, 3, 4, 10]
    @test sort(collect(Set([b.arc.from.number for
                b in collect(get_components(Branch,sys))]))) == [1, 2, 3, 4]
    @test sort(collect(Set([b.arc.to.number for
                b in collect(get_components(Branch,sys))]))) == [2, 3, 4, 10]

    # TODO: add test for loadzones testing MAPPING_BUSNUMBER2INDEX

end

@testset "Test unique bus numbers" begin
    number = 100
    bus1 = Bus(;
               number=number,
               name="bus100",
               bustype=PSY.PV::BusType,
               angle=1.0,
               voltage=1.0,
               voltagelimits=(min=-1.0, max=1.0),
               basevoltage=1.0
              )
    bus2 = Bus(;
               number=number,
               name="bus101",
               bustype=PSY.PV::BusType,
               angle=1.0,
               voltage=1.0,
               voltagelimits=(min=-1.0, max=1.0),
               basevoltage=1.0
              )

    add_component!(sys, bus1)
    @test_throws ArgumentError add_component!(sys, bus2)
end


text:
#=
Direct estimation of the probability of for a diffusion with X(0) = x₀ satisfying
X(T) ∈ B for the Muller potential.
=#

using Distributed
# nw = 4; # number of workers
# addprocs(nw);

using Statistics
using HypothesisTests
using Printf
@everywhere include("muller_setup.jl")
opts = MDOptions(n_iters=nΔt)

Random.seed!(100)
@everywhere n_trials = 10^4;
n_in_set = @sync @distributed (+) for j in 1:n_trials
    X = copy(x₀);
    sample_trajectory!(X, sampler, options=opts);
    f(X)
end
CI = confint(BinomialTest(n_in_set, n_trials));
@printf("95%% CI = (%g,%g)\n", CI[1], CI[2])


text:
"""
    loaddata(file)

Loads methylation data.

`file` is a `RData` file containing a matrix named `betas.dasen` of size `(n x m)` for `n` samples and `m` probes.

Sample names and probe names taken from colnames and rownames respectively.

"""
function loaddata(file)

    R"""
        load($file)
        samples <- colnames(betas.dasen)
        probes  <- rownames(betas.dasen)
        beta    <- betas.dasen
    """

    @rget samples probes beta

    samples, probes, beta

end

"""
    loadsamplemeta(file)

Loads sample meta, a tab or commma separated file that must contain the following first threecolumns (in any order):

| Sample    | Sex   | Age   |  Filter  |
|-----------|-------|-------|----------|
|  Sample_1 | F     | 1.0   |   1      |
|  Sample_2 | M     | 2.0   |   0      |
|  Sample_3 | M     | 10.0  |   1      |

- `Sample` unique string identifier
- `Sex`  {M, F} or missing
- `Age` floating point age
- `Filter` optional column of `{1, 0}`s to exclude samples, when `Filter = 0` sample is excluded.

`Sample` should match `colnames` of RData file.

"""
loadsamplemeta(file) = loadmeta(file, [:Sample, :Sex, :Age], "sample")

"""
    loadprobemeta(file)

Loads probe meta, a tab or commma separated file that must contain the following columns (in any order):

| Probe      |
|------------|
| cg00000029 |
| cg00000103 |
| cg00000109 |

Can be used to override the `rownames` in `RData` file.

"""
loadprobemeta(file) = loadmeta(file, [:Probe], "probe")

"""
    loadmeta(file, mandatorycols, label)

    Helper function used to load metafiles
"""
function loadmeta(file, mandatorycols, label)
    meta = CSV.read(file, DataFrame)

    if length(intersect(mandatorycols, propertynames(meta))) != length(mandatorycols)
        missingcols = setdist(mandatorycols, propertynames(meta))
        colstring = join(names(meta), ", ")
        mcolstring = join(missingcols, ", ")
        error("Loading $label meta. Found cols: $colstring\nMissing: $mcolstring")
    end
    meta
end

"""
    loadall(samplemetafile, betafile)

Loads samplefmetafile and `RData` file `betafile` ensures samples details correspond and filters samples
"""
function loadall(samplemetafile, betafile, verbose=true)

    verbose && println("[GPM]\tLoaded sample meta file ", samplemetafile)
    verbose && println("[GPM]\tLoaded probe data file ", betafile)
    meta = loadsamplemeta(samplemetafile)
    samples, probes, beta = loaddata(betafile)

    if meta.Sample != samples
        error("Sample meta does not correspond to beta colnames")
    end

    verbose && println("[GPM]\tLoaded meta information for ", size(meta, 1), " samples.")
    verbose && println("[GPM]\tLoaded probe data on  ", size(beta, 1), " probes and ", size(beta, 2), " samples.")

    ### Filter samples
    if :Filter ∈ propertynames(meta)
        ind = meta.Filter .== 1

        if verbose
            println("[GPM]\tExcluding ", sum(.!ind), " samples:")
            display(meta[.!ind, :])
        end
        meta = meta[ind, :]
        beta = beta[:, ind]

    end

    ### Converts row major R to column major Julia
    meta, probes, Matrix(beta')
end


text:
"""
	match_patches(img::AbstractArray{Float64},
			refIndex::Array{CartesianIndex{2},2},
			patchSize::CartesianIndex{2},
			searchWindow::CartesianIndex{2},
			nMatch::Int64,
			distF::Function)

Full-search block matching algorithm for BM3DDenoise
"""
function match_patches(
	img::AbstractArray{Float64,2},
	refIndex::Array{CartesianIndex{2},2},
	patchSize::CartesianIndex{2},
	searchWindow::CartesianIndex{2},
	nMatch::Int64,
	distF::Function,
)

	refIndexH, refIndexW = size(refIndex)

	matches = Array{CartesianIndex{2}}(undef, refIndexH, refIndexW, nMatch)

	dist = Array{Float64}(undef, refIndexH, refIndexW, nMatch + 1)# dist[i, j, end] is the largest item in current dist[i, j, 1:end-1]
	dist .= typemax(Float64)

	maxDistMatch = Array{UInt}(undef, refIndexH, refIndexW)# maxDistMatch[i, j] is the position in matches[i, j, :] corresponding to the largest distance block
	maxDistMatch .= 1

	lockPool = Array{ReentrantLock}(undef, refIndexH, refIndexW)
	for i in 1:length(lockPool)
		lockPool[i] = ReentrantLock()
	end

	@inbounds Threads.@threads for j in 1:refIndexW
		for i in 1:refIndexH
			p = refIndex[i, j]
			iE = min(refIndexH, i + searchWindow[1])
			jE = min(refIndexW, j + searchWindow[2])
			imgRef = @view img[p:(p + patchSize - CartesianIndex(1, 1))]

			# In order to skip the current block and save unnecessary judgment, it is split into three parts.
			for sj in (j + 1):jE, si in (i + 1):iE
				@views weight = distF(
					imgRef,
					img[refIndex[si, sj]:(refIndex[
						si,
						sj,
					] + patchSize - CartesianIndex(1, 1))],
				)
				lock(lockPool[i, j]) do
					if weight < dist[i, j, nMatch + 1]
						maxDistImg = maxDistMatch[i, j]
						matches[i, j, maxDistImg] = CartesianIndex(si, sj)
						dist[i, j, maxDistImg] = weight
						dist[i, j, nMatch + 1], maxDistMatch[i, j] =
							findmax(@views dist[i, j, 1:nMatch])
					end
				end

				lock(lockPool[si, sj]) do
					if weight < dist[si, sj, nMatch + 1]
						maxDistImg = maxDistMatch[si, sj]
						matches[si, sj, maxDistImg] = CartesianIndex(i, j)
						dist[si, sj, maxDistImg] = weight
						dist[si, sj, nMatch + 1], maxDistMatch[si, sj] =
							findmax(@views dist[si, sj, 1:nMatch])
					end
				end
			end
			for si in (i + 1):iE
				@views weight = distF(
					imgRef,
					img[refIndex[si, j]:(refIndex[si, j] + patchSize - CartesianIndex(1, 1))],
				)
				lock(lockPool[i, j]) do
					if weight < dist[i, j, nMatch + 1]
						maxDistImg = maxDistMatch[i, j]
						matches[i, j, maxDistImg] = CartesianIndex(si, j)
						dist[i, j, maxDistImg] = weight
						dist[i, j, nMatch + 1], maxDistMatch[i, j] =
							findmax(@views dist[i, j, 1:nMatch])
					end
				end

				lock(lockPool[si, j]) do
					if weight < dist[si, j, nMatch + 1]
						maxDistImg = maxDistMatch[si, j]
						matches[si, j, maxDistImg] = CartesianIndex(i, j)
						dist[si, j, maxDistImg] = weight
						dist[si, j, nMatch + 1], maxDistMatch[si, j] =
							findmax(@views dist[si, j, 1:nMatch])
					end
				end
			end
			for sj in (j + 1):jE
				@views weight = distF(
					imgRef,
					img[refIndex[i, sj]:(refIndex[i, sj] + patchSize - CartesianIndex(1, 1))],
				)
				lock(lockPool[i, j]) do
					if weight < dist[i, j, nMatch + 1]
						maxDistImg = maxDistMatch[i, j]
						matches[i, j, maxDistImg] = CartesianIndex(i, sj)
						dist[i, j, maxDistImg] = weight
						dist[i, j, nMatch + 1], maxDistMatch[i, j] =
							findmax(@views dist[i, j, 1:nMatch])
					end
				end

				lock(lockPool[i, sj]) do
					if weight < dist[i, sj, nMatch + 1]
						maxDistImg = maxDistMatch[i, sj]
						matches[i, sj, maxDistImg] = CartesianIndex(i, j)
						dist[i, sj, maxDistImg] = weight
						dist[i, sj, nMatch + 1], maxDistMatch[i, sj] =
							findmax(@views dist[i, sj, 1:nMatch])
					end
				end
			end
		end
	end
	matches
end

# For color images
function match_patches(
	img::Array{Float64,3},
	refIndex::Array{CartesianIndex{2},2},
	patchSize::CartesianIndex{2},
	searchWindow::CartesianIndex{2},
	nMatch::Int64,
	distF::Function,
)
	match_patches(view(img, :, :, 1), refIndex, patchSize, searchWindow, nMatch, distF)
end


text:
using Test
using DFTK: PlaneWaveBasis, Model, G_vectors, index_G_vectors
using LinearAlgebra

include("testcases.jl")

function test_pw_cutoffs(testcase, Ecut, fft_size)
    model = Model(testcase.lattice, n_electrons=testcase.n_electrons)
    pw = PlaneWaveBasis(model, Ecut, testcase.kcoords, testcase.ksymops; fft_size=fft_size)

    for (ik, kpt) in enumerate(pw.kpoints)
        for G in G_vectors(kpt)
            @test sum(abs2, model.recip_lattice * (kpt.coordinate + G)) ≤ 2 * Ecut
        end
    end
end

@testset "PlaneWaveBasis: Check struct construction" begin
    Ecut = 3
    fft_size = [15, 15, 15]
    model = Model(silicon.lattice, n_electrons=silicon.n_electrons)
    pw = PlaneWaveBasis(model, Ecut, silicon.kcoords, silicon.ksymops; fft_size=fft_size)

    @test pw.model.lattice == silicon.lattice
    @test pw.model.recip_lattice ≈ 2π * inv(silicon.lattice)
    @test pw.model.unit_cell_volume ≈ det(silicon.lattice)
    @test pw.model.recip_cell_volume ≈ (2π)^3 * det(inv(silicon.lattice))

    @test pw.Ecut == 3
    @test pw.fft_size == Tuple(fft_size)

    g_start = -ceil.(Int, (Vec3(pw.fft_size) .- 1) ./ 2)
    g_stop  = floor.(Int, (Vec3(pw.fft_size) .- 1) ./ 2)
    g_all = vec(collect(G_vectors(pw)))

    for (ik, kpt) in enumerate(pw.kpoints)
        kpt = pw.kpoints[ik]
        @test kpt.coordinate == silicon.kcoords[pw.krange_thisproc[ik]]

        for (ig, G) in enumerate(G_vectors(kpt))
            @test g_start <= G <= g_stop
        end
        @test g_all[kpt.mapping] == G_vectors(kpt)
    end
    @test pw.kweights == ([1, 8, 6, 12] / 27)[pw.krange_thisproc]
end

@testset "PlaneWaveBasis: Energy cutoff is respected" begin
    test_pw_cutoffs(silicon, 4.0, [15, 15, 15])
    test_pw_cutoffs(silicon, 3.0, [15, 13, 13])
    test_pw_cutoffs(silicon, 4.0, [11, 13, 11])
end

@testset "PlaneWaveBasis: Check cubic basis and cubic index" begin
    Ecut = 3
    fft_size = [15, 15, 15]
    model = Model(silicon.lattice, n_electrons=silicon.n_electrons)
    pw = PlaneWaveBasis(model, Ecut, silicon.kcoords, silicon.ksymops; fft_size=fft_size)
    g_all = collect(G_vectors(pw))

    for i in 1:15, j in 1:15, k in 1:15
        @test index_G_vectors(pw, g_all[i, j, k]) == CartesianIndex(i, j, k)
    end
    @test index_G_vectors(pw, [15, 1, 1]) === nothing
    @test index_G_vectors(pw, [-15, 1, 1]) === nothing
end

@testset "PlaneWaveBasis: Check index for kpoints" begin
    Ecut = 3
    fft_size = [7, 9, 11]
    model = Model(silicon.lattice, n_electrons=silicon.n_electrons)
    basis = PlaneWaveBasis(model, Ecut, silicon.kcoords, silicon.ksymops; fft_size=fft_size)
    g_all = collect(G_vectors(basis))

    for kpt in basis.kpoints
        for (iball, ifull) in enumerate(kpt.mapping)
            @test index_G_vectors(basis, kpt, g_all[ifull]) == iball
        end

        if kpt.coordinate == [1/3, 1/3, 0]
            @test index_G_vectors(basis, kpt, [-2, -3, -1]) == 62
        else
            @test index_G_vectors(basis, kpt, [-2, -3, -1]) === nothing
        end
        @test index_G_vectors(basis, kpt, [15, 1, 1]) === nothing
        @test index_G_vectors(basis, kpt, [-15, 1, 1]) === nothing
    end
end


text:
include("Balances.jl")
using Sundials;

# ----------------------------------------------------------------------------------- #
# Copyright (c) 2016 Varnerlab
# School of Chemical Engineering Purdue University
# W. Lafayette IN 46907 USA

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #
function SolveBalances(TSTART,TSTOP,Ts,data_dictionary)
# ----------------------------------------------------------------------------------- #
# SolveBalances.jl was generated using the Kwatee code generation system.
# SolveBalances: Solves model equations from TSTART to TSTOP given parameters in data_dictionary.
# Username: jeffreyvarner
# Type: GRN-JULIA
# Version: 1.0
# Generation timestamp: 02-26-2016 14:31:22
#
# Input arguments:
# TSTART  - Time start
# TSTOP  - Time stop
# Ts - Time step
# data_dictionary  - Data dictionary instance (holds model parameters)
#
# Return arguments:
# TSIM - Simulation time vector
# X - Simulation state array (NTIME x NSPECIES)
# ----------------------------------------------------------------------------------- #

# Get required stuff from DataFile struct -
TSIM = collect(TSTART:Ts:TSTOP);
initial_condition_vector = data_dictionary["INITIAL_CONDITION_ARRAY"];

# Call the ODE solver -
fbalances(t,y,ydot) = Balances(t,y,ydot,data_dictionary);
X = Sundials.cvode(fbalances,initial_condition_vector,TSIM,reltol=1e-4,abstol=1e-6);

return (TSIM,X);
end


text:
include("bomex_model.jl")

function main()
    # add a command line argument to specify the kind of surface flux
    # TODO: this will move to the future namelist functionality
    bomex_args = ArgParseSettings(autofix_names = true)
    add_arg_group!(bomex_args, "BOMEX")
    @add_arg_table! bomex_args begin
        "--surface-flux"
        help = "specify surface flux for energy and moisture"
        metavar = "prescribed|bulk"
        arg_type = String
        default = "prescribed"
        "--moisture-model"
        help = "specify cloud condensate model"
        metavar = "equilibrium|nonequilibrium"
        arg_type = String
        default = "equilibrium"
    end

    cl_args =
        ClimateMachine.init(parse_clargs = true, custom_clargs = bomex_args)

    surface_flux = cl_args["surface_flux"]
    moisture_model = cl_args["moisture_model"]

    FT = Float64
    config_type = AtmosLESConfigType

    # DG polynomial order
    N = 4
    # Domain resolution and size
    Δh = FT(100)
    Δv = FT(40)

    resolution = (Δh, Δh, Δv)

    # Prescribe domain parameters
    xmax = FT(6400)
    ymax = FT(6400)
    zmax = FT(3000)

    t0 = FT(0)

    # For a full-run, please set the timeend to 3600*6 seconds
    # and change the values in ConservationCheck
    # For the test we set this to == 20 minutes
    timeend = FT(1200)
    #timeend = FT(3600 * 6)
    CFLmax = FT(0.35)

    # Choose default IMEX solver
    ode_solver_type = ClimateMachine.IMEXSolverType()

    model = bomex_model(
        FT,
        config_type,
        zmax,
        surface_flux,
        moisture_model = moisture_model,
    )
    ics = model.problem.init_state_prognostic
    # Assemble configuration
    driver_config = ClimateMachine.AtmosLESConfiguration(
        "BOMEX",
        N,
        resolution,
        xmax,
        ymax,
        zmax,
        param_set,
        ics;
        solver_type = ode_solver_type,
        model = model,
    )

    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        init_on_cpu = true,
        Courant_number = CFLmax,
        CFL_direction = HorizontalDirection(),
    )
    dgn_config = config_diagnostics(driver_config)

    if moisture_model == "equilibrium"
        filter_vars = ("moisture.ρq_tot",)
    elseif moisture_model == "nonequilibrium"
        filter_vars = ("moisture.ρq_tot", "moisture.ρq_liq", "moisture.ρq_ice")
    end

    cbtmarfilter = GenericCallbacks.EveryXSimulationSteps(1) do
        Filters.apply!(
            solver_config.Q,
            filter_vars,
            solver_config.dg.grid,
            TMARFilter(),
        )
        nothing
    end

    check_cons = (
        ClimateMachine.ConservationCheck("ρ", "3000steps", FT(0.0001)),
        ClimateMachine.ConservationCheck("ρe", "3000steps", FT(0.0025)),
    )

    result = ClimateMachine.invoke!(
        solver_config;
        user_callbacks = (cbtmarfilter,),
        diagnostics_config = dgn_config,
        check_cons = check_cons,
        check_euclidean_distance = true,
    )
end

main()

max_forks_repo_name: jm-c/CLIMA


text:
using IterTools
using LinearAlgebra
using Distributions
using OrdinaryDiffEq

export get_max_num, master_operator, steady_state_masterOP, dynamics_masterOP,
		dynamics_masterOP_diag

"""
	get_max_num(x0::Vector{Float64}, precision::Float64)

	Calculates the vector of maximum number of species allowed where the
	Poisson distribution is below a certain probability (precision)

	Mean of Poisson distribution for each species is the initial condition x0

	Returns a vector of max_num

"""
function get_max_num(x0::Vector{Float64}, precision::Float64)

	num_species = length(x0)
	max_num = zeros(Int, num_species)
	for i in 1:num_species
		x,n = 1,0
		while abs(x) > precision
			x = 1-cdf(Poisson(x0[i]), n)
			n += 1
		end
		max_num[i] = n
	end
	return max_num

end
"""
	master_operator(p::Parameters, max_num; α=1.0)

	Calculates the master operator for the reaction network
	where the state space will be truncated at a maximum number
	of species `max_num`

	Input
	p		Parameters struct for the reaction network
	max_num	number of species where the state space is truncated
	α		Plefka expansion parameter, default is α=1.0

	Returns the master operator and state space in a tuple

"""
function master_operator(p::Parameters, max_num::Vector{Int}; α=1.0)

	num_species = length(p.k[1])	# number of species
	params = p.k	# reaction rates
	s_i = p.s_i		# Stoichiometric products
	r_i = p.r_i		# Stoichiometric reactants

	# state space
	x = [collect(0:max_num[j]) for j in 1:num_species]
	state_space = reduce(vcat, collect(Iterators.product(x...)))

	# master operator
	len_state_space = prod(max_num .+ 1)
	master = zeros(len_state_space, len_state_space)

	for s in state_space

		# get index from s (assume only one appearence)
		idx_s = findfirst(x->x==s, state_space)

		for j in 1:num_species

			# Annihiliation
			t = collect(s)
			t[j] += 1
			t = tuple(t...)
			if t[j] <= max_num[j]
				idx_t = findfirst(x->x==t, state_space)		# index of t
				master[idx_s, idx_t] += params[2][j]*state_space[idx_t][j]
			end

			# Creation
			t = collect(s)
			t[j] -= 1
			t = tuple(t...)
			if t[j] >= 0
				idx_t = findfirst(x->x==t, state_space) 	# index of t
				master[idx_s, idx_t] += params[1][j]
			end

		end

		# Interaction
		for β in 1:length(params[3])

			new_state = tuple([(s[m] - s_i[β,m] + r_i[β,m]) for m in 1:num_species]...)

			if all(state_space[1] .<= new_state) && all(state_space[end] .>= new_state)
				idx = findfirst(x->x==new_state, state_space)
				r = α*params[3][β]
				for n in 1:num_species
					for p in 0:r_i[β,n]-1
						r *= (state_space[idx][n] - p)
					end
				end
				master[idx_s, idx] += r
			end

		end

	end

	# sum along columns must be zero
	master[diagind(master)] = -sum(master, dims=1)

	return (master, state_space)

end

function master_operator(p::Parameters, max_num::Int; α=1.0)
	max_num_vec = max_num .* ones(Int, length(p.k[1]))
	return master_operator(p::Parameters, max_num_vec; α=α)
end

"""
	steady_state_masterOP(master, state_space)

	Calculates the steady state of the master operator and state space

	Returns the steady state

"""
function steady_state_masterOP(master, state_space, E::Eigen)

	num_species = length(state_space[1])	# number of species

	# index of largest eigenvalue
	max_eig_idx = argmax(real.(E.values))

	# sum of eigenvectors of the highest eigenvalue
	sum_eigvecs = sum(abs.(E.vectors[:,max_eig_idx]))

	x0 = zeros(num_species)		# returned steady state

	for i in 1:length(state_space)
		for j in 1:num_species
			x0[j] += state_space[i][j] * abs(E.vectors[i,max_eig_idx]) / sum_eigvecs
		end
	end

	return x0

end

function steady_state_masterOP(master, state_space)
	E = eigen(master, sortby=nothing)		# eigenvalues and eigenvectors
	return steady_state_masterOP(master, state_space, E)
end

"""
	calc_mean_masterOP(p, state_space)

	Calculates the mean concentration with the probability vector `p`
	and state space

	Returns the calculated mean

"""
function calc_mean_masterOP(p::Vector{Float64}, state_space)

	num_species = length(state_space[1])	# number of species
	x = zeros(num_species)					# returned mean

	# sum of probabilities, should be == 1
	sum_p = sum(abs.(p))

	for i in 1:length(state_space)
		for j in 1:num_species
			x[j] += state_space[i][j] * abs(p[i]) / sum_p
		end
	end

	return x

end

"""
	initial_distr_state_space(state_space, x0)

	Calculates the initial Poisson distribution with means `x0`
	for the different species.

	Returns the initial probability distribution vector

"""
function initial_distr_state_space(state_space, x0::Vector{Float64})

	num_species = length(state_space[1])	# number of species
	p0 = ones(length(state_space))			# initial prob. distr.

	for i in 1:length(state_space)
		for j in 1:num_species
			p0[i] *= pdf(Poisson(x0[j]), state_space[i][j])
		end
	end

	return p0./sum(p0)

end

"""
	dynamics_masterOP(master, state_space, tspan, x0, method)

	Runs the dynamics for the master operator using an euler step integrator.

	Input
	master		master operator
	state_space	truncated state space
	tspan		time grid
	x0			initial condition for copy numbers
	method		integration algorithm. If not specified: default is Tsit5()

	You can find a list of available solvers at https://diffeq.sciml.ai/stable/solvers/ode_solve/

	Returns a Result struct with time_grid and mean copy numbers

"""
function dynamics_masterOP(master, state_space, tspan::Vector{Float64}, x0::Vector{Float64}, method; rtn_prob=false)

	# Differential equation
	function f!(du,u,M,t)
		du .= M*u
	end

	# Solve the differential equation
	p0 = initial_distr_state_space(state_space, x0)
    P = ode_solver(f!, master, tspan, p0, method; savedt=true)

    # Convert probability distributions to mean copy numbers
	y = zeros(length(x0), length(tspan))
	for t in 1:length(tspan)
		y[:,t] = calc_mean_masterOP(P.data[:,t], state_space)
	end

	if rtn_prob
		return (P.data, Result(tspan, y))
	else
		return Result(tspan, y)
	end
end

function dynamics_masterOP(master, state_space, tspan::Vector{Float64}, x0::Vector{Float64})
	println("No integration algorithm specified! Using Tsit5() as default!")
	dynamics_masterOP(master, state_space, tspan, x0, Tsit5)
end

"""
	dynamics_masterOP_diag(master, state_space, E::Eigen, tspan, x0)

	Runs the dynamics for the master operator by diagonalizing the master
	operator and directly calculates the matrix exponential.

	Solves the master equation   dP/dt = M*P   by the solution

		P(t) = exp(M*t) * P0

	Input
	master		master operator
	state_space	truncated state space
	E			Eigen values and vectors, calculate by `eigen(master)`
	tspan		time grid
	x0			initial condition for copy numbers

	Returns a Result struct with time_grid and mean copy numbers

"""
function dynamics_masterOP_diag(master, state_space, E::Eigen, tspan, x0; rtn_prob=false)

	er = E.vectors
	el = inv(er)
	d = E.values #real.(E.values)

	p0 = initial_distr_state_space(state_space, x0)
	y = zeros(length(x0), length(tspan))
	p = zeros(length(state_space), length(tspan))

	elp0 = el*p0

	for j in 1:length(tspan)
		for i in 1:length(state_space)
			@inbounds p[i,j] = abs(sum(er[i,:].*exp.(d*tspan[j]).*elp0))
		end
		p[:,j] ./= sum(p[:,j])
		y[:,j] = CRN.calc_mean_masterOP(p[:,j], state_space)
	end

	if rtn_prob
		return (p, Result(tspan, y))
	else
		return Result(tspan, y)
	end

end

function dynamics_masterOP_diag(master, state_space, tspan, x0)
	E = eigen(master, sortby=nothing)		# eigenvalues and eigenvectors
	return dynamics_masterOP_diag(master, state_space, E, tspan, x0)
end

max_stars_repo_name: mkurjahn/CRN

max_issues_repo_name: mkurjahn/CRN

max_forks_repo_name: mkurjahn/CRN


text:
"""Basic interface for spline functions
"""

using StatsModels
using StatsBase

const SplineTermContext = Any

abstract type AbstractSplineTerm <: AbstractTerm end

shorthand(::Type{<:AbstractSplineTerm}) = "spline"
spline_term(v::AbstractSplineTerm) = throw("term not implemented")
spline_knots(v::AbstractSplineTerm) = throw("knots not implemented")
spline_deg(v::AbstractSplineTerm) = throw("deg not implemented")

function Base.show(io::IO, t::T) where T<:AbstractSplineTerm
    print(io, "$(shorthand(T))($(spline_term(t)), $(spline_knots(t)))")
end

function StatsBase.coefnames(t::T) where T <: AbstractSplineTerm
    c = coefnames(spline_term(t))
    k = 1:width(t)

    map(k) do knot
        "$(shorthand(T))($(c), $(t.knots.n))[$knot]"
    end
end


text:
using BinDeps

# Heavily taken from Sundials.jl build.jl file.
@BinDeps.setup

fannlib = library_dependency("libdoublefann")
uri = "https://github.com/libfann/fann/archive/master.zip"

provides(Sources, URI(uri), fannlib, unpacked_dir="fann-master")

prefix = joinpath(BinDeps.depsdir(fannlib),"usr")
srcdir = joinpath(BinDeps.depsdir(fannlib),"src", "fann-master")

provides(SimpleBuild,
    (@build_steps begin
        GetSources(fannlib)
        @build_steps begin
            ChangeDirectory(srcdir)
            `cmake -DCMAKE_INSTALL_PREFIX=$prefix .`
            `make install`
        end
    end), fannlib)

@BinDeps.install [:libdoublefann => :libfann]


text:
module Issue50
using Test
using Configurations

@option struct OptionA
    name::Maybe{String} = nothing
    int::Int = 1
end

@option struct OptionB
    name::Maybe{OptionA} = nothing
    int::Int = 1
end

@testset "to_dict style" begin
    option = OptionB()
    @test to_dict(option, TOMLStyle) ==
          to_dict(option; include_defaults=true, exclude_nothing=true)
    @test to_dict(option, YAMLStyle) ==
          to_dict(option; include_defaults=true, exclude_nothing=false)
    @test to_dict(option, JSONStyle) ==
          to_dict(option; include_defaults=true, exclude_nothing=false)
end

@testset "#50" begin
    @test from_dict(Issue50.OptionB, Dict{String,Any}("name" => nothing)) == OptionB()
end

end


text:
dev = device()

pool = memory_pool(dev)

pool2 = CuMemoryPool(dev)
@test pool2 != pool
memory_pool!(dev, pool2)
@test pool2 == memory_pool(dev)
@test pool2 != default_memory_pool(dev)

memory_pool!(dev, pool)
@test pool == memory_pool(dev)


text:
using Quadrature, Cuba, Cubature, Zygote, FiniteDiff, ForwardDiff
using Test

### One Dimensional

f(x,p) = sum(sin.(p[1] * x))
lb = 1.0
ub = 3.0
p = 2.0
prob = QuadratureProblem(f,lb,ub,p)
sol = solve(prob,QuadGKJL(),reltol=1e-3,abstol=1e-3)

function testf(lb,ub,p)
    prob = QuadratureProblem(f,lb,ub,p)
    sin(solve(prob,QuadGKJL(),reltol=1e-3,abstol=1e-3)[1])
end
dlb1,dub1,dp1 = Zygote.gradient(testf,lb,ub,p)
dlb2 = FiniteDiff.finite_difference_derivative(lb->testf(lb,ub,p),lb)
dub2 = FiniteDiff.finite_difference_derivative(ub->testf(lb,ub,p),ub)
dp2 = FiniteDiff.finite_difference_derivative(p->testf(lb,ub,p),p)

# dlb3 = ForwardDiff.derivative(lb->testf(lb,ub,p),lb)
# dub3 = ForwardDiff.derivative(ub->testf(lb,ub,p),ub)
dp3 = ForwardDiff.derivative(p->testf(lb,ub,p),p)

#@test dlb1 ≈ dlb3
#@test dub1 ≈ dub3
@test dp1 ≈ dp3

### N-dimensional

f(x,p) = sum(sin.(x .* p))
lb = ones(2)
ub = 3ones(2)
p = [1.5,2.0]
prob = QuadratureProblem(f,lb,ub,p)
sol = solve(prob,CubaCuhre(),reltol=1e-3,abstol=1e-3)

function testf(lb,ub,p)
    prob = QuadratureProblem(f,lb,ub,p)
    sin(solve(prob,CubaCuhre(),reltol=1e-6,abstol=1e-6)[1])
end

function testf2(lb,ub,p)
    prob = QuadratureProblem(f,lb,ub,p)
    sin(solve(prob,HCubatureJL(),reltol=1e-6,abstol=1e-6)[1])
end

dlb1,dub1,dp1 = Zygote.gradient(testf,lb,ub,p)
dlb2 = FiniteDiff.finite_difference_gradient(lb->testf(lb,ub,p),lb)
dub2 = FiniteDiff.finite_difference_gradient(ub->testf(lb,ub,p),ub)
dp2 = FiniteDiff.finite_difference_gradient(p->testf(lb,ub,p),p)

@test_broken dlb1 ≈ dlb2
@test_broken dub1 ≈ dub2
@test dp1 ≈ dp2

# dlb3 = ForwardDiff.gradient(lb->testf(lb,ub,p),lb)
# dub3 = ForwardDiff.gradient(ub->testf(lb,ub,p),ub)
dp3 = ForwardDiff.gradient(p->testf2(lb,ub,p),p)

#@test dlb1 ≈ dlb3
#@test dub1 ≈ dub3
@test dp1 ≈ dp3

# dlb4 = ForwardDiff.gradient(lb->testf(lb,ub,p),lb)
# dub4 = ForwardDiff.gradient(ub->testf(lb,ub,p),ub)
dp4 = ForwardDiff.gradient(p->testf(lb,ub,p),p)

#@test dlb1 ≈ dlb4
#@test dub1 ≈ dub4
@test dp1 ≈ dp4

### N-dimensional N-out

f(x,p) = sin.(x .* p)
lb = ones(2)
ub = 3ones(2)
p = [1.5,2.0]
prob = QuadratureProblem(f,lb,ub,p,nout=2)
sol = solve(prob,CubaCuhre(),reltol=1e-3,abstol=1e-3)

function testf(lb,ub,p)
    prob = QuadratureProblem(f,lb,ub,p,nout=2)
    sin(sum(solve(prob,CubaCuhre(),reltol=1e-6,abstol=1e-6)))
end

function testf2(lb,ub,p)
    prob = QuadratureProblem(f,lb,ub,p,nout=2)
    sin(sum(solve(prob,HCubatureJL(),reltol=1e-6,abstol=1e-6)))
end

dlb1,dub1,dp1 = Zygote.gradient(testf,lb,ub,p)
dlb2 = FiniteDiff.finite_difference_gradient(lb->testf(lb,ub,p),lb)
dub2 = FiniteDiff.finite_difference_gradient(ub->testf(lb,ub,p),ub)
dp2 = FiniteDiff.finite_difference_gradient(p->testf(lb,ub,p),p)

@test_broken dlb1 ≈ dlb2
@test_broken dub1 ≈ dub2
@test dp1 ≈ dp2

# dlb3 = ForwardDiff.gradient(lb->testf(lb,ub,p),lb)
# dub3 = ForwardDiff.gradient(ub->testf(lb,ub,p),ub)
dp3 = ForwardDiff.gradient(p->testf2(lb,ub,p),p)

#@test dlb1 ≈ dlb3
#@test dub1 ≈ dub3
@test dp1 ≈ dp3

# dlb4 = ForwardDiff.gradient(lb->testf(lb,ub,p),lb)
# dub4 = ForwardDiff.gradient(ub->testf(lb,ub,p),ub)
dp4 = ForwardDiff.gradient(p->testf(lb,ub,p),p)

#@test dlb1 ≈ dlb4
#@test dub1 ≈ dub4
@test dp1 ≈ dp4

### Batch Single dim
f(x,p) = x*p[1].+p[2]*p[3]

lb =1.0
ub = 3.0
p = [2.0, 3.0, 4.0]
prob = QuadratureProblem(f,lb,ub,p)

function testf3(lb,ub,p; f=f)
    prob = QuadratureProblem(f,lb,ub,p, batch = 10, nout=1)
    solve(prob, CubatureJLh(); reltol=1e-3,abstol=1e-3)[1]
end

dp1 = ForwardDiff.gradient(p->testf3(lb,ub,p),p)
dp2 = Zygote.gradient(p->testf3(lb,ub,p),p)[1]
dp3 = FiniteDiff.finite_difference_gradient(p->testf3(lb,ub,p),p)

@test dp1 ≈ dp3 #passes
@test dp2 ≈ dp3 #passes

### Batch single dim, nout
f(x,p) = (x*p[1].+p[2]*p[3]).*[1;2]

lb =1.0
ub = 3.0
p = [2.0, 3.0, 4.0]
prob = QuadratureProblem(f,lb,ub,p)

function testf3(lb,ub,p; f=f)
    prob = QuadratureProblem(f,lb,ub,p, batch = 10, nout=2)
    sum(solve(prob, CubatureJLh(); reltol=1e-3,abstol=1e-3))
end

dp1 = ForwardDiff.gradient(p->testf3(lb,ub,p),p)
# dp2 = Zygote.gradient(p->testf3(lb,ub,p),p)[1]
dp3 = FiniteDiff.finite_difference_gradient(p->testf3(lb,ub,p),p)

@test dp1 ≈ dp3 #passes
# @test dp2 ≈ dp3 #passes

### Batch multi dim
f(x,p) = x[1,:]*p[1].+p[2]*p[3]

lb =[1.0,1.0]
ub = [3.0,3.0]
p = [2.0, 3.0, 4.0]
prob = QuadratureProblem(f,lb,ub,p)

function testf3(lb,ub,p; f=f)
    prob = QuadratureProblem(f,lb,ub,p, batch = 10, nout=1)
    solve(prob, CubatureJLh(); reltol=1e-3,abstol=1e-3)[1]
end

dp1 = ForwardDiff.gradient(p->testf3(lb,ub,p),p)
dp2 = Zygote.gradient(p->testf3(lb,ub,p),p)[1]
dp3 = FiniteDiff.finite_difference_gradient(p->testf3(lb,ub,p),p)

@test dp1 ≈ dp3
@test dp2 ≈ dp3

### Batch multi dim nout
f(x,p) = x*p[1].+p[2]*p[3]

lb =[1.0,1.0]
ub = [3.0,3.0]
p = [2.0, 3.0, 4.0]
prob = QuadratureProblem(f,lb,ub,p)

function testf3(lb,ub,p; f=f)
    prob = QuadratureProblem(f,lb,ub,p, batch = 10, nout=2)
    sum(solve(prob, CubatureJLh(); reltol=1e-3,abstol=1e-3))
end

dp1 = ForwardDiff.gradient(p->testf3(lb,ub,p),p)
# dp2 = Zygote.gradient(p->testf3(lb,ub,p),p)[1]
dp3 = FiniteDiff.finite_difference_gradient(p->testf3(lb,ub,p),p)

@test dp1 ≈ dp3
# @test dp2 ≈ dp3

## iip Batch mulit dim
function g(dx, x,p)
    dx.= sum(x*p[1].+p[2]*p[3], dims=1)
end

lb =[1.0,1.0]
ub = [3.0,3.0]
p = [2.0, 3.0, 4.0]
prob = QuadratureProblem(g,lb,ub,p)

function testf3(lb,ub,p; f=g)
    prob = QuadratureProblem(f,lb,ub,p, batch = 10, nout=1)
    solve(prob, CubatureJLh(); reltol=1e-3,abstol=1e-3)[1]
end

testf3(lb,ub,p)

dp1 = ForwardDiff.gradient(p->testf3(lb,ub,p),p)
dp2 = Zygote.gradient(p->testf3(lb,ub,p),p)[1]
dp3 = FiniteDiff.finite_difference_gradient(p->testf3(lb,ub,p),p)

@test dp1 ≈ dp3
@test dp2 ≈ dp3

## iip Batch mulit dim nout
function g(dx, x,p)
    dx.= x*p[1].+p[2]*p[3]
end

lb =[1.0,1.0]
ub = [3.0,3.0]
p = [2.0, 3.0, 4.0]
prob = QuadratureProblem(g,lb,ub,p)

function testf3(lb,ub,p; f=g)
    prob = QuadratureProblem(f,lb,ub,p, batch = 10, nout=2)
    sum(solve(prob, CubatureJLh(); reltol=1e-3,abstol=1e-3))
end

dp1 = ForwardDiff.gradient(p->testf3(lb,ub,p),p)
# dp2 = Zygote.gradient(p->testf3(lb,ub,p),p)[1]
dp3 = FiniteDiff.finite_difference_gradient(p->testf3(lb,ub,p),p)

@test dp1 ≈ dp3
# @test dp2 ≈ dp3


text:
using FFTW
using Random
using Distributions
using LsqFit

function fourierFiltering(pw, h)
    N= 2^pw
    N2 = N ÷ 2
    β = 2*h+1
    f = zeros(Complex{Float64}, N2)

    rad = [(i^(-β/2)) for i in 1:N2]
    rad = rad .* rand(Normal(0.0,1.0), N2)
    ϕ = 2π .* rand(N2)
    exp_iϕ = cos.(ϕ) .+ sin.(ϕ) .* im
    f[1:N2] = rad.*exp_iϕ
#     f[N2] = real(f[N2])
    g = irfft(f, 2*N2-1)  # irfft is nomalized by 1/N2
    zh = (g .- mean(g)).*N2^(1+h)
end

function fbmfft(pw, h)
    hurst=-1
    n=1
    z = fourierFiltering(pw,h)
    while !isapprox(hurst,h; atol=0.01) && n<10000
        z = fourierFiltering(pw,h)
        x, y = dfa(z; scal=3:0.25:(pw-2), offset=50)
        p = curve_fit((x,p)-> p[1] .+ p[2].*x, x,y,[0.0,0.0])

        hurst = p.param[2]-1
        n+=1
    end
    return hurst, z
end


text:
# Result tracking

import Base.Random: uuid4, UUID

type Result
  id::UUID
  value
  data
end

Result(id, value) = Result(id, value, Dict())

displayinline(r::Result) =
  Collapsible(span(strong("Result "), fade(string(r.id))),
              DOM.div([applydisplayinline(r.value), applydisplayinline(r.data)]))

const results = Dict{UUID,Result}()

handle("result.clear") do _, id
  delete!(results, UUID(id))
end

function register_result(result, dict...)
  id = uuid4()
  results[id] = Result(id, result, dict...)
end

# Raise on results

raise(obj::UUID, event, args...) =
  raise(global_client, :raise, @d(:id => string(obj),
                                  :event => event,
                                  :args => args))

raise(obj::Result, args...) = raise(obj.id, args...)

jscall(code) = raise(_currentresult_, :eval, code)

jsescapestring(s::AbstractString) = @> s replace("\\", "\\\\") replace("\"", "\\\"") replace("\n", "\\n") replace("'", "\\'")

# Current result

_currentresult_ = nothing

function withcurrentresult(f, r::Result)
  global _currentresult_ = r
  try
    f()
  finally
    _currentresult_ = nothing
  end
end

withcurrentresult(f, r::AbstractString) =
  withcurrentresult(f, get(results, UUID(r), nothing))

withcurrentresult(f, r) = nothing

# Evaluate Julia code

handle("eval.julia") do ed, data
  withcurrentresult(data["id"]) do
    include_string(data["code"])
  end
end

jlcall(code) =
  _currentresult_ == nothing ?
    """jlcall('$(htmlescape(code))')""" :
    """jlcall('$(_currentresult_.id)', '$(htmlescape(code))');"""


text:
container_spec =
[
    (:scene, Vector{VegaScene}, nothing),
    (:scales, Vector{VegaScale}, nothing),
    (:axes, Vector{VegaAxis}, nothing),
    (:legends, Vector{VegaLegend}, nothing),
    (:marks, Vector{VegaMark}, nothing),
]

primitivefactory(:VegaContainer, container_spec)

max_stars_repo_name: jpthor/vega


text:
using Roots, EngEconomics
n = 10

f(i) = (-3 + seriesPresentAmountFactor(i, n)) * ((4000 - 1000 * sinkingFundFactor(i, n)) / (1 - 3 * capitalRecoveryFactor(i, n))) + 1000 * presentWorthFactor(i, n) - 20000

i = find_zero(f, 0.05)

A = (4000 - 1000 * sinkingFundFactor(i, n)) / (1 - 3 * capitalRecoveryFactor(i, n))
P = 3 * A

g(x) = -P + A * seriesPresentAmountFactor(x, n) + 1000 * presentWorthFactor(x, n)
i_star = find_zero(g, 0.01)


text:
@testset "226.invert-binary-tree.jl" begin
    @test invert_tree!(TreeNode{Int}([4, 2, 7, 1, 3, 6, 9])) ==
        TreeNode{Int}([4, 7, 2, 9, 6, 3, 1])
    @test invert_tree!(TreeNode{Int}([2, 1, 3])) == TreeNode{Int}([2, 3, 1])
    @test invert_tree!(nothing) === nothing
end


text:
"Arguments are equal in distribution"
function equaldist end