text:
const NUCLEOTIDES = ['C','G','A','T']

function count_nucleotides(strand::AbstractString)
    counts = zeros(Int,length(NUCLEOTIDES))
    l = Dict(v=>i for (i,v) in enumerate(NUCLEOTIDES))
    try
        for s in strand
            @inbounds counts[l[s]] += 1
        end
    catch(e)
        # If s is not a valid nucleotide
        isa(e,KeyError) && throw(DomainError(e))
        throw(e)
    end
    Dict(k=>counts[v] for (k,v) in l)
end

# This was my first attempt
# It doesn't throw domain errors and I couldn't work out how to do so without slowing it down a lot
#
# It's interesting that it's so fast even though it iterates through the list length(NUCLEOTIDES) times
#
# For large strings it's about 25% faster than the array based approach (thanks cmcaine for some optimisations)
function count_nucleotides_nosanitise(strand::AbstractString)
    Dict(n=>mapreduce(==(n),+,strand; init=0) for n in NUCLEOTIDES)
#   Dict(n=>sum(==(n),strand) for n in NUCLEOTIDES) # equivalent but less readable IMO
end

# ~#-real-cores speedup over the single-threaded one above (ta to cmcaine for pointing out codeunits for cheap conversion to array)
import ThreadsX
function count_nucleotides_threads(strand::AbstractString)
    Dict(Char(n)=>ThreadsX.mapreduce(==(n),+,codeunits(strand); init=0) for n in codeunits(join(NUCLEOTIDES)))
end

# GPU version
using CuArrays
# Getting data to the GPU is slow; this is approx 4x faster than the single threaded method
count_nucleotides_gpu(strand::String) = count_nucleotides_gpu(CuArray(codeunits(strand)))

# But once it's there we get a further 15x speedup (with a lowly 1080Ti)
count_nucleotides_gpu(strand::CuArray) = Dict(Char(n) => sum(==(n),strand) for n in codepoint.(NUCLEOTIDES))

# For testing
# using Random: randstring
# const s = randstring(NUCLEOTIDES,1_000_000_000);
# const c = codeunits(s);
# const d = CuArray(c);
# using BenchmarkTools
# @benchmark count_nucleotides(s)
# @benchmark count_nucleotides_threads(s)
# @benchmark count_nucleotides_gpu(d)


text:
using TheAlgorithms
using Test

using DataFrames
using DifferentialEquations
using GLM
using LinearAlgebra
using Plots
using Random

@testset "TheAlgorithms" begin

include("math.jl")
include("matrix.jl")
include("project-rosalind.jl")
include("searches.jl")
include("statistics.jl")

end

ext: jl

max_stars_repo_name: madisonestabrook/Julia

max_issues_repo_name: madisonestabrook/Julia

max_forks_repo_name: madisonestabrook/Julia


text:
struct EventObservations{S <: DiseaseStateSequence, M <: ILM}
  infection::Vector{Float64}
  removal::Union{Nothing, Vector{Float64}}
  seq::Union{Nothing, Vector{Union{Nothing, GeneticSeq}}}

  function EventObservations{S}(i::V,
                                r::V) where {
                                V <: Vector{Float64},
                                S <: Union{SEIR, SIR}}
    if length(i) != length(r)
      throw(ErrorException("Length of infection and removal times must be equal"))
    end
    return new{S, TNILM}(i, r, nothing)
  end

  function EventObservations{S}(i::V) where {
                                V <: Vector{Float64},
                                S <: Union{SEI, SI}}
    return new{S, TNILM}(i, nothing, nothing)
  end

  function EventObservations{S}(i::V,
                                r::V,
                                seq::VG) where {
                                V <: Vector{Float64},
                                G <: GeneticSeq,
                                VG <: Vector{Union{Nothing, G}},
                                S <: Union{SEIR, SIR}}
    if length(unique([length(i); length(r); length(seq)])) != 1
      throw(ErrorException("Length of infection times, removal times, and genetic sequence vectors must be equal"))
    end
    return new{S, PhyloILM}(i, r, seq)
  end

  function EventObservations{S}(i::V,
                                   s::VG) where {
                                   V <: Vector{Float64},
                                   G <: GeneticSeq,
                                   VG <: Vector{Union{Nothing, G}},
                                   S <: Union{SEI, SI}}
    if length(i) != length(s)
      throw(ErrorException("Length of infection time and genetic sequence vectors must be equal"))
    end
    return new{S, PhyloILM}(i, nothing, s)
  end
end

function EventObservations{S}(i::Array{Float64, 2}) where {
                              S <: Union{SEI, SI}}
  if size(i, 2) != 1
    throw(ErrorException("Invalid Array dimensions for observations of a $S model"))
  end
  return EventObservations{S}(i[:,1])
end

function EventObservations{S}(ir::Array{Float64, 2}) where {
                              S <: Union{SEIR, SIR}}
  if size(ir, 2) != 2
    throw(ErrorException("Invalid Array dimensions for observations of a $S model"))
  end
  return EventObservations{S}(ir[:, 1], ir[:, 2])
end

function EventObservations{S}(i::Array{Float64, 2},
                              s::VG) where {
                              G <: GeneticSeq,
                              VG <: Vector{Union{Nothing, G}},
                              S <: Union{SEI, SI}}
  if size(i, 2) != 1
    throw(ErrorException("Invalid Array dimensions for observations of a $S model"))
  end
  return EventObservations{S}(i[:,1], s)
end

function EventObservations{S}(ir::Array{Float64, 2},
                              s::VG) where {
                              G <: GeneticSeq,
                              VG <: Vector{Union{Nothing, G}},
                              S <: Union{SEIR, SIR}}
  if size(ir, 2) != 2
    throw(ErrorException("Invalid Array dimensions for observations of a $S model"))
  end
  return EventObservations{S}(ir[:, 1], ir[:, 2], s)
end

function individuals(x::EventObservations{S, M}) where {
                     S <: DiseaseStateSequence,
                     M <: ILM}
  return length(x.infection)
end

function Base.show(io::IO,
                   x::EventObservations{S, M}) where {
                   S <: DiseaseStateSequence,
                   M <: ILM}
  return print(io, "$S $M observations (n=$(individuals(x)))")
end

function Base.getindex(x::EventObservations{S, M},
                       new_state::DiseaseState) where {
                       S <: DiseaseStateSequence,
                       M <: ILM}
  if new_state == State_I
    return x.infection
  elseif (new_state == State_R) && (State_R ∈ S)
    return x.removal
  else
    error("Unrecognized indexing disease state")
  end
end

function Base.getindex(x::EventObservations{S, M},
                       states::DiseaseStates) where {
                       S <: DiseaseStateSequence,
                       M <: ILM}
  y = x[states[1]]
  for i = 2:length(states)
    y = hcat(y, x[states[i]])
  end
  return y
end

function Base.convert(::Type{Array{Float64, 2}}, x::EventObservations{S, M}) where {
                      S <: DiseaseStateSequence,
                      M <: ILM}
  states = State_R ∉ S ? [State_I] : [State_I, State_R]
  return x[states]
end

function Base.convert(::Type{Vector{Float64}}, x::EventObservations{S, M}) where {
                      S <: DiseaseStateSequence,
                      M <: ILM}
  states = State_R ∉ S ? [State_I] : [State_I, State_R]
  return x[states][:]
end

function Base.minimum(x::EventObservations{S, M}) where {
                      S <: DiseaseStateSequence,
                      M <: ILM}
  y = convert(Array{Float64, 2}, x)
  return minimum(y[.!isnan.(y)])
end

function Base.maximum(x::EventObservations{S, M}) where {
                      S <: DiseaseStateSequence,
                      M <: ILM}
  y = convert(Array{Float64, 2}, x)
  return maximum(y[.!isnan.(y)])
end


text:
export current_device, global_queue, synchronize, device_synchronize

"""
    current_device()::MtlDevice

Return the Metal GPU device associated with the current Julia task.

Since all M-series systems currently only externally show a single GPU, this function
effectively returns the only system GPU.
"""
function current_device()
    get!(task_local_storage(), :MtlDevice) do
        MtlDevice(1)
    end::MtlDevice
end

const global_queues = WeakKeyDict{MtlCommandQueue,Nothing}()

"""
    global_queue(dev::MtlDevice)::MtlCommandQueue

Return the Metal command queue associated with the current Julia thread.
"""
function global_queue(dev::MtlDevice)
    get!(task_local_storage(), (:MtlCommandQueue, dev)) do
        queue = MtlCommandQueue(dev)
        global_queues[queue] = nothing
        queue
    end::MtlCommandQueue
end

# TODO: Increase performance (currently ~15us)
"""
    synchronize(queue)

Wait for currently committed GPU work on this queue to finish.

Create a new MtlCommandBuffer from the global command queue, commit it to the queue,
and simply wait for it to be completed. Since command buffers *should* execute in a
First-In-First-Out manner, this synchronizes the GPU.
"""
function synchronize(queue::MtlCommandQueue=global_queue(current_device()))
    cmdbuf = MtlCommandBuffer(queue)
    commit!(cmdbuf)
    wait_completed(cmdbuf)
end

function device_synchronize()
    for queue in keys(global_queues)
        synchronize(queue)
    end
end


text:
__precompile__(true)

module LearnBase

# Only reexport required functions by default
import StatsBase: nobs, fit, fit!

# We temporary reexport issymmetric for smooth
# transition into 0.5
if VERSION >= v"0.5-"
    import Base.issymmetric
else
    const issymmetric = Base.issym
end

"""
Baseclass for any kind of cost. Notable examples for
costs are `Loss` and `Penalty`.
"""
abstract Cost

"""
Baseclass for all losses. A loss is some (possibly simplified)
function `L(features, targets, outputs)`, where `outputs` are the
result of some function `f(features)`.
"""
abstract Loss <: Cost

"""
A loss is considered **supervised**, if all the information needed
to compute `L(features, targets, outputs)` are contained in
`targets` and `outputs`, and thus allows for the simplification
`L(targets, outputs)`.
"""
abstract SupervisedLoss <: Loss

"""
A supervised loss, where the targets are in {-1, 1}, and which
can be simplified to `L(targets, outputs) = L(targets * outputs)`
is considered **margin-based**.
"""
abstract MarginLoss <: SupervisedLoss

"""
A supervised loss that can be simplified to
`L(targets, outputs) = L(targets - outputs)` is considered
distance-based.
"""
abstract DistanceLoss <: SupervisedLoss

"""
A loss is considered **usupervised**, if all the information needed
to compute `L(features, targets, outputs)` are contained in
`features` and `outputs`, and thus allows for the simplification
`L(features, outputs)`.
"""
abstract UnsupervisedLoss <: Loss

abstract Penalty <: Cost

function value end
function value! end
function meanvalue end
function sumvalue end
function meanderiv end
function sumderiv end
function deriv end
function deriv2 end
function deriv! end
function value_deriv end
function grad end
function grad! end
function addgrad! end
function value_grad end
function value_grad! end
function prox end
function prox! end

function isminimizable end
function isdifferentiable end
function istwicedifferentiable end
function isconvex end
function isstrictlyconvex end
function isstronglyconvex end
function isnemitski end
function isunivfishercons end
function isfishercons end
function islipschitzcont end
function islocallylipschitzcont end
function islipschitzcont_deriv end # maybe overkill
function isclipable end
function ismarginbased end
function isclasscalibrated end
function isdistancebased end

# fallback to supersets
isstrictlyconvex(any) = isstronglyconvex(any)
isconvex(any) = isstrictlyconvex(any)
islocallylipschitzcont(any) = islipschitzcont(any)

"""
Anything that takes an input and performs some kind
of function to produce an output. For example a linear
prediction function.
"""
abstract Transformation
abstract StochasticTransformation <: Transformation

function transform end
function transform! end

"""
Baseclass for any prediction model that can be minimized.
This means that an object of a subclass contains all the
information needed to compute its own current loss.
"""
abstract Minimizable <: Transformation

function getobs end
function getobs! end

function update end
function update! end
function learn end
function learn! end

import Base: AbstractSet

"A continuous range (inclusive) between a lo and a hi"
immutable IntervalSet{T<:Number} <: AbstractSet
    lo::T
    hi::T
end
function IntervalSet{A<:Number,B<:Number}(lo::A, hi::B)
    T = promote_type(A,B)
    IntervalSet{T}(convert(T,lo), convert(T,hi))
end
randtype(s::IntervalSet) = Float64
Base.rand{T<:Number}(s::IntervalSet{T}, dims::Integer...) = rand(dims...) * (s.hi - s.lo) + s.lo
Base.in{T<:Number}(x::Number, s::IntervalSet{T}) = s.lo <= x <= s.hi

"Set of discrete items"
immutable DiscreteSet{T<:AbstractArray} <: AbstractSet
    items::T
end
randtype(s::DiscreteSet) = eltype(s.items)
Base.rand(s::DiscreteSet, dims::Integer...) = rand(s.items, dims...)
Base.in(x, s::DiscreteSet) = x in s.items
Base.length(s::DiscreteSet) = length(s.items)
Base.getindex(s::DiscreteSet, i::Int) = s.items[i]

# operations on arrays of sets
randtype{S<:AbstractSet,N}(sets::AbstractArray{S,N}) = Array{promote_type(map(randtype, sets)...), N}
Base.rand{S<:AbstractSet}(sets::AbstractArray{S}) = eltype(randtype(sets))[rand(s) for s in sets]
function Base.rand{S<:AbstractSet}(sets::AbstractArray{S}, dim1::Integer, dims::Integer...)
    A = Array(randtype(sets), dim1, dims...)
    for i in eachindex(A)
        A[i] = rand(sets)
    end
    A
end
function Base.in{S<:AbstractSet}(xs::AbstractArray, sets::AbstractArray{S})
    size(xs) == size(sets) && all(map(in, xs, sets))
end

