text:
# Trust region method using the truncated conjugate gradient.

struct TrustRegion{T<:Real}
    η1::T
    η2::T
    γ1::T
    γ2::T
end

function trdefaults()
    return TrustRegion(0.01, 0.9, 0.5, 0.5)
end

mutable struct TrustRegionState
    k::Int64
    x::Vector{Float64}
    xcand::Vector{Float64}
    g::Vector{Float64}
    step::Vector{Float64}
    Δ::Float64
    ρ::Float64
    δ::Float64

    function TrustRegionState()
        state = new()
        state.δ = 1e-6
        return state
    end
end

function acceptcandidate!(state::TrustRegionState, tr::TrustRegion)
    if state.ρ >= tr.η1
        return true
    end
    return false
end

function updateradius!(state::TrustRegionState, tr::TrustRegion)
    if state.ρ >= tr.η2
        stepnorm = norm(state.step)
        state.Δ = min(1e20, max(4 * stepnorm, state.Δ))
    elseif state.ρ >= tr.η1
        state.Δ *= tr.γ2
    else
        state.Δ *= tr.γ1
    end
end

# Steihaug-Toint's method.
function tcg(H::Matrix, g::Vector, Δ::Float64)
    n::Int64 = length(g)
    s = zeros(n)
    normg0 = norm(g)
    v = g
    d = -v
    norm2d = gv = dot(g, v)
    norm2s::Float64 = 0.0
    sMd::Float64 = 0.0
    k::Int64 = 0
    while !stopcg(norm(g), normg0, k, n)
        Hd = H * d
        κ = dot(d, Hd)
        if κ <= 0
            σ = (-sMd + sqrt(sMd * sMd + norm2d * (Δ - dot(s, s)))) / norm2d
            s += σ * d
            break
        end
        α = gv / κ
        norm2s += α * (2 * sMd + α * norm2d)
        if norm2s >= Δ
            σ = (-sMd + sqrt(sMd * sMd + norm2d * (Δ - dot(s, s)))) / norm2d
            s += σ * d
            break
        end
        s += α * d
        g += α * Hd
        v = g
        newgv = dot(g, v)
        β = newgv / gv
        gv = newgv
        d = -v + β * d
        sMd = β * (sMd + α * norm2d)
        norm2d = gv + β * β * norm2d
        k += 1
    end
    return s
end

function stopcg(normg::Float64, normg0::Float64, k::Int64, kmax::Int64)
    χ::Float64 = 0.1
    θ::Float64 = 0.5
    if k == kmax || normg <= normg0 * min(χ, normg0^θ)
        return true
    end
    return false
end

function optimize(f::Function, g::Function, H::Function, step::Function,
                  x0::Vector{T}, approxh::Bool = false; kmax::Int64 = 5000) where {T<:Real}
    tr = trdefaults()
    state::TrustRegionState = TrustRegionState()
    state.k = 0
    state.x = x0
    state.g = g(state.x)
    state.Δ = 1.0
    state.δ *= state.δ
    n = length(state.x)
    d2fx = Matrix{Float64}(I, n, n)
    fx = f(state.x)
    if approxh
        y = zeros(n)
        gcand = zeros(n)
    else
        d2fx = H(state.x)
    end

    function model(s::Vector, g::Vector, H::Matrix)
        return dot(s, g) + 0.5 * dot(s, H * s)
    end

    while norm(state.g) > state.δ && state.k < kmax
        state.step = step(d2fx, state.g, state.Δ)
        state.xcand = state.x + state.step
        fcand = f(state.xcand)
        state.ρ = (fcand - fx) / model(state.step, state.g, d2fx)
        if approxh
            gcand = g(state.xcand)
            y = gcand - state.g
            d2fx = H(d2fx, y, state.step)
        end
        if acceptcandidate!(state, tr)
            state.x = copy(state.xcand)
            if !approxh
                state.g = g(state.x)
                d2fx = H(state.x)
            else
                state.g = copy(gcand)
            end
            fx = fcand
        end
        updateradius!(state, tr)
        state.k += 1
    end
    return state.x, state.k
end


text:
using Base.Test, DiffEqBase, SimpleTraits

f123(x,y) = 1
f123(::Val{:jac}, x) = 2
@test istrait(HasJac{typeof(f123)})

# @code_llvm istrait(HasJac{typeof(f123)})

g123(x,y) = 1
g123{T}(::Val{:jac}, x::T) = 2

h123(x,y) = 1
h123{T}(x::T) = 2

immutable T123 end
(::T123)(a) = 1
(::T123)(::Val{:jac}, a) = 1
t123 = T123()

immutable G123 end
(::G123)(a) = 1
(::G123)(b, a) = 1

@test has_jac(f123)
@test has_jac(g123)
@test !has_jac(h123)
@test has_jac(T123())
@test !has_jac(G123())

@traitfn tfn(f::::HasJac) = true
@traitfn tfn(f::::(!HasJac)) = false

@traitfn testing(f::::HasJac)   =2
@traitfn testing(f::::!(HasJac))=3.0

@inferred testing(f123)
@inferred testing(g123)
@inferred testing(h123)

@test tfn(f123)
@test tfn(g123)
@test !tfn(h123)
@test tfn(T123())
@test !tfn(G123())
@test !tfn(T123)  # NOTE this is inconsistent as has_jac(T123)==true
@test !tfn(G123)

using ParameterizedFunctions

type  LotkaVolterra <: AbstractParameterizedFunction{true}
          a::Float64
          b::Float64
end
(p::LotkaVolterra)(t,u,du) = begin
        du[1] = p.a * u[1] - p.b * u[1]*u[2]
        du[2] = -3 * u[2] + u[1]*u[2]
end
(p::LotkaVolterra)(::Val{:jac}, t, u, J) = 1

lv = LotkaVolterra(0.0,0.0)

@test has_jac(lv)
@test !has_syms(lv)
@test tfn(lv)

f = @ode_def_nohes LotkaVolterraTest begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a=>1.5 b=>1 c=3 d=1

@test has_syms(f)

ext: jl


text:
module ImageNet

import ..Metalhead
import ..Metalhead: ValidationImage, ValData, TestData, TrainData, ObjectClass

using Images

import ..valimgs, ..testimgs, ..labels

# Read ImageNet Metadata
const imagenet_val_labels_file = joinpath(@__DIR__, "..", "..", "datasets", "meta", "ILSVRC_val_labels.txt")
const imagenet_val_labels = String[strip(line) for line in eachline(imagenet_val_labels_file)]
const imagenet_labels = String[]
const synset_mapping = Dict{String, Int}()

const synset_mapping_file = joinpath(@__DIR__, "..", "..", "datasets", "meta", "ILSVRC_synset_mappings.txt")
for (idx, line) in enumerate(eachline(synset_mapping_file))
    synset = line[1:9]
    label = line[11:end]
    push!(imagenet_labels, label)
    synset_mapping[synset] = idx
end

struct ImageNet1k <: ObjectClass
    class::Int
end
labels(::Type{ImageNet1k}) = imagenet_labels

function Base.print(io::IO, class::ImageNet1k)
    print(io, imagenet_labels[class.class])
end

abstract type DataSet <: Metalhead.DataSet end

Base.size(v::ValData{<:DataSet}) = (50000,)
Base.size(v::TestData{<:DataSet}) = (100000,)

struct RawFS <: DataSet
    train_folder::String
    test_folder::String
    val_folder::String
end
function Base.show(io::IO, set::RawFS)
    println(io, "ImageNet DataSet on FileSystem")
    println(io, "      Train: $(set.train_folder)")
    println(io, "       Test: $(set.test_folder)")
    println(io, " Validation: $(set.val_folder)")
end

function RawFS(folder::String, layout=:none)
    if layout == :kaggle
        base = joinpath(folder, "Data", "CLS-LOC")
        return RawFS(
            joinpath(base, "train"),
            joinpath(base, "test"),
            joinpath(base, "val"),
        )
    else
        error("Unrecognized fs layout $(layout)")
    end
end

function Base.getindex(val::ValData{RawFS}, i::Int)
    fname = joinpath(val.set.val_folder, "ILSVRC2012_val_$(lpad(string(i), 8, '0')).JPEG")
    label = synset_mapping[imagenet_val_labels[i]][1]
    ValidationImage(DataSet, i, fname, ImageNet1k(label))
end
Base.getindex(val::TestData{RawFS}, i::Int) = joinpath(val.set.test_folder, "ILSVRC2012_test_$(lpad(string(i), 8, '0')).JPEG")

end


text:
"""
The Trigger supertype corresponds to the trigger component in DAQ.
"""
abstract type Trigger end

"""
The TrapFilter trigger is a dummy simulation of the trapezoidal filter
    algorithm to produce a trigger if the waveform amplitude passes
    a certain threshold.

This struct is currently mutable because in case of NoiseFromData modelling,
    we need to calculate the threshold post-factum after analysing the noise
    levels in data.
"""
@with_kw mutable struct TrapFilter <: Trigger
    "Lenghs of the three trap filter windows in number of samples"
    window_lengths::NTuple{3, Int} = (250,250,250)

    "Trigger threshold in keV"
    threshold_keV::typeof(1.0*energy_unit) = 0u"keV"

    "Trigger threshold after conversion"
    threshold::Real = 0

end

"""
    TrapFilter(sim_conf)

PropDict -> TrapFilter

Construct a TrapFilter instance based on simulation configuration given in <sim_conf>.
"""
function TrapFilter(sim_conf::PropDict)
    TrapFilter(
        window_lengths = (sim_conf.setup.trigger.window_lengths[1],
            sim_conf.setup.trigger.window_lengths[2], sim_conf.setup.trigger.window_lengths[3]),
        threshold_keV = T(sim_conf.setup.trigger.threshold)u"keV"
    )
end

"""
    Trigger(sim_conf)

PropDict -> <Trigger>

Construct a Trigger supertype instance based on settings given in <sim_conf>.
The returned type depends on the given settings.

Currently only TrapFilter type is implemented.

"""
function Trigger(sim_conf::PropDict)
    if sim_conf.setup.trigger.type == "trapezoidal"
        TrapFilter(sim_conf)
    else
        @info "Trigger type $(sim_config.setup.trigger.type) not implemented!\n
        Available type: trapezoidal"
    end
end

# -------------------------------------------------------------------

"""
    trap_filter(wf_values, sample_idx, trap_filter)

AbstractVector, Int, TrapFilter -> Real, Bool

Simulate the output of the trapezoidal filter algorithm
    applied to <wf_values> when the first sliding window starts
    at <sample_idx>.
The trapezoidal algorithm parameters are contained in <trap_filter.
The returned values are the online energy estimate corresponding to the
    difference between the 1st and the 3rd sliding window, and a Bool
    correspondong to whether there is a trigger issued or not.
"""
function trap_filter(wf_values::AbstractVector, sample_idx::Int, trap_filter::TrapFilter)
    # first window
    r1 = sample_idx : sample_idx + trap_filter.window_lengths[1]
    # third window
    r2 = sample_idx + trap_filter.window_lengths[1] + trap_filter.window_lengths[2] : sample_idx + sum(trap_filter.window_lengths)
    # difference between the third and first window
    r = mean(wf_values[r2]) - mean(wf_values[r1])
    r, r >= trap_filter.threshold
end

function trap_filter(wf::RDWaveform, sample_idx::Int, trigger::TrapFilter)
    trap_filter(wf.value, sample_idx, trigger)
end

"""
    simulate(wf, trap_filter)

RDWaveform, TrapFilter -> Int, Real

Simulate the results of applying <trap_filter> to <wf>
Returns the index on which <wf> triggered, and the estimated
    online energy based on <trap_filter> output.
"""
function simulate(wf::RDWaveform, trigger::TrapFilter)
    trigger_window_length = sum(trigger.window_lengths)

    online_filter_output = zeros(T, length(wf.value) - trigger_window_length)
    t0_idx::Int = 0
    trig = false

    # replace by a while loop?
    for i in eachindex(online_filter_output)
        online_filter_output[i], trig = trap_filter(wf, i, trigger)
        if trig && t0_idx == 0
            t0_idx = i + trigger.window_lengths[1] + trigger.window_lengths[2]
        end
    end

    t0_idx, maximum(online_filter_output)
end


text:
function oclformdetailslist(md)
    fl = findelements(md, :FormDef)
    df = DataFrame(OID = String[], Name = String[], Repeating = String[],
    ParentFormOID = String[], SectionLabel = String[], SectionTitle = String[], SectionSubtitle = String[], SectionInstructions = String[], SectionPageNumber = String[])
    for i in fl
        el = findelement(i, :FormDetails)
        if !(isnothing(el))
            pfoid  = attribute(el, "ParentFormOID")
            secdet = findelement(el, :SectionDetails)
            sec    = findelements(secdet, :Section)
            for s in sec
                seclab = attribute(s, "SectionLabel")
                sectit = attribute(s, "SectionTitle")
                secsubt = attribute(s, "SectionSubtitle")
                secinst = attribute(s, "SectionInstructions")
                secpn = attribute(s, "SectionPageNumber")
                push!(df, (attribute(i, "OID"), attribute(i, "Name"), attribute(i, "Repeating"), pfoid, seclab, sectit, secsubt, secinst, secpn))
            end
        else
            pfoid = ""
            seclab = ""
            sectit = ""
            secsubt = ""
            secinst = ""
            secpn = ""
            push!(df, (attribute(i, "OID"), attribute(i, "Name"), attribute(i, "Repeating"), pfoid, seclab, sectit, secsubt, secinst, secpn))
        end
    end
    df
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "LinearTurboFold"
# invented version number because there aren't any releases yet and we
# are installing from a git commit
version = v"0.0.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/LinearFold/LinearTurboFold/archive/3e58e713e2743ec57aad21d557289c43571e97f0.tar.gz",
                  "58b02b83cb71675af92a14b18a0bf64a1171e554cfaa5e914038d97c6589fd30"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/LinearTurboFold-*/