"Groups several heterogenous sets. Used mainly for proper dispatch."
immutable TupleSet{T<:Tuple} <: AbstractSet
    sets::T
end
TupleSet(sets::AbstractSet...) = TupleSet(sets)

# rand can return arrays or tuples, but defaults to arrays
randtype(sets::TupleSet, ::Type{Vector}) = Vector{promote_type(map(randtype, sets.sets)...)}
Base.rand(sets::TupleSet, ::Type{Vector}) = eltype(randtype(sets, Vector))[rand(s) for s in sets.sets]
randtype(sets::TupleSet, ::Type{Tuple}) = Tuple{map(randtype, sets.sets)...}
Base.rand(sets::TupleSet, ::Type{Tuple}) = map(rand, sets.sets)
function Base.rand{OT}(sets::TupleSet, ::Type{OT}, dim1::Integer, dims::Integer...)
    A = Array(randtype(sets, OT), dim1, dims...)
    for i in eachindex(A)
        A[i] = rand(sets, OT)
    end
    A
end

randtype(sets::TupleSet) = randtype(sets, Vector)
Base.rand(sets::TupleSet, dims::Integer...) = rand(sets, Vector, dims...)
Base.in(x, sets::TupleSet) = all(map(in, x, sets.sets))

"Returns an AbstractSet representing valid input values"
function inputdomain end

"Returns an AbstractSet representing valid output/target values"
function targetdomain end

export

    # Types
    Cost,
        Loss,
            SupervisedLoss,
                MarginLoss,
                DistanceLoss,
            UnsupervisedLoss,
        Penalty,

    Transformation,
        StochasticTransformation,

    Minimizable,

    AbstractSet,
        IntervalSet,
        DiscreteSet,
        TupleSet,

    # Functions
    getobs,
    getobs!,

    learn,
    learn!,
    update,
    update!,
    transform,
    transform!,
    value,
    value!,
    meanvalue,
    sumvalue,
    meanderiv,
    sumderiv,
    deriv,
    deriv!,
    grad,
    grad!,
    addgrad!,
    deriv2,
    value_deriv,
    value_deriv!,
    value_grad,
    value_grad!,
    prox,
    prox!,
    inputdomain,
    targetdomain,

    isminimizable,
    isdifferentiable,
    istwicedifferentiable,
    isconvex,
    isstrictlyconvex,
    isstronglyconvex,
    isnemitski,
    isunivfishercons,
    isfishercons,
    islipschitzcont,
    islocallylipschitzcont,
    islipschitzcont_deriv,
    isclipable,
    ismarginbased,
    isclasscalibrated,
    isdistancebased,

    # Base
    issymmetric,

    # StatsBase
    fit,
    fit!,
    nobs

end # module


text:
module ClimaCommsMPI

using ClimaComms
import ClimaComms.Neighbor

using KernelAbstractions
using MPI

struct MPICommsContext <: ClimaComms.AbstractCommsContext
    mpicomm::MPI.Comm
    neighbors::Vector{ClimaComms.Neighbor}
    recv_reqs::Vector{MPI.Request}
    send_reqs::Vector{MPI.Request}

    function MPICommsContext(
        neighbors::Vector{ClimaComms.Neighbor};
        mpicomm = MPI.COMM_WORLD,
    )
        nneighbors = length(neighbors)
        new(
            mpicomm,
            neighbors,
            fill(MPI.REQUEST_NULL, nneighbors),
            fill(MPI.REQUEST_NULL, nneighbors),
        )
    end
end

function Neighbor(
    ::Type{MPICommsContext},
    pid,
    AT,
    FT,
    send_dims,
    recv_dims = send_dims,
)
    if pid > 0
        if AT <: Array || singlebuffered(CC)
            kind = ClimaComms.SingleRelayBuffer
        else
            kind = ClimaComms.DoubleRelayBuffer
        end
        send_buf = ClimaComms.RelayBuffer{FT}(AT, kind, send_dims...)
        recv_buf = ClimaComms.RelayBuffer{FT}(AT, kind, recv_dims...)
        B = ClimaComms.RelayBuffer
    else
        send_buf = recv_buf = missing
        B = Missing
    end
    ClimaComms.Neighbor{B}(pid, send_buf, recv_buf)
end

function ClimaComms.init(::Type{MPICommsContext})
    if !MPI.Initialized()
        MPI.Init()
        atexit() do
            MPI.Finalize()
        end
    end

    return ClimaComms.mypid(MPICommsContext), ClimaComms.nprocs(MPICommsContext)
end

ClimaComms.mypid(::Type{MPICommsContext}) = MPI.Comm_rank(MPI.COMM_WORLD) + 1
ClimaComms.mypid(ctx::MPICommsContext) = MPI.Comm_rank(ctx.mpicomm) + 1
ClimaComms.iamroot(CC::Type{MPICommsContext}) = ClimaComms.mypid(CC) == 1
ClimaComms.nprocs(::Type{MPICommsContext}) = MPI.Comm_size(MPI.COMM_WORLD)
ClimaComms.nprocs(ctx::MPICommsContext) = MPI.Comm_size(ctx.mpicomm)
ClimaComms.singlebuffered(::Type{MPICommsContext}) = MPI.has_cuda()
ClimaComms.neighbors(ctx::MPICommsContext) = ctx.neighbors

function ClimaComms.start(ctx::MPICommsContext; dependencies = nothing)
    if !all(r -> r == MPI.REQUEST_NULL, ctx.recv_reqs)
        error("Must finish() before next start()")
    end

    progress = () -> iprobe_and_yield(ctx.mpicomm)
    nneighbors = length(ctx.neighbors)

    # start moving staged send data to transfer buffers
    events = ntuple(nneighbors) do n
        ClimaComms.prepare_transfer!(
            ctx.neighbors[n].send_buf,
            dependencies = dependencies,
            progress = progress,
        )
    end

    # post receives
    for n in 1:nneighbors
        tbuf = ClimaComms.get_transfer(ClimaComms.recv_buffer(ctx.neighbors[n]))
        ctx.recv_reqs[n] = MPI.Irecv!(
            tbuf,
            ClimaComms.pid(ctx.neighbors[n]) - 1,
            666,
            ctx.mpicomm,
        )
    end

    # wait for send data to reach the transfer buffers
    wait(CPU(), MultiEvent(events), progress)

    # post sends
    for n in 1:nneighbors
        tbuf = ClimaComms.get_transfer(ClimaComms.send_buffer(ctx.neighbors[n]))
        ctx.send_reqs[n] = MPI.Isend(
            tbuf,
            ClimaComms.pid(ctx.neighbors[n]) - 1,
            666,
            ctx.mpicomm,
        )
    end
end

function ClimaComms.progress(ctx::MPICommsContext)
    MPI.Iprobe(MPI.MPI_ANY_SOURCE, MPI.MPI_ANY_TAG, ctx.mpicomm)
end

function ClimaComms.finish(ctx::MPICommsContext; dependencies = nothing)
    if any(r -> r == MPI.REQUEST_NULL, ctx.recv_reqs)
        error("Must start() before finish()")
    end

    progress = () -> iprobe_and_yield(ctx.mpicomm)
    nneighbors = length(ctx.neighbors)

    # wait on previous receives
    MPI.Waitall!(ctx.recv_reqs)

    # move received data to stage buffers
    events = ntuple(nneighbors) do n
        ClimaComms.prepare_stage!(
            ClimaComms.recv_buffer(ctx.neighbors[n]);
            dependencies = dependencies,
            progress = progress,
        )
    end

    # ensure that sends have completed
    MPI.Waitall!(ctx.send_reqs)

    return MultiEvent(events)
end

function iprobe_and_yield(mpicomm)
    MPI.Iprobe(MPI.MPI_ANY_SOURCE, MPI.MPI_ANY_TAG, mpicomm)
    yield()
end

ClimaComms.barrier(ctx::MPICommsContext) = MPI.Barrier(ctx.mpicomm)

ClimaComms.reduce(ctx::MPICommsContext, val, op) =
    MPI.Reduce(val, op, 0, ctx.mpicomm)

ClimaComms.abort(ctx::MPICommsContext, status::Int) =
    MPI.Abort(ctx.mpicomm, status)

end # module


text:
using Documenter, JuDGE#, Literate

# for dir in ["examples"]
#     for file in sort(readdir(joinpath(@__DIR__, "src", dir)))
#         !endswith(file, ".jl") && continue
#         filename = joinpath(@__DIR__, "src", dir, file)
#         Literate.markdown(filename, dirname(filename); documenter=true)
#     end
# end
#
# const EXAMPLES = Any[
#     "examples/$(file)" for file  in filter(
#         f -> endswith(f, ".md"),
#         sort(readdir(joinpath(@__DIR__, "src", "examples")))
#     )
# ]

makedocs(sitename="JuDGE.jl: Julia Decomposition for Generalized Expansion",
        modules = [JuDGE],
        clean = true,
        format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"),
        pages = ["JuDGE" => "index.md",
                "Tutorials" => "tutorials.md",
                 "API Reference" => "api.md"])

deploydocs(repo = "github.com/reganbaucke/JuDGE.jl.git",devurl="docs")


text:
import SearchBall, ParticleSwarmOptimizer, DelimitedFiles, Profile, ProfileView
const sb = SearchBall
const pso = ParticleSwarmOptimizer
const df = DelimitedFiles

import Random
Random.seed!(1234);

myconfig = sb.get_config(["-f", "/home/tobias/doc/SearchBall/SearchBall.jl/configs/2dlinNN.config"])
function objective(x::AbstractVector{<:Number})
    # config = sb.get_config(["-f", "/home/tobias/doc/SearchBall/SearchBall.jl/configs/2dlinNN.config"]) # TODO: needs too long
    config = sb.Config(myconfig)
    strat = sb.get_strategy(config)
    strat.init_strategies(config, x=x, issimulation=true)
    myobjective = strat.objective
    dt = 1/30
    time_accu = 0.0
    while time_accu < 50.0
        sb.update_game!(config.ws, dt, config)
        time_accu += dt
        if config.global_communication["goal"]
            sb.free!(config)
            return time_accu + myobjective(config)
        end
        if config.global_communication["antigoal"] || config.global_communication["out"]
            sb.free!(config)
            return 70.0 + myobjective(config)
        end
    end
    sb.free!(config)
    return time_accu + myobjective(config)
end

objective_obj = pso.Objective(objective, [-7.2, 7.2], 12);

number_particle = 40;
neighbours = pso.LocalNeighbourhood(number_particle);
# neighbours = pso.HierachicalNeighbourhood(number_particle, 4);
# neighbours = pso.GlobalNeighbourhood(number_particle);

optimizer = pso.PSO(objective_obj, neighbours);

pso.optimize!(optimizer, 1);
println("everthing works!")
Profile.clear_malloc_data()
# @time pso.optimize!(optimizer, 300);
pso.optimize!(optimizer, 2900);
# pso.optimize!(optimizer, 600);
# Profile.@profile pso.optimize!(optimizer, 300);
# using Juno
# Juno.profiler()

# params, score = pso.getoptimum(optimizer);
# @show params, score
# df.writedlm("run01.txt", params)

# oldparams = vec(df.readdlm("run01.txt"))
# objective(oldparams)
# @time objective(oldparams)
# Profile.clear_malloc_data()
# Profile.clear()
# Profile.@profile objective(oldparams)

# ProfileView.view(C=true)


text:
include( "uses.jl")
include( "runner_libs.jl")
include( "table_libs.jl")
include( "dash_libs.jl")

const PREAMBLE = """
...

"""

function one_mt_box( name :: String, prefix :: String )
	dbc_row([
		dbc_col( dbc_label(name) ),
		dbc_col(
			dcc_radioitems(
			options =[
				Dict("label" => "Keep as-is", "value" => "$(prefix)_nc"),
				Dict("label" => "Abolish Completely", "value" => "$(prefix)_abolish"),
				Dict("label" => "Keep Housing Benefits only", "value" => "$(prefix)_asincome")
			],
			value="$(prefix)_nc") # radio
		) # col
	]) # row
end

function bi_ben_fieldset()
	bi = html_fieldset([
		html_legend( "Treatment of Existing Benefits"),
		one_mt_box( "Means Tested Benefits", "mtbs" ),

		dbc_row([
			dbc_col(
                dbc_label("Make BI Taxable?"; html_for="bi_tax")
			),
			dbc_col(
				dbc_checkbox( id="bi_tax", value=false)
			)

		]
		),
		dbc_row([
			dbc_col(
                dbc_label("Retain Sickness Benefits?"; html_for="sickben")
			),
			dbc_col(

			)
		])
    ])
	bi
end

"""
Create the block of sliders and radios on the LHS
"""
function get_input_block()
	return dbc_form([
		dbc_row([
			dbc_col( it_fieldset() ),
			dbc_col( bi_ben_fieldset() ),
			dbc_col( ben_fieldset() )
		])
		dbc_row([
			dbc_col([
				submit_button()
			]) # col
		]) # row
	]) # form
end

app = dash(
	external_stylesheets=[dbc_themes.UNITED],
	url_base_pathname="/basicincome/" )
# BOOTSTRAP|SIMPLEX|MINTY|COSMO|SANDSTONE|UNITED|SLATE|SOLAR|UNITED|

app.layout = dbc_container(fluid=true, className="p-5") do
	html_title( "Basic Income")
	html_h1("BI for Scotland"),
	dcc_markdown( PREAMBLE ),
	get_input_block(),
	dcc_loading(
		id="model_running",
		type="default",
		children = [
			html_div(
				id="loading-output-1"
			) # div
		] # children
	), # loading
	html_div( id="output-block" ),
	html_div( dcc_markdown( ENDNOTES ))
end # layout

function do_output( br, hr, tr, pa, ni_prim, ni_sec, cb, pen, uct, ucs, wtcb, scp )
	results = nothing
	sys = deepcopy( BASE_STATE.sys )
	# 21.15, 179.60, 10.0

	# 20, 41, 46, 12_570, 12, 13.8, 21.15, 179.60, 55, 324.84, 2_005, 10.0

	if (br != 20) || (hr !=41)||(tr !=46)||(uct != 55 )||(cb != 21.15)||(pen!= 179.60)||(scp!=10)||(pa!=12_570)||(ni_prim!=12)||(ni_sec!=13.8)||(ucs!=324.84)||(wtcb!=2_005)
		br /= 100.0
		hr /= 100.0
		tr /= 100.0
		uct /= 100.0
		pa /= WEEKS_PER_YEAR
		wtcb /= WEEKS_PER_YEAR
		ucs /= WEEKS_PER_MONTH

		ni_prim /= 100.0
		ni_sec /= 100.0

		if br == 0
			sys.it.non_savings_rates[1:3] .= 0.0
		else
			bincr = br-sys.it.non_savings_rates[2]
			sys.it.non_savings_rates[1:3] .+= bincr
			sys.it.non_savings_rates[1] = max(0, sys.it.non_savings_rates[1])
		end
		sys.it.non_savings_rates[4] = hr
		sys.it.non_savings_rates[5] = tr
		sys.it.personal_allowance = pa
		sys.uc.taper = uct
		sys.lmt.working_tax_credit.basic = wtcb

		if ucs == 0
			sys.uc.age_25_and_over = 0.0
			sys.uc.age_18_24 = 0.0
			sys.uc.couple_both_under_25 = 0.0
			sys.uc.couple_oldest_25_plus = 0.0
		else
			ucsd = ucs - sys.uc.age_25_and_over
			# move main uc allows equally, as in covid uplift
			sys.uc.age_25_and_over = max(0.0, ucsd+sys.uc.age_25_and_over)
			sys.uc.age_18_24 = max(0.0, ucsd+sys.uc.age_18_24)
			sys.uc.couple_both_under_25 = max(0.0, sys.uc.couple_both_under_25+ucsd)
			sys.uc.couple_oldest_25_plus = max(0.0,sys.uc.couple_oldest_25_plus+ucsd)
		end

		sys.nmt_bens.child_benefit.first_child = cb
		sys.nmt_bens.pensions.new_state_pension = pen
		sys.scottish_child_payment.amount = scp
		sys.ni.primary_class_1_rates[3] = ni_prim
		sys.ni.secondary_class_1_rates[2:3] .= ni_sec

		sys.ubi.abolished = true
		sys.ubi.adult_amount = 4_800.0
		sys.ubi.child_amount= 3_000.0
		sys.ubi.universal_pension = 8_780.0
		sys.ubi.adult_age = 17
		sys.ubi.retirement_age = 66

		#=
		sys.uc.abolished :
		sys.uc.other_income
		incomes     = LEGACY_MT_INCOME
        hb_incomes  = LEGACY_HB_INCOME
        pc_incomes  = LEGACY_PC_INCOME
        sc_incomes  = LEGACY_SAVINGS_CREDIT_INCOME
		=#

		results = do_run( sys )
	else
		results = (
			results = BASE_STATE.results,
			summary = BASE_STATE.summary,
			gain_lose = BASE_STATE.gain_lose )
	end

	println("sys.it.non_savings_rates $(sys.it.non_savings_rates)")
	println("BASE_STATE.sys.it.non_savings_rates $(BASE_STATE.sys.it.non_savings_rates)")
	return make_output_table(results,sys)
end

function no_nothings( things ...)::Bool
	for thing in things
		if isnothing( thing )
			return false
		end
	end
	return true
end

callback!(
    app,
    Output("model_running",  "children"),
	Output("output-block", "children"),
	Input("submit-button", "n_clicks"),

	State( "basic_rate", "value"),
	State( "higher_rate", "value"),
	State( "top_rate", "value"),
	State( "pa", "value"),

	State( "ni_prim", "value"),
	State( "ni_sec", "value"),

	State( "cb", "value"),
	State( "pen", "value"),
	State( "uctaper", "value"),
	State( "ucs", "value"),
	State( "wtcb", "value"),
	State( "scp", "value")

	) do n_clicks, basic_rate, higher_rate, top_rate, pa,
	     ni_prim, ni_sec,
		 cb, pen, uctaper, ucs, wtcb, scp

	println( "n_clicks = $n_clicks")
	# will return 'nothing' if something is out-of-range or not a number, or if no clicks on submit
	if no_nothings( n_clicks, basic_rate, higher_rate, top_rate, pa, ni_prim, ni_sec, cb, pen, uctaper, ucs, wtcb, scp )
		println( "running the live calc version")
		return [nothing, do_output( basic_rate, higher_rate, top_rate, pa, ni_prim, ni_sec, cb, pen, uctaper, ucs, wtcb, scp )]
	end
	println( "doing the do-nothing version")
	[nothing, do_output( 20, 41, 46, 12_570, 12, 13.8, 21.15, 179.60, 55, 324.84, 2_005, 10.0 )]
end

run_server( app, "0.0.0.0", 8053; debug=true )


text:
import Test: Test, record, finish
using Test: AbstractTestSet, DefaultTestSet, Result, Pass, Fail, Error
using Test: get_testset_depth, get_testset

struct InternalTestSet <: AbstractTestSet
    default_ts::DefaultTestSet
end

InternalTestSet(desc::AbstractString; verbose::Bool = false) =
    InternalTestSet(DefaultTestSet(desc, verbose = verbose))

record(ts::InternalTestSet, t::Result) = record(ts.default_ts, t)

function record(ts::InternalTestSet, t::Fail)
    push!(ts.default_ts.results, t)
    t
end

finish(ts::InternalTestSet) = finish(ts.default_ts)