make -j${nproc} CC=${CXX}
for b in bin/*; do
    install -Dvm 755 "${b}" "${bindir}/$(basename "${b}")${exeext}"
done

# needed parameter files
install -Dvm 755 src/data_tables/fam_hmm_pars.dat "${prefix}/share/LinearTurboFold/src/data_tables/fam_hmm_pars.dat"

install_license LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("linearturbofold", :linearturbofold),
    FileProduct("share/LinearTurboFold/src/data_tables/fam_hmm_pars.dat", :fam_hmm_pars),
]

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

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


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)
module Divergences

using StatsFuns
import Distances: evaluate, gradient, PreMetric
import Calculus: hessian

abstract Divergence <: PreMetric

immutable CressieRead <: Divergence
    α::Float64
    function CressieRead(α::Float64)
        @assert isempty(findin(α, [-1, 0])) "CressieRead is defined for all α != {-1,0}"
        new(α)
    end
end

CressieRead(α::Int64) = CressieRead(float(α))

immutable ChiSquared  <: Divergence end
immutable KullbackLeibler  <: Divergence end

immutable ReverseKullbackLeibler <: Divergence end

immutable ModifiedKullbackLeibler <: Divergence
    ϑ::Float64
    d::Divergence
    m::NTuple{4, Float64}
    function ModifiedKullbackLeibler(ϑ::Real)
        @assert ϑ > 0 "ModifiedKullbackLeibler is defined for ϑ > 0"
        uϑ = 1.0 + ϑ
        d  = KullbackLeibler()
        f0 = evaluate(d, [uϑ])
        f1 = gradient(d, uϑ)
        f2 = hessian(d, uϑ)
        new(float(ϑ), d, (f0, f1, f2, uϑ))
    end
end

immutable FullyModifiedKullbackLeibler <: Divergence
    φ::Float64
    ϑ::Float64
    d::Divergence
    m::NTuple{8, Float64}
    function FullyModifiedKullbackLeibler(φ::Real, ϑ::Real)
        @assert ϑ > 0 "FullyModifiedKullbackLeibler is defined for ϑ > 0"
        @assert φ > 0 && φ < 1.0 "FullyModifiedKullbackLeibler is defined for ϕ ∈ (0,1)"
        uϑ = 1.0 + ϑ
        d  = KullbackLeibler()
        f0 = evaluate(d, [uϑ])
        f1 = gradient(d, uϑ)
        f2 = hessian(d, uϑ)
        uφ  = float(φ)
        g0  = evaluate(d, [uφ])
        g1  = gradient(d, uφ)
        g2  = hessian(d, uφ)
        new(float(φ), float(ϑ), d, (f0, f1, f2, uϑ, g0, g1, g2, uφ))
    end
end

immutable ModifiedReverseKullbackLeibler <: Divergence
    ϑ::Float64
    d::Divergence
    m::NTuple{4, Float64}
    function ModifiedReverseKullbackLeibler(ϑ::Real)
        @assert ϑ > 0 "ModifiedReverseKullbackLeibler is defined for ϑ > 0"
        uϑ = 1.0 + ϑ
        d  = ReverseKullbackLeibler()
        f0 = evaluate(d, [uϑ])
        f1 = gradient(d, uϑ)
        f2 = hessian(d, uϑ)
        new(float(ϑ), d, (f0, f1, f2, uϑ))
    end
end

immutable FullyModifiedReverseKullbackLeibler <: Divergence
    φ::Float64
    ϑ::Float64
    d::Divergence
    m::NTuple{8, Float64}
    function FullyModifiedReverseKullbackLeibler(φ::Real, ϑ::Real)
        @assert ϑ > 0 "ModifiedReverseKullbackLeibler is defined for ϑ > 0"
        @assert φ > 0 && φ < 1.0 "ModifiedReverseKullbackLeibler is defined for φ ∈ (0,1)"
        d   = ReverseKullbackLeibler()
        uϑ  = 1.0 + ϑ
        f0  = evaluate(d, [uϑ])
        f1  = gradient(d, uϑ)
        f2  = hessian(d, uϑ)
        uφ  = float(φ)
        g0  = evaluate(d, [uφ])
        g1  = gradient(d, uφ)
        g2  = hessian(d, uφ)
        new(float(φ), float(ϑ), d, (f0, f1, f2, uϑ, g0, g1, g2, uφ))
    end
end

immutable ModifiedCressieRead <: Divergence
    α::Float64
    ϑ::Float64
    d::Divergence
    m::NTuple{4, Float64}
    function ModifiedCressieRead(α::Real, ϑ::Real)
        @assert isempty(findin(α, [-1, 0])) "ModifiedCressieRead is defined for all α! = {-1,0}."
        @assert ϑ > 0 "ModifiedCressieRead is defined for ϑ > 0"
        uϑ = 1.0 + ϑ
        d  = CressieRead(α)
        f0 = evaluate(d, [uϑ])
        f1 = gradient(d, uϑ)
        f2 = hessian(d, uϑ)
        new(float(α), float(ϑ), d, (f0, f1, f2, uϑ))
    end
end

immutable FullyModifiedCressieRead <: Divergence
    α::Float64
    φ::Float64
    ϑ::Float64
    d::Divergence
    m::NTuple{8, Float64}
    function FullyModifiedCressieRead(α::Real, φ::Real, ϑ::Real)
        @assert isempty(findin(α, [-1, 0])) "ModifiedCressieRead is defined for all α != {-1,0}"
        @assert ϑ > 0 "FullyModifiedCressieRead is defined for ϑ > 0"
        @assert φ > 0 && φ < 1.0 "FullyModifiedCressieRead is defined for φ ∈ (0, 1)"
        uϑ = 1.0 + ϑ
        d  = CressieRead(α)
        f0 = evaluate(d, [uϑ])
        f1 = gradient(d, uϑ)
        f2 = hessian(d, uϑ)
        uφ = float(φ)
        g0 = evaluate(d, [uφ])
        g1 = gradient(d, uφ)
        g2 = hessian(d, uφ)
        new(float(α), float(φ), float(ϑ),  d, (f0, f1, f2, uϑ, g0, g1, g2, uφ))
    end
end

typealias KL KullbackLeibler
typealias MKL ModifiedKullbackLeibler
typealias FMKL FullyModifiedKullbackLeibler

typealias RKL ReverseKullbackLeibler
typealias MRKL ModifiedReverseKullbackLeibler
typealias FMRKL FullyModifiedReverseKullbackLeibler

typealias CR CressieRead
typealias MCR ModifiedCressieRead
typealias FMCR FullyModifiedCressieRead

HD()  = CressieRead(-1/2)

include("common.jl")
include("cressieread.jl")
include("modified_cressieread.jl")
include("kl.jl")
include("reversekl.jl")
include("chisq.jl")

export
    Divergence,
    # KL
    KullbackLeibler,
    ModifiedKullbackLeibler,
    FullyModifiedKullbackLeibler,
    # RKL
    ReverseKullbackLeibler,
    ModifiedReverseKullbackLeibler,
    FullyModifiedReverseKullbackLeibler,
    # CR
    CressieRead,
    ModifiedCressieRead,
    FullyModifiedCressieRead,
    # Abbr.
    KL,
    MKL,
    FMKL,
    RKL,
    MRKL,
    FMRKL,
    CR,
    MCR,
    FMCR,
    HD,
    ChiSquared,
    evaluate,
    gradient!,
    hessian!,
    gradient,
    hessian

end # module


text:
############################################################################################
# Import External PackagesJK
using Test
using Random: Random, AbstractRNG, seed!
using ArgCheck: ArgCheck, @argcheck
using UnPack: UnPack, @unpack
using Distributions

############################################################################################
# Import Baytes Packages
using BaytesCore, ModelWrappers, BaytesFilters

############################################################################################
# Include Files
include("testhelper/TestHelper.jl");

############################################################################################
# Run Tests
@testset "All tests" begin
    # General
    include("test-memory.jl")
    # PF specific
    include("particlefilter/test-construction.jl")
    include("particlefilter/test-estimate.jl")
    include("particlefilter/test-types.jl")
    include("particlefilter/test-ancestors.jl")
    include("particlefilter/test-Ndata.jl")
end


text:
# # Sensitivity Analysis of SVM

#md # [![](https://img.shields.io/badge/show-github-579ACA.svg)](@__REPO_ROOT_URL__/docs/src/examples/sensitivity-analysis-svm.jl)

# This notebook illustrates sensitivity analysis of data points in a [Support Vector Machine](https://en.wikipedia.org/wiki/Support-vector_machine) (inspired from [@matbesancon](http://github.com/matbesancon)'s [SimpleSVMs](http://github.com/matbesancon/SimpleSVMs.jl).)

# For reference, Section 10.1 of https://online.stat.psu.edu/stat508/book/export/html/792 gives an intuitive explanation of what it means to have a sensitive hyperplane or data point. The general form of the SVM training problem is given below (with $\ell_2$ regularization):

# ```math
# \begin{split}
# \begin{array} {ll}
# \mbox{minimize} & \lambda||w||^2 + \sum_{i=1}^{N} \xi_{i} \\
# \mbox{s.t.} & \xi_{i} \ge 0 \quad \quad i=1..N  \\
#             & y_{i} (w^T X_{i} + b) \ge 1 - \xi_{i} \quad i=1..N  \\
# \end{array}
# \end{split}
# ```
# where
# - `X`, `y` are the `N` data points
# - `w` is the support vector
# - `b` determines the offset `b/||w||` of the hyperplane with normal `w`
# - `ξ` is the soft-margin loss
# - `λ` is the $\ell_2$ regularization.
#
# This tutorial uses the following packages

using JuMP     # The mathematical programming modelling language
import DiffOpt # JuMP extension for differentiable optimization
import Ipopt   # Optimization solver that handles quadratic programs
import Plots   # Graphing tool
import LinearAlgebra: dot, norm
import Random

# ## Define and solve the SVM

# Construct two clusters of data points.

N = 100
D = 2

Random.seed!(62)
X = vcat(randn(N ÷ 2, D), randn(N ÷ 2, D) .+ [2.0, 2.0]')
y = append!(ones(N ÷ 2), -ones(N ÷ 2))
λ = 0.05;

# Let's initialize a special model that can understand sensitivities

model = Model(() -> DiffOpt.diff_optimizer(Ipopt.Optimizer))
MOI.set(model, MOI.Silent(), true)

# Add the variables

@variable(model, ξ[1:N])
@variable(model, w[1:D])
@variable(model, b);

# Add the constraints.

@constraint(
    model,
    [i in 1:N],
    ξ[i] >= 0
);
@constraint(
    model,
    cons[i in 1:N],
    y[i] * (dot(X[i, :], w) + b) >= 1 - ξ[i]
);

# Define the objective and solve

@objective(
    model,
    Min,
    λ * dot(w, w) + sum(ξ),
)

optimize!(model)

# We can visualize the separating hyperplane.

loss = objective_value(model)

wv = value.(w)

bv = value(b)

svm_x = [-2.0, 4.0] # arbitrary points
svm_y = (-bv .- wv[1] * svm_x )/wv[2]

p = Plots.scatter(X[:,1], X[:,2], color = [yi > 0 ? :red : :blue for yi in y], label = "")
Plots.plot!(p, svm_x, svm_y, label = "loss = $(round(loss, digits=2))", width=3)

# ## Gradient of hyperplane wrt the data point coordinates

# Now that we've solved the SVM, we can compute the sensitivity of optimal
# values -- the separating hyperplane in our case -- with respect to
# perturbations of the problem data -- the data points -- using DiffOpt.

# How does a change in coordinates of the data points, `X`,
# affects the position of the hyperplane?
# This is achieved by finding gradients of `w` and `b` with respect to `X[i]`.

# Begin differentiating the model.
# analogous to varying θ in the expression:
# ```math
# y_{i} (w^T (X_{i} + \theta) + b) \ge 1 - \xi_{i}
# ```
∇ = zeros(N)
for i in 1:N
    for j in 1:N
        if i == j
            ## we consider identical perturbations on all x_i coordinates
            MOI.set(
                model,
                DiffOpt.ForwardInConstraint(),
                cons[j],
                y[j] * sum(w),
            )
        else
            MOI.set(model, DiffOpt.ForwardInConstraint(), cons[j], 0.0 * sum(w))
        end
    end
    DiffOpt.forward(model)
    dw = MOI.get.(
        model,
        DiffOpt.ForwardOutVariablePrimal(),
        w,
    )
    db = MOI.get(
        model,
        DiffOpt.ForwardOutVariablePrimal(),
        b,
    )
    ∇[i] = norm(dw) + norm(db)
end

# We can visualize the separating hyperplane sensitivity with respect to the data points.
# Note that all the small numbers were converted into 1/10 of the
# largest value to show all the points of the set.

p3 = Plots.scatter(
    X[:,1], X[:,2],
    color = [yi > 0 ? :red : :blue for yi in y], label = "",
    markersize = 2 * (max.(1.8∇, 0.2 * maximum(∇))),
)
Plots.yaxis!(p3, (-2, 4.5))
Plots.plot!(p3, svm_x, svm_y, label = "", width=3)
Plots.title!("Sensitivity of the separator to data point variations")


text:
import ReverseDiff
using ReverseDiff: SpecialInstruction, istracked, increment_deriv!, deriv, TrackedArray
using ReverseDiff: track, InstructionTape, TrackedReal, seed!, unseed!, reverse_pass!, record!
import ReverseDiff: special_reverse_exec!

mutable struct GFBackpropTraceState
    trace::GFTrace
    score::TrackedReal
    tape::InstructionTape
    visitor::AddressVisitor
    params::Dict{Symbol,Any}
    selection::AddressSet
    tracked_choices::HomogenousTrie{Any,TrackedReal}
    value_trie::DynamicAssignment
    gradient_trie::DynamicAssignment
end

function GFBackpropTraceState(trace, selection, params, tape)
    score = track(0., tape)
    visitor = AddressVisitor()
    tracked_choices = HomogenousTrie{Any,TrackedReal}()
    value_trie = DynamicAssignment()
    gradient_trie = DynamicAssignment()
    GFBackpropTraceState(trace, score, tape, visitor, params,
                       selection, tracked_choices, value_trie, gradient_trie)
end

get_args_change(state::GFBackpropTraceState) = nothing
get_addr_change(state::GFBackpropTraceState, addr) = nothing
set_ret_change!(state::GFBackpropTraceState, value) = begin end

function fill_gradient_trie!(gradient_trie::DynamicAssignment,
                             tracked_trie::HomogenousTrie{Any,TrackedReal})
    for (key, tracked) in get_leaf_nodes(tracked_trie)
        set_leaf_node!(gradient_trie, key, deriv(tracked))
    end
    # NOTE: there should be no address collision between these primitive
    # choices and the generator invocations, as enforced by the visitor
    for (key, node) in get_internal_nodes(tracked_trie)
        @assert !has_leaf_node(gradient_trie, key) && !has_internal_node(gradient_trie, key)
        gradient_trie_node = DynamicAssignment()
        fill_gradient_trie!(gradient_trie_node, node)
        set_internal_node!(gradient_trie, key, gradient_trie_node)
    end
end

function fill_value_trie!(value_trie::DynamicAssignment,
                          tracked_trie::HomogenousTrie{Any,TrackedReal})
    for (key, tracked) in get_leaf_nodes(tracked_trie)
        set_leaf_node!(value_trie, key, ReverseDiff.value(tracked))
    end
    # NOTE: there should be no address collision between these primitive
    # choices and the generator invocations, as enforced by the visitor
    for (key, node) in get_internal_nodes(tracked_trie)
        @assert !has_leaf_node(value_trie, key) && !has_internal_node(value_trie, key)
        value_trie_node = DynamicAssignment()
        fill_value_trie!(value_trie_node, node)
        set_internal_node!(value_trie, key, value_trie_node)
    end
end

function addr(state::GFBackpropTraceState, dist::Distribution{T}, args, addr) where {T}
    visit!(state.visitor, addr)
    call::CallRecord = get_primitive_call(state.trace, addr)
    retval::T = call.retval
    if has_internal_node(state.selection, addr)
        error("Got internal node but expected leaf node in selection at $addr")
    end
    if has_leaf_node(state.selection, addr)
        tracked_retval = track(retval, state.tape)
        set_leaf_node!(state.tracked_choices, addr, tracked_retval)
        score_tracked = logpdf(dist, tracked_retval, args...)
        state.score += score_tracked
        return tracked_retval
    else
        state.score += logpdf(dist, retval, args...)
        return retval
    end
end

struct BackpropTraceRecord
    generator::Generator
    subtrace::Any
    selection::AddressSet
    value_trie::DynamicAssignment
    gradient_trie::DynamicAssignment
    addr::Any
end

function addr(state::GFBackpropTraceState, gen::Generator{T}, args, addr, args_change) where {T}
    visit!(state.visitor, addr)
    if has_leaf_node(state.selection, addr)
        error("Cannot select a whole subtrace, tried to select $addr")
    end
    trace = get_subtrace(state.trace, addr)
    call::CallRecord = get_call_record(trace) # use the return value recorded in the trace
    retval::T = call.retval
    if accepts_output_grad(gen)
        retval_maybe_tracked = track(retval, state.tape)
        if !istracked(retval_maybe_tracked)
            error("Could not track return value at address $addr on AD tape.")
        end
    else
        retval_maybe_tracked = retval
    end
    # some of the args may be tracked (see special_reverse_exec!)
    # note: we still need to run backprop_params on gen, even if it does not
    # accept an output gradient, because it may make random choices.
    if has_internal_node(state.selection, addr)
        selection = get_internal_node(state.selection, addr)
    else
        selection = EmptyAddressSet()
    end
    record = BackpropTraceRecord(gen, trace, selection, state.value_trie,
                                 state.gradient_trie, addr)
    record!(state.tape, SpecialInstruction, record, (args...,), retval_maybe_tracked)
    retval_maybe_tracked
end

function backprop_trace(gf::GenFunction, trace::GFTrace, selection::AddressSet, retval_grad)
    tape = InstructionTape()
    state = GFBackpropTraceState(trace, selection, gf.params, tape)
    call = get_call_record(trace)
    args = call.args
    args_maybe_tracked = (map(maybe_track, args, gf.has_argument_grads, fill(tape, length(args)))...,)
    retval_maybe_tracked = exec(gf, state, args_maybe_tracked)
    if gf.accepts_output_grad
        if !istracked(retval_maybe_tracked)
            error("Output of $gf accepts gradient but return value was not tracked on AD tape.")
        end
        if retval_grad != nothing
            # it is not an error if retval_grad is nothing, it means the gradient is zero?
            deriv!(retval_maybe_tracked, retval_grad)
        end
    end
    seed!(state.score)
    reverse_pass!(tape)

    # fill trace gradient with gradients with respect to primitive random choices
    fill_gradient_trie!(state.gradient_trie, state.tracked_choices)
    fill_value_trie!(state.value_trie, state.tracked_choices)

    # return gradients with respect to inputs
    # NOTE: if a value isn't tracked the gradient is nothing
    input_grads::Tuple = (map((arg, has_grad) -> has_grad ? deriv(arg) : nothing,
                             args_maybe_tracked, gf.has_argument_grads)...,)

    (input_grads, state.value_trie, state.gradient_trie)
end

@noinline function special_reverse_exec!(instruction::SpecialInstruction{BackpropTraceRecord})
    record = instruction.func
    gen = record.generator
    args_maybe_tracked = instruction.input
    retval_maybe_tracked = instruction.output
    if accepts_output_grad(gen)
        @assert istracked(retval_maybe_tracked)
        retval_grad = deriv(retval_maybe_tracked)
    else
        retval_grad = nothing
    end
    (arg_grads, value_trie, grad_trie) = backprop_trace(
        gen, record.subtrace, record.selection, retval_grad)

    @assert !has_internal_node(record.gradient_trie, record.addr)
    @assert !has_leaf_node(record.gradient_trie, record.addr)
    set_internal_node!(record.gradient_trie, record.addr, grad_trie)
    set_internal_node!(record.value_trie, record.addr, value_trie)

    for (arg, grad, has_grad) in zip(args_maybe_tracked, arg_grads, has_argument_grads(gen))
        if has_grad && istracked(arg)
            increment_deriv!(arg, grad)
        end
    end
    nothing
end

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
function get_non_dominated_solutions_perm(population)
    ids = Int[1]
    n = length(population)

    for i in 2:n
        j = 1
        while j <= length(ids)
            jj = ids[j]
            relation = compare(population[i], population[jj])

            if relation == 2 # j dominates i
                break
            elseif relation == 1 # i dominates j
                deleteat!(ids, j)
                continue
            end

            j += 1
        end

        if j > length(ids)
            push!(ids, i)
        end

    end

    return ids
end

function get_non_dominated_solutions(population)

    mask = get_non_dominated_solutions_perm(population)
    return population[mask]

end

"""
    gen_ref_dirs(dimension, n_paritions)
Generates Das and Dennis's structured reference points. `dimension` could be
the number of objective functions in multi-objective functions.
"""
function gen_ref_dirs(dimension, n_paritions)
    return  gen_weights(dimension, n_paritions)
end

function gen_weights(a, b)
    nobj = a;
    H    = b;
    a    = zeros(nobj);
    d    = H;
    w    = [];
    produce_weight!(a, 1, d, H, nobj, w)
    return Array.(w)
end

function  produce_weight!(a, i, d, H, nobj, w)
    for k=0:d
        if i<nobj
            a[i] = k;
            d2   = d - k;
            produce_weight!(a, i+1, d2, H, nobj, w);
        else
            a[i] = d;
            push!(w, a/H)
            break;
        end
    end
end

"""
    ideal(points)

Computes the ideal point from a provided array of `Vector`s or a population or row vectors
in a `Matrix`.
"""
function ideal(points::Array{Vector{T}}) where T <: Real

    isempty(points) && isempty(points[1]) && return zeros(0)

    ideal = points[1]

    for point in points
        ideal = min.(ideal, point)
    end

    return ideal

end

"""
    nadir(points)

Computes the nadir point from a provided array of `Vector`s or a population or row vectors
in a `Matrix`.
"""
function nadir(points::Array{Vector{T}})  where T <: Real

    isempty(points) && isempty(points[1]) && return zeros(0)

    nadir = points[1]

    for point in points
        nadir = max.(nadir, point)
    end

    return nadir
end

function ideal(population::Array{xFgh_indiv})
    mask = sum_violations.(population) .== 0

    ideal(fval.(population[mask]))
end
ideal(A::Matrix) = ideal([A[i,:]  for i in 1:size(A,1)])

function nadir(population::Array{xFgh_indiv})
    mask = sum_violations.(population) .== 0

    nadir(fval.(population[mask]))
end

nadir(A::Matrix) = nadir([A[i,:]  for i in 1:size(A,1)])


text:
"""
    ngramize(lang, tokens, n)

Compute the ngrams of `tokens` of the order `n`.

# Example

```julia-repl
julia> ngramize(Languages.English(), ["To", "be", "or", "not", "to"], 3)
Dict{AbstractString,Int64} with 3 entries:
  "be or not" => 1
  "or not to" => 1
  "To be or"  => 1
```
"""
function ngramize(lang::S, words::Vector{T}, nlist::Integer...) where {S <: Language, T <: AbstractString}
    (length(nlist) == 1) && (first(nlist) == 1) && return onegramize(lang, words)

    n_words = length(words)

    tokens = Dict{AbstractString, Int}()

    for n in nlist
        for index in 1:(n_words - n + 1)
            token = join(words[index:(index + n - 1)], " ")
            tokens[token] = get(tokens, token, 0) + 1
        end
    end
    return tokens
end

ngramize(lang::Language, str::AbstractString, n::Integer) = ngramize(lang, tokenize(lang, str), n)

"""
    onegramize(lang, tokens)

Create the unigrams dict for input tokens.

# Example

```julia-repl
julia> onegramize(Languages.English(), ["To", "be", "or", "not", "to", "be"])
Dict{String,Int64} with 5 entries:
  "or"  => 1
  "not" => 1
  "to"  => 1
  "To"  => 1
  "be"  => 2
```
"""
function onegramize(lang::S, words::Vector{T}) where {S <: Language, T <: AbstractString}
    n_words = length(words)
    tokens = Dict{T, Int}()

    for word in words
        tokens[word] = get(tokens, word, 0) + 1
    end

    return tokens
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libunistring_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libunistring")
JLLWrappers.@generate_main_file("libunistring", UUID("6db05002-db9d-53dd-a359-17d4854bdf22"))
end  # module libunistring_jll


text:
# Scenario-dependent attributes #
# ========================== #
struct ScenarioDependentOptimizerAttribute <: MOI.AbstractOptimizerAttribute
    stage::Int
    scenario_index::Int
    attr::MOI.AbstractOptimizerAttribute
end

struct ScenarioDependentModelAttribute <: MOI.AbstractModelAttribute
    stage::Int
    scenario_index::Int
    attr::MOI.AbstractModelAttribute
end

struct ScenarioDependentVariableAttribute <: MOI.AbstractVariableAttribute
    stage::Int
    scenario_index::Int
    attr::MOI.AbstractVariableAttribute
end

struct ScenarioDependentConstraintAttribute <: MOI.AbstractConstraintAttribute
    stage::Int
    scenario_index::Int
    attr::MOI.AbstractConstraintAttribute
end

const AnyScenarioDependentAttribute = Union{ScenarioDependentOptimizerAttribute,
                                            ScenarioDependentModelAttribute,
                                            ScenarioDependentVariableAttribute,
                                            ScenarioDependentConstraintAttribute}

MOI.is_set_by_optimize(attr::AnyScenarioDependentAttribute) = MOI.is_set_by_optimize(attr.attr)

is_structure_independent(attr::MOI.AnyAttribute) = true
is_structure_independent(::Union{MOI.ObjectiveFunction,
                                 MOI.ObjectiveFunctionType,
                                 MOI.ObjectiveSense}) = false
is_structure_independent(attr::AnyScenarioDependentAttribute) = is_structure_independent(attr.attr)


text:
@interface a(x::I)::I where {I}
@interface b(f::Vector{I}, x::O) where {I, O}
@interface c(f::Vector{O}) where {O}
@interface d(f)

@impl begin
    function a(x::Int)::Int
        c = 1
        return x + c
    end
end

@impl begin
    struct MyB
        precision::AbstractFloat = 0.1
        dummy::Int = 2
    end

    function b(x::Vector{I}, y::Number) where I
        return x.+y.+precision/2
    end
end

@impl begin
    function b(x::Vector{I}, y::AbstractFloat) where I
        return x.+y
    end
end

@testset "Basic" begin
    @test a(0) == 1
    @test b([1.,2.,3.],2.) == [3.,4.,5.]
    @test b([1.,2.,3.],2.) == [3.,4.,5.]
    @test all(abs.(b(MyB(0.1, 2),[1.,2.,3.], 2.) .- [3.,4.,5.]) .< 0.1)
    println(methods(a))
    @test list_impls(b) == [MyB, B1]
    @test list_impls(b, Tuple{Vector, Int}) == [MyB]
    @test all(b(MyB(0.1, 2),[1.,2.,3.], 2.) == MyB(0.1,2)([1.,2.,3.],2.))
    println(get_method(MyB))
end

module A_
    using MultiInterface

    export four

    @interface one(a::Vector{<:T}, b::T) where {T}
    @interface two(a::Signed, b::AbstractFloat)
    @interface three(a, b::Int64)
    @interface four(a::Vector{<:T})::T where {T}
end

module B_
    using MultiInterface
    import ..A_
    import ..A_: one, two

    @impl begin
        struct One_ end
        function one(a::Vector{<:AbstractFloat}, b::AbstractFloat)
            1.0
        end
        function one(a::Vector{<:Integer}, b::Integer)
            1
        end
        function one(a::Vector{Float32}, b::AbstractFloat)
            2.0
        end

    end

    @impl begin
        struct One_NoPolicy end
        function one(a::Vector{Symbol}, b::Symbol)
            :a
        end
    end

    @impl begin
        struct Two_ end

        function two(a::Int64, b::Float32)
            (typeof(a), typeof(b))
        end

        function two(a::Int32, b::Float32)
            (typeof(a), typeof(b))
        end

        function two(a::Int64, b::Float64)
            (typeof(a), typeof(b))
        end

        function two(::Signed, ::AbstractFloat)
            (Signed, AbstractFloat)
        end
    end

    @impl begin
        struct Three_ end

        function A_.three(a, b::Int64)
            b
        end
    end

    error = nothing

    try
        @impl begin
            struct Two_Error end

            function two(a::Symbol, b::Float64)
                (typeof(a), typeof(b))
            end
        end
    catch e
        global error = e
    end

    @impl begin
        struct Four_ end

        function A_.four(a::Vector{Float64})::Float64
            isempty(a) ? 0.0 : a[1]
        end
    end
end

using .A_
import .B_

import MultiInterface: get_imp

struct Policy_One <: Policy end
get_imp(policy::Policy_One, ::Type{A_.One}, args...) = B_.One_()
get_imp(policy::Policy_One, ::Type{A_.Two}, args...) = B_.Two_()

@testset "module testing" begin
    @test A_.one(Float64[1.0], 1.0) == 1.0
    @test A_.one(Int[1], 1) == 1
    @test A_.one([:a], :a) == :a

    @test A_.two(-1, 1.0) == (Int64, Float64)
    @test A_.two(1, 1.0f0) == (Int64, Float32)
    @test A_.two(convert(Int32, 1), 1.0f0) == (Int32, Float32)
    @test A_.two(convert(Int32, 1), 1.0) == (Signed, AbstractFloat)

    @test A_.three(:a, 1) == 1
    @test four([1.0,2.0]) == 1.0
    @test four(Float64[]) == 0.0

    # TODO these no longer error
    # @test B_.error isa ErrorException
    # @test startswith(B_.error.msg,"Invalid definition of two")

    with_policy(Policy_One()) do
        @test A_.one(Float64[1.0], 1.0) == 1.0
        @test A_.one(Float32[1.0], 1.0) == 2.0
        @test A_.one(Int[1], 1) == 1
        @test_throws MethodError A_.one([:a], :a)

        @test A_.two(-1, 1.0) == (Int64, Float64)
        @test A_.two(1, 1.0f0) == (Int64, Float32)
        @test A_.two(convert(Int32, 1), 1.0f0) == (Int32, Float32)
        @test A_.two(convert(Int32, 1), 1.0) == (Signed, AbstractFloat)
    end
end


text:
#3rd party
import Base.convert
using MAT
using ForwardDiff
using Optim
using Polynomials
using JLD
using PyPlot

# path_functions="/home/genis/model_fitting_julia/scripts/functions/"
path_functions="/home/genis/fit_behaviour_julia/scripts/functions/"

include(path_functions*"parameter_functions.jl")
include(path_functions*"functions_julia.jl")
include(path_functions*"utils.jl")

model_num="6"

Nstim=50
# path_data="/home/genis/model_fitting_julia/synthetic_data/"
path_data="/home/genis/fit_behaviour_julia/synthetic_data/"
# path="/home/genis/model_fitting_julia/"
# filename="synthetic_data_DW_Nstim"*string(Nstim)*"_5sigmas_test_largeParam.jld"
# filename="synthetic_data_DW_Nstim"*string(Nstim)*"_5sigmas_urgency_SLogisticTransform.jld"
filename="synthetic_data_DW_Nstim"*string(Nstim)*"_model"*model_num*"_traces.jld"

data=load(path_data*filename)

# save_dict=Dict(["d"=>d,"stim"=>stim,"param"=>param,"sigma_a"=>sigma_a,"Tframe"=>T,"PR"=>PR_sim,"kernel"=>kernel_sigma, "PR_sim_test"=>PR_sim_test, "stim_test"=>stim_test,"d_test"=>d_test])
# save(path_save*filename,save_dict)
stim=data["stim"]
x_traces_sim=data["x_trace"]
Tframe=data["Tframe"]
param=data["param"]
Nstim,Ntraces,NT=size(x_traces_sim)
model_num="6"
# args=[ ("c2",2),"c4",("st",4),("x0",2),"bias_d","hbias_d","sigma_a"]
# #   c21 c22  c4  st1  st2 st3  st4   x0  bias_d sigma
# ub=[10., 10.,10., 10.,10., 10., 10.,10, 10, 10.,10,10.]
# #   c21  c22  c4    st1  st2  st3  st4    x0  bias_d  sigma
# lb=[-10.,-10,-10 ,0.,  0., -10, -10.,-10,-10,-10,-10., 0]
lambda=0
reg=[1,1,1,1,1,1,1,1,1,1,1,1]
model_f=Dict(:c2_f=>c2_urgency,:c4_f=>c4_const,:c6_f=>c6_zeros,:bias_d_f=>bias_d_const,
:hbias_d_f=>history_bias_d_const,:st_f=>st_logistic_bias2,:x0_f=>initial_condition_bias_hbias,
:ro_f=>read_out_perfect,:lambda=>lambda,:reg=>reg)

# param=Dict("c2"=>[1,0],"c4"=>[0.5],"x0"=>[0.1,0.05], "st"=>[1.,3,0.1,0.2],"bias_d"=>[0.1],
# "hbias_d"=>[0.2],"sigma_a"=>0.3)

fig,axes=subplots(5,2)
i_stim=[1,2,3,4,5]
T=Tframe*length(stim[1,:])
t=0:T/(length(x_traces_sim[1,1,:])-1):T
for i in 1:length(i_stim)
    figure()
    x,bins=pdf_x_time_1stim(stim[i_stim[i],:],[0.],Tframe,param,model_f)
    imshow(transpose(x),origin="lower",extent=[0,T,bins[1],bins[end]],aspect="auto",vmax=0.1,cmap="hot")

    for itrial in 1:Ntraces
        plot(t,x_traces_sim[i_stim[i],itrial,:],"-",color="darkcyan")
    end

    figure()
    plot(stim[i_stim[i],:])

end
show()

max_stars_repo_name: emmafdz/FinalDegreeProject

max_issues_repo_name: emmafdz/FinalDegreeProject

max_forks_repo_name: emmafdz/FinalDegreeProject


text:
module IntsWithInfinity

import Base: hash, string, show, convert, promote_rule,
    zero, one, signbit, copysign, flipsign, abs,
    typemin, typemax, widen,
    iseven, isodd, isinf, isfinite,
    round,trunc,
    isequal, isless, (==),(!=),(<),(<=),(>=),(>),
    (~),(+),(-),(*),(/),(//),(\),(%),(^),
    widemul,div,fld,cld,rem,mod,divrem,fldmod,
    gcd, lcm, gcdx, isqrt, factorial, binomial

export InfInt8, InfInt16, InfInt32, InfInt64, InfInt128,
    iszero,nonzero,isone,isneg,nonneg,ispos,nonpos

using Genera

include("type/InfInts.jl")

end # IntsWithInfinity


text:
module NewtonsMethod

using ForwardDiff, LinearAlgebra

function newtonroot(f, f′; x_0, tol=1E-7, maxiter=1000)
    # setup the algorithm
    x_old = x_0
    normdiff = Inf
    iter = 1
    while normdiff > tol && iter <= maxiter
        x_new = x_old - ( f(x_old) / f′(x_old) ) # use Newton's method
        normdiff = norm(x_new - x_old)
        x_old = x_new
        iter = iter + 1
    end
    if iter <= maxiter
        return (value = x_old, normdiff=normdiff, iter=iter) # A named tuple
    else
        println("No unit root found with tolerance level $tol and Maximum Number of Iterations $maxiter")
        return nothing
    end
end

function newtonroot(f; x_0, tol=1E-7, maxiter=1000)
    # setup the algorithm
    x_old = x_0
    normdiff = Inf
    iter = 1
    D(f) = x -> ForwardDiff.derivative(f, x)
    f′ = D(f) # use auto diff
    while normdiff > tol && iter <= maxiter
        x_new = x_old - ( f(x_old) / f′(x_old) ) # use Newton's method
        normdiff = norm(x_new - x_old)
        x_old = x_new
        iter = iter + 1
    end
    if iter <= maxiter
        return (value = x_old, normdiff=normdiff, iter=iter) # A named tuple
    else
        println("No unit root found with tolerance level $tol and Maximum Number of Iterations $maxiter")
        return nothing
    end
end

export newtonroot

end


text:
# Parameters:
#	N::Int
#	- the problem size
#	FFT::Array{Array{Float64, 2}, 1}
#	- a Fast Fourier Transform of size N
#	- should be the output of sn_fft() or sn_fft_sp()
#	YOR::Array{Array{Array{SparseMatrixCSC, 1}, 1}, 1}
#	- output1 from yor()
#	PT::Array{Array{Array{Int, 1}, 1}, 1}
#	- output2 from yor()
# Return Values:
#	SNF::Array{Float64, 1}
#	- the function over Sn that corresponds to FFT
function sn_ifft(N::Int, FFT::Array{Array{Float64, 2}, 1},  YOR::Array{Array{Array{SparseMatrixCSC, 1}, 1}, 1}, PT::Array{Array{Array{Int, 1}, 1}, 1})
	np = nprocs()
	if np == 1 || N < 10
		SNF = Array(Float64, factorial(N))
		compute_ifft(N, SNF, FFT, YOR, PT, Counter(1))
		return SNF
	else
		pSNFA = Array(Array{Float64, 1}, N)
		RR_FFT = Array(RemoteRef, np)
		RR_YOR = Array(RemoteRef, np)
		RR_PT = Array(RemoteRef, np)
		for p = 1:np
			if p != myid()
				RR_FFT[p] = RemoteRef(p)
				put!(RR_FFT[p], FFT)
				RR_YOR[p] = RemoteRef(p)
				put!(RR_YOR[p], YOR)
				RR_PT[p] = RemoteRef(p)
				put!(RR_PT[p], PT)
			end
		end
		i = 1
		nextidx() = (idx = i; i += 1; idx)
		@sync begin
			for p = 1:np
				if p != myid()
					@async begin
						while true
							idx = nextidx()
							if idx > N
								break
							end
							pSNFA[idx] = remotecall_fetch(p, compute_sifft_remote, N, idx, RR_FFT[p], RR_YOR[p], RR_PT[p])
						end
					end
				end
			end
		end
		SNF = Array(Float64, factorial(N))
		BS = factorial(N - 1)
		i = 1
		for n = 1:N
			pSNF = pSNFA[n]
			for si = 1:BS
				SNF[i] = pSNF[si]
				i += 1
			end
		end
		return SNF
	end
end

# Parameters:
#	N::Int
#	- the problem size
#	SNF::Array{Float64, 1}
#	- the function over Sn that is being calculated
#	FFT::Array{Array{Float64, 2}, 1}
#	- a Fast Fourier Transform of size N
#	- should be the output of sn_fft() or sn_fft_sp()
#	YOR::Array{Array{Array{SparseMatrixCSC, 1}, 1}, 1}
#	- output1 from yor()
#	PT::Array{Array{Array{Int, 1}, 1}, 1}
#	- output2 from yor()
#	C::Counter
#	- a wrapper for an Int that is used to coordinate the recursive defining of the elements of SNF
# Return Values:
#	None
#	- once the recursion is finished, SNF is full
function compute_ifft(N::Int, SNF::Array{Float64, 1}, FFT::Array{Array{Float64, 2}, 1}, YOR::Array{Array{Array{SparseMatrixCSC, 1}, 1}, 1}, PT::Array{Array{Array{Int, 1}, 1}, 1}, C::Counter)
	if N == 2
		YORn = YOR[2]
		sFFT = 0.5 * (full(YORn[1][1]) * FFT[1] + full(YORn[2][1]) * FFT[2])
		SNF[C.N] = sFFT[1,1]
		C.N += 1
		sFFT = 0.5 * (FFT[1] + FFT[2])
		SNF[C.N] = sFFT[1,1]
		C.N += 1
	else
		YORn = YOR[N]
		NPn = length(YORn)
		YORd = YOR[N - 1]
		NPd = length(YORd)
		PTn = PT[N]
		sFFT = Array(Array{Float64, 2}, NPd)
		for n = 1:N
			for p = 1:NPd
				Dim = size(YORd[p][1], 1)
				sFFT[p] = zeros(Float64, Dim, Dim)
			end
			for p = 1:NPn
				update_sfft!(N, n, sFFT, FFT[p], YORn[p], YORd, PTn[p])
			end
			compute_ifft(N - 1, SNF, sFFT, YOR, PT, C)
		end
	end
end

# Parameters:
#	N::Int
#	- the problem size
#	n::Int
#	- determines which left-sided coset we are defining the FFT for
#	sFFT::Array{Array{Float64, 2}, 1}
#	- the FFT of the nth left-sided coset that we are defining
#	FFTp::Array{Float64, 2}
#	- the pth component of the FFT of size N
#	YORnp::Array{SparseMatrixCSC, 1}, 1}
#	- Youngs Orthogonal Representations for the pth Partition of N
#	YORd::Array{Array{SparseMatrixCSC, 1}, 1}
#	- Youngs Orthogonal Reprsentations for the Partitions of N - 1
#	PTnp::Array{Array{Array{Int, 1}, 1}, 1}
#	- The decomposition indices for the pth Partition of N
# Return Values:
#	None
#	- sFFT is fully defined once all of the components of FFT have been used
function update_sfft!(N::Int, n::Int, sFFT::Array{Array{Float64, 2}, 1}, FFTp::Array{Float64, 2}, YORnp::Array{SparseMatrixCSC, 1}, YORd::Array{Array{SparseMatrixCSC, 1}, 1}, PTnp::Array{Int, 1})
	Dim = size(YORnp[1], 1)
	M = eye(Float64, Dim, Dim)
	for ccn = n:(N - 1)
		M = M * YORnp[ccn]
	end
	M = transpose(M)
	M = M * FFTp
	lb = 1
	for d = 1:length(PTnp)
		index = PTnp[d]
		sDim = size(YORd[index][1], 1)
		ub = lb + sDim - 1
		sFFT[index] += (Dim / (sDim * N)) * M[lb:ub, lb:ub]
		lb += sDim
	end
end

# Parameters:
#	N::Int
#	- the problem size
#	n::Int
#	- determines which left-sided coset we are defining the function for
#	FFT::Array{Array{Float64, 2}, 1}
#	- a Fast Fourier Transform of size N
#	- should be the output of sn_fft() or sn_fft_sp()
#	YOR::Array{Array{Array{SparseMatrixCSC, 1}, 1}, 1}
#	- output1 from yor()
#	PT::Array{Array{Array{Int, 1}, 1}, 1}
#	- output2 from yor()
# Return Values:
#	pSNF::Array{Float64, 1}
#	- the function on the nth left-sided coset of Sn
function compute_sifft(N::Int, n::Int, FFT::Array{Array{Float64, 2}, 1}, YOR::Array{Array{Array{SparseMatrixCSC, 1}, 1}, 1},  PT::Array{Array{Array{Int, 1}, 1}, 1})
	pSNF = Array(Float64, factorial(N - 1))
	C = Counter(1)
	YORn = YOR[N]
	NPn = length(YORn)
	YORd = YOR[N - 1]
	NPd = length(YORd)
	PTn = PT[N]
	sFFT = Array(Array{Float64, 2}, NPd)
	for p = 1:NPd
		Dim = size(YORd[p][1], 1)
		sFFT[p] = zeros(Float64, Dim, Dim)
	end
	for p = 1:NPn
		update_sfft!(N, n, sFFT, FFT[p], YORn[p], YORd, PTn[p])
	end
	compute_ifft(N - 1, pSNF, sFFT, YOR, PT, C)
	return pSNF
end

function compute_sifft_remote(N::Int, n::Int, RR_FFT::RemoteRef, RR_YOR::RemoteRef, RR_PT::RemoteRef)
	FFT = fetch(RR_FFT)
	YOR = fetch(RR_YOR)
	PT = fetch(RR_PT)
	pSNF = compute_sifft(N, n, FFT, YOR, PT)
	return pSNF
end


text:
# Busjle et. al. 2009
# ===================

function _buslje09(aln, usegap, clusters, lambda, apc)
    mi = estimateincolumns(aln, ResidueCount{Float64, 2, usegap}, MutualInformation{Float64}(),
                           AdditiveSmoothing{Float64}(lambda), clusters, false)
    if apc
        APC!(mi)
    end
    mi
end

"""
This function takes a MSA or a file and a `Format` as first arguments.
Calculates a Z score and a corrected MI/MIp as described on **Busjle et. al. 2009**

Argument, type, default value and descriptions:

```
  - lambda      Float64   0.05    Low count value
  - clustering  Bool      true    Sequence clustering (Hobohm I)
  - threshold             62      Percent identity threshold for clustering
  - maxgap      Float64   0.5     Maximum fraction of gaps in positions included in calculation
  - apc         Bool      true    Use APC correction (MIp)
  - usegap      Bool      false   Use gaps on statistics
  - samples     Int       100     Number of samples for Z-score
  - fixedgaps   Bool      true    Fix gaps positions for the random samples
```

This function returns:

```
  - Z score
  - MI or MIp
```
"""
function buslje09(aln::AbstractMatrix{Residue}; lambda::Float64=0.05,
                  clustering::Bool=true, threshold=62,
                  maxgap::Float64=0.5, apc::Bool=true, samples::Int=100,
                  usegap::Bool=false, fixedgaps::Bool=true)
    used = gapfraction(aln,1) .<= maxgap
    col_labels = columnlabels(aln)[used]
    ncol = ncolumns(aln)
    aln = filtercolumns(aln, used)
    clusters = clustering ? hobohmI(aln, threshold) : NoClustering()
    mi = _buslje09(aln, usegap, clusters, lambda, apc)
    labels!(mi, col_labels)
    if samples > 0
        rand_mi = Array(typeof(mi), samples)
        for ns in 1:samples
            fixedgaps ? shuffle_residues_sequencewise!(aln) : shuffle_sequencewise!(aln)
            rand_mi[ns] = _buslje09(aln, usegap, clusters, lambda, apc)
        end
        zmi = zscore(rand_mi, mi)
        labels!(zmi, col_labels)
        return(zmi, mi)
    else
        return(zeros(mi), mi)
    end
end

function buslje09{T <: Format}(filename::AbstractString, format::Type{T}; kargs...)
    aln = read(filename, T, AnnotatedMultipleSequenceAlignment, generatemapping=true)
    buslje09(aln; kargs...)
end

# MIToS BLMI: Blosum MI
# =====================

function _BLMI(aln, clusters, alpha, beta, apc, lambda::Float64=zero(Float64))
    mi = estimateincolumns(aln, Float64, ResidueProbability{Float64, 2, false},
                           alpha, beta,
                           MutualInformation{Float64}(),
                           AdditiveSmoothing(lambda), clusters, false, 0.0)
    if apc
        APC!(mi)
    end
    mi
end

"""
This function takes a MSA or a file and a `Format` as first arguments.
Calculates a Z score (ZBLMI) and a corrected MI/MIp as described on **Busjle et. al. 2009** but using using BLOSUM62 pseudo frequencies instead of a fixed pseudocount.

Argument, type, default value and descriptions:

```
  - beta        Float64   8.512   β for BLOSUM62 pseudo frequencies
  - lambda      Float64   0.0     Low count value
  - threshold             62      Percent identity threshold for sequence clustering (Hobohm I)
  - maxgap      Float64   0.5     Maximum fraction of gaps in positions included in calculation
  - apc         Bool      true    Use APC correction (MIp)
  - samples     Int       50      Number of samples for Z-score
  - fixedgaps   Bool      true    Fix gaps positions for the random samples
```

This function returns:

```
  - Z score (ZBLMI)
  - MI or MIp using BLOSUM62 pseudo frequencies (BLMI/BLMIp)
```
"""
function BLMI(aln::AbstractMatrix{Residue}; beta::Float64=8.512, threshold=62,
              maxgap::Float64=0.5, apc::Bool=true, samples::Int=50,
              fixedgaps::Bool=true, lambda::Float64=zero(Float64))
    used = gapfraction(aln,1) .<= maxgap
    col_labels = columnlabels(aln)[used]
    ncol = ncolumns(aln)
    aln = filtercolumns(aln, used)
    clusters = hobohmI(aln, threshold)
    numbercl = nclusters(clusters)
    mi = _BLMI(aln, clusters, numbercl, beta, apc, lambda)
    labels!(mi, col_labels)
    if samples > 0
        rand_mi = Array(typeof(mi), samples)
        for ns in 1:samples
            fixedgaps ? shuffle_residues_sequencewise!(aln) : shuffle_sequencewise!(aln)
            rand_mi[ns] = _BLMI(aln, clusters, numbercl, beta, apc, lambda)
        end
        zmi = zscore(rand_mi, mi)
        labels!(zmi, col_labels)
        return(zmi, mi)
    else
        return(zeros(mi), mi)
    end
end

function BLMI{T <: Format}(filename::AbstractString, format::Type{T}; kargs...)
    aln = read(filename, T, AnnotatedMultipleSequenceAlignment, generatemapping=true)
    BLMI(aln; kargs...)
end

# MIToS Pairwise Gap Percentage
# =============================

"""
This function takes a MSA or a file and a `Format` as first arguments.
Calculates the percentage of gaps on columns pairs (union and intersection) using sequence clustering (Hobohm I).

Argument, type, default value and descriptions:

  - clustering  Bool      true    Sequence clustering (Hobohm I)
  - threshold             62      Percent identity threshold for sequence clustering (Hobohm I)

This function returns:

  - pairwise gap percentage (union)
  - pairwise gap percentage (intersection)
"""
function pairwisegapfraction(aln::AbstractMatrix{Residue}; clustering::Bool=true, threshold=62)
    col_labels = columnlabels(aln)
    clusters = clustering ? hobohmI(aln, threshold) : NoClustering()
    gu = estimateincolumns(aln, ResidueCount{Float64, 2, true}, GapUnionPercentage{Float64}(), zero(AdditiveSmoothing{Float64}), clusters, true)
    gi = estimateincolumns(aln, ResidueCount{Float64, 2, true}, GapIntersectionPercentage{Float64}(), zero(AdditiveSmoothing{Float64}), clusters, true)
    labels!(gu, col_labels)
    labels!(gi, col_labels)
    gu, gi
end

function pairwisegapfraction{T <: Format}(filename::AbstractString, format::Type{T}; kargs...)
    aln = read(filename, T, AnnotatedMultipleSequenceAlignment, generatemapping=true)
    pairwisegapfraction(aln; kargs...)
end


text:
@testset "minimum" begin

    # Call function and determine threshold on some synthetic data.
    samples = map(x->gaussian_mixture(x, P₁ = 0.5, μ₁ = 64, σ₁ = 15, P₂ = 0.5, μ₂ = 192, σ₂ = 15), 0:255)
    edges = 0:255
    t = find_threshold(MinimumIntermodes(), samples, edges)
    @test t == 128

    # Call function and determine threshold on cameraman image for which we know
    # what the correct threshold should be (based on an implementation in
    # another project).
    img = testimage("cameraman")
    edges, counts = build_histogram(img,  256)
    t = find_threshold(MinimumIntermodes(), counts[1:end], edges)
    @test t == 0.296875

end


text:
@testset "CUDNN" begin

if !isdefined(CuArrays, :CUDNN)
@warn "Not testing CUDNN"
else
using CuArrays.CUDNN
@info "Testing CUDNN $(CUDNN.version())"

@testset "NNlib" begin
  using NNlib
  using NNlib: ∇conv_data, ∇conv_filter,
               maxpool, meanpool, ∇maxpool, ∇meanpool,
               softmax, ∇softmax, logsoftmax, ∇logsoftmax

  @test testf(NNlib.conv, rand(Float64, 10, 10, 3, 1), rand(Float64, 2, 2, 3, 4))
  @test testf(∇conv_data, rand(Float64, 9, 9, 4, 1), rand(Float64, 2, 2, 3, 4))
  @test testf(∇conv_filter, rand(Float64, 9, 9, 4, 1), rand(Float64, 10, 10, 3, 1))
  @test testf(CuArrays.CUDNN.∇conv_bias!, cu(rand(Float64, 1, 1, 10, 1)), cu(rand(Float64, 10, 10, 10, 1)))

  @test testf(NNlib.conv, rand(Float64, 10, 10, 3, 1), rand(Float64, 2, 2, 3, 4); dilation=2)
  @test testf(∇conv_data, rand(Float64, 8, 8, 4, 1), rand(Float64, 2, 2, 3, 4); dilation=2)
  @test testf(∇conv_filter, rand(Float64, 8, 8, 4, 1), rand(Float64, 10, 10, 3, 1); dilation=2)

  @test testf(NNlib.crosscor, rand(Float64, 10, 10, 3, 1), rand(Float64, 2, 2, 3, 4))
  @test testf(∇conv_data, rand(Float64, 9, 9, 4, 1), rand(Float64, 2, 2, 3, 4); flipkernel=1)
  @test testf(∇conv_filter, rand(Float64, 9, 9, 4, 1), rand(Float64, 10, 10, 3, 1); flipkernel=1)

  @test_nowarn NNlib.conv!(cu(zeros(Float64, 9, 9, 3, 1)), cu(rand(Float64, 10, 10, 1, 1)), cu(rand(Float64, 2, 2, 1, 3)), algo=1)
  @test_nowarn NNlib.∇conv_data!(cu(zeros(Float64, 10, 10, 1, 1)), cu(ones(Float64, 9, 9, 3, 1)), cu(rand(Float64, 2, 2, 1, 3)), algo=1)
  @test_nowarn NNlib.∇conv_filter!(cu(zeros(Float64, 2, 2, 1, 3)), cu(ones(Float64, 9, 9, 3, 1)), cu(rand(Float64, 10, 10, 1, 1)), algo=1)

  @test testf(NNlib.conv, rand(Float64, 10, 10, 10, 3, 1), rand(Float64, 2, 2, 2, 3, 4))
  @test testf(∇conv_data, rand(Float64, 9, 9, 9, 4, 1),  rand(Float64, 2, 2, 2, 3, 4))
  @test testf(∇conv_filter, rand(Float64, 9, 9, 9, 4, 1), rand(Float64, 10, 10, 10, 3, 1))

  @test testf(NNlib.conv, rand(Float64, 10, 10, 10, 3, 1), rand(Float64, 2, 2, 2, 3, 4); dilation=2)
  @test testf(∇conv_data, rand(Float64, 8, 8, 8, 4, 1), rand(Float64, 2, 2, 2, 3, 4); dilation=2)
  @test testf(∇conv_filter, rand(Float64, 8, 8, 8, 4, 1), rand(Float64, 10, 10, 10, 3, 1); dilation=2)

  @test testf(NNlib.crosscor, rand(Float64, 10, 10, 10, 3, 1), rand(Float64, 2, 2, 2, 3, 4))
  @test testf(∇conv_data, rand(Float64, 9, 9, 9, 4, 1), rand(Float64, 2, 2, 2, 3, 4); flipkernel=1)
  @test testf(∇conv_filter, rand(Float64, 9, 9, 9, 4, 1), rand(Float64, 10, 10, 10, 3, 1); flipkernel=1)

  @test testf(x -> maxpool(x, (2,2)), rand(Float64, 10, 10, 3, 1))
  @test testf(x -> meanpool(x, (2,2)), rand(Float64, 10, 10, 3, 1))
  @test testf((x, dy) -> ∇maxpool(dy, maxpool(x, (2,2)), x, (2,2)), rand(Float64, 10, 10, 3, 1), rand(Float64, 5, 5, 3, 1))
  @test testf((x, dy) -> ∇meanpool(dy, meanpool(x, (2,2)), x, (2,2)), rand(Float64, 10, 10, 3, 1), rand(Float64, 5, 5, 3, 1))

  @test testf(x -> maxpool(x, (2,2,2)), rand(Float64, 10, 10, 10, 3, 1))
  @test testf(x -> meanpool(x, (2,2,2)), rand(Float64, 10, 10, 10, 3, 1))
  @test testf((x, dy) -> ∇maxpool(dy, maxpool(x, (2,2,2)), x, (2,2,2)), rand(Float64, 10, 10, 10, 3, 1), rand(Float64, 5, 5, 5, 3, 1))
  @test testf((x, dy) -> ∇meanpool(dy, meanpool(x, (2,2,2)), x, (2,2,2)), rand(Float64, 10, 10, 10, 3, 1), rand(Float64, 5, 5, 5, 3, 1))

  for dims in [(5,5), (5,)]
    @test testf(softmax, rand(Float64, dims))
    @test testf(∇softmax, rand(Float64, dims), rand(Float64, dims))
    @test testf(logsoftmax, rand(Float64, dims))
    @test testf(∇logsoftmax, rand(Float64, dims), rand(Float64, dims))
  end
end

@testset "Activations and Other Ops" begin
  @test testf(CuArrays.CUDNN.cudnnAddTensor, cu(rand(Float64, 10, 10, 3, 1)), cu(rand(Float64, 10, 10, 3, 1)))
  @test testf(CuArrays.CUDNN.cudnnActivationForward, cu(rand(Float64, 10, 10, 3, 1)), cu(rand(Float64, 10, 10, 3, 1)))
  @test testf(CuArrays.CUDNN.cudnnActivationBackward, cu(rand(Float64, 10, 10, 3, 1)), cu(rand(Float64, 10, 10, 3, 1)), cu(rand(Float64, 10, 10, 3, 1)), cu(rand(Float64, 10, 10, 3, 1)))
end

end

end


text:
# Adapted from:
# SOSDEMO8 --- Bounds in Probability
# Section 3.8 of SOSTOOLS User's Manual

@testset "SOSDEMO8 with $(factory.optimizer_constructor)" for factory in sdp_factories
    @polyvar x

    # The probability adds up to one.
    μ0 = 1
    # Mean
    μ1  = 1
    # Variance
    σ = 1/2
    # E(x^2)
    μ2 = σ^2+μ1^2
    # Support of the random variable
    R = [0,5]

    # Event whose probability we want to bound
    E = [4,5]

    m = SOSModel(factory)

    @variable m a
    @variable m b
    @variable m c

    P = a + b*x + c*x^2

    # Nonnegative on the support
    K = @set 0 <= x && x <= 5
    cons = @constraint(m, P >= 0, domain = K)

    # Greater than one on the event
    @constraint(m, P >= 1, domain = (@set 4 <= x && x <= 5))

    # The bound
    bnd =  a * μ0 + b * μ1 + c * μ2

    @objective m Min bnd

    JuMP.optimize!(m)
    @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT

    @test isapprox(JuMP.objective_value(m), 1/37, rtol=1e-5)
    @test isapprox(JuMP.value(P), ((12/37)x-11/37)^2, rtol=1e-3)

    @test sos_decomposition(cons, K, 0.0001) isa SOSDecompositionWithDomain
end


text:
function normalize(c::Char, isarabic::Bool, encoder::AbstractEncoder)
    ch = isarabic ? Symbol(c) : encoder.decode[Symbol(c)]
    if !isarabic
        return string(encoder.encode[SP_DEDIAC_MAPPING[ch]])
    else
        return string(SP_DEDIAC_MAPPING[ch])
    end
end

"""
    normalize(s::String)

Normalize a Arabic or Buckwalter `String` characters.

# Examples
```julia-repl
julia> data = QuranData()
julia> crps, tnzl = load(data)
julia> crpsdata = table(crps);
julia> tnzldata = table(tnzl);
julia> normalize(dediac(verses(crpsdata[1][1])[1]))
"bsm Allh AlrHmAn AlrHym"
```
"""
function normalize(s::String)
    trans = Transliterator()
    isarabic = in(Symbol(s[1]), collect(keys(trans.encode))) ? true : false
    word = ""
    for c in s
        if c === ' '
            word *= " "
            continue
        end
        isnormalize = !isarabic ? in(trans.decode[Symbol(c)], SP_DEDIAC_KEYS) : in(Symbol(c), SP_DEDIAC_KEYS)
        if isnormalize
            word *= normalize(c, isarabic, trans)
        else
            word *= c
        end
    end
    word = normalize(word, :tatweel)
    return word
end

"""
    normalize(s::String, chars::Array{Symbol,1})

Normalize a specific Arabic or Buckwalter `String` character/s (`chars`).

# Examples
```julia-repl
julia> data = QuranData()
julia> crps, tnzl = load(data)
julia> crpsdata = table(crps);
julia> tnzldata = table(tnzl);
julia> normalize(arabic(verses(crpsdata[1][1])[1]), [:alif_khanjareeya, :hamzat_wasl])
"بِسْمِ اللَّهِ الرَّحْمَانِ الرَّحِيمِ"
```
"""
function normalize(s::String, chars::Array{Symbol,1})
    for char in chars
        s = normalize(s, char)
    end
    return s
end

"""
    normalize(s::String, char::Symbol)

Normalize a specific Arabic or Buckwalter `String` character (`chars`).

# Examples
```julia-repl
julia> data = QuranData()
julia> crps, tnzl = load(data)
julia> crpsdata = table(crps);
julia> tnzldata = table(tnzl);
julia> normalize(arabic(verses(crpsdata[1][1])[1]), :hamzat_wasl)
"بِسْمِ اللَّهِ الرَّحْمَٰنِ الرَّحِيمِ"
```
"""
function normalize(s::String, char::Symbol)
    trans = Transliterator()
    isarabic = in(Symbol(s[1]), collect(keys(trans.encode))) ? true : false

    s = isarabic ? s : arabic(s)
    if char === :tatweel
        word = replace(s, Regex(string(Char(0x0640)[1])) => s"")
    elseif char === :alif_maddah
        word = replace(s, Regex(string(Char(0x0622)[1])) => string(Char(0x0627)))
        word = replace(word, Regex(string(Char(0x0653)[1])) => s"")
    elseif char === :alif_hamza_above
        word = replace(s, Regex(string(Char(0x0623)[1])) => string(Char(0x0627)))
    elseif char === :alif_khanjareeya
        word = replace(s, Regex(string(Char(0x0670)[1])) => string(Char(0x0627)))
    elseif char === :hamzat_wasl
        word = replace(s, Regex(string(Char(0x0671)[1])) => string(Char(0x0627)))
    elseif char === :alif_hamza_below
        word = replace(s, Regex(string(Char(0x0625)[1])) => string(Char(0x0627)))
    elseif char === :waw_hamza_above
        word = replace(s, Regex(string(Char(0x0624)[1])) => string(Char(0x0648)))
    elseif char === :ya_hamza_above
        word = replace(s, Regex(string(Char(0x0626)[1])) => string(Char(0x064A)))
    elseif char === :alif_maksura
        word = replace(s, Regex(string(Char(0x0649)[1])) => string(Char(0x064A)))
    elseif char === :ta_marbuta
        word = replace(s, Regex(string(Char(0x0629)[1])) => string(Char(0x0647)))
    else
        throw(DomainError(char, "Character not found"))
    end
    return isarabic ? word : encode(word)
end


text:
using RadeonProRender, GeometryBasics, Colors
const RPR = RadeonProRender

function translationmatrix(t::Vec{3,T}) where {T}
    T0, T1 = zero(T), one(T)
    return Mat{4}(T1, T0, T0, T0, T0, T1, T0, T0, T0, T0, T1, T0, t[1], t[2], t[3], T1)
end

isdefined(Main, :context) && RPR.release(context)

context = RPR.Context()
scene = RPR.Scene(context)
matsys = RPR.MaterialSystem(context, 0)

camera = RPR.Camera(context)
lookat!(camera, Vec3f(7, 0, 5), Vec3f(2, 0, 2), Vec3f(0, 0, 1))
RPR.rprCameraSetFocalLength(camera, 45.0)
set!(scene, camera)
set!(context, scene)

env_light = RadeonProRender.EnvironmentLight(context)
image_path = joinpath(@__DIR__, "studio026.exr")
img = RPR.Image(context, image_path)
set!(env_light, img)
setintensityscale!(env_light, 1.1)
push!(scene, env_light)

light = RPR.PointLight(context)
transform!(light, translationmatrix(Vec3f(2, 0, 5)))
f = 20
RPR.setradiantpower!(light, 500 / f, 641 / f, 630 / f)
push!(scene, light)

function add_shape!(scene, context, matsys, mesh; material=RPR.RPR_MATERIAL_NODE_DIFFUSE,
                    color=colorant"green", roughness=0.01)
    rpr_mesh = RadeonProRender.Shape(context, mesh)
    push!(scene, rpr_mesh)
    m = RPR.MaterialNode(matsys, material)
    set!(m, RPR.RPR_MATERIAL_INPUT_COLOR, color)
    set!(m, RPR.RPR_MATERIAL_INPUT_ROUGHNESS, roughness, roughness, roughness, roughness)
    set!(rpr_mesh, m)
    return rpr_mesh, m
end

mesh, mat = add_shape!(scene, context, matsys, Rect3f(Vec3f(-50, -50, -1), Vec3f(50*2, 50*2, 1));
                       color=colorant"white")

x, y = collect(-8:0.5:8), collect(-8:0.5:8)
z = [sinc(√(X^2 + Y^2) / π) for X ∈ x, Y ∈ y]

M, N = size(z)
points = map(x-> x ./ Point3f(4, 4, 1.0) .+ Point3f(0, 0, 1.0), vec(Point3f.(x, y', z)))
# Connect the vetices with faces, as one would use for a 2D Rectangle
# grid with M,N grid points
faces = decompose(LineFace{GLIndex}, Tesselation(Rect2(0, 0, 1, 1), (M, N)))
indices = RPR.rpr_int[]
for elem in faces
    push!(indices, first(elem)-1, first(elem)-1, last(elem)-1, last(elem)-1)
end
curve = RPR.Curve(context, points, indices, [0.01], [Vec2f(0.1)], [length(faces)])

m = RPR.MaterialNode(matsys, RPR.RPR_MATERIAL_NODE_DIFFUSE)
set!(m, RPR.RPR_MATERIAL_INPUT_COLOR, colorant"red")
set!(m, RPR.RPR_MATERIAL_INPUT_ROUGHNESS, 0.0, 0.0, 0.0, 0.0)
set!(curve, m)
push!(scene, curve)

fb_size = (800, 600)
frame_buffer = RPR.FrameBuffer(context, RGBA, fb_size)
frame_bufferSolved = RPR.FrameBuffer(context, RGBA, fb_size)
set!(context, RPR.RPR_AOV_COLOR, frame_buffer)

set_standard_tonemapping!(context)

begin
    clear!(frame_buffer)
    RPR.rprContextSetParameterByKey1u(context, RPR.RPR_CONTEXT_ITERATIONS, 200)
    @time RPR.render(context)
    RPR.rprContextResolveFrameBuffer(context, frame_buffer, frame_bufferSolved, false)
    RPR.save(frame_bufferSolved, "curves.png")
end


text:
struct RibbonGeometry <: Gadfly.GeometryElement
    default_statistic::Gadfly.StatisticElement
    fill::Bool
    tag::Symbol
end
RibbonGeometry(default_statistic=Gadfly.Stat.identity(); fill=true, tag=empty_tag) =
        RibbonGeometry(default_statistic, fill, tag)

"""
    Geom.ribbon[(; fill=true)]

Draw a ribbon at the positions in `x` bounded above and below by `ymax` and
`ymin`, respectively.  Optionally draw multiple ribbons by grouping with `color` and `alpha` (for `fill=true`),
 or `color` and `linestyle` (for `fill=false`).
"""
const ribbon = RibbonGeometry

default_statistic(geom::RibbonGeometry) = geom.default_statistic

element_aesthetics(::RibbonGeometry) = [:x, :ymin, :ymax, :color, :linestyle, :alpha]

function render(geom::RibbonGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.ribbon", aes, :x, :ymin, :ymax)
    Gadfly.assert_aesthetics_equal_length("Geom.ribbon", aes, element_aesthetics(geom)...)

    default_aes = Gadfly.Aesthetics()
    default_aes.linestyle = fill(1, length(aes.x))
    default_aes.color = fill(theme.default_color, length(aes.x))
    default_aes.alpha = fill(1, length(aes.x))
    aes = inherit(aes, default_aes)

    aes_x, aes_ymin, aes_ymax, aes_color, aes_linestyle, aes_alpha =
         concretize(aes.x, aes.ymin, aes.ymax, aes.color, aes.linestyle, aes.alpha)
    XT, CT, LST, AT = eltype(aes_x), eltype(aes_color), eltype(aes_linestyle), eltype(aes_alpha)
    YT = eltype(aes_ymin)
    groups = collect((Tuple{CT, LST, AT}), zip(aes_color, aes_linestyle, aes_alpha))
    ug = unique(groups)

    V = Vector{Tuple{XT, YT}}
    K = Tuple{CT, LST, AT}

    max_points = Dict{K, V}(g=>V[] for g in ug)
    for (x, y, c, ls, a) in zip(aes_x, aes_ymax, aes_color, aes_linestyle, aes_alpha)
        push!(max_points[(c,ls,a)], (x, y))
    end

    min_points = Dict{K, V}(g=>V[] for g in ug)
    for (x, y, c, ls, a) in zip(aes_x, aes_ymin, aes_color, aes_linestyle, aes_alpha)
        push!(min_points[(c,ls,a)], (x, y))
    end

    for k in keys(max_points)
        sort!(max_points[k], by=first)
        sort!(min_points[k], by=first, rev=true)
    end

    kys = keys(max_points)
    polys = [collect(Tuple{XT, YT}, Iterators.flatten((min_points[k], max_points[k]))) for k in kys]
    lines = [collect(Tuple{XT, YT}, Iterators.flatten((min_points[k], [(last(min_points[k])[1], oneunit(YT)*NaN)], max_points[k]))) for k in kys]

    n = length(kys)
    colors = Vector{Union{Colorant, String}}(undef, n)
    linestyles = Vector{Vector{Measure}}(undef, n)
    alphas = Vector{Float64}(undef, n)
    alpha_discrete  = AT <: Int

    for (i, (c,ls,a)) in enumerate(kys)
        colors[i] = theme.lowlight_color(c)
        linestyles[i] = Gadfly.get_stroke_vector(theme.line_style[ls])
        alphas[i] = alpha_discrete ? theme.alphas[a] : a
    end

    ctx = context()

    geom.fill ? compose!(ctx, Compose.polygon(polys, geom.tag), fill(colors), fillopacity(alphas)) :
        compose!(ctx, Compose.line(lines, geom.tag), fill(nothing), stroke(colors), strokedash(linestyles))

    return compose!(ctx, svgclass("geometry"), linewidth(theme.line_width))
end


text:
using FractionalCalculus
using Test

@testset "Test Direct Fractional Integral" begin
    @test isapprox(fracint(x->x^4, 0.5, 0, 1.23, 1e-8, RLDirect())[1], 1.163931646862977; atol = 0.00001)
    @test isapprox(fracint(x->x^4, 0.25, 0, 1.23, 1e-8, RLDirect())[1], 1.64294134724409; atol = 0.00001)
    @test isapprox(fracint(x->x^4, 0.15, 0, 1.23, 1e-8, RLDirect())[1], 1.878978485456182; atol = 0.00001)
    @test isapprox(fracint(x->sin(x), 0.5, 0, 0.97, 1e-8, RLDirect())[1], 0.6442601778843754; atol = 0.00001)
    @test isapprox(fracint(x->sin(x), 0.25, 0, 0.97, 1e-8, RLDirect())[1], 0.7448298999485293; atol = 0.00001)
    @test isapprox(fracint(x->sin(x), 0.15, 0, 0.97, 1e-8, RLDirect())[1], 0.7802555196986999; atol = 0.00001)
end

@testset "Test Piecewise Fractional Integral" begin
    @test isapprox(fracint(x->x^4, 0.5, 1.23, 1e-4, RLPiecewise()), 1.163931646862977; atol = 1e-5)
    @test isapprox(fracint(x->x^4, 0.25, 1.23, 1e-4, RLPiecewise()), 1.64294134724409; atol = 1e-5)
    @test isapprox(fracint(sin, 0.5, 0.97, 1e-4, RLPiecewise()), 0.6442601778843754; atol = 1e-5)

    result = fracint(x->x^4, 0.5, [1.23, 2, 3], 1e-4, RLPiecewise())
    @test isapprox(result[1], 1.163931646862977; atol = 1e-3)
    @test isapprox(result[2], 10.375035632494745; atol = 1e-3)
    @test isapprox(result[3], 64.3280200256207; atol = 1e-3)
end

@testset "Test RLInt_Approx" begin
    @test isapprox(fracint(x->x^4, 0.5, 1.23, 1e-4, RLIntApprox()), 1.163931646862977; atol = 1e-4)

    result = fracint(x->x^4, 0.5, [1.23, 2, 3], 1e-4, RLIntApprox())
    @test isapprox(result[1], 1.163931646862977; atol = 1e-3)
    @test isapprox(result[2], 10.375035632494745; atol = 1e-3)
    @test isapprox(result[3], 64.3280200256207; atol = 1e-3)
end

@testset "Test RLInt_Matrix" begin
    @test isapprox(fracint(x->x, 0.5, 1, 0.5, RLIntMatrix()), [0; 0.3535533905932738; 0.8838834764831844]; atol=1e-3)
end

@testset "Test RL_LinearInterp" begin
    @test isapprox(fracint(x->x^4, 0.5, 1.23, 1e-4, RLLinearInterp()), 1.163931646862977; atol = 1e-4)

    result = fracint(x->x^4, 0.5, [1.23, 2, 3], 1e-4, RLLinearInterp())
    @test isapprox(result[1], 1.163931646862977; atol = 1e-3)
    @test isapprox(result[2], 10.375035632494745; atol = 1e-3)
    @test isapprox(result[3], 64.3280200256207; atol = 1e-3)
end

@testset "Test Simpson Fractional Integral" begin
    @test isapprox(fracint(x->x^4, 0.5, 1.23, 1e-7, RLIntSimpson()), 1.163931646862977; atol = 1e-2)

    result = fracint(x->x^4, 0.5, [1.23, 2, 3], 1e-7, RLIntSimpson())
    @test isapprox(result[1], 1.163931646862977; atol = 1e-2)
    @test isapprox(result[2], 10.375035632494745; atol = 1e-2)
    @test isapprox(result[3], 64.3280200256207; atol = 1)
end

@testset "Test Trapezoidal Fractional Integral" begin
    @test isapprox(fracint(x->x^4, 0.5, 1.23, 1e-7, RLIntTrapezoidal()), 1.163931646862977; atol = 1e-4)
#=
    result = fracint(x->x^4, 0.5, [1.23, 2, 3], 1e-7, RLInt_Trapezoidal())
    @test isapprox(result[1], 1.163931646862977; atol = 1e-3)
    @test isapprox(result[2], 10.375035632494745; atol = 1e-3)
    @test isapprox(result[3], 64.3280200256207; atol = 1e-3)
    =#
end

@testset "Test Rectangular Fractional Integral" begin
    @test isapprox(fracint(x->x^4, 0.5, 1.23, 1e-6, RLIntRectangular()), 1.163931646862977; atol = 1e-3)

    @test isapprox(fracint(x->x^4, 0.5, [1.23, 2, 3], 0.001, RLIntRectangular()), [1.1617493910169032; 10.363130001411157; 64.27898282619807]; atol=1e-2)
end

@testset "Test Cubic Spline Interpolation Fractional Integral" begin
    @test isapprox(fracint(x->x^4, 0.5, 1.23, 0.001, RLIntCubicSplineInterp()), 1.163931646862977; atol = 1e-5)
    @test isapprox(fracint(sin, 0.5, 0.97, 0.01, RLIntCubicSplineInterp()), 0.6442601778843754; atol = 1e-6)

    result = fracint(x->x^4, 0.5, [1.23, 2, 3], 0.01, RLIntCubicSplineInterp())
    @test isapprox(result[1], 1.163931646862977; atol = 1e-3)
    @test isapprox(result[2], 10.375035632494745; atol = 1e-3)
    @test isapprox(result[3], 64.3280200256207; atol = 1e-3)
end

@testset "Test fractional integral macros" begin
    @test isapprox(@fracint(x->x^4, 0.5, 1.23), 1.163931646862977; atol = 1e-4)
    @test isapprox(@semifracint(x->x^4, 1.23), 1.163931646862977; atol = 1e-4)
end

@testset "Test Hadamard integral" begin
    (a, b)=fracint(sin, 0.5, 1, 2, 1/3, 0.1, HadamardMat())

    @test isapprox(a, [ 1.0
    1.023136962522694
    1.077552576065
    1.1626221795557758
    1.2764357600635179
    1.414213562373095
    1.5668630279525204
    1.7202493081322234
    1.85605792647593
    1.954772501883529
    2.0]; atol=1e-3)
    @test isapprox(b, [Inf
    3.2545536740283265
    1.9422879792123904
    1.5034729488252023
    1.291625449283718
    1.1492602172225017
    1.0157975575699352
    0.8713525896641822
    0.725484956837982
    0.6074878993825767
    0.5500616274393197]; atol=1e-3)
end


text:
using HTTP

#WebSocket Examples
@async HTTP.WebSockets.listen("127.0.0.1", UInt16(8081)) do ws
    while !eof(ws)
        data = readavailable(ws)
        write(ws, data[end:-1:begin])
    end
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# Truncated normal distribution

TruncatedNormal(mu::Float64, sigma::Float64, a::Float64, b::Float64) =
    Truncated(Normal(mu, sigma), a, b)

TruncatedNormal(mu::Real, sigma::Real, a::Real, b::Real) =
    TruncatedNormal(Float64(mu), Float64(sigma), Float64(a), Float64(b))

### statistics

minimum{T <: Real}(d::Truncated{Normal{T},Continuous}) = d.lower
maximum{T <: Real}(d::Truncated{Normal{T},Continuous}) = d.upper

function mode{T <: Real}(d::Truncated{Normal{T},Continuous})
    μ = mean(d.untruncated)
    d.upper < μ ? d.upper :
    d.lower > μ ? d.lower : μ
end

modes{T <: Real}(d::Truncated{Normal{T},Continuous}) = [mode(d)]

function mean{T <: Real}(d::Truncated{Normal{T},Continuous})
    d0 = d.untruncated
    μ = mean(d0)
    σ = std(d0)
    a = (d.lower - μ) / σ
    b = (d.upper - μ) / σ
    μ + ((normpdf(a) - normpdf(b)) / d.tp) * σ
end

function var{T <: Real}(d::Truncated{Normal{T},Continuous})
    d0 = d.untruncated
    μ = mean(d0)
    σ = std(d0)
    a = (d.lower - μ) / σ
    b = (d.upper - μ) / σ
    z = d.tp
    φa = normpdf(a)
    φb = normpdf(b)
    aφa = isinf(a) ? 0.0 : a * φa
    bφb = isinf(b) ? 0.0 : b * φb
    t1 = (aφa - bφb) / z
    t2 = abs2((φa - φb) / z)
    abs2(σ) * (1 + t1 - t2)
end

function entropy{T <: Real}(d::Truncated{Normal{T},Continuous})
    d0 = d.untruncated
    z = d.tp
    μ = mean(d0)
    σ = std(d0)
    a = (d.lower - μ) / σ
    b = (d.upper - μ) / σ
    aφa = isinf(a) ? 0.0 : a * normpdf(a)
    bφb = isinf(b) ? 0.0 : b * normpdf(b)
    0.5 * (log2π + 1.) + log(σ * z) + (aφa - bφb) / (2.0 * z)
end

### sampling

## Benchmarks doesn't seem to show that this specialized
## sampler is faster than the generic quantile-based method

# function rand(d::Truncated{Normal,Continuous})
#     d0 = d.untruncated
#     μ = mean(d0)
#     σ = std(d0)
#     a = (d.lower - μ) / σ
#     b = (d.upper - μ) / σ
#     z = randnt(a, b, d.tp)
#     return μ + σ * z
# end

# Rejection sampler based on algorithm from Robert (1995)
#
#  - Available at http://arxiv.org/abs/0907.4010
#
# function randnt(lb::Float64, ub::Float64, tp::Float64)
#     r::Float64
#     if tp > 0.3   # has considerable chance of falling in [lb, ub]
#         r = randn()
#         while r < lb || r > ub
#             r = randn()
#         end
#         return r

#     else
#         span = ub - lb
#         if lb > 0 && span > 2.0 / (lb + sqrt(lb^2 + 4.0)) * exp((lb^2 - lb * sqrt(lb^2 + 4.0)) / 4.0)
#             a = (lb + sqrt(lb^2 + 4.0))/2.0
#             while true
#                 r = rand(Exponential(1.0 / a)) + lb
#                 u = rand()
#                 if u < exp(-0.5 * (r - a)^2) && r < ub
#                     return r
#                 end
#             end
#         elseif ub < 0 && ub - lb > 2.0 / (-ub + sqrt(ub^2 + 4.0)) * exp((ub^2 + ub * sqrt(ub^2 + 4.0)) / 4.0)
#             a = (-ub + sqrt(ub^2 + 4.0)) / 2.0
#             while true
#                 r = rand(Exponential(1.0 / a)) - ub
#                 u = rand()
#                 if u < exp(-0.5 * (r - a)^2) && r < -lb
#                     return -r
#                 end
#             end
#         else
#             while true
#                 r = lb + rand() * (ub - lb)
#                 u = rand()
#                 if lb > 0
#                     rho = exp((lb^2 - r^2) * 0.5)
#                 elseif ub < 0
#                     rho = exp((ub^2 - r^2) * 0.5)
#                 else
#                     rho = exp(-r^2 * 0.5)
#                 end
#                 if u < rho
#                     return r
#                 end
#             end
#         end
#     end
# end


text:
using PyCall
py"""import validate_candidate"""
iter = [t.iteration for t in trace]
data = [ trace[i+1,1,1].metadata["pop"] for i in iter ]
evo_loss
model = Chain(Dense(d, 15, relu), Dense(15, nclasses))
@info "MLP" loss=loss(data, evomodel) accuracy = accuracy(data, evomodel)

function loss(model)
    σee = model[1]
    pee = model[2]
    σei = model[3]
    pei = model[4]
    P1, C1 = make_net(Ne, Ni, σee = σee, pee = pee, σei = σei, pei = pei)#,a=a)
    E1, I1 = P1
    SNN.monitor([E1, I1], [:fire])
    sim_length = 500
    @inbounds for t = 1:sim_length*ms
        E1.I = vec([11.5 for i = 1:sim_length])#vec(E_stim[t,:])#[i]#3randn(Ne)
        SNN.sim!(P1, C1, 1ms)
    end

    spkd_found = get_trains(P1[1])
    println("Ground Truth \n")
    SNN.raster([E]) |> display
    println("Best Candidate \n")

    SNN.raster([E1]) |> display

    error = raster_difference(spkd_ground, spkd_found)
    error
end
#loss(model) = (x,y)->logitcrossentropy(model(x), y)
#loss(model,x,y) = loss(model)(x, y)
#loss(xy, model) = loss(model)(hcat(map(first,xy)...), hcat(map(last,xy)...))

opt = ADAM(1e-4)
evalcb = Flux.throttle(() -> @show(loss(data, model), accuracy(data, model)), 5)
for i in 1:500
    Flux.train!(loss(model), params(model), data, opt, cb = evalcb)
end
f_nn = (x1,x2) -> Tracker.data(mod([x1;x2]))[1]
#
plot(x1,x2,f_nn,seriestype=:surface,color=:blues)
plot!(X[1,:],X[2,:],Y',seriestype=:scatter,markersize=MS1,marker=M1,markercolor=MC1,label="data")
plot!(xlabel="Weight [kg]", ylabel="Age [years]", zlabel="Blood fat content")
plot!(title="Humans: blood fat vs. weight, age",camera=(-60,20))


text:
module TestBinaryLGBM

using MLJBase
using Test
using Random: seed!

import CategoricalArrays
import LightGBM

## CLASSIFIER -- shamelessly copied from MLJModels/test/XGBoost.jl

model = LightGBM.MLJInterface.LGBMClassifier(objective="binary", num_iterations=100)

# test binary case:
N = 2
seed!(0)

X       = (x1=rand(1000), x2=rand(1000), x3=rand(1000))
ycat    = string.(mod.(round.(Int, X.x1 * 10), N)) |> MLJBase.categorical
weights = Float64.(MLJBase.int(ycat)) # just use the 1's/2's directly as multipliers

y = MLJBase.identity.(ycat) # make plain Vector with categ. elements (actually not sure what this is doing)

# fit once, without weights
train, test              = MLJBase.partition(MLJBase.eachindex(y), 0.6)
fitresult, cache, report = MLJBase.fit(model, 0, MLJBase.selectrows(X, train), y[train];)
yhat                     = MLJBase.mode.(MLJBase.predict(model, fitresult, MLJBase.selectrows(X, test)))

# fit again with weights
fitresult, cache, report = MLJBase.fit(model, 0, MLJBase.selectrows(X, train), y[train], weights[train])
yhat_with_weights        = MLJBase.mode.(MLJBase.predict(model, fitresult, MLJBase.selectrows(X, test)))

misclassification_rate   = sum(yhat .!= y[test])/length(test)

# Well, although XGBoost gets misclassification below 0.01, LightGBM can't do it with the default settings ...
# It gets to exactly 0.01...
@test misclassification_rate < 0.015

# All we can really say about fitting with/without weights for this example is that the solutions shouldn't be identical
@test yhat_with_weights != yhat

# Cache contains iterations counts history
@test cache isa NamedTuple
@test cache.num_boostings_done == [100]

@test isa(report, NamedTuple)

expected_return_type = Tuple{
    LightGBM.LGBMClassification,
    CategoricalArrays.CategoricalArray,
    LightGBM.MLJInterface.LGBMClassifier,
}

@test isa(fitresult, expected_return_type)

# here we test the bit where we try to fit a nonbinary using a binary estimator -- we expect this to throw
crazyN = 5 # testing nonbinary
ycrazy = string.(mod.(round.(Int, X.x1 * 10), crazyN)) |> MLJBase.categorical
@test_throws ArgumentError MLJBase.fit(model, 0, X, ycrazy)

end # module
true


text:
ispath("test") && cd("test")
include("test_func_defs.jl")

testplots() = testplot("No frills           1")
function testplots(pref::String)
    function fna(n)
        if Test.get_testset() isa Test.FallbackTestSet
            strip(pref) * " temp " * string(n)
        else
            strip(pref) * Test.get_testset().description * " " *string(n)
        end
    end
    function _png(n)
        title!(fna(n))
        isfile(fna(n)) && println("Writing over $(fna(n))")
        png(fna(n))
        # uncomment for interactive use:
        #current()
    end
    _png() = _png("")
    @testset verbose = true "-1 One series" begin
        @testset "-1-1 One axis" begin
            @testset "-1-1-1 Unitless" begin
                plot(vr)        ;_png()
                plot(f_r_r)      ;_png(1)
                plot(f_r_r, xlims = (-2, 2))           ;_png(2)
            @test true;end
            @testset "-1-1-2 Quantity" begin
                plot(vq, guide = "Units are postfixed:")           ; _png()
                plot(f_r_q)           ; _png(1)
                plot(x-> f_q_q(x∙s))           ;_png(2)
                plot(f_q_q, xlims = (-2, 2)s)           ;_png(3)
            @test true;end
        end
        @testset "-1-2 Two axis" begin
            @testset "-1-2-1 Unitless" begin
                plot(s1x, s1y)           ; _png()
                plot(s1x, f_r_r) ; _png(1)
                plot(f_r_r, s1x) ;_png(2)
                p = plot(s1x, s1y, (x,y ) -> f_r_r(x - y)) # f creates a Surface, not shown
               _png(3)
                @test p[1][1][:x][1] == s1x[1]
                @test p[1][1][:y][1]== s1y[1]
            end
            @testset "-1-2-2 Quantity" begin
                plot(s1x∙m, s1y∙s, xguide = "Units are postfixed:", yguide ="y axis") ; _png()
                plot(s1x∙s, f_q_q) ; _png(1)
                plot(f_q_q, s1x∙s) ;_png(2)
                # Surfaces are type-checked to be Float64. User can create the surfaces on their own.
                @test_throws TypeError plot(s1x∙s, s1y∙s, (x,y ) -> f_q_q(x - y))
            end
        end
        @testset "-1-3 Three axis" begin
            @testset "-1-3-1 Unitless" begin
                plot(s1x, s1y, s1z) ; _png()
                p = plot(f_r_r, s1y, s1z) # f(y)
                _png(1)
                @test p[1][1][:x][1] == f_r_r(s1y[1])
                p = plot(s1x, f_r_r, s1z) # f(x)
               _png(2)
                @test p[1][1][:y][1] == f_r_r(s1x[1])
            end
            @testset "-1-3-2 Quantity" begin
                plot(s1x∙m, s1y∙s, s1z∙s) ; _png()
                p = plot(f_q_q, s1y∙s, s1z∙s) # f(y)
                _png(1)
                @test p[1][1][:x][1]N == f_q_q(s1y[1]∙s)
               _png(2)
                p = plot(s1x∙s, f_q_q, s1z∙s)    # f(x)
               _png(3)
                @test p[1][1][:y][1]N == f_q_q(s1x[1]∙s)
            end
        end
    end
    @testset verbose = true "-2 Two series" begin
        @testset "-2-1 One axis" begin
            @testset "-2-1-1 Unitless" begin
                plot(mxr) ; _png()
                # Lenient syntax
                plot([f_r_r  x->f_r_r(2x)], xlims =(-2,2)) ; _png(1)
                # Same plot
                plot([f_r_r, x->f_r_r(2x)], xlims =(-2,2)) ;_png(2)

            @test true;end
            @testset "-2-1-2 Quantity" begin
                plot(mxq) ;_png()
                # Lenient syntax
                plot([f_q_q  x->f_q_q(2x)], xlims =(-2,2)s, ylims = (-10,5)N) ; _png(1)
                # Same plot
                plot([f_q_q, x->f_q_q(2x)], xlims =(-2,2)s, ylims = (-10,5)N) ;_png(2)

            @test true;end
        end
        @testset "-2-2 Two axis" begin
            @testset "-2-2-1 Unitless" begin
                plot(mxr, myr) ; _png()
                plot(mxr, [f_r_r, x-> 0.8f_r_r(1.25x)]) ;_png(2)
                plot([f_r_r , x-> 0.8f_r_r(1.25x)], myr) ;_png(3)
                p = plot(hcat(1:0.1:2, 5:0.1:6), [f_r_r, x-> 0.8f_r_r(1.25x)])
               _png(3)
                @test p[1][1][:y][1] == f_r_r(1)
                @test p[1][2][:y][1] == 0.8f_r_r(1.25 * 5)
            end
            @testset "-2-2-2 Quantity" begin
                plot(mxq, myq) ; _png()
                plot(myq, [f_q_q, x-> 0.8f_q_q(1.25x)]) ; _png(1)
                plot([f_q_q, x-> 0.8f_q_q(1.25x)], myq) ;_png(2)
                # We can't dispatch on this without also dispatching on unitless function mappings, which we won't do.
                # An alternative would be type piracy on function expand_extrema!, which would be so bad.
                @test_throws DimensionError p = plot(hcat(1:0.1:2, 5:0.1:6), [f_r_q, x-> 0.8f_r_q(1.25x)], xlims = (0,7))
            end
            @testset "-2-2-3 Quantity mixed units" begin
                plot(mxqm, myqm) ; _png()
                plot(myq, [f_q_q, x-> 0.8f_q_q(1.25x)∙m]) ; _png(1)
                plot([f_q_q, x-> 0.8f_q_q(1.25x)∙m], myq) ;_png(2)
                plot([x->f_q_q(x∙s/m), x-> 0.8f_q_q(1.25x∙s/m²)∙m], mxqm) ;_png(3)
                plot([x->f_q_q(x∙s/m)  x-> 0.8f_q_q(1.25x∙s/m²)∙m], mxqm,
                    xguide= "Load",
                    yguide = "Response",
                    xlims = (0, 1.1)m,
                    ylims = (0,10)m,
                    framestyle = :origin) ;_png(4)
            end
        end
        @testset "-2-3 Three axis" begin
            @testset "-2-3-1 Unitless" begin
                plot(s1x, s1y, s1z) ; _png()
                p = plot(f_r_r, s1y, s1z) # f(y)
                _png(1)
                @test p[1][1][:x][1] == f_r_r(s1y[1])
                p = plot(s1x, f_r_r, s1z) # f(x)
               _png(2)
                @test p[1][1][:y][1] == f_r_r(s1x[1])
            end
            @testset "-2-3-2 Quantity" begin
                plot(s1x∙m, s1y∙s, s1z∙s) ; _png()
                p = plot(f_q_q, s1y∙s, s1z∙s) # f(y)
                _png(1)
                @test p[1][1][:x][1]N == f_q_q(s1y[1]∙s)
                p = plot(s1x∙s, f_q_q, s1z∙s)    # f(x)
               _png(2)
                @test p[1][1][:y][1]N == f_q_q(s1x[1]∙s)
                # Surfaces are type-checked to be Float64. User can create the surfaces on their own.
                @test_throws TypeError plot(s1x∙s, s1y∙s, (x, y) -> f_q_q(x - y))
            end
        end
    end
end

testplots_ribbons() = testplots_ribbon("Ribbon           1")
function testplots_ribbons(pref::String)
    function fna(n)
        if Test.get_testset() isa Test.FallbackTestSet
            strip(pref) * " temp " * string(n)
        else
            strip(pref) * Test.get_testset().description * " " *string(n)
        end
    end
    function _png(n)
        title!(fna(n))
        isfile(fna(n)) && println("Writing over $(fna(n))")
        png(fna(n))
        # uncomment for interactive use:
        current()
    end
    _png() = _png("")

    @testset verbose = true "-1 One series" begin
        @testset "-1-1 One axis" begin
            @testset "-1-1-1 Unitless" begin
                plot(vr, ribbon = 2)        ; _png()
                plot(f_r_r, ribbon = 2)      ; _png(1)
                plot(f_r_r, ribbon = x-> 20 + 2f_r_r(2x / 3)) ; _png(2)
                plot(vr, ribbon = -0.5:0.5:5)  ; _png(3)
                plot(f_r_r, ribbon = -0.5:0.5:5)  ; _png(4)
                plot(f_r_r, ribbon = (5:1:50, 5:-0.1:0.5))  ; _png(5)
            @test true;end
            @testset "-1-1-2 Quantity" begin
                plot(vq, ribbon = 2s)           ; _png()
                plot(f_r_q, ribbon = 2N)           ; _png(1)
                plot(f_r_q, ribbon = x-> 20N + 2f_r_q(2x / 3)); _png(2)
                plot(vq, ribbon = (-0.5:0.5:5)s)  ; _png(3)
                plot(f_r_q, ribbon = (-0.5:0.5:5)N)  ; _png(4)
                plot(f_r_q, ribbon = (5:1:50, 5:-0.1:0.5)N)  ; _png(5)
            @test true;end
        end
        @testset "-1-2 Two axis" begin
            @testset "-1-2-1 Unitless" begin
                plot(s1x, s1y, ribbon = 2)   ; _png()
                plot(s1x, f_r_r, ribbon = 2) ; _png(1)
                plot(f_r_r, s1x, ribbon = x-> 20 + 2f_r_r(2x / 3)) ; _png(2)
                plot(s1x, s1y, ribbon = -0.5:1.5:5) ; _png(3)
                plot(s1x, f_r_r, ribbon = -0.5:1.5:5) ; _png(4)
                plot(f_r_r, s1x, ribbon = (5:10:50, 5:-1:0.5)) ; _png(5)
                plot(xlims =(-5, 5), f_r_r, ribbon = x-> 1.25 * f_r_r(x)); _png(6)
            @test true;end
            @testset "-1-2-2 Quantity" begin
                plot(s1x∙m, s1y∙m, ribbon = 2.0m)   ; _png()
                plot(s1x∙s, f_q_q, ribbon = 2N) ; _png(1)
                plot(f_q_q, s1x∙s, ribbon = x-> 20N + 2f_q_q(2x / 3)) ; _png(2)
                plot(s1x∙s, s1y∙N, ribbon = (-0.5:1.5:5)N) ; _png(3)
                plot(s1x∙s, f_q_q, ribbon = (-0.5:1.5:5)N) ; _png(4)
                plot(f_q_q, s1x∙s, ribbon = (5:10:50, 5:-1:0.5)N) ; _png(5)
                plot(xlims =(-5s, 5s), f_q_q, ribbon = x-> 1.25 * f_q_q(x)); _png(6)
            @test true;end
        end
    end
    @testset verbose = true "-2 Two series" begin
        @testset "-2-1 One axis" begin
            @testset "-2-1-1 Unitless" begin
                plot(mxr, ribbon = [2 1]) ; _png()
                plot(mxr, ribbon = [2 1; 1 0.5 ; 0.5 0.25; 0.03 0.1]) ; _png(1)
                plot([f_r_r  x-> f_r_r(0.5x)], ribbon = [-1.5:0.15:1.15  -3.0:0.30:2.3 ]) ; _png(2)
                plot([f_r_r  x-> f_r_r(0.5x)], ribbon = ([-15:1.5:11.5  -3.0:0.30:2.3 ], [-1.5:0.15:1.15  -15:1.5:11.5  ])) ; _png(3)
                plot([f_r_r  x-> f_r_r(0.5x)], ribbon = x-> 0.2f_r_r(x)) ; _png(4)
                plot([f_r_r, x-> f_r_r(0.5x)], ribbon = [x-> 0.1f_r_r(x), x-> 0.5f_r_r(x)]) ; _png(5)
            @test true;end
            @testset "-2-1-2 Quantity" begin
                plot(mxq, ribbon = [2.0 1.0]m) ; _png()
                plot(mxq, ribbon = [2 1; 1 0.5 ; 0.5 0.25; 0.03 0.1]m) ; _png(1)
                # With unknown abscissa, it is harder to implement unit conversions without
                # also affecting other types of plots.
            @test true;end
        end
        @testset "-2-2 Two axis" begin
            @testset "-2-2-1 Unitless" begin
                plot(mxr, myr, ribbon = [2 1]) ; _png()
                plot(mxr, [f_r_r, x-> 0.8f_r_r(1.25x)], ribbon = [2 1; 1 0.5 ; 0.5 0.25; 0.03 0.1]) ; _png(1)
                plot([f_r_r , x-> 0.8f_r_r(1.25x)], myr, ribbon = [1.5:-0.55:-1.15  15:-5.5:-11.5 ]) ; _png(2)
                p = plot(hcat(1:0.1:3, 3:0.1:5),
                    [f_r_r, x-> 0.8f_r_r(1.25x)],
                    ribbon = ([s_1_down  s_2_down], [s_1_up  s_2_up]),
                    labels = ["s_1" "s_2"])
                _png(3)
                @test p[1][1][:y][1] == f_r_r(1)
                @test p[1][2][:y][1] == 0.8f_r_r(1.25 * 3)
                plot( xlims =(-5, 5), [f_r_r, x-> 0.5f_r_r(1.5x)],
                        ribbon = [x-> f_r_r(0.25x), x-> f_r_r(0.5x)]); _png(4)
            end
            @testset "-2-2-2 Quantity" begin
                plot(mxq, myq, ribbon = [2 1]s) ; _png()
                plot(mxq, [x->f_q_q(x∙s/m), x-> 0.8f_q_q(1.25x∙s/m)], ribbon = [2 1; 1 0.5 ; 0.5 0.25; 0.03 0.1]N) ; _png(1)
                plot([x-> f_q_q(x) , x-> 0.8f_q_q(1.25x)], myq, ribbon = [1.5:-0.55:-1.15  15:-5.5:-11.5 ]N); _png(2)
                p = plot(hcat(1:0.1:3, 3:0.1:5)s,
                    [f_q_q, x-> 0.8f_q_q(1.25x)],
                    ribbon = ([s_1_down  s_2_down]N, [s_1_up  s_2_up]N),
                    labels = ["s_1" "s_2"])
                _png(3)
                @test p[1][1][:y][1]N == f_q_q(1s)
                @test p[1][2][:y][1]N == 0.8f_q_q(1.25 * 3s)
                # We can't dispatch on this without also dispatching on unitless function mappings, which we won't do.
                @test_throws DimensionError p = plot(hcat(1:0.1:2, 5:0.1:6), [f_r_q, x-> 0.8f_r_q(1.25x)], xlims = (0,7))
                plot( xlims =(-5s, 5s), [f_q_q, x-> 0.5f_q_q(1.5x)],
                     ribbon = [x-> f_q_q(0.25x), x-> f_q_q(0.5x)]); _png(4)
            end
            @testset "-2-2-3 Quantity mixed units" begin
                plot(mxqm, myqm, ribbon = [2s 1s²], xguide = "Several units, postfixed:", yguide = "Other units, postfixed:") ; _png()
                plot(myq, [f_q_q, x-> 0.8f_q_q(1.25x)∙m], ribbon = [2N 1N∙m]) ; _png(1)
                plot([f_q_q, x-> 0.8f_q_q(1.25x)∙m], myq, ribbon = [2N 1N∙m]) ; _png(2)
                p = plot(hcat(1:0.1:3, 3:0.1:5)s,
                    [f_q_q, x-> 0.8f_q_q(1.25x)],
                    ribbon = ([s_1_down  s_2_down]N, [s_1_up  s_2_up]N),
                    labels = ["s_1" "s_2"])
                @test p[1][1][:y][1]N == f_q_q(1s)
                @test p[1][2][:y][1]N == 0.8f_q_q(1.25 * 3s)
                plot(myq, [f_q_q, x-> 0.8f_q_q(1.25x)∙m],  ribbon = [2N 1N∙m]) ; _png(2)
                plot([x->f_q_q(x∙s/m), x-> 0.8f_q_q(1.25x∙s/m²)s], mxqm,
                    ribbon = ([s_1_down∙N s_2_down∙N∙s], [s_1_up∙N  s_2_up∙N∙s])) ;_png(3)
                plot( xlims =(-5s, 5s), [f_q_q, x-> 0.5f_q_q(1.5x)∙m],
                    ribbon = [x-> f_q_q(0.25x), x-> 0.5∙f_q_q(0.5x)∙m],
                    yguide = "Load", xguide = "Time"); _png(4)
                plot( xlims =(-5s, 5s),
                     [f_q_q  x-> 0.5f_q_q(1.5x)∙m],
                    ribbon = ([s_1_down∙N  s_2_down∙N∙m], [s_1_up∙N   s_2_up∙N∙m]),
                    yguide = "Load", xguide = "Time"); _png(5)
            end
        end
    end
end

@testset verbose = true "No frills           1" begin
    testplots(Test.get_testset().description)
end
@testset verbose = true "Ribbons             2" begin
    testplots_ribbons(Test.get_testset().description)
end
# Further tests: Set additional parameters with Plots.default(...)