text:
int_rules_4_3_10 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.3.10*(c+d*x)^m*(a+b*tan)^n =#
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * tan(~(e') + ~(k') * pi + ~(f') * Complex(0, ~fz) * ~x), ~x) => (-im * (~c + ~d * ~x) ^ (~m + 1)) / (~d * (~m + 1)) + 2 * im * Antiderivative(((~c + ~d * ~x) ^ ~m * ℯ ^ (-2 * im * ~k * pi) * ℯ ^ (2 * (-im * ~e + ~f * ~fz * ~x))) / (1 + ℯ ^ (-2 * im * ~k * pi) * ℯ ^ (2 * (-im * ~e + ~f * ~fz * ~x))), ~x) <-- FreeQ([~c, ~d, ~e, ~f, ~fz], ~x) && (IntegerQ(4 * ~k) && IGtQ(~m, 0))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * tan(~(e') + ~(k') * pi + ~(f') * ~x), ~x) => (im * (~c + ~d * ~x) ^ (~m + 1)) / (~d * (~m + 1)) - 2 * im * Antiderivative(((~c + ~d * ~x) ^ ~m * ℯ ^ (2 * im * ~k * pi) * ℯ ^ (2 * im * (~e + ~f * ~x))) / (1 + ℯ ^ (2 * im * ~k * pi) * ℯ ^ (2 * im * (~e + ~f * ~x))), ~x) <-- FreeQ([~c, ~d, ~e, ~f], ~x) && (IntegerQ(4 * ~k) && IGtQ(~m, 0))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * tan(~(e') + ~(f') * Complex(0, ~fz) * ~x), ~x) => (-im * (~c + ~d * ~x) ^ (~m + 1)) / (~d * (~m + 1)) + 2 * im * Antiderivative(((~c + ~d * ~x) ^ ~m * ℯ ^ (2 * (-im * ~e + ~f * ~fz * ~x))) / (1 + ℯ ^ (2 * (-im * ~e + ~f * ~fz * ~x))), ~x) <-- FreeQ([~c, ~d, ~e, ~f, ~fz], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * tan(~(e') + ~(f') * ~x), ~x) => (im * (~c + ~d * ~x) ^ (~m + 1)) / (~d * (~m + 1)) - 2 * im * Antiderivative(((~c + ~d * ~x) ^ ~m * ℯ ^ (2 * im * (~e + ~f * ~x))) / (1 + ℯ ^ (2 * im * (~e + ~f * ~x))), ~x) <-- FreeQ([~c, ~d, ~e, ~f], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((~b * (~c + ~d * ~x) ^ ~m * (~b * tan(~e + ~f * ~x)) ^ (~n - 1)) / (~f * (~n - 1)) - ((~b * ~d * ~m) / (~f * (~n - 1))) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * (~b * tan(~e + ~f * ~x)) ^ (~n - 1), ~x)) - (~b) ^ 2 * Antiderivative((~c + ~d * ~x) ^ ~m * (~b * tan(~e + ~f * ~x)) ^ (~n - 2), ~x) <-- FreeQ([~b, ~c, ~d, ~e, ~f], ~x) && (GtQ(~n, 1) && GtQ(~m, 0))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (((~c + ~d * ~x) ^ ~m * (~b * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~b * ~f * (~n + 1)) - ((~d * ~m) / (~b * ~f * (~n + 1))) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * (~b * tan(~e + ~f * ~x)) ^ (~n + 1), ~x)) - (1 / (~b) ^ 2) * Antiderivative((~c + ~d * ~x) ^ ~m * (~b * tan(~e + ~f * ~x)) ^ (~n + 2), ~x) <-- FreeQ([~b, ~c, ~d, ~e, ~f], ~x) && (LtQ(~n, -1) && GtQ(~m, 0))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * (~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand((~c + ~d * ~x) ^ ~m, (~a + ~b * tan(~e + ~f * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && (IGtQ(~m, 0) && IGtQ(~n, 0))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') / (~a + ~(b') * tan(~(e') + ~(f') * ~x)), ~x) => ((~c + ~d * ~x) ^ (~m + 1) / (2 * ~a * ~d * (~m + 1)) - (~a * (~c + ~d * ~x) ^ ~m) / (2 * ~b * ~f * (~a + ~b * tan(~e + ~f * ~x)))) + ((~a * ~d * ~m) / (2 * ~b * ~f)) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) / (~a + ~b * tan(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && GtQ(~m, 0))
    @apply_utils Antiderivative(1 / ((~(c') + ~(d') * ~x) ^ 2 * (~a + ~(b') * tan(~(e') + ~(f') * ~x))), ~x) => (-1 / (~d * (~c + ~d * ~x) * (~a + ~b * tan(~e + ~f * ~x))) + (~f / (~b * ~d)) * Antiderivative(cos(2 * ~e + 2 * ~f * ~x) / (~c + ~d * ~x), ~x)) - (~f / (~a * ~d)) * Antiderivative(sin(2 * ~e + 2 * ~f * ~x) / (~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && EqQ((~a) ^ 2 + (~b) ^ 2, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~m / (~a + ~(b') * tan(~(e') + ~(f') * ~x)), ~x) => (~f * (~c + ~d * ~x) ^ (~m + 2)) / (~b * (~d) ^ 2 * (~m + 1) * (~m + 2)) + (~c + ~d * ~x) ^ (~m + 1) / (~d * (~m + 1) * (~a + ~b * tan(~e + ~f * ~x))) + ((2 * ~b * ~f) / (~a * ~d * (~m + 1))) * Antiderivative((~c + ~d * ~x) ^ (~m + 1) / (~a + ~b * tan(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && (LtQ(~m, -1) && NeQ(~m, -2)))
    #= Antiderivative(((~c')+(~d')*(~x))^(~m)/((~a)+(~b')*tan((~e')+(~f')*(~x))),~x) :=  (c+d*x)^(m+1)/(d*(m+1)*(a+b*tan(e+f*x))) + f/(b*d*(m+1))*Antiderivative((c+d*x)^(m+1),x) + 2*b*f/(a*d*(m+1))*Antiderivative((c+d*x)^(m+1)/(a+b*tan(e+f*x)),x) <--  FreeQ([a,b,c,d,e,f],x) && EqQ(a^2+b^2,0) && LtQ(m,-1) =#
    @apply_utils Antiderivative(1 / ((~(c') + ~(d') * ~x) * (~a + ~(b') * tan(~(e') + ~(f') * ~x))), ~x) => log(~c + ~d * ~x) / (2 * ~a * ~d) + (1 / (2 * ~a)) * Antiderivative(cos(2 * ~e + 2 * ~f * ~x) / (~c + ~d * ~x), ~x) + (1 / (2 * ~b)) * Antiderivative(sin(2 * ~e + 2 * ~f * ~x) / (~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && EqQ((~a) ^ 2 + (~b) ^ 2, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~m / (~a + ~(b') * tan(~(e') + ~(f') * ~x)), ~x) => (~c + ~d * ~x) ^ (~m + 1) / (2 * ~a * ~d * (~m + 1)) + (1 / (2 * ~a)) * Antiderivative((~c + ~d * ~x) ^ ~m * ℯ ^ (((2 * ~a) / ~b) * (~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && Not(IntegerQ(~m)))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~m * (~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => Antiderivative(ExpandIntegrand((~c + ~d * ~x) ^ ~m, (1 / (2 * ~a) + cos(2 * ~e + 2 * ~f * ~x) / (2 * ~a) + sin(2 * ~e + 2 * ~f * ~x) / (2 * ~b)) ^ -(~n), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && (ILtQ(~m, 0) && ILtQ(~n, 0)))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~m * (~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => Antiderivative(ExpandIntegrand((~c + ~d * ~x) ^ ~m, (1 / (2 * ~a) + ℯ ^ (((2 * ~a) / ~b) * (~e + ~f * ~x)) / (2 * ~a)) ^ -(~n), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && ILtQ(~n, 0))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * (~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => With([u = IntHide((~a + ~b * tan(~e + ~f * ~x)) ^ ~n, ~x)], Dist((~c + ~d * ~x) ^ ~m, u, ~x) - ~d * ~m * Antiderivative(Dist((~c + ~d * ~x) ^ (~m - 1), u, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && (ILtQ(~n, -1) && GtQ(~m, 0)))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') / (~a + ~(b') * tan(~(e') + ~(k') * pi + ~(f') * ~x)), ~x) => (~c + ~d * ~x) ^ (~m + 1) / (~d * (~m + 1) * (~a + im * ~b)) + 2 * im * ~b * Antiderivative(((~c + ~d * ~x) ^ ~m * ℯ ^ (2 * im * ~k * pi) * ℯ ^ Simp(2 * im * (~e + ~f * ~x), ~x)) / ((~a + im * ~b) ^ 2 + ((~a) ^ 2 + (~b) ^ 2) * ℯ ^ (2 * im * ~k * pi) * ℯ ^ Simp(2 * im * (~e + ~f * ~x), ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IntegerQ(4 * ~k) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && IGtQ(~m, 0)))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') / (~a + ~(b') * tan(~(e') + ~(f') * ~x)), ~x) => (~c + ~d * ~x) ^ (~m + 1) / (~d * (~m + 1) * (~a + im * ~b)) + 2 * im * ~b * Antiderivative(((~c + ~d * ~x) ^ ~m * ℯ ^ Simp(2 * im * (~e + ~f * ~x), ~x)) / ((~a + im * ~b) ^ 2 + ((~a) ^ 2 + (~b) ^ 2) * ℯ ^ Simp(2 * im * (~e + ~f * ~x), ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && IGtQ(~m, 0))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) / (~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ 2, ~x) => (-((~c + ~d * ~x) ^ 2) / (2 * ~d * ((~a) ^ 2 + (~b) ^ 2)) - (~b * (~c + ~d * ~x)) / (~f * ((~a) ^ 2 + (~b) ^ 2) * (~a + ~b * tan(~e + ~f * ~x)))) + (1 / (~f * ((~a) ^ 2 + (~b) ^ 2))) * Antiderivative((~b * ~d + 2 * ~a * ~c * ~f + 2 * ~a * ~d * ~f * ~x) / (~a + ~b * tan(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && NeQ((~a) ^ 2 + (~b) ^ 2, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * (~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => Antiderivative(ExpandIntegrand((~c + ~d * ~x) ^ ~m, (1 / (~a - im * ~b) - (2 * im * ~b) / ((~a) ^ 2 + (~b) ^ 2 + (~a - im * ~b) ^ 2 * ℯ ^ (2 * im * (~e + ~f * ~x)))) ^ -(~n), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (ILtQ(~n, 0) && IGtQ(~m, 0)))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) * sqrt(~a + ~(b') * tan(~(e') + ~(f') * ~x)), ~x) => (-(sqrt(2)) * ~b * (~c + ~d * ~x) * atanh(sqrt(~a + ~b * tan(~e + ~f * ~x)) / (sqrt(2) * Rt(~a, 2)))) / (Rt(~a, 2) * ~f) + ((sqrt(2) * ~b * ~d) / (Rt(~a, 2) * ~f)) * Antiderivative(atanh(sqrt(~a + ~b * tan(~e + ~f * ~x)) / (sqrt(2) * Rt(~a, 2))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && EqQ((~a) ^ 2 + (~b) ^ 2, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) * sqrt(~(a') + ~(b') * tan(~(e') + ~(f') * ~x)), ~x) => (((-im * Rt(~a - im * ~b, 2) * (~c + ~d * ~x)) / ~f) * atanh(sqrt(~a + ~b * tan(~e + ~f * ~x)) / Rt(~a - im * ~b, 2)) + ((im * Rt(~a + im * ~b, 2) * (~c + ~d * ~x)) / ~f) * atanh(sqrt(~a + ~b * tan(~e + ~f * ~x)) / Rt(~a + im * ~b, 2)) + ((im * ~d * Rt(~a - im * ~b, 2)) / ~f) * Antiderivative(atanh(sqrt(~a + ~b * tan(~e + ~f * ~x)) / Rt(~a - im * ~b, 2)), ~x)) - ((im * ~d * Rt(~a + im * ~b, 2)) / ~f) * Antiderivative(atanh(sqrt(~a + ~b * tan(~e + ~f * ~x)) / Rt(~a + im * ~b, 2)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && NeQ((~a) ^ 2 + (~b) ^ 2, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) / sqrt(~a + ~(b') * tan(~(e') + ~(f') * ~x)), ~x) => (1 / (2 * ~a)) * Antiderivative((~c + ~d * ~x) * sqrt(~a + ~b * tan(~e + ~f * ~x)), ~x) + (~a / 2) * Antiderivative(((~c + ~d * ~x) * sec(~e + ~f * ~x) ^ 2) / (~a + ~b * tan(~e + ~f * ~x)) ^ (3 / 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && EqQ((~a) ^ 2 + (~b) ^ 2, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) / sqrt(~(a') + ~(b') * tan(~(e') + ~(f') * ~x)), ~x) => (((-im * (~c + ~d * ~x)) / (~f * Rt(~a - im * ~b, 2))) * atanh(sqrt(~a + ~b * tan(~e + ~f * ~x)) / Rt(~a - im * ~b, 2)) + ((im * (~c + ~d * ~x)) / (~f * Rt(~a + im * ~b, 2))) * atanh(sqrt(~a + ~b * tan(~e + ~f * ~x)) / Rt(~a + im * ~b, 2)) + ((im * ~d) / (~f * Rt(~a - im * ~b, 2))) * Antiderivative(atanh(sqrt(~a + ~b * tan(~e + ~f * ~x)) / Rt(~a - im * ~b, 2)), ~x)) - ((im * ~d) / (~f * Rt(~a + im * ~b, 2))) * Antiderivative(atanh(sqrt(~a + ~b * tan(~e + ~f * ~x)) / Rt(~a + im * ~b, 2)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && NeQ((~a) ^ 2 + (~b) ^ 2, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * tan(~(e') + ~(f') * ~x) ^ ~(n'), ~x) => If(MatchQ(~f, ~(f1') * Complex(0, ~j)), If(MatchQ(~e, ~(e1') + pi / 2), im ^ ~n * Unintegrable((~c + ~d * ~x) ^ ~m * coth(-im * (~e - pi / 2) - im * ~f * ~x) ^ ~n, ~x), im ^ ~n * Unintegrable((~c + ~d * ~x) ^ ~m * tanh(-im * ~e - im * ~f * ~x) ^ ~n, ~x)), If(MatchQ(~e, ~(e1') + pi / 2), (-1) ^ ~n * Unintegrable((~c + ~d * ~x) ^ ~m * cot((~e - pi / 2) + ~f * ~x) ^ ~n, ~x), Unintegrable((~c + ~d * ~x) ^ ~m * tan(~e + ~f * ~x) ^ ~n, ~x))) <-- FreeQ([~c, ~d, ~e, ~f, ~m, ~n], ~x) && IntegerQ(~n)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * (~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => Unintegrable((~c + ~d * ~x) ^ ~m * (~a + ~b * tan(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x)
    @apply_utils Antiderivative((~u) ^ ~(m') * (~(a') + ~(b') * tan(~v)) ^ ~(n'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * (~a + ~b * tan(ExpandToSum(~v, ~x))) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~m, ~n], ~x) && (LinearQ([~u, ~v], ~x) && Not(LinearMatchQ([~u, ~v], ~x)))
    @apply_utils Antiderivative((~u) ^ ~(m') * (~(a') + ~(b') * cot(~v)) ^ ~(n'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * (~a + ~b * cot(ExpandToSum(~v, ~x))) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~m, ~n], ~x) && (LinearQ([~u, ~v], ~x) && Not(LinearMatchQ([~u, ~v], ~x)))
end

max_stars_repo_path: src/rules/4 Trig functions/4.3 Tangent/4.3.10 (c+d x)^m (a+b tan)^n.jl


text:
"""
    Line(f::File,no::Int)

Create a reference to line `no` in file `f`.
"""
mutable struct Line <: Reference
    file::Union{File,Nothing}
    ln::Union{Int,Nothing}
    function Line(f::File,no::Int)
        @assert length(f.data) >= no
        @assert no > 0
        l = new(f,no)
        push!(f.references,l)
        return l
    end
end

==(a::Line,b::Line) = a.ln == b.ln && a.file === b.file

@alive l value(l::Line) = l.file.data[l.ln]

function show(io::IO,l::Line)
    if l.file == nothing
        print(io,"*Deleted*")
    else
        print(io,"Line($(l.ln)) #$(value(l)|> truncate)")
    end
end

"""
    next(l::Line,i::Int=1)

Returns the `i`-th line after `l`.
"""
@alive l next(l::Line,i::Int=1) = Line(l.file,l.ln+i)
"""
    prev(l::Line,i::Int=1)

Returns the `i`-th line before `l`.
"""
@alive l prev(l::Line,i::Int=1) = Line(l.file,l.ln-i)

+(l::Line,i::Int) = next(l,i)
-(l::Line,i::Int) = prev(l,i)
+(i::Int,l::Line) = next(l,i)
-(i::Int,l::Line) = prev(l,i)

"""
    copy(l::Line)

Create a copy of line `l`.
"""
@alive l copy(l::Line) = Line(l.file,l.ln)

#=
Delete functions
=#
"""
    delete!(l::Line)

Delete line `l` from the document it's in.
"""
@alive l function delete!(l::Line)
    deletefromreferences!(l.file,l)
    deleteat!(l.file.data,l.ln)
    notify_delete(l.file,l)
    destroy!(l)
    return nothing
end

@alive l function destroy!(l::Line)
    l.file = nothing
    l.ln = nothing
end
isdestroyed(l::Line) = l.ln == nothing && l.file == nothing

function notify_delete(self::Line,del::Line)
    if self.ln > del.ln
        self.ln -= 1
    elseif self.ln == del.ln
        deletefromreferences!(self.file,self)
        destroy!(self)
    end
end

#=
Insert functions
=#
"""
    insert!(l::Line,str::String)

Insert `str` before line `l`.
"""
@alive l function insert!(l::Line,str::String)
    notify_insert(l.file,l,1)
    insert!(l.file.data,l.ln,str)
    l.ln += 1
    return Line(l.file,l.ln-1)
end

"""
    insert!(l::Line,from::Reference)

Insert the value of reference `from` before line `l`.
"""
@alive l from insert!(l::Line,from::T) where T<:Reference = insert!(l,value(from))

function notify_insert(self::Line,ins::Line,lines::Int)
    if self.ln >= ins.ln
        self.ln += lines
    end
end

#=
Append functions
=#

"""
    append!(l::Line,str::String)

Append `str` after line `l`.
"""
@alive l function append!(l::Line,str::String)
    notify_append(l.file,l,1)
    if l.ln < length(l.file.data)
        insert!(l.file.data,l.ln+1,str)
    else
        push!(l.file.data,str)
    end
    return Line(l.file,l.ln+1)
end

"""
    append!(l::Line,from::Reference)

Append the value of reference `from` after line `l`.
"""
@alive l from append!(l::Line,from::T) where T<:Reference = append!(l,value(from))

function notify_append(self::Line,app::Line,lines::Int)
    if self.ln > app.ln
        self.ln += lines
    end
end

#=
Replace functions
=#

"""
    replace!(l::Line,str::String)

Replace the value of line `l` with `str`
"""
@alive l function replace!(l::Line,str::String)
    l.file.data[l.ln] = str
    return l
end

"""
    replace!(l::Line,from::Reference)

Replace line `l` with the value of `from`
"""
@alive l from replace!(l::Line,from::T) where T<:Reference = replace!(l,value(from))

function notify_replace(self::Line,rep::Line,lines::Int)
    if self.ln > rep.ln
        self.ln -= lines
    end
end

###############################################################################
#                                     Move                                    #
###############################################################################

movebeforeline(ref::Line) = ref.ln
moveafterline(ref::Line) = ref.ln

"""
    move!(self::Line,to::T)

Move `self` to `to`. This is a combination of `insert!(to,self)` and `delete!(self)`.
"""
@alive self to function move!(self::Line,to::T) where T<:Reference
    # Insert duplicate
    dup = insert!(to,self)
    # Delete original but keep reference
    notify_delete(self.file,self)
    deleteat!(self.file.data,self.ln)
    # Move pointer
    self.ln = dup.ln
    # Delete reference for dup
    deletefromreferences!(dup.file,dup)
    return self
end

"""
    moveafter!(self::Line,to::T)

Move `self` to the position after `to`. This is a combination of `append!(to,self)` and `delete!(self)`.
"""
@alive self to function moveafter!(self::Line,to::T) where T<:Reference
    # Insert duplicate
    dup = append!(to,self)
    # Delete original but keep reference
    notify_delete(self.file,self)
    deleteat!(self.file.data,self.ln)
    # Move pointer
    self.ln = dup.ln
    # Delete reference for dup
    deletefromreferences!(dup.file,dup)
    return self
end


text:
"""
    rglob(filepattern, dir = pwd(), depth = 4)

Recursive glob up to `depth` layers deep within `dir`.
"""
function rglob(filepattern = "*", dir = pwd(), depth = 4)
    patterns = [repeat("*/", i) * filepattern for i in 0:(depth - 1)]

    return vcat([glob(pattern, dir) for pattern in patterns]...)
end

"""
    FileDataset([loadfn = FileIO.load,] paths)
    FileDataset([loadfn = FileIO.load,] dir, pattern = "*", depth = 4)

Wrap a set of file `paths` as a dataset (traversed in the same order as `paths`).
Alternatively, specify a `dir` and collect all paths that match a glob `pattern`
(recursively globbing by `depth`). The glob order determines the traversal order.
"""
struct FileDataset{F, T<:AbstractString} <: AbstractDataContainer
    loadfn::F
    paths::Vector{T}
end

FileDataset(paths) = FileDataset(FileIO.load, paths)
FileDataset(loadfn,
            dir::AbstractString,
            pattern::AbstractString = "*",
            depth = 4) = FileDataset(loadfn, rglob(pattern, string(dir), depth))
FileDataset(dir::AbstractString, pattern::AbstractString = "*", depth = 4) =
    FileDataset(FileIO.load, dir, pattern, depth)

Base.getindex(dataset::FileDataset, i::Integer) = dataset.loadfn(dataset.paths[i])
Base.getindex(dataset::FileDataset, is::AbstractVector) = map(Base.Fix1(getobs, dataset), is)
Base.length(dataset::FileDataset) = length(dataset.paths)


text:
using LightXML, TimeZones, DataFrames

export parse_run, isrun
export parse_laps

mutable struct Trackpoint
    dt::ZonedDateTime
    distance::Float64
    altitude::Float64
    hr::Float64
    speed::Float64
    pos::Tuple{Float64, Float64}
end

function parse_trackpoint(tp::XMLElement)
    @assert name(tp) == "Trackpoint"
    dt = ZonedDateTime(content(find_element(tp, "Time")))
    distance = parse(Float64, content(find_element(tp, "DistanceMeters")))
    altitude = parse(Float64, content(find_element(tp, "AltitudeMeters")))
    hr_el = find_element(tp, "HeartRateBpm")
    hr = NaN
    if !isnothing(hr_el)
        hr = parse(Float64, content(hr_el))
    end
    ext_el = find_element(tp, "Extensions")
    n3_el = find_element(ext_el, "TPX")
    speed_el = find_element(n3_el, "Speed")
    speed = parse(Float64, content(speed_el))
    pos_el = find_element(tp, "Position")
    lat = parse(Float64, content(find_element(pos_el, "LatitudeDegrees")))
    long = parse(Float64, content(find_element(pos_el, "LongitudeDegrees")))

    return Trackpoint(dt, distance, altitude, hr, speed, (lat, long))
end

function count_points(activity::XMLElement)
    @assert name(activity) == "Activity"

    n = 0

    laps = activity["Lap"]

    for lap in laps

        tracks = lap["Track"]

        for track in tracks
            n += length(track["Trackpoint"])
        end
    end

    return n
end

function get_activity(path::String)
    xdoc = parse_file(path)

    activities_el = root(xdoc)["Activities"][1]

    activity_els = activities_el["Activity"]

    @assert length(activity_els) == 1 #TODO: handle activities

    return activity_els[1]
end

function parse_laps(path::String)
    activity = get_activity(path)

    laps = activity["Lap"]
    n = length(laps)

    df = DataFrame()
    df.seconds = Vector{Float64}(undef, n)
    df.meters = Vector{Float64}(undef, n)
    for i = 1:length(laps)

        lap = laps[i]
        dm = find_element(lap, "DistanceMeters")
        df.meters[i] = parse(Float64, content(dm))

        tsec = find_element(lap, "TotalTimeSeconds")
        df.seconds[i] = parse(Float64, content(tsec))
    end

    return df
end

function parse_run(path::String)

    activity = get_activity(path)

    @assert attribute(activity, "Sport") == "Running" # TODO: handle other sports

    n = count_points(activity)

    df = DataFrame()
    df.time = Vector{ZonedDateTime}(undef, n)
    df.distance = zeros(n)
    df.hr = zeros(n)
    df.speed = zeros(n)
    df.altitude = zeros(n)
    df.pos = Vector{Tuple{Float64, Float64}}(undef, n)

    counter = 0
    for lap in activity["Lap"]
        for track in lap["Track"]
            for tp_el in track["Trackpoint"]
                counter += 1
                tp = parse_trackpoint(tp_el)
                if isnan(tp.hr)
                    tp.hr = df.hr[counter - 1] # TODO: handle missing data better
                end

                df.time[counter] = tp.dt
                df.distance[counter] = tp.distance
                df.hr[counter] = tp.hr
                df.speed[counter] = tp.speed
                df.altitude[counter] = tp.altitude
                df.pos[counter] = tp.pos

            end
        end
    end

    # Activity ID is same as start time
    start_time = ZonedDateTime(content(find_element(activity, "Id")))

    rs = RunSummary(n, start_time)
    init_speed = df.speed[1]

    #sometimes initial speed is zero
    if init_speed  == 0.0
        ind = 2
        while init_speed == 0.0
            init_speed = df.speed[ind]
            ind += 1
        end

    end

    rs.time[1] = round(df.distance[1] / init_speed) #should I just assume first tp is start???

    rs.dist[1] = df.distance[1]
    rs.alt[1] = 0.0
    rs.hr[1] = df.hr[1]

    for i = 2:n

        rs.time[i] = (df.time[i] - df.time[i - 1]).value / 1e3 # time stored in milliseconds
        rs.alt[i] = df.altitude[i] - df.altitude[i - 1]
        rs.dist[i] = geo_dist(df.pos[i], df.pos[i - 1])
        rs.hr[i] = df.hr[i]
    end

    return rs
end

const E_RADIUS = 6.3781370e6
const P_RADIUS = 6.3567523e6
const DEG_TO_RAD = 2.0 * π / 360.0

function geo_dist(p1::Tuple{Float64, Float64}, p2::Tuple{Float64, Float64})
    dlat = (p1[1] - p2[1]) * DEG_TO_RAD
    dlong = (p1[2]- p2[2]) * DEG_TO_RAD
    mlat = 0.5 * (p1[1] + p2[1]) * DEG_TO_RAD

    cmlat = cos(mlat)

    x = 1.0 / sqrt(1.0 / E_RADIUS^2 + tan(mlat)^2 / P_RADIUS^2) #distance from center of ellipse
    radius = x / cmlat

    return radius * sqrt(dlat^2 + (cmlat * dlong)^2)
end

function isrun(path::String)
    activity = get_activity(path)

    return attribute(activity, "Sport") == "Running"
end


text:
abstract type AbstractServiceFormulation end

mutable struct ServiceModel{D <: PSY.Service, B <: AbstractServiceFormulation}
    service_type::Type{D}
    formulation::Type{B}
end


text:
# # Cohen-Bergstresser model
#
# This example considers the Cohen-Bergstresser model[^CB1966],
# reproducing the results of the original paper. This model is particularly
# simple since its linear nature allows one to get away without any
# self-consistent field calculation.
#
# [^CB1966]: M. L. Cohen and T. K. Bergstresser Phys. Rev. **141**, 789 (1966) DOI [10.1103/PhysRev.141.789](https://doi.org/10.1103/PhysRev.141.789)

# We build the lattice using the tabulated lattice constant from the original paper,
# stored in DFTK:
using DFTK

Si = ElementCohenBergstresser(:Si)
lattice = Si.lattice_constant / 2 .* [[0 1 1.]; [1 0 1.]; [1 1 0.]]
atoms = [Si => [ones(3)/8, -ones(3)/8]];

# Next we build the rather simple model and discretise it with moderate `Ecut`:
Ecut = 10.0
model = Model(lattice; atoms=atoms, terms=[Kinetic(), AtomicLocal()])
basis = PlaneWaveBasis(model, Ecut, kgrid=(1, 1, 1));

# We diagonalise at the Gamma point to find a Fermi level ...
ham = Hamiltonian(basis)
eigres = diagonalize_all_kblocks(DFTK.lobpcg_hyper, ham, 6)
εF = find_fermi_level(basis, eigres.λ)

# ... and compute and plot 8 bands:
using Plots

n_bands = 8
ρ0 = guess_density(basis)  # Just dummy, has no meaning in this model
p = plot_bandstructure(basis, ρ0, n_bands, εF=εF, kline_density=10)
ylims!(p, (-5, 6))


text:
include("distribution_types.jl")
include("distributions.jl")
include("samples.jl")


text:
mutable struct Envelope
    fs::Float64
    attack::Float64 #attack in seconds
    decay::Float64  #decay in seconds
    sustain::Float64 #sustain
    sustainTime::Float64 # sustain duration in seconds
    release::Float64 #release in seconds
end
function getAttack(env)
    size = env.fs * env.attack
    return 0.0 : inv(size) : 1
end
function getDecay(env)
    size = env.fs * env.decay
    return 1:inv(size)*-(1-env.sustain):env.sustain
end
function getRelease(env)
    size =env.fs * env.release
    return env.sustain : inv(size)*-env.sustain : 0
end
function getSustain(env)
    size = round(env.fs * env.sustainTime)
    return zeros(Int(floor(size))).+env.sustain
end

function applyEnvelope(samples::AbstractVector{T}, env)where
    {T<:Real}
    returnVector = T[] #signal to return
    size = length(samples) #length of the original signal
    #get vectors that represent attack, decay of a given envelope
    attack = getAttack(env)
    decay = getDecay(env)
    release = getRelease(env)
    decayAttack = length(attack)+length(decay)
    sustain = getSustain(env)
    sustainLength = length(sustain)
    #appending attack part of the envelope
    append!(returnVector, attack.*samples[1:length(attack)] )
    #appending decay
    append!(returnVector, decay.*samples[1+length(attack):decayAttack] )
    #appending sustain
    append!(returnVector, sustain .* samples[decayAttack+1:sustainLength+decayAttack])
    #appending release
    append!(returnVector, release.*samples[sustainLength+decayAttack+1:sustainLength+decayAttack+length(release)])
    return returnVector
end

max_stars_repo_name: justinbroce/julia-audio-stuffs

max_issues_repo_name: justinbroce/julia-audio-stuffs

max_forks_repo_name: justinbroce/julia-audio-stuffs


text:
using CodecZlib, Requires

abstract type ArchiveFormat{S} end

decompressor_stream(::Type{ArchiveFormat{:Gzip}}, io::IOStream) = GzipDecompressorStream(io)

# Handles the case when we are loading a supported archive format but the codec is not loaded
function decompressor_stream(::Type{ArchiveFormat{S}}, io::IOStream) where S
    throw(CodecMissingError(S))
end

# called by __init__() to enable support for optional codecs
function define_optional_decompressor_streams()
    @require CodecBzip2="523fee87-0ab8-5b00-afb7-3ecf72e48cfd" begin
        decompressor_stream(::Type{ArchiveFormat{:Bzip2}}, io::IOStream) = CodecBzip2.Bzip2DecompressorStream(io)
    end
    @require CodecXz="ba30903b-d9e8-5048-a5ec-d1f5b0d4b47b" begin
        decompressor_stream(::Type{ArchiveFormat{:Xz}}, io::IOStream) = CodecXz.XzDecompressorStream(io)
    end
end

# Magic numbers that identify the supported compression formats
const MAGIC_GZIP = b"\x1F\x8B"
const MAGIC_BZIP2 = b"BZh"
const MAGIC_XZ = b"\xFD7zXZ\x00"

function decompress(io)
    # Read the first 6 bytes to obtain the magic number if there is one.
    buffer = zeros(UInt8, 6)
    readbytes!(io, buffer)
    seekstart(io)

    # Check for any of the gzip, bzip2 or xz magic numbers
    if buffer[1:2] == MAGIC_GZIP
        return decompressor_stream(ArchiveFormat{:Gzip}, io)
    elseif buffer[1:3] == MAGIC_BZIP2
        return decompressor_stream(ArchiveFormat{:Bzip2}, io)
    elseif buffer == MAGIC_XZ
        return decompressor_stream(ArchiveFormat{:Xz}, io)
    end

    # If none of the magic numbers match, we assume the file is not compressed

    return io
end


text:
module BinaryLambdaCalculus

include("coding.jl")
include("grygiel_lescanne.jl")

end


text:
abstract type TranslationModel <:EoSModel end
function init_model(model::TranslationModel,components,userlocations,verbose)
    return model
end

function init_model(model::Type{<:TranslationModel},components,userlocations,verbose)
    verbose && @info("""Now creating mixing model:
    $model""")
    return model(components;userlocations,verbose)
end

has_sites(::Type{<:TranslationModel})=false

include("NoTranslation.jl")
include("Rackett.jl")
include("Peneloux.jl")
include("MT.jl")


text:
## header to provide surface and context
using Cairo
c = CairoRGBSurface(256,256);
cr = CairoContext(c);

save(cr);
set_source_rgb(cr,0.8,0.8,0.8);    # light gray
rectangle(cr,0.0,0.0,256.0,256.0); # background
fill(cr);
restore(cr);

save(cr);

## api example, use set_line_type

dashes_by_name = [
    "solid",
    "dot",
    "dotdashed",
    "longdashed",
    "dash",
    "dotdotdashed",
    "dotdotdotdashed"]

set_line_width(cr, 3.0);

for k=1:7
    save(cr)
    set_line_type(cr, dashes_by_name[k])
    move_to(cr,16.0,16.0+(k*32.0))
    line_to(cr,240.0,16.0+(k*32.0))
    stroke(cr)
    restore(cr)
end

## mark picture with current date
restore(cr);
move_to(cr,0.0,12.0);
set_source_rgb(cr, 0,0,0);
show_text(cr,Libc.strftime(time()));
write_to_png(c,"sample_set_dash.png");


text:
# The following block allows reading Union types written prior to v0.2
const LEGACY_H5TYPE_UNION = VariableLengthDatatype(H5TYPE_DATATYPE)

function jlconvert(::ReadRepresentation{Union, Vlen{DataTypeODR()}}, f::JLDFile,
                   ptr::Ptr, header_offset::RelOffset)
    v = Union{jlconvert(ReadRepresentation{DataType,Vlen{DataTypeODR()}}(), f, ptr, NULL_REFERENCE)...}
    track_weakref!(f, header_offset, v)
    v
end

constructrr(::JLDFile, ::Type{T}, dt::VariableLengthDatatype, ::Vector{ReadAttribute}) where {T<:Union} =
    dt == LEGACY_H5TYPE_UNION ? (ReadRepresentation{Union,Vlen{DataTypeODR()}}(), true) :
                         throw(UnsupportedFeatureException())

# The following definition is needed to correctly load Strings written
# with JLD2 with versions v0.1.12 - v0.3.1
function read_array(f::JLDFile, dataspace::ReadDataspace,
                    rr::FixedLengthString{String}, data_length::Int,
                    filter_id::UInt16, header_offset::RelOffset,
                    attributes::Union{Vector{ReadAttribute},Nothing})
    rrv = ReadRepresentation{UInt8,odr(UInt8)}()
    v = read_array(f, dataspace, rrv, data_length, filter_id, NULL_REFERENCE, attributes)
    String(v)
end


text:
using BlackBoxOptim
using Distributions

abstract NaturalEvolutionStrategyOpt <: PopulationOptimizer

type SeparableNESOpt <: NaturalEvolutionStrategyOpt
  d::Integer                      # Number of dimensions
  lambda::Integer                 # Number of samples to take per iteration
  mu::Array{Float64,2}            # Average position for sampling each position
  sigma::Array{Float64,2}         # Average std deviation for sampling in position
  distr::Distribution             # Distribution we sample the step sizes from
  mu_learnrate::Float64
  sigma_learnrate::Float64
  last_s::Array{Float64,2}        # The s values sampled in the last call to ask
  population::Array{Float64,2}    # The last sampled values, now being evaluated
  utilities::Array{Float64,2}     # The fitness shaping utility vector

  SeparableNESOpt(searchSpace; lambda = false, mu_learnrate = 1.0,
    sigma_learnrate = false) = begin

    numDimensions = length(mins(searchSpace))

    mu = rand(numDimensions, 1)
    sigma = ones(numDimensions, 1)
    distr = Normal(0, 1)

    lambda = lambda || convert(Int64, 4 + ceil(log(3*numDimensions)))
    sigma_learnrate = sigma_learnrate || calc_sigma_learnrate_for_snes(numDimensions)

    new(numDimensions, lambda, mu, sigma, distr,
      mu_learnrate, sigma_learnrate,
      eye(numDimensions), eye(numDimensions),
      # Most modern NES papers use log rather than linear fitness shaping.
      fitness_shaping_utilities_log(lambda))
  end
end

# We use a different ordering of the dimensions than other optimizers, so transpose.
population(o::NaturalEvolutionStrategyOpt) = o.population'

NES_DefaultOptions = {
  "lambda" => false,          # If false it will be set based on the number of dimensions
  "mu_learnrate" => 1.0,
  "sigma_learnrate" => false, # If false it will be set based on the number of dimensions
}

function separable_nes(parameters)
  params = merge(NES_DefaultOptions, parameters)
  SeparableNESOpt(params[:SearchSpace];
    lambda = params["lambda"],
    mu_learnrate = params["mu_learnrate"],
    sigma_learnrate = params["sigma_learnrate"])
end

calc_sigma_learnrate_for_snes(d) = (3 + log(d)) / (5 * sqrt(d))
calc_sigma_learnrate_for_nes(d) = (9 + 3 * log(d)) / (5 * d * sqrt(d))

# Get a set of new individuals to be ranked based on fitness.
function ask(snes::SeparableNESOpt)
  # Sample from N(0, 1)
  snes.last_s = s = randn(snes.d, snes.lambda)

  # Add in the mu and sigma's...
  #sampled_solutions = repmat(snes.mu, 1, l) + repmat(snes.sigma, 1, l) .* s
  # Quicker version as dimensions increase:
  snes.population = sampled_solutions = broadcast(+, snes.mu, broadcast(*, snes.sigma, s))

  # The rest of BlackBoxOptim still uses row-major order so transpose the
  # individuals before returning them.
  mix_with_indices( sampled_solutions', 1:snes.lambda )
end

function mix_with_indices(candidates, indices = false)
  if indices == false
    indices = 1:size(candidates,1)
  end
  ary = Any[]
  for i in indices
    push!(ary, (candidates[i,:], i))
  end
  ary
end

# Tell the sNES the ranking of a set of candidates.
function tell!(snes::SeparableNESOpt, rankedCandidates)
  u = assign_weights(rankedCandidates, snes.utilities)

  # Calc gradient
  gradient_mu = snes.last_s * u
  sq_s_minus1 = snes.last_s.^2 - 1
  gradient_sigma = sq_s_minus1 * u

  # Update the mean and sigma vectors based on the gradient
  snes.mu = snes.mu + snes.mu_learnrate * (snes.sigma .* gradient_mu)
  snes.sigma = snes.sigma .* exp(snes.sigma_learnrate / 2 * gradient_sigma)

  # There is no notion of how many was better in NES so return 0
  0
end

function assign_weights(rankedCandidates, u)
  n = length(rankedCandidates)
  # We must reorder the samples according to the order of the fitness in
  # rankedCandidates!!! Or we must reorder the utilities accordingly. The latter
  # is the preferred method and we can use the indices in rankedCandidates to
  # accomplish it.
  u_ordered = zeros(n, 1)
  for(i in 1:n)
    u_ordered[rankedCandidates[i][2]] = u[i]
  end

  u_ordered
end

# xNES is nice but scales badly with increasing dimension.
type XNESOpt <: NaturalEvolutionStrategyOpt
  d::Integer                      # Number of dimensions
  lambda::Integer                 # Number of samples to take per iteration
  utilities::Array{Float64,2}     # Fitness utility to give to each rank
  x_learnrate::Float64
  a_learnrate::Float64
  A::Array{Float64,2}
  expA::Array{Float64,2}
  population::Array{Float64,2}    # The last sampled values, now being evaluated
  x::Array{Float64,2}             # The current incumbent (aka most likely value, mu etc)
  Z::Array{Float64,2}

  XNESOpt(searchSpace; lambda = false) = begin
    d = numdims(searchSpace)
    lambda = lambda || convert(Int64, 4 + 3*floor(log(d)))
    x_learnrate = 1
    a_learnrate = 0.5 * minimum([1.0 / d, 0.25])
    x = rand_individual(searchSpace)

    new(d, lambda, fitness_shaping_utilities_log(lambda),
      x_learnrate, a_learnrate,
      zeros(d, d), zeros(d, d), zeros(d, lambda), x, zeros(d, d))
  end
end

function xnes(parameters)
  params = merge(NES_DefaultOptions, parameters)
  XNESOpt(params[:SearchSpace]; lambda = params["lambda"])
end

function ask(xnes::XNESOpt)
  xnes.expA = expm(xnes.A)
  xnes.Z = randn(xnes.d, xnes.lambda)
  xnes.population = broadcast(+, xnes.x, (xnes.expA * xnes.Z))

  # The rest of BlackBoxOptim still uses row-major order so transpose the
  # individuals before returning them.
  mix_with_indices( xnes.population', 1:xnes.lambda )
end

function tell!(xnes::XNESOpt, rankedCandidates)
  u = assign_weights(rankedCandidates, xnes.utilities)

  G = (repmat(u', xnes.d, 1) .* xnes.Z) * xnes.Z' - eye(xnes.d)
  dx = xnes.x_learnrate * xnes.expA * (xnes.Z*u)
  dA = xnes.a_learnrate * G

  xnes.x += dx
  xnes.A += dA

  0
end

# Calculate the fitness shaping utilities vector using the log method.
function fitness_shaping_utilities_log(n)
  u = maximum(hcat(zeros(n, 1), log(n / 2 + 1.0) - log(1:n)), 2)
  u / sum(u)
end

# Calculate the fitness shaping utilities vector using the steps method.
function fitness_shaping_utilities_linear(n)
  # Second half has zero utility.
  treshold = convert(Int64, floor(n/2))
  second_half = zeros(n - treshold, 1)

  # While first half's utility decreases in linear steps
  step_size = 1 / treshold
  first_half = linspace(1.0, step_size, treshold)

  # But the utilities should sum to 1 so we normalize, then return
  u = vcat(first_half, second_half)
  u / sum(u)
end


text:
# Helper function for changing initial value and time type of the following DDEs

"""
    remake_dde_constant_u0_tType(prob::DDEProblem, u₀, tType)

Create a new delay differential problem by replacing the initial state of problem `prob`
with `u0` and setting the type of time `t` to `tType`.

This function makes special assumptions about the structure of `prob` and is intended to be
used for all problems of name `prob_dde_constant_*` in `DDEProblemLibrary`. The functions of
these delay differential equation problems with constant delays are purposefully
implemented such that they work for arbitrary types of state `u` and time `t`, and hence in
particular for number types with units. The element type of `u` is saved as parameter `p` to
ensure that the return type of the history functions is correct (which are functions of `p`
but not of `u`).
"""
function remake_dde_constant_u0_tType(prob::DDEProblem, u₀, tType)
  remake(prob; u0 = u₀, tspan = tType.(prob.tspan), p = eltype(u₀),
         constant_lags = tType.(prob.constant_lags))
end

# History functions

function h_dde_constant_ip(p, t; idxs = 1)
  idxs == 1 || error("history function is only implemented for the first component")
  t ≤ zero(t) || error("history function is only implemented for t ≤ 0")

  p === nothing ? zero(t) : zero(p)
end

function h_dde_constant_oop(p, t)
  t ≤ zero(t) || error("history function is only implemented for t ≤ 0")

  p === nothing ? [zero(t)] : [zero(p)]
end

function h_dde_constant_scalar(p, t)
  t ≤ zero(t) || error("history function is only implemented for t ≤ 0")

  p === nothing ? zero(t) : zero(p)
end

# Single constant delay

## Short time span

### In-place function

@doc raw"""
    prob_dde_constant_1delay_ip

Delay differential equation

```math
u'(t) = -u(t - 1)
```

for ``t \in [0, 1]`` with history function ``\phi(t) = 0`` if ``t < 0`` and ``\phi(0) = 1``.

# Solution

The analytical solution for ``t \in [0, 10]`` can be obtained by the method of steps and
is provided in this implementation.
"""
prob_dde_constant_1delay_ip

function f_dde_constant_1delay_ip!(du, u, h, p, t)
  e = oneunit(t)
  du[1] = - h(p, t - e; idxs = 1) * inv(e)

  nothing
end

function fanalytic_dde_constant_1delay(u₀, ::Union{typeof(h_dde_constant_ip),
                                                   typeof(h_dde_constant_oop),
                                                   typeof(h_dde_constant_scalar)},
                                       p, t)
  z = t * inv(oneunit(t))
  0 ≤ z ≤ 10 || error("analytical solution is only implemented for t ∈ [0, 10]")

  if z < 1
    copy(u₀)
  else
    if z < 2
      c = @evalpoly(z, 2, -1)
    elseif z < 3
      c = @evalpoly(z, 4, -3, 1//2)
    elseif z < 4
      c = @evalpoly(z, 17//2, -15//2, 2, -1//6)
    elseif z < 5
      c = @evalpoly(z, 115//6, -109//6, 6, -5//6, 1//24)
    elseif z < 6
      c = @evalpoly(z, 1085//24, -1061//24, 197//12, -35//12, 1//4, -1//120)
    elseif z < 7
      c = @evalpoly(z, 13201//120, -13081//120, 521//12, -107//12, 1, -7//120, 1//720)
    elseif z < 8
      c = @evalpoly(z, 39371//144, -39227//144, 27227//240, -3685//144, 487//144,
                    -21//80, 1//90, -1//5040)
    elseif z < 9
      c = @evalpoly(z, 1158379//1680, -1156699//1680, 212753//720, -51193//720,
                    1511//144, -701//720, 1//18, -1//560, 1//40320)
    else
      c = @evalpoly(z, 23615939//13440, -23602499//13440, 7761511//10080,
                    -279533//1440, 89269//2880, -1873//576, 323//1440, -11//1120,
                    1//4032, -1//362880)
    end

    c .* u₀
  end
end

const prob_dde_constant_1delay_ip =
  DDEProblem(DDEFunction(f_dde_constant_1delay_ip!;
                         analytic = fanalytic_dde_constant_1delay),
             [1.0], h_dde_constant_ip, (0.0, 10.0), typeof(1.0);
             constant_lags = [1])

### Out-of-place function

"""
    prob_dde_constant_1delay_oop

Same delay differential equation as [`prob_dde_constant_1delay_ip`](@ref), but purposefully
implemented with an out-of-place function.
"""
prob_dde_constant_1delay_oop

function f_dde_constant_1delay_oop(u, h, p, t)
  e = oneunit(t)
  .- h(p, t - e) .* inv(e)
end

const prob_dde_constant_1delay_oop =
  DDEProblem(DDEFunction(f_dde_constant_1delay_oop;
                         analytic = fanalytic_dde_constant_1delay),
             [1.0], h_dde_constant_oop, (0.0, 10.0), typeof(1.0);
             constant_lags = [1])

### Scalar function

"""
    prob_dde_constant_1delay_scalar

Same delay differential equation as [`prob_dde_constant_1delay_ip`](@ref), but purposefully
implemented with a scalar function.
"""
prob_dde_constant_1delay_scalar

function f_dde_constant_1delay_scalar(u, h, p, t)
  e = oneunit(t)
  - h(p, t - e) * inv(e)
end

const prob_dde_constant_1delay_scalar =
  DDEProblem(DDEFunction(f_dde_constant_1delay_scalar;
                         analytic = fanalytic_dde_constant_1delay),
             1.0, h_dde_constant_scalar, (0.0, 10.0), typeof(1.0);
             constant_lags = [1])

## Long time span

### In-place function

@doc raw"""
    prob_dde_constant_1delay_long_ip

Delay differential equation

```math
u'(t) = u(t) - u(t - 1/5)
```

for ``t \in [0, 100]`` with history function ``\phi(t) = 0`` if ``t < 0`` and
``\phi(0) = 1``.
"""
prob_dde_constant_1delay_long_ip

function f_dde_constant_1delay_long_ip!(du, u, h, p, t)
  T = typeof(t)
  du[1] = (u[1] - h(p, t - T(1/5); idxs = 1)) * inv(oneunit(t))

  nothing
end

const prob_dde_constant_1delay_long_ip =
  DDEProblem(f_dde_constant_1delay_long_ip!, [1.0], h_dde_constant_ip, (0.0, 100.0),
             typeof(1.0);
             constant_lags = [1/5])

### Out-of-place function

"""
    prob_dde_constant_1delay_long_oop

Same delay differential equation as [`prob_dde_constant_1delay_long_ip`](@ref), but
purposefully implemented with an out-of-place function.
"""
prob_dde_constant_1delay_long_oop

function f_dde_constant_1delay_long_oop(u, h, p, t)
  T = typeof(t)
  (u .- h(p, t - T(1/5))) .* inv(oneunit(t))
end

const prob_dde_constant_1delay_long_oop =
  DDEProblem(f_dde_constant_1delay_long_oop, [1.0], h_dde_constant_oop, (0.0, 100.0),
             typeof(1.0);
             constant_lags = [1/5])

### Scalar function

"""
    prob_dde_constant_1delay_long_scalar

Same delay differential equation as [`prob_dde_constant_1delay_long_ip`](@ref), but
purposefully implemented with a scalar function.
"""
prob_dde_constant_1delay_long_scalar

function f_dde_constant_1delay_long_scalar(u, h, p, t)
  T = typeof(t)
  (u - h(p, t - T(1/5))) * inv(oneunit(t))
end

const prob_dde_constant_1delay_long_scalar =
  DDEProblem(f_dde_constant_1delay_long_scalar, 1.0, h_dde_constant_scalar, (0.0, 100.0),
             typeof(1.0);
             constant_lags = [1/5])

# Two constant delays

## Short time span

### In-place function

@doc raw"""
    prob_dde_constant_2delays_ip

Delay differential equation

```math
u'(t) = -u(t - 1/3) - u(t - 1/5)
```

for ``t \in [0, 1]`` with history function ``\phi(t) = 0`` if ``t < 0`` and ``\phi(0) = 1``.

# Solution

The analytical solution for ``t \in [0, 10]`` can be obtained by the method of steps and
is provided in this implementation.
"""
prob_dde_constant_2delays_ip

function f_dde_constant_2delays_ip!(du, u, h, p, t)
  T = typeof(t)
  du[1] = - (h(p, t - T(1/3); idxs = 1) + h(p, t - T(1/5); idxs = 1)) * inv(oneunit(t))

  nothing
end

function fanalytic_dde_constant_2delays(u₀, ::Union{typeof(h_dde_constant_ip),
                                                    typeof(h_dde_constant_oop),
                                                    typeof(h_dde_constant_scalar)},
                                        p, t)
  z = t * inv(oneunit(t))
  0 ≤ z ≤ 1 || error("analytical solution is only implemented for t ∈ [0, 1]")

  if z < 1/5
    copy(u₀)
  else
    if z < 1/3
      c = @evalpoly(z, 6//5, -1)
    elseif z < 2/5
      c = @evalpoly(z, 23//15, -2)
    elseif z < 8/15
      c = @evalpoly(z, 121//75, -12//5, 1//2)
    elseif z < 3/5
      c = @evalpoly(z, 427//225, -52//15, 3//2)
    elseif z < 2/3
      c = @evalpoly(z, 4351//2250, -547//150, 9//5, -1//6)
    elseif z < 11/15
      c = @evalpoly(z, 539//250, -647//150, 23//10, -1//6)
    elseif z < 4/5
      c = @evalpoly(z, 7942//3375, -128//25, 17//5, -2//3)
    elseif z < 13/15
      c = @evalpoly(z, 39998//16875, -1952//375, 89//25, -4//5, 1//24)
    elseif z < 14/15
      c = @evalpoly(z, 10109//3750, -1583//250, 243//50, -13//10, 1//24)
    else
      c = @evalpoly(z, 171449//60750, -139199//20250, 2579//450, -173//90, 5//24)
    end

    c .* u₀
  end
end

const prob_dde_constant_2delays_ip =
  DDEProblem(DDEFunction(f_dde_constant_2delays_ip!;
                         analytic = fanalytic_dde_constant_2delays),
             [1.0], h_dde_constant_ip, (0.0, 1.0), typeof(1.0);
             constant_lags = [1/3, 1/5])

### Out-of-place function

"""
    prob_dde_constant_2delays_oop

Same delay differential equation as [`prob_dde_constant_2delays_ip`](@ref), but purposefully
implemented with an out-of-place function.
"""
prob_dde_constant_2delays_oop

function f_dde_constant_2delays_oop(u, h, p, t)
  T = typeof(t)
  .- (h(p, t - T(1/3)) .+  h(p, t - T(1/5))) .* inv(oneunit(t))
end

const prob_dde_constant_2delays_oop =
  DDEProblem(DDEFunction(f_dde_constant_2delays_oop;
                         analytic = fanalytic_dde_constant_2delays),
             [1.0], h_dde_constant_oop, (0.0, 1.0), typeof(1.0);
             constant_lags = [1/3, 1/5])

### Scalar function

"""
    prob_dde_constant_2delays_scalar

Same delay differential equation as [`prob_dde_constant_2delays_ip`](@ref), but purposefully
implemented with a scalar function.
"""
prob_dde_constant_2delays_scalar

function f_dde_constant_2delays_scalar(u, h, p, t)
  T = typeof(t)
  - (h(p, t - T(1/3)) +  h(p, t - T(1/5))) * inv(oneunit(t))
end

const prob_dde_constant_2delays_scalar =
  DDEProblem(DDEFunction(f_dde_constant_2delays_scalar;
                         analytic = fanalytic_dde_constant_2delays),
             1.0, h_dde_constant_scalar, (0.0, 1.0), typeof(1.0);
             constant_lags = [1/3, 1/5])

## Long time span

### In-place function

@doc raw"""
    prob_dde_constant_2delays_long_ip

Delay differential equation

```math
u'(t) = - u(t - 1/3) - u(t - 1/5)
```

for ``t \in [0, 100]`` with history function ``\phi(t) = 0`` if ``t < 0`` and
``\phi(0) = 1``.
"""
prob_dde_constant_2delays_long_ip

function f_dde_constant_2delays_long_ip!(du, u, h, p, t)
  T = typeof(t)
  du[1] = - (h(p, t - T(1/3); idxs = 1) + h(p, t - T(1/5); idxs = 1)) * inv(oneunit(t))

  nothing
end

const prob_dde_constant_2delays_long_ip =
  DDEProblem(f_dde_constant_2delays_long_ip!, [1.0], h_dde_constant_ip, (0.0, 100.0),
             typeof(1.0);
             constant_lags = [1/3, 1/5])

### Out-of-place function

"""
    prob_dde_constant_2delays_long_oop

Same delay differential equation as [`prob_dde_constant_2delays_long_ip`](@ref), but
purposefully implemented with an out-of-place function.
"""
prob_dde_constant_2delays_long_oop

function f_dde_constant_2delays_long_oop(u, h, p, t)
  T = typeof(t)
  .- (h(p, t - T(1/3)) .+ h(p, t - T(1/5))) .* inv(oneunit(t))
end

const prob_dde_constant_2delays_long_oop =
  DDEProblem(f_dde_constant_2delays_long_oop, [1.0], h_dde_constant_oop, (0.0, 100.0),
             typeof(1.0);
             constant_lags = [1/3, 1/5])

#### Scalar function

"""
    prob_dde_constant_2delays_long_scalar

Same delay differential equation as [`prob_dde_constant_2delays_long_ip`](@ref), but
purposefully implemented with a scalar function.
"""
prob_dde_constant_2delays_long_scalar

function f_dde_constant_2delays_long_scalar(u, h, p, t)
  T = typeof(t)
  - (h(p, t - T(1/3)) + h(p, t - T(1/5))) * inv(oneunit(t))
end

const prob_dde_constant_2delays_long_scalar =
  DDEProblem(f_dde_constant_2delays_long_scalar, 1.0, h_dde_constant_scalar, (0.0, 100.0),
             typeof(1.0);
             constant_lags = [1/3, 1/5])


text:
# # Using UMLS Concepts with MeSH

#md # [![nbviewer](https://img.shields.io/badge/jupyter_notebook-nbviewer-orange.svg)](http://nbviewer.jupyter.org/github/bcbi/BioMedQuery.jl/tree/master/docs/src/notebooks/2_umls_map_and_filter.ipynb)

# The Medical Subject Headings (MeSH) terms returned from a PubMed search can be further analyzed
# by mapping them to Unified Medical Language System (UMLS) concepts, as well as
# filtering the MeSH Terms by concepts.

# For both mapping MeSH to UMLS Concepts and filtering MeSH by concept, the following backends are supported:
# * MySQL
# * SQLite
# * DataFrames

# ### Set Up
using SQLite
using MySQL
using BioMedQuery.DBUtils
using BioMedQuery.Processes
using BioServices.UMLS
using BioMedQuery.PubMed
using DataFrames

# Credentials are environment variables (e.g set in your .juliarc.jl)
umls_user = ENV["UMLS_USER"];
umls_pswd = ENV["UMLS_PSSWD"];
email = ""; # Only needed if you want to contact NCBI with inqueries
search_term = """(obesity[MeSH Major Topic]) AND ("2010"[Date - Publication] : "2012"[Date - Publication])""";
umls_concept = "Disease or Syndrome";
max_articles = 5;
results_dir = ".";
verbose = true;

results_dir = ".";

# ## MySQL

# ### Map Medical Subject Headings (MeSH) to UMLS

# This example demonstrates the typical workflow to populate a MESH2UMLS database
# table relating all concepts associated with all MeSH terms in the input database.

# *Note: this example reuses the MySQL DB from the PubMed Search and Save example.*
#
# Create MySQL DB connection
host = "127.0.0.1";
mysql_usr = "root";
mysql_pswd = "";
dbname = "pubmed_obesity_2010_2012";

const mysql_conn = DBUtils.init_mysql_database(host, mysql_usr, mysql_pswd, dbname) # hide
PubMed.create_tables!(mysql_conn) # hide
Processes.pubmed_search_and_save!(email, search_term, max_articles, mysql_conn, verbose) # hide
sleep(1) # hide
MySQL.disconnect(mysql_conn) #hide

db_mysql = MySQL.connect(host, mysql_usr, mysql_pswd, db = dbname);

# Map MeSH to UMLS
@time map_mesh_to_umls_async!(db_mysql, umls_user, umls_pswd; append_results=false, timeout=3);

# #### Explore the output table

db_query(db_mysql, "SELECT * FROM mesh2umls")

# ### Filtering MeSH terms by UMLS concept

# Getting the descriptor to index dictionary and the occurence matrix
@time labels2ind, occur = umls_semantic_occurrences(db_mysql, umls_concept);

# Descriptor to Index Dictionary
labels2ind

# Output Data Matrix
Matrix(occur)

MySQL.disconnect(db_mysql) #src

# ## SQLite

# This example demonstrates the typical workflow to populate a MESH2UMLS database
# table relating all concepts associated with all MeSH terms in the input database.

# *Note: this example reuses the SQLite DB from the PubMed Search and Save example.*
#
# Create SQLite DB connection
db_path = "$(results_dir)/pubmed_obesity_2010_2012.db";
db_sqlite = SQLite.DB(db_path);

if isfile(db_path) # hide
    rm(db_path) # hide
end # hide
db_sqlite = SQLite.DB(db_path); # hide
PubMed.create_tables!(db_sqlite); # hide
Processes.pubmed_search_and_save!(email, search_term, max_articles, db_sqlite, false) # hide
sleep(1) # hide

# ### Map MeSH to UMLS
@time map_mesh_to_umls_async!(db_sqlite, umls_user, umls_pswd; append_results=false, timeout=3);

# Explore the output table
db_query(db_sqlite, "SELECT * FROM mesh2umls;")

# ### Filtering MeSH terms by UMLS concept

# Getting the descriptor to index dictionary and occurence matrix
@time labels2ind, occur = umls_semantic_occurrences(db_sqlite, umls_concept);

# Descriptor to Index Dictionary
labels2ind

# Output Data Matrix
Matrix(occur)

# ## DataFrames

# This example demonstrates the typical workflow to create a MeSH to UMLS map as a DataFrame
# relating all concepts associated with all MeSH terms in the input dataframe.

# Get the articles (same as example in PubMed Search and Parse)
dfs = Processes.pubmed_search_and_parse(email, search_term, max_articles, verbose)

# Map MeSH to UMLS and explore the output table
@time res = map_mesh_to_umls_async(dfs["mesh_desc"], umls_user, umls_pswd)

# Getting the descriptor to index dictionary and occurence matrix
@time labels2ind, occur = umls_semantic_occurrences(dfs, res, umls_concept);

# Descriptor to Index Dictionary
labels2ind

# Output Data Matrix
Matrix(occur)


text:
"""
    TextDatasetRecipe(tablefile; catcols, contcols, kwargs...])

Recipe for loading a `TextDataset`. `tablefile` is the path of a file that can
be read as a table. `catcols` and `contcols` indicate the categorical and
continuous columns of the text table. If they are not given, they are detected
automatically.
"""


text:
"""
    f_c(m,Z,B)

Compute the (unsigned) cyclotron frequency of a charged species in MHz.
If `m` does not have units, it is assumed to be in kg.
If `B` does not have Unitful units, it is assumed to be in Tesla.
"""
f_c(m,Z,B) = abs(Z*qe*check_or_assume(B,u"Gauss")/check_or_assume(m,u"kg") /(2π)) |> u"MHz"

"""
    f_pe(m,Z,n)

Compute the plasma frequency of a charged species in MHz.
If `m` does not have units, it is assumed to be in kg.
If `n` does not have Unitful units, it is assumed to be in m^-3.
"""
f_p(m,Z,n) = sqrt(check_or_assume(n,u"m^-3")*Z^2*qe^2/(check_or_assume(m,u"kg")*ε0)) /(2π) |> u"MHz"

"""
    rL(m,Z,B,vperp)

Compute the gyro radius of a charged species in m.
If `m` does not have units, it is assumed to be in kg.
If `B` does not have Unitful units, it is assumed to be in Tesla.
If `vperp` does not have Unitful units, it is assumed to be in m/s.
"""
r_L(m,Z,B,vperp) = abs(check_or_assume(m,u"kg")*check_or_assume(vperp,u"m/s")/(Z*qe*check_or_assume(B,u"Gauss"))) |> u"m"

"""
    cth(m,Z,T)

Compute the thermal speed of a charged species in m/s.
This ignores numerical factors and gammae.
"""
c_th(m,Z,T) = sqrt(Z*force_energy_units!(check_or_assume(T,u"eV"))/check_or_assume(m,u"kg")) |> u"m/s"

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

"""
    f_ce(B)
    f_ce(P::QuasineutralPlasma)

Compute the (unsigned) cyclotron frequency of electrons in GHz.
If `B` does not have Unitful units, it is assumed to be in Tesla.
"""
f_ce(B) = f_c(me,1,B) |> u"GHz"
f_ce(P::QuasineutralPlasma) = f_ce(P.B)

"""
    f_ci(mi,Zi,B)
    f_ci(P::QuasineutralPlasma)

Compute the (unsigned) cyclotron frequency of ions of mass `mi` in kHz.
If `B` does not have Unitful units, it is assumed to be in Tesla.
"""
f_ci(mi,Zi,B) = f_c(mi,Zi,B) |> u"kHz"
f_ci(P::QuasineutralPlasma) = f_ci(P.mi,P.Zi,P.B)

"""
    f_pe(ne)
    f_pe(P::QuasineutralPlasma)

Compute the plasma frequency of electrons in GHz.
If `ne` does not have Unitful units, it is assumed to be in m^-3.
"""
f_pe(ne) = f_p(me,1,ne) |> u"GHz"
f_pe(P::QuasineutralPlasma) = f_pe(P.np)

"""
    f_pi(mi,Z,ni)
    f_pi(P::QuasineutralPlasma)

Compute the plasma frequency of ions of mass `mi` in MHz.
If `ni` does not have Unitful units, it is assumed to be in m^-3.
"""
f_pi(mi,Zi,ni) = f_p(mi,Zi,ni) |> u"kHz"
f_pi(P::QuasineutralPlasma) = f_pi(P.mi,P.Zi,P.np)

"""
    r_Le(B,vperpe)
    r_Le(P::QuasineutralPlasma)

Compute the electron gyro radius in mm.
If `B` does not have Unitful units, it is assumed to be in Tesla.
If `vperpe` does not have Unitful units, it is assumed to be in m/s.
"""
r_Le(B,vperpe) = r_L(me,1,B,vperpe) |> u"mm"
r_Le(P::QuasineutralPlasma) = r_Le(P.B,sqrt(P.Te/me))

"""
    r_Lsi(mi,Z,B,Te)
    r_Lsi(P::QuasineutralPlasma)

Compute the sonic ion gyro radius based on the sound speed, in cm.
If `mi` does not have units, it is assumed to be in kg.
If `B` does not have Unitful units, it is assumed to be in Tesla.
If `Te` does not have Unitful units, it is assumed to be in eV.
"""
r_Lsi(mi,Zi,B,Te) = r_L(mi,Zi,B,sqrt(force_energy_units!(check_or_assume(Te,u"eV"))/check_or_assume(mi,u"kg"))) |> u"cm"
r_Lsi(P::QuasineutralPlasma) = r_Lsi(P.mi,P.Zi,P.B,P.Te)

"""
    λ_De(ne,Te)
    λ_De(P::QuasineutralPlasma)

Compute the Debye length of the plasm in mm.
If ne does not have Unitful units, it is assumed to be in m^-3.
If Te does not have Unitful units, it is assumed to be in eV.
"""
λ_De(ne,Te) = sqrt(ε0*force_energy_units!(check_or_assume(Te,u"eV"))/(check_or_assume(ne,u"m^-3")*qe^2)) |> u"mm"
λ_De(P::QuasineutralPlasma) = λ_De(P.np,P.Te)

"""
    c_s(mi,Zi,Te)
    c_s(P::QuasineutralPlasma)

Compute the sound speed of ions in m/s.
This ignores numerical factors and gammae.
"""
c_s(mi,Zi,Te) = c_th(mi,Zi,Te) # Just an alias for cth
c_s(P::QuasineutralPlasma) = c_s(P.mi,P.Zi,P.Te)

"""
    c_the(Te)
    c_the(P::QuasineutralPlasma)

Compute the thermal speed of electrons in km/s.
This ignores numerical factors and gammae.
"""
c_the(Te) = c_th(me,1,Te) |> u"km/s"
c_the(P::QuasineutralPlasma) = c_the(P.Te)

"""
    CoulombLog(ne,Te)

Approximate the Coulomb logarithm.
"""
CoulombLog(np,Te) = 9+0.5*log(1e18u"m^-3"/np * Te^3/1u"eV^3") # approximation

"""
    ν_ei(Zi,np,Te)
    ν_ei(P::QuasineutralPlasma)

Electron-ion momentum collision frequency in MHz.
"""
function ν_ei(Zi,np,Te)
    np = check_or_assume(np,u"m^-3")
    Te = force_energy_units!(check_or_assume(Te,u"eV"))
    return sqrt(2) * np * Zi^2 * qe^4 * CoulombLog(np,Te) /
        (12π^(3/2) * ε0^2 * sqrt(me) * Te^(3/2)) |> u"MHz"
end
ν_ei(P::QuasineutralPlasma) = ν_ei(P.Zi,P.np,P.Te)

"""
    ν_ie(mi,Zi,np,Te)
    ν_ie(P::QuasineutralPlasma)

Ion-electron momentum collision frequency in MHz.
"""
ν_ie(mi,Zi,np,Te) = me/mi*ν_ei(Zi,np,Te)
ν_ie(P::QuasineutralPlasma) = ν_ie(P.mi,P.Zi,P.np,P.Te)

"""
    ν_ee(ne,Te)
    ν_ee(P::QuasineutralPlasma)

Electron-electron momentum collision frequency in MHz.
"""
function ν_ee(ne,Te)
    ne = check_or_assume(ne,u"m^-3")
    Te = force_energy_units!(check_or_assume(Te,u"eV"))
    return ne * qe^4 * CoulombLog(ne,Te) /
        (12π^(3/2) * ε0^2 * sqrt(me) * Te^(3/2)) |> u"MHz"
end
ν_ee(P::QuasineutralPlasma) = ν_ee(P.np,P.Te)

"""
    χ(Zi,np,Te,B)
    χ(P::QuasineutralPlasma)

Hall parameter based on the electron-ion collision frequency only.
"""
χ(Zi,np,Te,B) = f_ce(B)/ν_ei(Zi,np,Te) |> Unitful.NoUnits
χ(P::QuasineutralPlasma) = χ(P.Zi,P.np,P.Te,P.B)


text:
using Test
using BenchmarkTools

include("../src/Parsers.jl")

# Generates data for the Lorenz system.
ds = LorenzSystem()
data = datagen(ds)

# Generates data for a lotka volterra model
ds = LotkaVolterra()
data = datagen(ds)

# Generates data for a first order linear ODE
ds = ODE1()
data = datagen(ds)

# Generates data for a first order non linear ODE
ds = ODE3()
data = datagen(ds)

# Reads the data for the Helicopter Model
ds = HelicopterData("src/data/heli.csv")
data = datagen(ds)


text:
struct TrackedReal{T <: Real} <: Real
    value::T
    record::Record{T}
end

struct TrackedComplex{T <: Complex} <: Number
    value::T
    record::Record{T}
end

const TrackedNumber = Union{TrackedReal, TrackedComplex}

track(x::Real, record::Record=leaf(x)) = TrackedReal(x, record)
track(x::Complex, record::Record=leaf(x)) = TrackedComplex(x, record)

is_tracked(::Union{TrackedReal, TrackedComplex}) = true

function Base.show(io::IO, mime::MIME"text/plain", x::Union{TrackedReal, TrackedComplex})
    show(io, mime, x.value)
    print(io, " (tracked)")
end


text:
abstract type Inflow end

discharge(i::Inflow, t) =
    error("Please implement a discharge function for type ", typeof(i))

struct ConstantInflow <: Inflow
    Q::Float64
end

discharge(i::ConstantInflow, t) = i.Q

struct DiscreteHydrograph <: Inflow
    n::Int
    t::Vector{Float64}
    Q::Vector{Float64}
end

discharge(i::DiscreteHydrograph, t) = _linear_interpolation(i.n, i.t, i.Q, t)


text:
import Test: @test, @testset

@testset "Full 3D propagation" begin
using Luna
Luna.set_fftw_mode(:estimate)
import FFTW
import Luna.PhysData: wlfreq
import LinearAlgebra: norm

gas = :Ar
pres = 1

τ = 30e-15
λ0 = 800e-9

w0 = 500e-6
energy = 1e-9
L = 2

R = 4e-3
N = 128

grid = Grid.RealGrid(L, 800e-9, (400e-9, 2000e-9), 0.2e-12)
xygrid = Grid.FreeGrid(R, N)

x = xygrid.x
y = xygrid.y
energyfun, energyfunω = Fields.energyfuncs(grid, xygrid)

dens0 = PhysData.density(gas, pres)
densityfun(z) = dens0

responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)

linop = LinearOps.make_const_linop(grid, xygrid, PhysData.ref_index_fun(gas, pres))
normfun = NonlinearRHS.const_norm_free(grid, xygrid, PhysData.ref_index_fun(gas, pres))

inputs = Fields.GaussGaussField(λ0=λ0, τfwhm=τ, energy=energy, w0=w0)

Eω, transform, FT = Luna.setup(grid, xygrid, densityfun, normfun, responses, inputs)

output = Output.MemoryOutput(0, grid.zmax, 21)

Luna.run(Eω, grid, linop, transform, FT, output, max_dz=Inf, init_dz=1e-1)

Eout = output.data["Eω"] # (ω, ky, kx, z)

ω0idx = argmin(abs.(grid.ω .- wlfreq(λ0)))
λ0 = 2π*PhysData.c/grid.ω[ω0idx]
w1 = w0*sqrt(1+(L*λ0/(π*w0^2))^2)
Iω0_analytic = Maths.gauss.(xygrid.x, w1/2) # analytical solution (in paraxial approx)

Eω0yx = FFTW.ifft(Eout[ω0idx, :, :, end], (1, 2))
Iω0yx = abs2.(Eω0yx)
Iω0y = Maths.normbymax(dropdims(sum(Iω0yx, dims=2), dims=2))
Iω0x = Maths.normbymax(dropdims(sum(Iω0yx, dims=1), dims=1))

@test maximum(abs.(Iω0x .- Iω0_analytic)/norm(Iω0x)) < 5e-5
@test maximum(abs.(Iω0y .- Iω0_analytic)/norm(Iω0y)) < 5e-5

end


text:
using Test, Glob, MarketObservables

tests = glob("test_*")

for test_file ∈ tests
  include(test_file)
end


text:
"""
    load_scfres(filename)

Load back an `scfres`, which has previously been stored with `save_scfres`.
Note the warning in `save_scfres`.
"""
function load_scfres end

"""
    save_scfres(filename, scfres)

Save an `scfres` obtained from `self_consistent_field` to a JLD2 file.

!!! warning "No compatibility guarantees"
    No guarantees are made with respect to this function at this point.
    It may change incompatibly between DFTK versions or stop working / be removed
    in the future.
"""
function save_scfres end


text:
"""
Standard basic quantum circuit simulator building blocks.

This is a component package for [Yao.jl](https://github.com/QuantumBFS/Yao.jl). It contains the abstract definitions and basic implementation of Yao's circuit building blocks.
"""
module YaoBlocks

using YaoBase
using LinearAlgebra

import YaoBase:
    apply!,
    apply_back!,
    chcontent,
    chsubblocks,
    content,
    dispatch!,
    expect,
    fidelity,
    focus!,
    getiparams,
    iparams_eltype,
    iscommute,
    isreflexive,
    isunitary,
    mat,
    mat_back!,
    niparams,
    nqubits,
    nqudits,
    nlevel,
    occupied_locs,
    operator_fidelity,
    parameters,
    parameters_eltype,
    print_block,
    render_params,
    setiparams!,
    subblocks,
    ishermitian,
    nparameters

export AbstractBlock,
    AbstractContainer,
    CompositeBlock,
    LocationConflictError,
    PrimitiveBlock,
    QubitMismatchError,
    RemoveMeasured,
    ResetTo,
    TagBlock,
    apply!,
    apply,
    apply_back!,
    chcontent,
    chsubblocks,
    content,
    dispatch!,
    dispatch,
    expect,
    getiparams,
    iparams_eltype,
    iscommute,
    isreflexive,
    isunitary,
    mat,
    mat_back!,
    niparams,
    nqubits,
    nqudits,
    nlevel,
    occupied_locs,
    operator_fidelity,
    parameters,
    parameters_eltype,
    print_block,
    render_params,
    setiparams!,
    setiparams,
    subblocks,
    ishermitian,
    nparameters

export applymatrix, cache_key

include("utils.jl")
# include("traits.jl")

include("abstract_block.jl")

# concrete blocks
include("routines.jl")
include("primitive/primitive.jl")
include("composite/composite.jl")

include("algebra.jl")
include("blocktools.jl")
include("measure_ops.jl")
include("layout.jl")
include("treeutils/treeutils.jl")

include("autodiff/autodiff.jl")
export AD, Optimise

include("deprecations.jl")

end # YaoBlocks


text:
using ShaderToy, FileIO, GLAbstraction

shadertoy("clouds.frag", Dict{Symbol, Any}(
	:iChannel0 => Texture(load("noise.png").data, x_repeat=:repeat, minfilter=:linear),
))