text:
# Bogumił Kamiński, 2021

# Codes for chapter 2

# Code for listing 2.1

1
true
"Hello world!"
0.1
[1, 2, 3]

# Code for listing 2.2

typeof(1)
typeof(true)
typeof("Hello world!")
typeof(0.1)
typeof([1, 2, 3])

# Code for showing bit representation of numbers

bitstring(1)
bitstring(1.0)
bitstring(Int8(1))

# Code showing to what Int alias expands

Int

# Code for checking if value is of some type

[1, 2, 3] isa Vector{Int}
[1, 2, 3] isa Array{Int64, 1}

# Code for section 2.2

x = 1
y = [1, 2, 3]

x = 1
x
typeof(x)
x = 0.1
x
typeof(x)

Kamiński = 1
x₁ = 0.5
ε = 0.0001

?₁
?ε

# Code for listing 2.3

x = -7
if x > 0
    println("positive")
elseif x < 0
    println("negative")
elseif x == 0
    println("zero")
else
    println("unexpected condition")
end

# Code showing that logical condition must be Bool

x = -7
if x
    println("condition was true")
end

# Code showing comparisons against NaN

NaN > 0
NaN >= 0
NaN < 0
NaN <= 0
NaN == 0

NaN != 0
NaN != NaN

# Code showing that floating point arithmetic is only approximate

0.1 + 0.2 == 0.3

0.1 + 0.2

isapprox(0.1 + 0.2, 0.3)

0.1 + 0.2 ≈ 0.3

# Code showing combining conditions

x = -7
x > 0 && x < 10
x < 0 || log(x) > 10

x = -7
log(x)

# Code showing typical one-line conditional execution expressions

x > 0 && println(x)

if x > 0
    println(x)
end

x > 0 || println(x)

if !(x > 0)
    println(x)
end

x = -7
x < 0 && println(x^2)
iseven(x) || println("x is odd")

x = -7
if x < 0
    println(x^2)
end
if !iseven(x)
    println("x is odd")
end

x = -7
if x < 0 && x^2
    println("inside if")
end

# Code showing ternary operator

x > 0 ? sqrt(x) : sqrt(-x)

if x > 0
    sqrt(x)
else
    sqrt(-x)
end

x = -7
x > 0 ? println("x is positive") : println("x is not positive")

# Code showing that conditional statements return a value

x = -4.0
y = if x > 0
        sqrt(x)
    else
        sqrt(-x)
    end
y

x = 9.0
y = x > 0 ? sqrt(x) : sqrt(-x)
y

# Code from listing 2.4

for i in [1, 2, 3]
    println(i, " is ", isodd(i) ? "odd" : "even")
end

# Code from listing 2.5

i = 1
while i < 4
    println(i, " is ", isodd(i) ? "odd" : "even")
    global i += 1
end

# Code showing break and continue keywords

i = 0
while true
    global i += 1
    i > 6 && break
    isodd(i) && continue
    println(i, " is even")
end

# Code from listing 2.6

x = -7
x < 0 && begin
    println(x)
    x += 1
    println(x)
    2 * x
end
x > 0 ? (println(x); x) : (x += 1; println(x); x)

# Code from section 2.3.4

x = [8, 3, 1, 5, 7]
k = 1

y = sort(x)

for i in 1:k
    y[i] = y[k + 1]
    y[end - i + 1] = y[end - k]
end
y

s = 0
for v in y
    global s += v
end
s
s / length(y)

# Code from listing 2.7

function times_two(x)
    return 2 * x
end
times_two(10)

# Code from listing 2.8

function compose(x, y=10; a, b=10)
    return x, y, a, b
end
compose(1, 2; a=3, b=4)
compose(1, 2; a=3)
compose(1; a=3)
compose(1)
compose(; a=3)

# Code from listing 2.9

times_two(x) = 2 * x
compose(x, y=10; a, b=10) = x, y, a, b

# Code showing the use of map function

map(times_two, [1, 2, 3])

# Code from listing 2.10

map(x -> 2 * x, [1, 2, 3])

# Code showing sum taking a function as a first argument

sum(x -> x ^ 2, [1, 2, 3])

# Code showing do-end syntax

sum([1, 2, 3]) do x
    println("processing ", x)
    return x ^ 2
end

# Code showing the difference between sort and sort!

x = [5, 1, 3, 2]
sort(x)
x
sort!(x)
x

# Code showing a simple implementation of winsorized_mean function

function winsorized_mean(x, k)
    y = sort(x)
    for i in 1:k
        y[i] = y[k + 1]
        y[end - i + 1] = y[end - k]
    end
    s = 0
    for v in y
        s += v
    end
    return s / length(y)
end
winsorized_mean([8, 3, 1, 5, 7], 1)

# Code from section 2.5

function fun1()
    x = 1
    return x + 1
end
fun1()
x

function fun2()
    if true
        x = 10
    end
    return x
end
fun2()

function fun3()
    x = 0
    for i in [1, 2, 3]
        if i == 2
            x = 2
        end
    end
    return x
end
fun3()

function fun4()
    for i in [1, 2, 3]
        if i == 2
            x = 2
        end
    end
    return x
end
fun4()

function fun5()
    for i in [1, 2, 3]
        if i == 1
            x = 1
        else
            x += 1
        end
        println(x)
    end
end
fun5()

function fun6()
    x = 0
    for i in [1, 2, 3]
        if i == 1
            x = 1
        else
            x += 1
        end
        println(x)
    end
end
fun6()

max_stars_repo_name: Mo-Gul/JuliaForDataAnalysis

max_issues_repo_name: Mo-Gul/JuliaForDataAnalysis

max_forks_repo_name: Mo-Gul/JuliaForDataAnalysis


text:
#=
  From BLOG Swift example birthday.blog
  https://github.com/lileicc/swift/blob/master/example/birthday.blog
  Result
  """
  Distribution of values for exists Person x (exists Person y ((!(x = y) & (Birthday(x) = Birthday(y)))))
	true	0.5446927374301685
	false	0.4553072625698328
  """

=#
using Turing, StatsPlots, DataFrames, MCMCChains
include("jl_utils.jl")

@model function birthday(num_people=missing)
    # var day = _.range(365);
    num_days = 365

    if num_people === missing
        num_people ~ truncated(Poisson(20),1,num_days)
    end

    # Note: I call it zbirthdays so it's displayed last
    zbirthdays = Vector{Int}(undef, num_people)
    for d in 1:num_people
       zbirthdays[d] ~ DiscreteUniform(1,num_days)
    end

    #  is there a duplicate of birthdays?
    # hakank: I'm not sure if the is the best (or even good) way
    #         to model this...
    b = make_hash(zbirthdays) # Collect the birthdays
    s = sum(values(b) .> 1)
    if s == 0
        t ~ Bernoulli(0)
    else
        t ~ Bernoulli(1)
    end

end

model = birthday(23)
# model = birthday()

num_chns = 4
println("Prior distribution")
chns_prior = sample(model, Prior(), MCMCThreads(), 1000, num_chns)
display(chns_prior)
# display(plot(chns_prior))

println("\nPosterior distribution")
num_chns = 4

# chns = sample(model, MH(), MCMCThreads(), 10_000, num_chns)
# chns = sample(model, MH(), MCMCThreads(), 30_000, num_chns)

chns = sample(model, IS(), MCMCThreads(), 10_000, num_chns) # ORIG

# chns = sample(model, PG(20), MCMCThreads(), 1000, num_chns)

# chns = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chns)

## LoadError: MethodError: Cannot `convert` an object of type Float64 to an object of type Symbol
# chns = sample(model, Gibbs(NUTS(0.65,5,1,1,0,:lambda_1,:lambda_2), PG(20,:tau) ), 1000)
# chns = sample(model, Gibbs(HMC(0.1,5,:lambda_1,:lambda_2), PG(20,:tau) ), MCMCThreads(), 1000, num_chns)

# Error
# chns = sample(model, Gibbs(HMCDA(0.15,0.65,0.3,:lambda_1,:lambda_2), PG(20,:tau) ), MCMCThreads(), 1000, num_chns) # Error...

# Skipping this since it contains a lot of uninteresting zbirthdays[b]
# display(chns)
println("t:")
display(chns[:t])

# plot(chns)

if :num_people in chns.name_map.parameters
    println("num_people:")
    display(chns[:num_people])
end

# df = DataFrame(chns)
# display(df)

ext: jl

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using Waterfall
import Luxor
import Base
include(joinpath(WATERFALL_DIR,"bin","io.jl"))
# include(joinpath(WATERFALL_DIR,"src","utils","correlate.jl"))
include(joinpath(WATERFALL_DIR,"src","utils","draw.jl"))
# include(joinpath(WATERFALL_DIR,"src","Figure","Axis.jl"))
# include(joinpath(WATERFALL_DIR,"src","Figure","Plot.jl"))
# include(joinpath(WATERFALL_DIR,"src","Cascade","Vertical.jl"))
# include(joinpath(WATERFALL_DIR,"src","Cascade","Parallel.jl"))

filename(p::Plot) = filename(p.cascade)

function filename(x::Cascade{T}) where T <: Geometry
    path = joinpath(WATERFALL_DIR,"fig","perm",lowercase(string(T)))
    !isdir(path) && mkpath(path)
    order = [Printf.@sprintf("-%02d", perm) for perm in x.permutation]
    file = lowercase(string(T, order..., ".png"))
    path = joinpath(path, file)
    println("Saving $path")
    return path
end

nstep = 3
nsample = 1
nperm = 8
ncor = 1000
N = 13
permute = true
perms = [[collect(1:N)]; Waterfall.random_permutation(1:N, min(factorial(N),nperm))]

perm = perms[2]
# args = (quantile=1.0,)

# for T in [Parallel,Horizontal,Vertical]
for T in [Vertical]
    for perm in perms
        x = Cascade(df; permutation=perm, interactivity=(0.01,0.3), nsample=nsample, ncor=ncor, kwargs...)

        pdata = Plot(copy(x); ylabel="Efficiency (%)")
        p = convert(Plot{T}, pdata, 1.0; permute=true)

        Luxor.@png begin
            Luxor.fontsize(18)
            Luxor.setmatrix([1 0 0 1 LEFT_BORDER TOP_BORDER])

            draw(p; distribution=distribution, samples=nsample)
        end WIDTH+LEFT_BORDER+RIGHT_BORDER HEIGHT+TOP_BORDER+BOTTOM_BORDER filename(p)
    end
end

# # _convert(::Type{T}, x, args...; kwargs...) where T<:Geometry = _convert!(T, copy(x), args...; kwargs...)

# # function offdiag_rowonal(v::Vector, ii)
# # end

# # function Base.convert(::Type{Matrix}, lst) where T<:Real
# #     return _convert(::Type{Matrix}, lst, 2)
# # end

# # function permute!(x::Cascade)
# # end

# # xv = convert(Cascade{Vertical}, x)
# #


text:
instance, device = init(; device_extensions = ["VK_KHR_synchronization2"])
frag_shader = resource("decorations.frag")

@testset "Shader cache" begin
    spec = ShaderSpecification(frag_shader, GLSL)
    cache = ShaderCache(device)
    find_shader!(ShaderCache(device), spec) # trigger JIT compilation
    t = @elapsed find_shader!(cache, spec)
    @test t > 0.01
    t = @elapsed find_shader!(cache, spec)
    @test t < 1e-5
end

@testset "Descriptors" begin
    da = DescriptorAllocator(device)
    spec = ShaderSpecification(frag_shader, GLSL)
    cache = ShaderCache(device)
    shader = Rhyolite.find_shader!(cache, spec)
    dset_layout_cache = DescriptorSetLayoutCache(device)
    sets = allocate_descriptor_sets!(da, dset_layout_cache, [shader])
    @test length(da.pools) == 1
    pool_state = first(da.pools)
    @test pool_state.allocated == dictionary([DESCRIPTOR_TYPE_STORAGE_IMAGE => 1])

    # check resources are reused by the descriptor allocator
    sets2 = allocate_descriptor_sets!(da, dset_layout_cache, [shader])
    @test pool_state.allocated == dictionary([DESCRIPTOR_TYPE_STORAGE_IMAGE => 2])

    # check that new sets were allocated
    @test handle.(sets) ≠ handle.(sets2)

    # check resources get cleaned up
    free_descriptor_sets!(da, [sets; sets2])

    @test pool_state.allocated == dictionary([DESCRIPTOR_TYPE_STORAGE_IMAGE => 0])
end


text:
using WebIO
using Test

@testset "bundling javascript" begin
    # TODO: if we want to be picky, we can even test that
    # bundlejs produces exactly the same output stored in the
    # current bundle. But that requires pinning down exact
    # versions of all the dependencies, which we don't currently do.
    WebIO.bundlejs()
end


text:
export arrowClass, generalClass, getTConsHead
export InstRec, GlobalTC, LocalTC, show_hints
const TCState = Base._return_type(mk_tcstate, (Vector{HMT}, ))
Gensym = Symbol

"""
datatype for type class instances
"""
mutable struct InstRec
    t :: HMT             # instance type
    ln :: LineNumberNode # where the instance defined
    gensym :: ExprImpl   # the actual gensym(IR level) to access this instance
    isPruned :: Bool     # is the instance type pruned?
end
@typedIO InstRec(HMT, LineNumberNode, ExprImpl, Bool)

const InstResolCtx = LinkedList{InstRec}

function toVec(t::Type{LinkedList{E}}, x::Cons{E}) where E
    [toVec(E, x.head), toVec(t, x.tail)]
end

function toVec(::Type{LinkedList{E}}, ::Nil{E}) where E
    nothing
end

function fromVec(::Type{LinkedList{E}}, x::Nothing) where E
    Some(nil(E))
end

function fromVec(t::Type{LinkedList{E}}, x::Vector) where E
    length(x) !== 2 && return nothing
    head = fromVec(E, x[1])
    head === nothing && return nothing
    tail = fromVec(t, x[2])
    tail === nothing && return nothing
    Some(cons(head.value, tail.value))
end

InstRec(t::HMT, ln::LineNumberNode, gensym::Symbol, isPruned::Bool) = InstRec(t, ln, IR.EVar(gensym), isPruned)

const arrowClass = Nom(:arrow_class)
const generalClass = Nom(:general_class)

getTConsHead(h::HMT) =
    @match h begin
        Nom(_) => h
        App(t, _) ||
        Forall(_, t) => getTConsHead(t)
        Arrow(Implicit(_), t) => getTConsHead(t)
        Arrow(_, _) => arrowClass
        _ => generalClass
    end

struct GlobalTC
    tcstate :: TCState
    count :: Ref{UInt}
    globalImplicits :: Dict{Nom, Vector{InstRec}}
    globalImplicitDeltas :: Dict{Nom, Int32}
    queries :: Vector{Pair{String, HMT}}
    moduleName :: Symbol
end

Base.show(io::IO, ::GlobalTC) = Base.show(io, "<GlobalTC>")

Base.empty(::Type{GlobalTC}, moduleName::Symbol = :main) =
    GlobalTC(
        mk_tcstate(HMT[]),
        Ref(UInt(0)),
        Dict{Nom, Vector{InstRec}}(),
        Dict{Nom, UInt32}(),
        Pair{String, HMT}[],
        moduleName
    )

function show_hints(globalTC::GlobalTC)
    prune = globalTC.tcstate.prune
    for (k, v) in globalTC.queries
        println(k, ": ",  prune(v))
    end
end
struct LocalTC
    typeEnv :: Store{Symbol, HMT}
    symmap :: Store{Symbol, Gensym}
    localImplicits :: InstResolCtx
    ln :: LineNumberNode
end

Base.empty(::Type{LocalTC}) =
    LocalTC(
        Store{Symbol, HMT}(),
        Store{Symbol, Symbol}(),
        nil(InstRec),
        LineNumberNode(1, :unknown))

const αβ = UInt8[('a':'z')...]

symgen(g::GlobalTC) =
    let i = g.count[] + 1
        g.count[] = i
        chars = UInt8[]
        while i != 0
            j = i%26
            push!(chars, αβ[j + 1])
            i = div(i, 26)
        end
        Symbol(:|, g.moduleName, :|, Symbol(chars))
    end


text:
function sayhello()
   println("hi I am worker number $(myid()), I live on $(gethostname())")
end

max_stars_repo_name: yutec/FightBrand

max_issues_repo_name: yutec/FightBrand

max_forks_repo_name: yutec/FightBrand


text:
"""
    notebook()

Lorem ipsum dolor sit amet.
"""
function notebook(args::Dict)
  notebook()
end

function notebook()

  target_name = bump()

  Base.run(`jupyter notebook --notebook-dir=./lib/notebooks`)

end


text:
__precompile__(false)
module JuliaVenv

using Pkg

const venv_depot = joinpath(@__DIR__, "depot")

const orig_DEPOT_PATH = Ref{Vector{String}}()
const orig_exename = Ref{String}()
const orig_bindir = Ref{String}()

const current_exename = Ref{String}()

function __init__()
    ### Set `DEPOT_PATH[1] = HERE/depot`
    orig_DEPOT_PATH[] = copy(Base.DEPOT_PATH)

    # Default `Base.DEPOT_PATH[1] = "~/.julia"` (`Base.init_depot_path`)
    default_depot1 = joinpath(homedir(), ".julia")

    # `DEPOT_PATH` without "~/.julia":
    depots = filter(!isequal(default_depot1), Base.DEPOT_PATH)
    # TODO: Maybe add an option to "overlay" venv depot on top of the
    #       standard depot including ~/.julia?
    # TODO: Handle `haskey(ENV, "JULIA_DEPOT_PATH")`?

    if isempty(depots) || depots[1] != venv_depot
        pushfirst!(depots, venv_depot)
    end
    append!(empty!(Base.DEPOT_PATH), depots)

    ### Prepare for `monkey_patching_julia_cmd`:
    orig_bindir[] = Sys.BINDIR
    orig_exename[] = Base.julia_exename()
    current_exename[] = orig_exename[]

    @eval JuliaVenv Base.julia_exename() = current_exename[]
end

"""
    monkey_patching_julia_cmd(f)

Do `f` while `Base.julia_cmd` is patched to return the command using
`julia-venv`.
"""
function monkey_patching_julia_cmd(f)
    try
        current_exename[] = "julia-venv"
        @eval Sys BINDIR = $(dirname(ENV["JULIA_VENV_PYTHON"]))

        @info "Monkey patched julia_cmd() = $(Base.julia_cmd())"
        return f()
    finally
        current_exename[] = orig_exename[]
        @eval Sys BINDIR = $(orig_bindir[])
    end
end

"""
    precompile_package(pkg::Union{String, Vector{String}})
    precompile_package(packages::AbstractVector{Union{AbstractString, Symbol}})

Precompile package(s) using `julia-venv`.
"""
precompile_package(package::AbstractString) = precompile_package([package])

function precompile_package(packages::AbstractVector{<: AbstractString})
    pkgids = Base.identify_package.(packages)
    not_found = findall(isequal(nothing), pkgids)
    if !isempty(not_found)
        error("Package(s) not fond: $(packages[not_found])")
    end
    monkey_patching_julia_cmd() do
        for id in pkgids
            Base.compilecache(id)
        end
    end
end

"""
    add(pkg::Union{String, Vector{String}})
    add(pkg::Union{PackageSpec, Vector{PackageSpec}})

Like `Pkg.add` but precompile packages using `precompile_package`, so
that it's safe to import them afterward in `julia-venv`.
"""
add(package::AbstractString) = add([package])
add(packages::AbstractVector{<: AbstractString}) = add(PackageSpec.(packages))
add(package::PackageSpec) = add([package])

function add(packages::AbstractVector{PackageSpec})
    Pkg.add(copy(packages))  # Pkg.add mutates `packages`
    precompile_package([p.name for p in packages])
end

end  # module

max_stars_repo_name: tkf/julia-venv

max_issues_repo_name: tkf/julia-venv

max_forks_repo_name: tkf/julia-venv


text:
module Router # Bukdu

import ..Bukdu: Naming, Deps, Routing, DirectRequest, request_handler

"""
    Router.call

To get the object of the action without going through HTTP server.
Results are named tuple `(got= , resp= , route= )`.
"""
function call
end

"""
    Router.call(verb, path::String)
"""
function call(verb, path::String)
    method = Naming.verb_name(verb)
    req = Deps.Request(method, path)
    call(req)
end

"""
    Router.call(req::Deps.Request)
"""
function call(req::Deps.Request)
    route = Routing.handle(req)
    dreq = DirectRequest(req)
    request_handler(route, dreq) # (got=, resp=)
end

module Helpers # Bukdu.Router

import ...Bukdu: ApplicationController, Routing

struct URLPathError <: Exception
    msg
end

function url_path(verb, C::Type{<:ApplicationController}, action)
    path = Routing.route_path(verb, C, action)
    if path isa Nothing
        URLPathError(string("failed to get the path: ", join([verb, C, action], ' ')))
    else
        path
    end
end

end # module Bukdu.Router.Helpers

end # module Bukdu.Router

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
using Random, JSON, Statistics

noFunc() = nothing

include("utils/GeneticProgramming/nodes.jl")
include("utils/GeneticProgramming/nodeSetsParser.jl")
include("utils/miscFunctions.jl")
include("tools/individual.jl")
include("tools/experimentSummary.jl")

"
Auxiliary structure that keeps some generic information for GP experiments.
It is the abstract type from which every other GP experiment info inherits.

See also: [`ExperimentInfo`](@ref), [`CGPInfo`](@ref), [`STGPInfo`](@ref),
[`GEPInfo`](@ref), [`GEInfo`](@ref)
"
abstract type GPExperimentInfo end
# abstract type

"""
    getVariables(gpExperimentInfo::GPExperimentInfo)

Gets the identifiers of the variables of GP problem.

# Arguments
- `gpExperimentInfo::GPExperimentInfo`: information about the GP experiment.

# Returns
Array with variables' identifiers as `strings`.
"""
getVariables(gpExperimentInfo::GPExperimentInfo) = gpExperimentInfo._variables
# function

"
Auxiliary structure that keeps some generic information for experiments.

# Fields
- `_individualType::DataType`: type of the individual that determines which
    methods are going to be used.
- `_GPExperimentInfo::GPExperimentInfo`: if the problem to solve is GP, this field
    will contain its information.
- `_randomSeed::UInt`: seed used for rand functions. By default, it will be 1234.
- `_rng::Random.AbstractRNG`: random number generator used by the system. By
    default, it will be set to Random.GLOBAL_RNG.
- `_algorithm::Function`: genetic algorithm used for the experiment. By default,
    it will be set to the classic genetic algorithm (see [`basicExperiment`](@ref)).
- `_algorithmArgs::Array{Any}`: arguments for the genetic algorithm, if needed.

See also: [`GeneticJulia`](@ref)
"
mutable struct ExperimentInfo
    _individualType::DataType
    _GPExperimentInfo::GPExperimentInfo
    _experimentSummary::ExperimentSummary
    _randomSeed::UInt
    _rng::Random.AbstractRNG
    _algorithm::Function
    _algorithmArgs::Array{Any}

    ExperimentInfo() = new()
end # struct

include("tools/stopCondition.jl")
include("tools/evaluator.jl")
include("tools/initialization.jl")
include("tools/selector.jl")
include("tools/crossover.jl")
include("tools/mutation.jl")
include("tools/replacement.jl")

"
Main structure of the system, gathers all the tools necessary to run a
population-based algorithm. By default, the main structure for `code users` is
called `GenJ`.

# Fields
- `_experimentInfo::ExperimentInfo`: auxiliary structure that keeps some generic
    information for experiments.
- `_stopCondition::StopCondition`: tool that defines the stop conditions and
    checks them throughout the experiment.
- `_population::Array{Individual}`: array that stores the population as Individual.
- `_evaluator::Evaluator`: defines the methods and structure for evaluate individuals.
- `_generator::Generator`: defines the methods and structure for generate individuals.
- `_selector::SelectionOperator`: defines the methods and structure for select
    \"parents\" among the population.
- `_crossoverOp::CrossoverOperator`: defines the methods and structure for cross
    \"parents\" and obtain the offspring.
- `_mutationOp::MutationOperator`: defines the methods and structure for mutate
    an offspring's individual.
- `_replacementOp::ReplacementOperator`: defines the methods and structure for
    replace the old population with the offspring.
"
mutable struct GeneticJuliaStructure
    _experimentInfo::ExperimentInfo
    _stopCondition::StopCondition
    _population::Array{Individual}
    _evaluator::Evaluator
    _generator::Generator
    _selector::SelectionOperator
    _crossoverOp::CrossoverOperator
    _mutationOp::MutationOperator
    _replacementOp::ReplacementOperator

    GeneticJuliaStructure() = new()
end # struct

# Global structure generation for code user transparency.
GenJ = GeneticJuliaStructure()

# Due to optional options, it is created here.
GenJ._experimentInfo = ExperimentInfo()
GenJ._experimentInfo._randomSeed = 1234
GenJ._experimentInfo._rng = Random.GLOBAL_RNG

include("utils/GeneticProgramming/includes.jl")
include("utils/GeneticAlgorithms/includes.jl")
include("utils/algorithms.jl")
include("utils/defaultSettings.jl")
include("utils/error.jl")
include("utils/parser.jl")
include("utils/verboseFunctions.jl")
include("interface/includes.jl")

"""
    getIndividualType(experimentInfo::ExperimentInfo)

Gets the individual type of the main structure.
"""
getIndividualType(experimentInfo::ExperimentInfo) = experimentInfo._individualType
 # function

"""
    getIndividualType(genj::GeneticJulia)

Gets the individual type of the main structure.
"""
getIndividualType(genj::GeneticJuliaStructure) = getIndividualType(genj._experimentInfo)
 # function

"""
    setIndividualType(type::DataType; genj::GeneticJulia = GenJ)

Sets the type of the individuals' representation for the experiment.

# Arguments
- `type::DataType`: data type of the individuals of the experiment. Should be
    one of the implemented ones (see [`IndividualRep`](@ref)). The user could
    provide their own type, as long as they also provide their own methods for
    population generation, crossover and mutation.

# Keyword Arguments
- `genj::GeneticJulia = GenJ`: the main structure. For `code user`: don't modify
    the default value unless multiple experiments are going to be run.

# Examples
```jldoctest
julia> setIndividualType(BinaryGenotype)
BinaryGenotype
```
"""
function setIndividualType(type::DataType; genj::GeneticJuliaStructure = GenJ)

    if findfirst(x->type <: x, [BinaryGenotype, IntegerGenotype, CGPGenotype,
                                          STGPGenotype, GEPGenotype, GEGenotype]) == nothing
        @warn "The individual type provided does not match any of the predefined
               ones. Personalized methods for generation, crossover and mutation
               shall be provided"
    end

    genj._experimentInfo._individualType = type
end # function

"""
    setAlgorithm(algorithm::Function, varArgs...; genj::GeneticJulia = GenJ)

Sets the genetic algorithm that is going to be run for the experiment.

# Arguments
- `algorithm::Function`: genetic algorithm that is going to run the experiment.
    By default, it is set to the classic one, which follows this steps:
    - Generation of the population.
    - Fitness evaluation.
    - Parent selection.
    - Crossover.
    - Mutation.
    - Fitness evaluation.
    - Replacement of the population.
- `varArgs`: arguments for the algorithm, if needed.

# Keyword Arguments
- `genj::GeneticJulia = GenJ`: the main structure. For `code user`: don't modify
    the default value unless multiple experiments are going to be run.

# Examples
```jldoctest
julia> setAlgorithm(basicExperiment)
```
"""
function setAlgorithm(algorithm::Function, varArgs...; genj::GeneticJuliaStructure = GenJ)
    genj._experimentInfo._algorithm = algorithm
    varArgs = [i for i in varArgs]
    genj._experimentInfo._algorithmArgs = varArgs
end # function

"""
    setRandomSeed(seed::UInt; genj::GeneticJulia = GenJ)

Sets the random number generator as a MersenneTwister generator initialized with
the specified seed.

# Arguments
- `seed::Integer`: seed to initialize the random number generator with.

# Keyword Arguments
- `genj::GeneticJulia = GenJ`: the main structure. For `code user`: don't modify
    the default value unless multiple experiments are going to be run.

# Examples
```jldoctest
julia> setRandomSeed(2)
```
"""
function setRandomSeed(seed::Integer; genj::GeneticJuliaStructure = GenJ)
    if seed < 0
        error("RandomSeed ($seed) must be a positive integer")
    end
    genj._experimentInfo._randomSeed = convert(UInt, seed)
    genj._experimentInfo._rng = Random.MersenneTwister(seed)
end # function

"""
    getRNG(experimentInfo::ExperimentInfo)

Gets the individual type of the main structure.
"""
getRNG(experimentInfo::ExperimentInfo) = experimentInfo._rng
 # function

"""
    getRNG(genj::GeneticJulia)

Gets the random number generator of the main structure.
"""
getRNG(genj::GeneticJuliaStructure) = getRNG(genj._experimentInfo)
# function

"""
    setRNG(rng::Random.AbstractRNG; genj::GeneticJulia = GenJ)

Sets the random number generator for the system to run the experiment.

# Arguments
- `seed::Integer`: seed to initialize the random number generator with.

# Keyword Arguments
- `genj::GeneticJulia = GenJ`: the main structure. For `code user`: don't modify
    the default value unless multiple experiments are going to be run.

# Examples
```jldoctest
julia> rng = Random.MersenneTwister()

julia> setRNG(rng)
```
"""
function setRNG(rng::Random.AbstractRNG; genj::GeneticJuliaStructure = GenJ)
    genj._experimentInfo._rng = rng
end # function

"""
    runGenJ(jsonFile::String; verbose::Bool = true, outputFile::String = "")

Runs the experiments specified in the configuration file.

# Arguments
- `configFile::String`: the configuration file (in JSON format) with all the information about the
    experiments that are going to be run.

# Keyword Arguments
- `verbose::Bool = true`: specifies wether information about the experiment is going
    to be displayed or not.
- `outputFile::String = ""`: specifies the output stream in which the information
    displayed by verbose is going to be dumped.

# Returns
The array of experiments run.

For guidance about how to create the information file about the experiments
visit our GitHub repository at:
[https://github.com/SergioGmezM/Evolutionary\\_Computation\\_Julia](https://github.com/SergioGmezM/Evolutionary_Computation_Julia)
"""
function runGenJ(configFile::String; verbose::Bool = true, outputFile::String = "")
    experiments = generateMainStructure(configFile)
    nExperiments = length(experiments)

    if !verbose
        verbose = outputFile != ""
    end

    for i=1:nExperiments
        # Sets the default settings according to the individual type provided in the experiment
        setDefaultSettings(genj=experiments[i], i=i)
        # Checks that all has been set
        checkParametersDefined(experiments[i], i)
        # Checks that all has been set correctly
        checkParametersCorrect(experiments[i], i)
    end

    for i=1:nExperiments
        println("Experiment $i \n==========================")
        # Prints the information about the experiment
        verbose && printExperimentInfo(outputFile, experiments[i])
        # Runs the experiment
        experiments[i]._experimentInfo._algorithm(experiments[i]._experimentInfo._algorithmArgs..., genj=experiments[i])
    end

    return experiments
end # function

"""
    runGenJ(genj::GeneticJulia = GenJ; verbose::Bool = true, outputFile::String = "")

Runs a single experiment for `code users`.

# Arguments
- `genj::GeneticJulia = GenJ`: the main structure. For `code user`: don't modify
    the default value unless multiple experiments are going to be run.

# Keyword Arguments
- `verbose::Bool = true`: specifies wether information about the experiment is going
    to be displayed or not.
- `outputFile::String = ""`: specifies the output stream in which the information
    displayed by verbose is going to be dumped.

For guidance about how to create the information file about the experiments
visit our GitHub repository at:
[https://github.com/SergioGmezM/Evolutionary\\_Computation\\_Julia](https://github.com/SergioGmezM/Evolutionary_Computation_Julia)
"""
function runGenJ(genj::GeneticJuliaStructure = GenJ; verbose::Bool = true, outputFile::String = "")

    if !verbose
        verbose = outputFile != ""
    end

    # Sets the default settings according to the individual type provided in the experiment
    setDefaultSettings(genj=genj)
    # Checks that all has been set
    checkParametersDefined(genj)
    # Checks that all has been set correctly
    checkParametersCorrect(genj)
    # Prints the information about the experiment
    verbose && printExperimentInfo(outputFile, genj)
    # Runs the experiment
    genj._experimentInfo._algorithm(genj._experimentInfo._algorithmArgs..., genj=GenJ)
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
# ***** KEEP *****
# The paper begins with an example in which swapping the first two (plotted) steps makes a
# high impact on the values of these two steps. In this file, we examine the impact of
# plotting ONLY the first two steps -- which could be any two of N total steps.
using Waterfall
include(joinpath(WATERFALL_DIR,"src","includes.jl"))
include(joinpath(WATERFALL_DIR,"bin","io.jl"))

localkwargs = (
    nsample = 1000,
    ncor = true,
    correlate = true,
    permute = true,
)

function idk(df, swaps; kwargs...)
    cascades = Dict(swap => define_from(Cascade{Data}, df; permutation=swap, kwargs...)
        for swap in swaps)
    return Dict(swap => collect_value(cascade) for (swap,cascade) in cascades)
end

# Find all combinations of four steps AND swap the first two elements.
N = size(df,1)-2
xsrt = collect(Combinatorics.combinations(1:N,4))
xrev = swapat!.(copy.(xsrt),1,2)

vsrt = idk(df, xsrt; localkwargs..., kwargs...)
vrev = idk(df, xrev; localkwargs..., kwargs...)

# Summarize in DataFrame. What's the greatest difference in MEAN FINAL VALUE between [1,2]
dfsrt = DataFrames.DataFrame(srt=xsrt, rev=xrev, vsrt_2=0.0, vsrt_3=0.0, vsrt_end=0.0)
[dfsrt[ii,:vsrt_2] = Statistics.mean(vsrt[dfsrt[ii,:srt]][2,:]) for ii in 1:size(dfsrt,1)]
[dfsrt[ii,:vsrt_3] = Statistics.mean(vsrt[dfsrt[ii,:srt]][3,:]) for ii in 1:size(dfsrt,1)]
[dfsrt[ii,:vsrt_end] = Statistics.mean(vsrt[dfsrt[ii,:srt]][end,:]) for ii in 1:size(dfsrt,1)]

dfrev = DataFrames.DataFrame(rev=xrev, vrev_2=0.0, vrev_3=0.0, vrev_end=0.0)
[dfrev[ii,:vrev_2] = Statistics.mean(vrev[dfrev[ii,:rev]][2,:]) for ii in 1:size(dfrev,1)]
[dfrev[ii,:vrev_3] = Statistics.mean(vrev[dfrev[ii,:rev]][3,:]) for ii in 1:size(dfrev,1)]
[dfrev[ii,:vrev_end] = Statistics.mean(vrev[dfrev[ii,:rev]][end,:]) for ii in 1:size(dfrev,1)]

dfsum = DataFrames.innerjoin(dfsrt, dfrev, on=:rev)
dfsum[!,:diff_2] .= abs.(dfsum[:,:vsrt_2] .- dfsum[:,:vrev_2])
dfsum[!,:diff_3] .= abs.(dfsum[:,:vsrt_3] .- dfsum[:,:vrev_3])
dfsum[!,:diff_end] .= abs.(dfsum[:,:vsrt_end] .- dfsum[:,:vrev_end])

dfsum[!,:comp] .= dfsum[:,:diff_2] .* dfsum[:,:diff_3] .* dfsum[:,:diff_end]
DataFrames.select!(dfsum, [:srt,:rev,:comp,:diff_2,:diff_3,:diff_end])

# ii = collect(1:N)
# swaps = collect(Combinatorics.combinations(1:N,2))
# swaps = swaps[getindex.(swaps,1).==1]

# perms = Dict(swap => swapat!(collect(1:N), swap...) for swap in swaps)

# nsample = 1000

# v = Dict(swap => update_stop!(collect_value(
#     define_from(Cascade{Data}, df[[1:N+1;size(df,1)],:];
#         nsample=nsample,
#         ncor=false,
#         correlate=false,
#         permute=true,
#         permutation=perms[swap],
#         kwargs...,
#         interactivity=(0.1,0.2),
#     ))) for swap in swaps)

# # Summary DataFrame of all of the pairwise swaps that could be made.
# N = length(swaps)
# df = DataFrames.DataFrame(swap=swaps,)
# df[!,:cumulative] .= [cumulative_v(v[df[ii,:swap]]) for ii in 1:N]
# df[!,:max] .= [Statistics.maximum(df[ii,:cumulative][2:end-1,1]) for ii in 1:N]
# df[!,:min] .= [Statistics.minimum(df[ii,:cumulative][2:end-1,1]) for ii in 1:N]
# df[!,:diff] .= df[:,:max] .- df[:,:min]

# df[!,:endmax] .= [Statistics.maximum(df[ii,:cumulative][end-2:end-1,:]) for ii in 1:N]
# df[!,:endmin] .= [Statistics.minimum(df[ii,:cumulative][end-2:end-1,:]) for ii in 1:N]
# df[!,:enddiff] .= df[:,:endmax] .- df[:,:endmin]

# #
# # swaps = collect(Combinatorics.combinations(1:N,2))


text:
mutable struct Gen
    GEN_BUS::Vector{Int64}  # "bus number"
    PG::Vector{Float64}     # "Pg, real power output (MW)"
    QG::Vector{Float64}     # "Qg, reactive power output (MVAr)"
    QMAX::Vector{Float64}   # "Qmax, maximum reactive power output at Pmin (MVAr)"
    QMIN::Vector{Float64}   # "Qmin, minimum reactive power output at Pmin (MVAr)"
    VG::Vector{Float64}     # "Vg, voltage magnitude setpoint (p.u.)"
    MBASE::Vector{Float64}  # "mBase, total MVA base of this machine, defaults to baseMVA"
    GEN_STATUS::Vector{Bool}    # "status, 1 - machine in service, 0 - machine out of service"
    PMAX::Vector{Float64}   # "Pmax, maximum real power output (MW)"
    PMIN::Vector{Float64}   # "Pmin, minimum real power output (MW)"
    PC1::Vector{Float64}    # "Pc1, lower real power output of PQ capability curve (MW)"
    PC2::Vector{Float64}    # "Pc2, upper real power output of PQ capability curve (MW)"
    QC1MIN::Vector{Float64} # "Qc1min, minimum reactive power output at Pc1 (MVAr)"
    QC1MAX::Vector{Float64} # "Qc1max, maximum reactive power output at Pc1 (MVAr)"
    QC2MIN::Vector{Float64} # "Qc2min, minimum reactive power output at Pc2 (MVAr)"
    QC2MAX::Vector{Float64} # "Qc2max, maximum reactive power output at Pc2 (MVAr)"
    RAMP_AGC::Vector{Float64}   # "ramp rate for load following/AGC (MW/min)"
    RAMP_10::Vector{Float64}    # "ramp rate for 10 minute reserves (MW)"
    RAMP_30::Vector{Float64}    # "ramp rate for 30 minute reserves (MW)"
    RAMP_Q::Vector{Float64}     # "ramp rate for reactive power (2 sec timescale) (MVAr/min)"
    APF::Vector{Float64}    # "area participation factor"
    MU_PMAX::Vector{Float64}    # "Kuhn-Tucker multiplier on upper Pg limit (u/MW)"
    MU_PMIN::Vector{Float64}    # "Kuhn-Tucker multiplier on lower Pg limit (u/MW)"
    MU_QMAX::Vector{Float64}    # "Kuhn-Tucker multiplier on upper Qg limit (u/MVAr)"
    MU_QMIN::Vector{Float64}    # "Kuhn-Tucker multiplier on lower Qg limit (u/MVAr)"
end

"""
Extract gen information from mpc Dict. Return instance of Gen.

IN:

* mpc. Instance of MPC dict (loaded, for instance, by loadcase())

OUT:

* Instance of Gen containing all gen data.
"""
function extract_gen(mpc)
    gen = mpc["gen"]
    nc = size(gen, 2)
    nf = length(fieldnames(Gen))
    blanks = fill([NaN], nf - nc)
    return Gen([gen[:, i] for i in 1:nc]..., blanks...)
end

"""
Shortcut for returning gen info for one of the cases listed by
`casenames()`.
"""
function extract_gen(cname::AbstractString)
    mpc = loadcase(cname; describe=false)
    extract_gen(mpc)
end


text:
# Romeo LTI systems
#
# System type hierarchy.

# TODO: Parallel and series type systems should be able to have multiple inner systems instead of just two.

# TODO: Consider introducing `PolynomialSystem` as a special case of `ParallelSystem`.

"""
    Base type representing all SISO LTI systems.
"""
abstract type SisoLtiSystem end

Base.:+(sys::SisoLtiSystem) = sys
Base.:-(sys::SisoLtiSystem) = -1 * sys

"""
    Zero-system

A system that blocks all input signals: ideal no-pass system.
"""
struct ZeroSys <: SisoLtiSystem end

"""
    Unit-system

A system that passes all input signals: ideal all-pass system.
"""
struct UnitSys <: SisoLtiSystem end

"""
    Differintegral.

Differentiator of order `α` (when α>0), integrator of order `-α` (when α<0), or
identity operator (when α=0).

Note that when the order `α` is actually zero, a special object of type `UnitSys`
is returned.
"""
struct Diff{T<:Number} <: SisoLtiSystem
    α::T

    Diff{T}(α) where {T<:Number} = (α ≠ 0) ? new(α) : UnitSys()
end

Diff(α::T) where {T<:Number} = Diff{T}(α)

Base.:(==)(first::Diff{<:Number}, second::Diff{<:Number}) = first.α == second.α

"""
    A SISO LTI system obtained by series connection of a scaler (pure gain) and another SISO LTI system.
"""
struct ScaledSystem{T<:Number} <: SisoLtiSystem
    k::T
    inner::SisoLtiSystem

    function ScaledSystem{T}(k, inner) where {T<:Number}
        if k == 1
            inner
        elseif k == 0 || inner == ZeroSys()
            ZeroSys()
        else
            new(k, inner)
        end # if
    end # function
end

ScaledSystem(k::T, inner::SisoLtiSystem) where {T<:Number} = ScaledSystem{T}(k, inner)
ScaledSystem(k::T, inner::ScaledSystem{<:Number}) where {T<:Number} = ScaledSystem{T}(k * inner.k, inner.inner)

Base.:(==)(first::ScaledSystem{<:Number}, second::ScaledSystem{<:Number}) =
    first.inner == second.inner && first.k == second.k

Base.:*(num::Number, sys::SisoLtiSystem) = ScaledSystem(num, sys)
Base.:*(sys::SisoLtiSystem, num::Number) = num * sys

"""
    A SISO LTI system obtained by series connection of two other SISO LTI systems.
"""
struct SeriesSystem <: SisoLtiSystem
    first::SisoLtiSystem
    second::SisoLtiSystem

    SeriesSystem(first::SisoLtiSystem, second::SisoLtiSystem) = new(first, second)
    SeriesSystem(first::SisoLtiSystem, second::UnitSys) = first
    SeriesSystem(first::UnitSys, second::SisoLtiSystem) = second
    SeriesSystem(first::SisoLtiSystem, second::ZeroSys) = ZeroSys()
    SeriesSystem(first::ZeroSys, second::SisoLtiSystem) = ZeroSys()
    SeriesSystem(first::SisoLtiSystem, second::ScaledSystem{<:Number}) = ScaledSystem(second.k, first * second.inner)
    SeriesSystem(first::ScaledSystem{<:Number}, second::SisoLtiSystem) = ScaledSystem(first.k, first.inner * second)
    SeriesSystem(first::ScaledSystem{<:Number}, second::ScaledSystem{<:Number}) = ScaledSystem(first.k * second.k, first.inner * second.inner)
    SeriesSystem(first::Diff{<:Number}, second::Diff{<:Number}) = Diff(first.α + second.α)
end

Base.:*(first::SisoLtiSystem, second::SisoLtiSystem) = SeriesSystem(first, second)

"""
    A SISO LTI system obtained by parallel connection of two other SISO LTI systems.
"""
struct ParallelSystem <: SisoLtiSystem
    first::SisoLtiSystem
    second::SisoLtiSystem

    ParallelSystem(first::SisoLtiSystem, second::SisoLtiSystem) = (first == second) ? 2 * first : new(first, second)
    ParallelSystem(first::SisoLtiSystem, second::ZeroSys) = first
    ParallelSystem(first::ZeroSys, second::SisoLtiSystem) = second
    ParallelSystem(first::SisoLtiSystem, second::ScaledSystem{<:Number}) = (second.inner == first) ? (second.k + 1) * first : new(first, second)
    ParallelSystem(first::ScaledSystem{<:Number}, second::SisoLtiSystem) = ParallelSystem(second, first)
    ParallelSystem(first::ScaledSystem{<:Number}, second::ScaledSystem{<:Number}) = (first.inner == second.inner) ? (first.k + second.k) * first.inner : new(first, second)
end

Base.:+(first::SisoLtiSystem, second::SisoLtiSystem) = ParallelSystem(first, second)
Base.:+(sys::SisoLtiSystem, num::Number) = sys + num * Diff(0)
Base.:+(num::Number, sys::SisoLtiSystem) = sys + num

Base.:-(first::SisoLtiSystem, second::SisoLtiSystem) = first + (-second)
Base.:-(sys::SisoLtiSystem, num::Number) = sys + (-num)
Base.:-(num::Number, sys::SisoLtiSystem) = num + (-sys)

"""
    A SISO LTI system obtained by taking a power of another system in the complex domain.
"""
struct PowerSystem{T<:Number} <: SisoLtiSystem
    α::T
    inner::SisoLtiSystem

    function PowerSystem{T}(α::T, inner::SisoLtiSystem) where {T<:Number}
        if α == 0
            UnitSys()
        elseif α == 1
            inner
        else
            new(α, inner)
        end # if
    end # function

    PowerSystem{T}(α::T, inner::ZeroSys) where {T<:Number} = (α ≠ 0) ? ZeroSys() : UnitSys()
    PowerSystem{T}(α::T, inner::UnitSys) where {T<:Number} = UnitSys()
    PowerSystem{T}(α::T, inner::Diff{<:Number}) where {T<:Number} = Diff(inner.α * α)
    PowerSystem{T}(α::T, inner::ScaledSystem{<:Number}) where {T<:Number} = ScaledSystem(inner.k^α, inner.inner^α)
end

PowerSystem(α::T, inner::SisoLtiSystem) where {T<:Number} = PowerSystem{T}(α, inner)

Base.:(==)(first::PowerSystem{<:Number}, second::PowerSystem{<:Number}) =
    first.inner == second.inner && first.α == second.α

Base.sqrt(sys::SisoLtiSystem) = PowerSystem(1 // 2, sys)
Base.:^(sys::SisoLtiSystem, n::Number) = PowerSystem(n, sys)

"""
    A SISO LTI system obtained by dividing two SISO LTI systems in the complex domain.
"""
struct RationalSystem <: SisoLtiSystem
    num::SisoLtiSystem
    den::SisoLtiSystem

    function RationalSystem(num::SisoLtiSystem, den::SisoLtiSystem)
        if num == den
            UnitSys()
        else
            new(num, den)
        end
    end

    RationalSystem(num::SisoLtiSystem, den::UnitSys) = num
    RationalSystem(num::SisoLtiSystem, den::ZeroSys) = @error "Cannot create rational systems with zero denominator."
    RationalSystem(num::SisoLtiSystem, den::ScaledSystem{<:Number}) = ScaledSystem(1 / den.k, num / den.inner)
    RationalSystem(num::SisoLtiSystem, den::PowerSystem{<:Number}) = (den.inner == num) ? PowerSystem(1 - den.α, den.inner) : new(num, den)

    RationalSystem(num::UnitSys, den::SisoLtiSystem) = PowerSystem(-1, den)
    RationalSystem(num::ZeroSys, den::SisoLtiSystem) = ZeroSys()

    RationalSystem(num::ScaledSystem{<:Number}, den::SisoLtiSystem) = ScaledSystem(num.k, num.inner / den)
    RationalSystem(num::ScaledSystem{<:Number}, den::ScaledSystem{<:Number}) = ScaledSystem(num.k / den.k, num.inner / den.inner)
    RationalSystem(num::Diff{<:Number}, den::Diff{<:Number}) = Diff(num.α - den.α)

    RationalSystem(num::PowerSystem{<:Number}, den::UnitSys) = num
    RationalSystem(num::UnitSys, den::PowerSystem{<:Number}) = PowerSystem(-den.α, den.inner)
    RationalSystem(num::PowerSystem{<:Number}, den::SisoLtiSystem) = (num.inner == den) ? PowerSystem(num.α - 1, num.inner) : new(num, den)
    RationalSystem(num::PowerSystem{<:Number}, den::PowerSystem{<:Number}) = (num.inner == den.inner) ? PowerSystem(num.α - den.α, num.inner) : new(num, den)

    # TODO: Improve construction of rational systems from other rational systems.
end

Base.:/(num::SisoLtiSystem, den::SisoLtiSystem) = RationalSystem(num, den)
Base.:/(sys::SisoLtiSystem, den::Number) = ScaledSystem(1 / den, sys)
Base.:/(num::Number, sys::SisoLtiSystem) = RationalSystem(num * Diff(0), sys)

Base.:\(den::SisoLtiSystem, num::SisoLtiSystem) = RationalSystem(num, den)
Base.:\(sys::SisoLtiSystem, num::Number) = RationalSystem(num * Diff(0), sys)
Base.:\(den::Number, sys::SisoLtiSystem) = ScaledSystem(1 / den, sys)

export SisoLtiSystem, Diff, ScaledSystem, ParallelSystem, SeriesSystem, RationalSystem, PowerSystem


text:
using Test
using DiffEqOperators, OrdinaryDiffEq

@testset "KdV equation (Single Solition)" begin
    N = 21
    Δx = 1/(N-1)
    r = 0.5

    # x = 10:Δx:30;
    x = -10:Δx:10;
    # ϕ(x,t) = (r/2)*sech.((sqrt(r)*(x-r*t)/2)-7).^2 # solution of the single forward moving wave
    ϕ(x,t) = (1/2)*sech.((x-t)/2).^2 # solution of the single forward moving wave

    u0 = ϕ(x,0);
    oriu = zeros(size(x));

    const du3 = zeros(size(x));
    const temp = zeros(size(x));

    # A = DerivativeOperator{Float64}(1,2,Δx,length(x),:Dirichlet0,:Dirichlet0);
    A = UpwindOperator{Float64}(1,3,Δx,length(x),true.|BitVector(undef,length(x)),
                                :Dirichlet0,:Dirichlet0);
    # C = DerivativeOperator{Float64}(3,2,Δx,length(x),:Dirichlet0,:Dirichlet0);
    C = UpwindOperator{Float64}(3,3,Δx,length(x),true.|BitVector(undef,length(x)),
                                :Dirichlet0,:Dirichlet0);

    function KdV(du, u, p, t)
       C(t,u,du3)
       A(t,u,du)
       @. temp = -0.5*u*du - 0.25*du3
       copyto!(du,temp)
    end

    single_solition = ODEProblem(KdV, u0, (0.,5.));
    soln = solve(single_solition,Tsit5(),abstol=1e-6,reltol=1e-6);

    #=
    for t in 0:0.5:5
        @show maximum(soln(t)-ϕ(x,t))
    end
    =#

    for t in 0:0.5:5
        @test_skip soln(t) ≈ ϕ(x,t) atol = 0.01;
    end
end

# Conduct interesting experiments by referring to
# http://lie.math.brocku.ca/~sanco/solitons/kdv_solitons.php
@testset "KdV equation (Double Solition)" begin
    N = 10
    Δx = 1/(N-1)

    x = -50:Δx:50;
    c1,c2 = 25,16
    # ϕ1(x,t) = 3*p1^2*sech(.5*(p1*(x+2+t))).^2+3*p2^2*sech(.5*(p2*(x+1+t))).^2;

    function ϕ(x,t)
        # t = t-10
        num1 = 2(c1-c2)*(c1*cosh.(√c2*(x-c2*t)/2).^2 + c2*sinh.(√c1*(x-c1*t)/2).^2)
        den11 = (√c1-√c2)*cosh.((√c1*(x-c1*t) + √c2*(x-c2*t))/2)
        den12 = (√c1+√c2)*cosh.((√c1*(x-c1*t) - √c2*(x-c2*t))/2)
        den1 = (den11+den12).^2
        return num1./den1
    end

    const du3 = zeros(x);
    const temp = zeros(x);

    A = UpwindOperator{Float64}(1,1,Δx,length(x),false.*BitVector(length(x)),
                                :Dirichlet0,:nothing);
    C = UpwindOperator{Float64}(3,1,Δx,length(x),false.*BitVector(length(x)),
                                :Dirichlet0,:nothing);

    function KdV(du, u, p, t)
       C(t,u,du3)
       A(t, u, du)
       @. temp = -6*u*du - du3
       copyto!(du,temp)
    end

    u0 = ϕ(x,0);
    double_solition = ODEProblem(KdV, u0, (0.,5.));
    soln = solve(double_solition,Tsit5());

    # The solution is a forward moving soliton wave with speed = 1
    for t in 0:0.1:9
        @test_skip soln(t) ≈ ϕ(x,t)
    end
end


text:
# This file includes code originally part of Julia.n
# Licensed under MIT License, see LICENSE.md

# The plan is to rewrite all of the functionality to not use the utf8proc library,
# and to use tables loaded up on initialization, as with StringLiterals.jl

# Currently, this depends on the following C calls:
# utf8proc_errmsg
# utf8proc_decompose
# utf8proc_reencode

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

function utf8proc_map(::Type{T}, str::MaybeSub{T}, options::Integer) where {T<:Str}
    nwords = ccall(:utf8proc_decompose, Int, (Ptr{UInt8}, Int, Ptr{UInt8}, Int, Cint),
                   str, sizeof(str), C_NULL, 0, options)
    nwords < 0 && utf8proc_error(nwords)
    buffer = Base.StringVector(nwords*4)
    nwords = ccall(:utf8proc_decompose, Int, (Ptr{UInt8}, Int, Ptr{UInt8}, Int, Cint),
                   str, sizeof(str), buffer, nwords, options)
    nwords < 0 && utf8proc_error(nwords)
    nbytes = ccall(:utf8proc_reencode, Int, (Ptr{UInt8}, Int, Cint), buffer, nwords, options)
    nbytes < 0 && utf8proc_error(nbytes)
    Str(cse(T), String(resize!(buffer, nbytes)))
end

utf8proc_map(str::MaybeSub{T}, options::Integer) where {T<:Str} =
    utf8proc_map(T, convert(UTF8Str, str), options)
utf8proc_map(str::MaybeSub{<:Str{UTF8CSE}}, options::Integer) =
    utf8proc_map(UTF8Str, str, options)

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

function _normalize(::Type{T}, str::AbstractString;
                   stable::Bool      = false,
                   compat::Bool      = false,
                   compose::Bool     = true,
                   decompose::Bool   = false,
                   stripignore::Bool = false,
                   rejectna::Bool    = false,
                   newline2ls::Bool  = false,
                   newline2ps::Bool  = false,
                   newline2lf::Bool  = false,
                   stripcc::Bool     = false,
                   casefold::Bool    = false,
                   lump::Bool        = false,
                   stripmark::Bool   = false,
                   ) where {T<:Str}
    (compose & decompose) && strerror(StrErrors.DECOMPOSE_COMPOSE)
    (!(compat | stripmark) & (compat | stripmark)) && strerror(StrErrors.COMPAT_STRIPMARK)
    newline2ls + newline2ps + newline2lf > 1 && strerror(StrErrors.NL_CONVERSION)
    flags =
        ifelse(stable,      Uni.STABLE, 0) |
        ifelse(compat,      Uni.COMPAT, 0) |
        ifelse(decompose,   Uni.DECOMPOSE, 0) |
        ifelse(compose,     Uni.COMPOSE, 0) |
        ifelse(stripignore, Uni.IGNORE, 0) |
        ifelse(rejectna,    Uni.REJECTNA, 0) |
        ifelse(newline2ls,  Uni.NLF2LS, 0) |
        ifelse(newline2ps,  Uni.NLF2PS, 0) |
        ifelse(newline2lf,  Uni.NLF2LF, 0) |
        ifelse(stripcc,     Uni.STRIPCC, 0) |
        ifelse(casefold,    Uni.CASEFOLD, 0) |
        ifelse(lump,        Uni.LUMP, 0) |
        ifelse(stripmark,   Uni.STRIPMARK, 0)
    T(utf8proc_map(str, flags))
end

normalize(str::T, options::Integer) where {T<:Str} = _normalize(T, UTF8Str(str), options)
normalize(str::Str{UTF8CSE}, options::Integer) = _normalize(UTF8Str, str, options)

normalize(str::Str, nf::Symbol) =
    utf8proc_map(str, nf == :NFC ? (Uni.STABLE | Uni.COMPOSE) :
                 nf == :NFD ? (Uni.STABLE | Uni.DECOMPOSE) :
                 nf == :NFKC ? (Uni.STABLE | Uni.COMPOSE | Uni.COMPAT) :
                 nf == :NFKD ? (Uni.STABLE | Uni.DECOMPOSE | Uni.COMPAT) :
                 strerror(StrErrors.NORMALIZE, nf))


text:
############ Historic Forex Ticks  ####################
"""
    historic_forex_ticks(opts::PolyOpts, from="AUD", to="USD", date="2020-10-14"; limit=100, kwargs...)

Get historic ticks for a forex currency pair.

# Arguments
 * opts::PolyOpts: The PolyOpts object used to configure the request.
 * from: The "from" symbol of the currency pair. Example: For USD/JPY the from would be USD.
 * to: The "to" symbol of the currency pair. Example: For USD/JPY the to would be JPY.
 * date: The date/day of the historic ticks to retrieve.
 * limit: TLimit the size of the response, max 10000. Defualt 100.
 * kwargs: Additional parameters to pass to the Poly API.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> historic_forex_ticks(opts, "AUD", "AUD", "2020-10-14")
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array
 * See https://polygon.io/docs/get_v1_historic_forex__from___to___date__anchor for documentation on response attributes and supported keyword arguments.
"""
function historic_forex_ticks(opts::PolyOpts, from="AUD", to="AUD", date="2020-10-14"; limit=100, kwargs...)
    params = Dict(
        "apiKey" => opts.api_key,
        "limit" => limit
        )

    merge!(params, Dict(kwargs))

    historic_forex_ticks_url = "$historic_forex_ticks_base_url/$from/$to/$date"
    return generate_output_from_url(YesSinkYesTicks(), historic_forex_ticks_url, params, opts.sink)
end

############ Real Time Currency Conversion  ####################
"""
    real_time_currency_conversion(opts::PolyOpts, from="AUD", to="USD"; amount=100, precision=2)

Get currency conversions using the latest market conversion rates.
Note than you can convert in both directions. For example USD to CAD or CAD to USD.

# Arguments
 * opts::PolyOpts: The PolyOpts object used to configure the request.
 * from: The "from" symbol of the pair.
 * to: The "to" symbol of the pair.
 * amount: The amount to convert.
 * precision: The decimal precision of the conversion. Defaults to 2 which is 2 decimal places accuracy.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> real_time_currency_conversion(opts, "AUD", "USD"; amount=100, precision=2)
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array
 * See https://polygon.io/docs/get_v1_conversion__from___to__anchor for documentation on response attributes and supported keyword arguments.
"""
function real_time_currency_conversion(opts::PolyOpts, from="AUD", to="USD"; amount=100, precision=2)
    params = Dict(
        "apiKey" => opts.api_key,
        "amount" => amount,
        "precision" => precision
        )
    real_time_currency_conversion_url = "$real_time_currency_conversion_base_url/$from/$to"
    return generate_output_from_url(YesSinkNoResults(), real_time_currency_conversion_url, params, opts.sink)
end

############ Last Quote Currency Pair  ####################
"""
    last_quote_currency_pair(opts::PolyOpts, from="AUD", to="USD")

Get the last quote tick for a forex currency pair.

# Arguments
 * opts::PolyOpts: The PolyOpts object used to configure the request.
 * from: The "from" symbol of the pair.
 * to: The "to" symbol of the pair.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> last_quote_currency_pair(opts, "AUD", "USD")
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array
 * See https://polygon.io/docs/get_v1_last_quote_currencies__from___to__anchor for documentation on response attributes and supported keyword arguments.
"""
function last_quote_currency_pair(opts::PolyOpts, from="AUD", to="USD")
    params = Dict(
        "apiKey" => opts.api_key
        )
    last_quote_currency_pair_url = "$last_quote_currency_pair_base_url/$from/$to"
    return generate_output_from_url(YesSinkNoResults(), last_quote_currency_pair_url, params, opts.sink)
end

############ Forex Grouped Daily Bars  ####################
"""
    forex_grouped_daily_bars(opts::PolyOpts, date="2020-10-14"; adjusted=true)

Get the daily open, high, low, and close (OHLC) for the entire forex markets.

# Arguments
 * opts::PolyOpts: The PolyOpts object used to configure the request.
 * date: The beginning date for the aggregate window.
 * adjusted: Whether or not the results are adjusted for splits.
    By default, results are adjusted. Set this to false to get results that are NOT adjusted for splits.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> forex_grouped_daily_bars(opts, "2020-10-14")
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array
 * See https://polygon.io/docs/get_v2_aggs_grouped_locale_global_market_fx__date__anchor for documentation on response attributes and supported keyword arguments.
"""
function forex_grouped_daily_bars(opts::PolyOpts, date="2020-10-14"; adjusted=true)
    params = Dict(
        "apiKey" => opts.api_key,
        "adjusted" => adjusted
        )
    forex_grouped_daily_bars_url = "$forex_grouped_daily_bars_base_url/$date"
    return generate_output_from_url(YesSinkYesResults(), forex_grouped_daily_bars_url, params, opts.sink)
end

############ Forex Previous Close  ####################
"""
    forex_previous_close(opts::PolyOpts, forexTicker="C:EURUSD"; adjusted=true)

Get the previous day's open, high, low, and close (OHLC) for the specified forex pair.

# Arguments
 * opts::PolyOpts: The PolyOpts object used to configure the request.
 * forexTicker: The ticker symbol of the currency pair.
 * adjusted: Whether or not the results are adjusted for splits.
    By default, results are adjusted. Set this to false to get results that are NOT adjusted for splits.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> forex_previous_close(opts, "C:EURUSD")
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array
 * See https://polygon.io/docs/get_v2_aggs_ticker__forexTicker__prev_anchor for documentation on response attributes and supported keyword arguments.
"""
function forex_previous_close(opts::PolyOpts, forexTicker="C:EURUSD"; adjusted=true)
    params = Dict(
        "apiKey" => opts.api_key,
        "adjusted" => adjusted
        )
    forex_previous_close_url = "$forex_previous_close_base_url/$forexTicker/prev"
    return generate_output_from_url(YesSinkYesResults(), forex_previous_close_url, params, opts.sink)
end

############ Forex Aggregates Bars  ####################
"""
    forex_aggregates_bars(opts::PolyOpts, forexTicker="C:EURUSD", multiplier=1, timespan="day", from="2020-10-14", to="2020-10-14";
                        adjusted=true, sort="asc", limit=120)

Get aggregate bars for a forex pair over a given date range in custom time window sizes.
For example, if timespan = ‘minute’ and multiplier = ‘5’ then 5-minute bars will be returned.

# Arguments
 * opts::PolyOpts: The PolyOpts object used to configure the request.
 * forexTicker: The ticker symbol of the currency pair.
 * multiplier: The size of the timespan multiplier.
 * timespan: The size of the time window.
 * from: The start of the aggregate time window.
 * to: The end of the aggregate time window.
 * adjusted: Whether or not the results are adjusted for splits.
    By default, results are adjusted. Set this to false to get results that are NOT adjusted for splits.
 * sort: Sort the results by timestamp. asc will return results in ascending order (oldest at the top),
    desc will return results in descending order (newest at the top).
 * limit: Limits the number of base aggregates queried to create the aggregate results. Max 50000 and Default 120.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> forex_aggregates_bars(opts, "C:EURUSD", "5", "minute", "2020-10-14", "2020-10-14")
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array
 * See https://polygon.io/docs/get_v2_aggs_ticker__forexTicker__range__multiplier___timespan___from___to__anchor for documentation on response attributes and supported keyword arguments.
"""
function forex_aggregates_bars(opts::PolyOpts, forexTicker="C:EURUSD", multiplier=1, timespan="day", from="2020-10-14", to="2020-10-14";
                               adjusted=true, sort="asc", limit=120)
    params = Dict(
        "apiKey" => opts.api_key,
        "adjusted" => adjusted,
        "sort" => sort,
        "limit" => limit
        )

    forex_aggregates_bars_url = "$forex_aggregates_bars_base_url/$forexTicker/range/$multiplier/$timespan/$from/$to"
    return generate_output_from_url(YesSinkYesResults(), forex_aggregates_bars_url, params, opts.sink)
end

############ Forex Snapshot Ticker  ####################
"""
    forex_snapshot_ticker(opts::PolyOpts, forexTicker="C:EURUSD")

Get the current minute, day, and previous day’s aggregate, as well as the last trade and quote for a single traded currency symbol.
Note: Snapshot data is cleared at 12am EST and gets populated as data is received from the exchanges.

# Arguments
 * opts::PolyOpts: The PolyOpts object used to configure the request.
 * forexTicker: The forex ticker.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> forex_snapshot_ticker(opts, "C:EURUSD")
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array
 * See https://polygon.io/docs/get_v2_snapshot_locale_global_markets_forex_tickers__ticker__anchor for documentation on response attributes and supported keyword arguments.
"""
function forex_snapshot_ticker(opts::PolyOpts, forexTicker="C:EURUSD")
    params = Dict(
        "apiKey" => opts.api_key
        )
    forex_snapshot_ticker_url = "$forex_snapshot_ticker_base_url/$forexTicker"
    return generate_output_from_url(YesSinkYesTicker(), forex_snapshot_ticker_url, params, opts.sink)
end

############ Forex Snapshot All Tickers  ####################
"""
    forex_snapshot_all_tickers(opts::PolyOpts; kwargs...)

Get the current minute, day, and previous day’s aggregate, as well as the last trade and quote for all traded forex symbols.
Note: Snapshot data is cleared at 12am EST and gets populated as data is received from the exchanges. This can happen as early as 4am EST.

# Arguments
 * opts::PolyOpts: The PolyOpts object used to configure the request.
 * kwargs: The keyword arguments to pass to the PolyOpts.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> forex_snapshot_all_tickers(opts)
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array
 * See https://polygon.io/docs/get_v2_snapshot_locale_global_markets_forex_tickers_anchor for documentation on response attributes and supported keyword arguments.
"""
function forex_snapshot_all_tickers(opts::PolyOpts; kwargs...)
    params = Dict(
        "apiKey" => opts.api_key
        )
    merge!(params, Dict(kwargs))
    forex_snapshot_all_tickers_url = "$forex_snapshot_all_tickers_base_url"
    return generate_output_from_url(YesSinkYesTickers(), forex_snapshot_all_tickers_url, params, opts.sink)
end

############ Forex Snapshot Gainers/Losers  ####################
"""
    forex_snapshot_gainers_losers(opts::PolyOpts, direction="gainers")

Get the current top 20 gainers or losers of the day in forex markets.
Top gainers are those tickers whose price has increased by the highest percentage since the previous day's close.
Top losers are those tickers whose price has decreased by the highest percentage since the previous day's close.
Note: Snapshot data is cleared at 12am EST and gets populated as data is received from the exchanges.

# Arguments
 * opts::PolyOpts: The PolyOpts object used to configure the request.
 * direction: The direction of the snapshot results to return. "gainers" or "losers".

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> forex_snapshot_gainers_losers(opts, "losers")
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array
 * See https://polygon.io/docs/get_v2_snapshot_locale_global_markets_forex__direction__anchor for documentation on response attributes and supported keyword arguments.
"""
function forex_snapshot_gainers_losers(opts::PolyOpts, direction="gainers")
    params = Dict(
        "apiKey" => opts.api_key
        )
    forex_snapshot_gainers_losers_url = "$forex_snapshot_gainers_losers_base_url/$direction"
    return generate_output_from_url(YesSinkYesTickers(), forex_snapshot_gainers_losers_url, params, opts.sink)
end


text:
# This file holds the definition of the functions pertaining to the HistoryEntry
#   type.

# The functions of the HistoryEntry type require no additional types.
requiredTypes = [ "historyEntry" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

function Base.:(==)( entryLeft::HistoryEntry, entryRight::HistoryEntry )
    return entryLeft.timestamp == entryRight.timestamp
end  # ==( entryLeft, entryRight )

function Base.:(==)( entry::HistoryEntry, num::T ) where T <: Real
    return entry.timestamp == num
end  # ==( entry, num )

function Base.:(==)( num::T, entry::HistoryEntry ) where T <: Real
    return num == entry.timestamp
end  # ==( num, entry )

# XXX The "<" operator must be defined by the Base.isless function. Otherwise
#   Julia generates an error.
function Base.isless( entryLeft::HistoryEntry, entryRight::HistoryEntry )
    return entryLeft.timestamp < entryRight.timestamp
end  # isless( entryLeft, entryRight )

function Base.isless( entry::HistoryEntry, num::T ) where T <: Real
    return entry.timestamp < num
end  # isless( entry, num )

function Base.isless( num::T, entry::HistoryEntry ) where T <: Real
    return num < entry.timestamp
end  # isless( num, entry )

function Base.show( io::IO, entry::HistoryEntry )
    print( io, "$(entry.timestamp): $(entry.newVal)")
end  # show( io, entry )


text:
## auxiliary

const lowest = v"0"
const highest = v"999"

# some version comparisons need to ignore part of the version number
strip_patch(ver) = VersionNumber(ver.major, ver.minor)
strip_minor(ver) = VersionNumber(ver.major)

## GCC compilers supported by the CUDA toolkit

# Source: CUDA/include/host_config.h or include/crt/host_config.h on 9.0+
const cuda_gcc_db = Dict(
    v"5.5"  => lowest:v"4.8",  # (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ > 8)) && #error
    v"6.0"  => lowest:v"4.8",  # (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ > 8)) && #error
    v"6.5"  => lowest:v"4.8",  # (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ > 8)) && #error
    v"7.0"  => lowest:v"4.9",  # (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ > 9)) && #error
    v"7.5"  => lowest:v"4.9",  # (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ > 9)) && #error
    v"8.0"  => lowest:v"6-",   # (__GNUC__ > 5)                                          && #error
    v"9.0"  => lowest:v"7-",   # (__GNUC__ > 6)                                          && #error
    v"9.1"  => lowest:v"7-",   # (__GNUC__ > 6)                                          && #error
    v"9.2"  => lowest:v"8-",   # (__GNUC__ > 7)                                          && #error
    v"10.0" => lowest:v"8-",   # (__GNUC__ > 7)                                          && #error
    v"10.1" => lowest:v"9-",   # (__GNUC__ > 8)                                          && #error
)

function gcc_supported(gcc::VersionNumber, toolkit::VersionNumber)
    gcc_range = get(cuda_gcc_db, strip_patch(toolkit)) do
        error("no support for CUDA $toolkit")
    end

    return in(strip_patch(gcc), gcc_range)
end

## MSVC compilers supported by the CUDA toolkit

# Source: CUDA/include/host_config.h or include/crt/host_config.h on 9.0+
#
# NOTE: plain integers, because from MSVC 14.1 on the number increases monotonically
#       https://blogs.msdn.microsoft.com/vcblog/2016/10/05/visual-c-compiler-version/
const cuda_msvc_db = Dict(
    v"5.5"     => 1500:1700,  # (_MSC_VER < 1500 || _MSC_VER > 1700) && #error
    v"6.0"     => 1500:1700,  # (_MSC_VER < 1500 || _MSC_VER > 1700) && #error
    v"6.5"     => 1600:1800,  # (_MSC_VER < 1600 || _MSC_VER > 1800) && #error
    v"7.0"     => 1600:1800,  # (_MSC_VER < 1600 || _MSC_VER > 1800) && #error
    v"7.5"     => 1600:1800,  # (_MSC_VER < 1600 || _MSC_VER > 1800) && #error
    v"8.0"     => 1600:1900,  # (_MSC_VER < 1600 || _MSC_VER > 1900) && #error
    v"9.0"     => 1600:1910,  # (_MSC_VER < 1600 || _MSC_VER > 1910) && #error
    v"9.0.176" => 1600:1911,  # (_MSC_VER < 1600 || _MSC_VER > 1911) && #error
    v"9.1"     => 1600:1911,  # (_MSC_VER < 1600 || _MSC_VER > 1911) && #error
    v"9.2"     => 1600:1913,  # (_MSC_VER < 1600 || _MSC_VER > 1913) && #error
    v"10.0"    => 1700:1920,  # (_MSC_VER < 1700 || _MSC_VER > 1920) && #error
    v"10.1"    => 1700:1930,  # (_MSC_VER < 1700 || _MSC_VER > 1930) && #error
)

function msvc_supported(msvc::VersionNumber, toolkit::VersionNumber)
    msvc_range = get(cuda_msvc_db, toolkit) do
        get(cuda_msvc_db, strip_patch(toolkit)) do
            error("no support for CUDA $toolkit")
        end
    end

    msvc_num = msvc.major * 100 + msvc.minor
    return in(msvc_num, msvc_range)
end

## devices supported by the CUDA toolkit

# Source:
# - https://en.wikipedia.org/wiki/CUDA#GPUs_supported
# - ptxas |& grep -A 10 '\--gpu-name'
const dev_cuda_db = Dict(
    v"1.0" => lowest:v"6.5",
    v"1.1" => lowest:v"6.5",
    v"1.2" => lowest:v"6.5",
    v"1.3" => lowest:v"6.5",
    v"2.0" => lowest:v"8.0",
    v"2.1" => lowest:v"8.0",
    v"3.0" => v"4.2":highest,
    v"3.2" => v"6.0":highest,
    v"3.5" => v"5.0":highest,
    v"3.7" => v"6.5":highest,
    v"5.0" => v"6.0":highest,
    v"5.2" => v"7.0":highest,
    v"5.3" => v"7.5":highest,
    v"6.0" => v"8.0":highest,
    v"6.1" => v"8.0":highest,
    v"6.2" => v"8.0":highest,
    v"7.0" => v"9.0":highest,
    v"7.2" => v"9.2":highest,
    v"7.5" => v"10.0":highest,
)

function devices_for_cuda(ver::VersionNumber)
    caps = Set{VersionNumber}()
    for (cap,r) in dev_cuda_db
        if strip_patch(ver) in r
            push!(caps, cap)
        end
    end
    return caps
end

## PTX ISAs supported by the CUDA toolkit

# Source:
# - PTX ISA document, Release History table
const isa_cuda_db = Dict(
    v"1.0" => v"1.0":highest,
    v"1.1" => v"1.1":highest,
    v"1.2" => v"2.0":highest,
    v"1.3" => v"2.1":highest,
    v"1.4" => v"2.2":highest,
    v"1.5" => v"2.2":highest,
    v"2.0" => v"3.0":highest,
    v"2.1" => v"3.1":highest,
    v"2.2" => v"3.2":highest,
    v"2.3" => v"4.2":highest,
    v"3.0" => v"4.1":highest,
    v"3.1" => v"5.0":highest,
    v"3.2" => v"5.5":highest,
    v"4.0" => v"6.0":highest,
    v"4.1" => v"6.5":highest,
    v"4.2" => v"7.0":highest,
    v"4.3" => v"7.5":highest,
    v"5.0" => v"8.0":highest,
    v"6.0" => v"9.0":highest,
    v"6.1" => v"9.1":highest,
    v"6.2" => v"9.2":highest,
    v"6.3" => v"10.0":highest,
    v"6.4" => v"10.1":highest,
)

function isas_for_cuda(ver::VersionNumber)
    caps = Set{VersionNumber}()
    for (cap,r) in isa_cuda_db
        if strip_patch(ver) in r
            push!(caps, cap)
        end
    end
    return caps
end

## devices supported by the LLVM NVPTX back-end

# Source: LLVM/lib/Target/NVPTX/NVPTX.td
const dev_llvm_db = Dict(
    v"2.0" => v"3.2":highest,
    v"2.1" => v"3.2":highest,
    v"3.0" => v"3.2":highest,
    v"3.2" => v"3.7":highest,
    v"3.5" => v"3.2":highest,
    v"3.7" => v"3.7":highest,
    v"5.0" => v"3.5":highest,
    v"5.2" => v"3.7":highest,
    v"5.3" => v"3.7":highest,
    v"6.0" => v"3.9":highest,
    v"6.1" => v"3.9":highest,
    v"6.2" => v"3.9":highest,
    v"7.0" => v"6.0":highest,
    v"7.2" => v"7.0":highest,
    v"7.5" => v"8.0":highest,
)

function devices_for_llvm(ver::VersionNumber)
    caps = Set{VersionNumber}()
    for (cap,r) in dev_llvm_db
        if strip_patch(ver) in r
            push!(caps, cap)
        end
    end
    return caps
end

## PTX ISAs supported by the LVM NVPTX back-end

# Source: LLVM/lib/Target/NVPTX/NVPTX.td
const isa_llvm_db = Dict(
    v"3.0" => v"3.2":v"3.5",
    v"3.1" => v"3.2":v"3.5",
    v"3.2" => v"3.5":highest,
    v"4.0" => v"3.5":highest,
    v"4.1" => v"3.7":highest,
    v"4.2" => v"3.7":highest,
    v"4.3" => v"3.9":highest,
    v"5.0" => v"3.9":highest,
    v"6.0" => v"6.0":highest,
    v"6.1" => v"7.0":highest,
    v"6.3" => v"8.0":highest,
)

function isas_for_llvm(ver::VersionNumber)
    caps = Set{VersionNumber}()
    for (cap,r) in isa_llvm_db
        if strip_patch(ver) in r
            push!(caps, cap)
        end
    end
    return caps
end

## other

shader(cap::VersionNumber) = "sm_$(cap.major)$(cap.minor)"


text:
"""
    regpoly(X,Y,m)

adalah fungsi untuk mencari persamaan polinomial `y=C[m]x^m + ... + C[2]x^2 + C[1]x`
dari variabel dependen `Y` dan independen `X` dengan pangkat tertinggi `m`.

# Example
```jldoctest
julia> x = [-3, 0, 2, 4];

julia> y = [ 3, 1, 1, 3];

julia> C = regpoly(x,y,2)
3-element Array{Float64,1}:
  0.17846247712019525
 -0.19249542403904824
  0.8505186089078707
```
return nilai koefisien `C = [C[m],...,C[2],C[1]]`.
"""
function regpoly(X,Y,m)
  F = zeros(length(X),m+1)
  for k = 1: m+1
    F[:,k]=X.^(k-1);
  end
  A = F'*F;
  B = F'*Y;
  C = A\B;
  C = reverse(C,dims=1)
end


text:
@time __precompile__()
module MySQLforQpcrAnalysis

import DataStructures.OrderedDict

## for MySQL versions outputing named tuples by default instead of data frames
## supposedly, "To get the results as a DataFrame, you can just do MySQL.query(conn, sql, DataFrame)",
## but there was some challenges that i don't remember exactly
# using Compat, NamedTuples
# nt1 = @NT(a=[1,2],b=[3,4])
# DataFrame(OrderedDict(Compat.pairs(nt1)))

using DataFrames, DataStructures, JSON, MySQL

# constants

const STEP_NAME_to_CALIB_DATA_KEY = OrderedDict(
    "calib_single" => OrderedDict(
        "Water" => "water",
        "Signal" => "channel_1"
    ),
    "calib_dual" => OrderedDict(
        "Water" => "water",
        "FAM" => "channel_1",
        "HEX" => "channel_2"
    ),
    "ot_single" => OrderedDict(
        "LED off" => "baseline",
        "LED on" => "excitation"
    ),
    "ot_dual" => OrderedDict(
        "Baseline" => "baseline",
        "Water" => "water",
        "FAM" => "channel_1",
        "HEX" => "channel_2",
    )
)

# functions

# top-level
function id2json(
    db_conn ::MySQL.MySQLHandle,
    action ::String; # calib_like, amplification, meltcurve, analyze
    guid ::String="",
    # 0 is not valid
    exp_id ::Integer=0,
    calib_id ::Integer=0,
    stage_id ::Integer=0,
    step_id ::Integer=0,
    ramp_id ::Integer=0,
    well_nums ::Vector=[],
    step_name_to_calib_data_key ::OrderedDict{String,OrderedDict{String,String}}=STEP_NAME_to_CALIB_DATA_KEY
    )

    # sanity check and resolving id if necessary
    if action == "calib_like" && guid == "" # for a calibration experiment not part of an `analyze` call, if `calib_id == 0` use `exp_id` to infer `calib_id`, at the end use `calib_id` as `exp_id` for further processing.
        if calib_id == 0
            if exp_id == 0
                error("One of exp_id and calib_id must be valid (not 0) when retrieving calibration data.")
            else
                exp2calib_qry = "SELECT id, calibration_id FROM experiments WHERE id = $exp_id"
                calib_id = MySQL.mysql_execute(db_conn, exp2calib_qry)[1][:calibration_id][1]
                println("calib_id not specified, inferred from exp_id as $calib_id.")
            end
        end
        exp_id = calib_id
    elseif exp_id == 0
        error("exp_id must be valid except when retrieving calibration information not used for \"analyze\".")
    end

    well_constraint = length(well_nums) == 0 ? "" : "AND well_num in ($(join(well_nums, ",")))"

    data = nothing

    # find steps

    step_ramp_qry = "SELECT DISTINCT
            steps.id AS step_id,
            steps.name AS step_name,
            steps.collect_data AS step_collect_data,
            ramps.id AS ramp_id,
            ramps.collect_data AS ramp_collect_data
        FROM experiments
            LEFT JOIN protocols ON experiments.experiment_definition_id = protocols.experiment_definition_id
            LEFT JOIN stages ON protocols.id = stages.protocol_id
            LEFT JOIN steps ON stages.id = steps.stage_id
            LEFT JOIN ramps ON steps.id = ramps.next_step_id
        WHERE
                experiments.id = $exp_id
            AND stages.stage_type <> \'meltcurve\'
        ORDER BY step_id
    "
    step_ramp_qry_res = MySQL.mysql_execute(db_conn, step_ramp_qry)[1]
    step_count = length(step_ramp_qry_res[:step_id])

    # collect all the steps, not just useful ones (e.g. Water, Signal, FAM, HEX)
    exp_steps = OrderedDict(
        begin
            step_name_tmp = step_ramp_qry_res[:step_name][i]
            step_id_tmp = step_ramp_qry_res[:step_id][i]
            step_name_adj = ismissing(step_name_tmp) ? "missing_$step_id_tmp" : step_name_tmp
            step_name_adj => step_id_tmp
        end
        for i in 1:step_count
    )
    println("exp_steps: $exp_steps")

    # process according to action

    if action == "analyze"
        if guid in ("optical_cal", "optical_test_single_channel", "optical_test_dual_channel")
            action = "calib_like"
        elseif guid == "thermal_consistency"
            action = "meltcurve"
            stage_id = 4
        elseif guid == "thermal_performance_diagnostic"
            data_qry = "
            SELECT
                    heat_block_zone_1_temp,
                    heat_block_zone_2_temp,
                    elapsed_time,
                    lid_temp
                FROM temperature_logs
                WHERE experiment_id = $exp_id
                ORDER BY elapsed_time
            "
        end # guid
    end # action == "analyze"

    # separate from `if action == "analyze"` because need to consider `action` updated during `if action == "analyze"`

    if action == "amplification"

        if step_id == 0 && ramp_id == 0
            println("step(s) where data is collected on the step or the ramp leading up to the step:")
            for i in 1:step_count
                step_collect_data = step_ramp_qry_res[:step_collect_data][i]
                ramp_collect_data = step_ramp_qry_res[:ramp_collect_data][i]
                if step_collect_data == 1 || ramp_collect_data == 1
                    step_id = step_ramp_qry_res[:step_id][i]
                    ramp_id = step_ramp_qry_res[:ramp_id][i]
                    step_name = step_ramp_qry_res[:step_name][i]
                    cyc_qry = "
                        SELECT DISTINCT cycle_num
                            FROM fluorescence_data
                            WHERE
                                    experiment_id = $exp_id
                                AND (step_id = $step_id
                                OR ramp_id = $ramp_id)
                                $well_constraint
                    "
                    cyc_res = MySQL.mysql_execute(db_conn, cyc_qry)[1]
                    num_cycles = length(cyc_res[:cycle_num])
                    #
                    info_clause = "info on the $(i)th step of this exp"
                    println("
                        $info_clause, start {
                        number of cycles - $num_cycles,
                        step_id - $step_id,
                        step_collect_data - $step_collect_data,
                        ramp_id - $ramp_id,
                        ramp_collect_data - $ramp_collect_data,
                        step_name - $step_name
                        } end, $info_clause
                    ")
                end # if step_collect_data
            end # for i
            error("step_id and ramp_id cannot be 0 at the same time for amplification. please consider using the information above to specify a step_id or ramp_id")

        elseif step_id != 0 && ramp_id != 0
            error("neither step_id or ramp_id is 0, it is unclear which one should be used. please specify only one as non-0 and the other as 0.")

        elseif step_id != 0 && ramp_id == 0
            step_ramp_constraint = "step_id = $step_id"

        elseif step_id == 0 && ramp_id != 0
            step_ramp_constraint = "ramp_id = $ramp_id"

        end # if step_id

        data_qry = "
        SELECT fluorescence_value, well_num, cycle_num, channel
            FROM fluorescence_data
            WHERE
                    experiment_id = $exp_id
                AND $step_ramp_constraint
                $well_constraint
            ORDER BY channel, well_num, cycle_num
        "

    elseif action == "meltcurve"
        if stage_id == 0
            stage_qry = "SELECT DISTINCT stages.id AS stage_id
                FROM experiments
                    LEFT JOIN protocols ON experiments.experiment_definition_id = protocols.experiment_definition_id
                    LEFT JOIN stages ON protocols.id = stages.protocol_id
                WHERE
                        experiments.id = $exp_id
                        AND stages.stage_type = \'meltcurve\'
            "
            println(MySQL.mysql_execute(db_conn, stage_qry)[1])
            error("stage_id must be valid for meltcurve. please consider the information above and specify a stage_id.")
        end
        data_qry = "
        SELECT well_num, temperature, fluorescence_value, channel
            FROM melt_curve_data
            WHERE
                    experiment_id = $exp_id
                AND stage_id = $stage_id
                $well_constraint
            ORDER BY channel, well_num, temperature
        "

    elseif action == "calib_like" # compute data
        # get data
        if guid == "optical_test_single_channel"
            step_config = "ot_single"
        elseif guid == "optical_test_dual_channel"
            step_config = "ot_dual"
        elseif "Signal" in keys(exp_steps)
            step_config = "calib_single"
        else
            step_config = "calib_dual"
        end # if guid
        data = OrderedDict(
            pair.second => begin
                step_id = exp_steps[pair.first]
                data_qry = "SELECT fluorescence_value, channel
                    FROM fluorescence_data
                    WHERE
                            experiment_id = $exp_id
                        AND step_id = $step_id
                        $well_constraint
                        AND cycle_num = 1
                	ORDER BY channel, well_num
                "
                step_data_raw = MySQL.mysql_execute(db_conn, data_qry)[1]
                if step_config in ("calib_single", "ot_single")
                    step_data = [step_data_raw[:fluorescence_value], nothing]
                elseif step_config in ("calib_dual", "ot_dual")
                    step_data = map(unique(step_data_raw[:channel])) do channel_num
                        step_data_raw[:fluorescence_value][step_data_raw[:channel] .== channel_num]
                    end # do channel_num
                end # if
                step_data # to return
            end # begin
            for pair in step_name_to_calib_data_key[step_config]
        )

    end # if action

    if data == nothing # data was not computed during `if action`. currently equivalent to `action != "calib_like"`
        data_qry_res = MySQL.mysql_execute(db_conn, data_qry)[1]
        data = OrderedDict(
            string(coln) => map(data_qry_res[coln]) do val
                ismissing(val) ? nothing : val
            end # do val
            for coln in data_qry_res.colindex.names # assuming qry_res is a DataFrame. need to change if `qry_res` is a NamedTuple
        )
    end # if data

    # return data # for testing
    return JSON.json(data)

end # id2json

end # module MySQLforQpcrAnalysis

#

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
@system BoundaryLayer(WeatherStub) begin
    w: leaf_width => 0.1 ~ preserve(u"m", parameter)

    sr: stomatal_ratio => begin
        # maize is an amphistomatous species, assume 1:1 (adaxial:abaxial) ratio.
        #sr = 1.0
        # switchgrass adaxial : abaxial (Awada 2002)
        # https://doi.org/10.4141/P01-031
        #sr = 1.28
        1.0
    end ~ preserve(parameter)
    scr(sr): sides_conductance_ratio => ((sr + 1)^2 / (sr^2 + 1)) ~ preserve

    # multiply by 1.4 for outdoor condition, Campbell and Norman (1998), p109
    ocr: outdoor_conductance_ratio => 1.4 ~ preserve

    u(u=wind): wind_velocity ~ track(u"m/s", min=0.1)
    # characteristic dimension of a leaf, leaf width in m
    d(w): characteristic_dimension => 0.72w ~ track(u"m")
    v: kinematic_viscosity_of_air_at_20 => 1.51e-5 ~ preserve(u"m^2/s", parameter)
    κ: thermal_diffusivity_of_air_at_20 => 21.5e-6 ~ preserve(u"m^2/s", parameter)
    Re(u, d, v): reynolds_number => u*d/v ~ track
    Nu(Re): nusselt_number => 0.60sqrt(Re) ~ track
    gh(κ, Nu, d, scr, ocr, P_air, Tk_air): boundary_layer_heat_conductance => begin
        g = κ * Nu / d
        # multiply by ratio to get the effective blc (per projected area basis), licor 6400 manual p 1-9
        g *= scr * ocr
        # including a multiplier for conversion from mm s-1 to mol m-2 s-1
        g * P_air / (u"R" * Tk_air)
    end ~ track(u"mmol/m^2/s")
    # 1.1 is the factor to convert from heat conductance to water vapor conductance, an average between still air and laminar flow (see Table 3.2, HG Jones 2014)
    gb(gh, P_air): boundary_layer_conductance => 0.147/0.135*gh / P_air ~ track(u"mol/m^2/s/bar" #= H2O =#)
end


text:
module DuckDB

using DBInterface
using WeakRefStrings
using Tables
using Base.Libc
using Dates
using DataFrames
using Tables
using UUIDs

export DBInterface, DuckDBException

include("helper.jl")
include("exceptions.jl")
include("ctypes.jl")
include("api.jl")
include("logical_type.jl")
include("value.jl")
include("validity_mask.jl")
include("vector.jl")
include("data_chunk.jl")
include("config.jl")
include("database.jl")
include("statement.jl")
include("result.jl")
include("transaction.jl")
include("ddl.jl")
include("appender.jl")
include("table_function.jl")
include("replacement_scan.jl")
include("data_frame_scan.jl")
include("old_interface.jl")

end # module


text:
mutable struct VIPData
    F::Array{JuMP.NonlinearExpression,1}
    var::Array{JuMP.Variable,1}
    relation::Dict{JuMP.Variable, JuMP.NonlinearExpression}
end

function VIPModel(solver=Ipopt.IpoptSolver(print_level=0))
    m = Model(solver=solver)
    m.ext[:VIP] = VIPData(Array{JuMP.NonlinearExpression}[], Array{JuMP.Variable}[], Dict{JuMP.Variable, JuMP.NonlinearExpression}() )
    return m
end

function getVIPData(m::Model)
    if haskey(m.ext, :VIP)
        return m.ext[:VIP]::VIPData
    else
        error("The 'getVIPData' function is only for VIP models as in VariationalInequality.jl")
    end
end

@eval const $(Symbol("@mapping")) = $(Symbol("@NLexpression"))

macro innerproduct(args...)
  if length(args) != 3
    error("3 arguments are required in @innerproduct(...)")
  end

  m = esc(args[1])
  F_name = esc(args[2])
  var = esc(args[3])

  code = quote
    innerproduct($m, $F_name, $var)
  end

  return code
end

################################################################################
# correspond interface
# The most basic/important one. All other interfaces call this method.
function innerproduct(m::Model, expression::JuMP.NonlinearExpression, variable::JuMP.Variable)
    data = getVIPData(m)
    data.relation[variable] = expression
end
# Alternative
function innerproduct(m::Model, variable::JuMP.Variable, expression::JuMP.NonlinearExpression)
    innerproduct(m, expression, variable)
end

# Do we need the below?
# function correspond(m::Model, variables::Array{JuMP.Variable,1}, expressions::Array{JuMP.NonlinearExpression,1})
#     correspond(m, expressions, variables)
# end
# function correspond(m::Model, expressions::Array{JuMP.NonlinearExpression,1}, variables::Array{JuMP.Variable,1})
#     @assert length(expressions) == length(variables)
#     for i in 1:length(variables)
#         correspond(m, expressions[i], variables[i])
#     end
# end

function innerproduct(m::Model, expressions::Array{JuMP.NonlinearExpression}, variables::Array{JuMP.Variable})
    expressions = collect(expressions)
    variables = collect(variables)
    @assert length(expressions) == length(variables)
    for i in 1:length(variables)
        innerproduct(m, expressions[i], variables[i])
    end
end
function innerproduct(m::Model, variables::Array{JuMP.Variable}, expressions::Array{JuMP.NonlinearExpression})
    innerproduct(m, expressions, variables)
end

# Do we need the below?
# function correspond(m::Model, expressions::JuMP.JuMPArray, variables::JuMP.JuMPArray)
#     variables = collect(variables.innerArray)
#     expressions = collect(expressions.innerArray)
#     @assert length(expressions) == length(variables)
#     for i in 1:length(variables)
#         correspond(m, expressions[i], variables[i])
#     end
# end
################################################################################

# function setVIP(m::Model, expressions::Array{JuMP.NonlinearExpression,1}, variables::Array{JuMP.Variable,1})
#     @assert length(expressions) == length(variables)
#
#     data = getVIPData(m)
#     data.F = expressions
#     data.var = variables
#
#     return m
# end

# function solveVIP(m::Model    ; step_size=0.01,
#                                 algorithm=:fixed_point,
#                                 tolerance=1e-10,
#                                 max_iter=1000            )
#
#     F = getVIPData(m).F
#     var = getVIPData(m).var
#     @assert length(F) == length(var)
#     dim = length(F)
#
#     x0 = zeros(dim)
#
#     return solveVIP(m, x0, step_size=step_size, algorithm=algorithm, tolerance=tolerance, max_iter=max_iter)
#
# end

# Copied from JuMP.jl
# internal method that doesn't print a warning if the value is NaN
_get_value(v::JuMP.Variable) = v.m.colVal[v.col]

function clear_values(m)
    relation = getVIPData(m).relation
    for variable in keys(relation)
        setvalue(variable, NaN)
    end
end

function initial_projection(m::JuMP.Model)
    relation = getVIPData(m).relation

    initial_values = Dict{JuMP.Variable, Float64}()
    for variable in keys(relation)
        val = _get_value(variable)
        if isnan(val)
            val = 0.0
        end
        initial_values[variable] = val
    end

    @objective(m, Min, sum( ( variable - initial_values[variable] )^2 for variable in keys(relation) ) )
    status = solve(m)
    @assert status==:Optimal
end

function gap_function(m)
    relation = getVIPData(m).relation
    var = get_variables(relation)
    y = get_current_x(relation)
    Fy = get_current_F(relation)

    @objective(m, Max,
        sum( Fy[j] * ( y[j] - var[j] ) for j in 1:length(var))
    )
    solve(m)
    return getobjectivevalue(m)
end

function save_solution(m)
    relation = getVIPData(m).relation
    var = get_variables(relation)
    x = get_current_x(relation)
    F = get_current_F(relation)

    solution = Dict(zip(var, x))
    F_value = Dict(zip(var, F))
    gap = gap_function(m)

    return solution, F_value, gap
end

function solveVIP(m::Model;    step_size=0.01,
                                algorithm=:fixed_point,
                                tolerance=1e-6,
                                max_iter=1000            )

    initial_projection(m)

    if algorithm==:fixed_point
        _fixed_point(m, step_size, tolerance, max_iter)
    elseif algorithm==:extra_gradient
        _extra_gradient(m, step_size, tolerance, max_iter)
    # elseif algorithm==:hyperplane
    #     _hyperplane(m, step_size, tolerance, max_iter)
    end

    # sol, Fval, gap = save_solution(m)
    return save_solution(m)
end


text:
export rule

@rule Gamma(:out, Marginalisation) (m_α::PointMass, m_θ::PointMass) = Gamma(mean(m_α), mean(m_θ))

@rule Gamma(:out, Marginalisation) (q_α::PointMass, q_θ::PointMass) = Gamma(mean(q_α), mean(q_θ))


text:
using OceanTurb

using OceanTurb: @add_standard_model_fields, AbstractParameters, AbstractModel

import OceanTurb.KPP: ∂B∂z

const nsol = 4
@solution U V T S

Base.@kwdef struct ConvectiveAdjustmentParameters{T} <: AbstractParameters
    κ   :: T = 1
end

struct ConvectiveAdjustmentModel{TS, G, T} <: AbstractModel{TS, G, T}
    @add_standard_model_fields
    parameters :: ConvectiveAdjustmentParameters{T}
    constants  :: Constants{T}
end

function ConvectiveAdjustmentModel(; N, L,
          grid = UniformGrid(N, L),
     constants = Constants(),
    parameters = Parameters(),
       stepper = :ForwardEuler,
           bcs = BoundaryConditions((ZeroFluxBoundaryConditions() for i=1:nsol)...)
    )

    solution = Solution((CellField(grid) for i=1:nsol)...)
    K = (U=KU, V=KV, T=KT, S=KS)
    R = (U=RU, V=RV, T=RT, S=RS)
    eqn = Equation(K=K, R=R)
    lhs = OceanTurb.build_lhs(solution)

    timestepper = Timestepper(stepper, eqn, solution, lhs)

    return ConvectiveAdjustmentModel(Clock(), grid, timestepper, solution, bcs, parameters, constants)
end

function temperature_gradient(T, i)
    return ∂z(T, i)
end

temperature_gradient(m, i) = temperature_gradient(m.solution.T, i)

KU() = 0
KT(∂T∂z, κ) = κ * (∂T∂z < 0)

KU(m, i) = KU()
KT(m, i) = KT(temperature_gradient(m, i), m.parameters.κ)

const KV = KU
const KS = KT

RU(m, i) =   m.constants.f * m.solution.V[i]
RV(m, i) = - m.constants.f * m.solution.U[i]
RT(m, i) = 0
RS(m, i) = 0


text:
using Test, LinearMaps, LinearAlgebra

# new type
struct SimpleFunctionMap <: LinearMap{Float64}
    f::Function
    N::Int
end
struct SimpleComplexFunctionMap <: LinearMap{Complex{Float64}}
    f::Function
    N::Int
end
Base.size(A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}) = (A.N, A.N)
Base.:(*)(A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}, v::Vector) = A.f(v)
LinearAlgebra.mul!(y::Vector, A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}, x::Vector) = copyto!(y, *(A, x))

@testset "composition" begin
    F = @inferred LinearMap(cumsum, y -> reverse(cumsum(reverse(x))), 10; ismutating=false)
    FC = @inferred LinearMap{ComplexF64}(cumsum, y -> reverse(cumsum(reverse(x))), 10; ismutating=false)
    A = 2 * rand(ComplexF64, (10, 10)) .- 1
    B = rand(size(A)...)
    M = @inferred 1 * LinearMap(A)
    N = @inferred LinearMap(B)
    v = rand(ComplexF64, 10)
    @test @inferred (F * F) * v == @inferred F * (F * v)
    @test @inferred (F * A) * v == @inferred F * (A * v)
    @test @inferred (A * F) * v == @inferred A * (F * v)
    @test @inferred A * (F * F) * v == @inferred A * (F * (F * v))
    @test @inferred (F * F) * (F * F) * v == @inferred F * (F * (F * (F * v)))
    @test @inferred Matrix(M * transpose(M)) ≈ A * transpose(A)
    @test @inferred !isposdef(M * transpose(M))
    @test @inferred isposdef(M * M')
    @test @inferred issymmetric(N * N')
    @test @inferred ishermitian(N * N')
    @test @inferred !issymmetric(M' * M)
    @test @inferred ishermitian(M' * M)
    @test @inferred issymmetric(F'F)
    @test @inferred ishermitian(F'F)
    @test @inferred !issymmetric(FC'FC)
    @test @inferred ishermitian(FC'FC)
    @test @inferred isposdef(transpose(F) * F * 3)
    @test @inferred isposdef(transpose(F) * 3 * F)
    @test @inferred !isposdef(-5*transpose(F) * F)
    @test @inferred isposdef((M * F)' * M * 4 * F)
    @test @inferred transpose(M * F) == @inferred transpose(F) * transpose(M)
    @test @inferred (4*((-3*M)*2)) == @inferred -12M*2
    @test @inferred (4*((3*(-M))*2)*(-5)) == @inferred -12M*(-10)
    L = @inferred 3 * F + 1im * A + F * M' * F
    LF = 3 * Matrix(F) + 1im * A + Matrix(F) * Matrix(M)' * Matrix(F)
    @test Array(L) ≈ LF
    R1 = rand(ComplexF64, 10, 10); L1 = LinearMap(R1)
    R2 = rand(ComplexF64, 10, 10); L2 = LinearMap(R2)
    R3 = rand(ComplexF64, 10, 10); L3 = LinearMap(R3)
    CompositeR = prod(R -> LinearMap(R), [R1, R2, R3])
    @test @inferred L1 * L2 * L3 == CompositeR
    @test @inferred transpose(CompositeR) == transpose(L3) * transpose(L2) * transpose(L1)
    @test @inferred adjoint(CompositeR) == L3' * L2' * L1'
    @test @inferred adjoint(adjoint((CompositeR))) == CompositeR
    @test transpose(transpose((CompositeR))) == CompositeR
    Lt = @inferred transpose(LinearMap(CompositeR))
    @test Lt * v ≈ transpose(R3) * transpose(R2) * transpose(R1) * v
    Lc = @inferred adjoint(LinearMap(CompositeR))
    @test Lc * v ≈ R3' * R2' * R1' * v

    # test inplace operations
    w = similar(v)
    mul!(w, L, v)
    @test w ≈ LF * v

    # test new type
    F = SimpleFunctionMap(cumsum, 10)
    FC = SimpleComplexFunctionMap(cumsum, 10)
    @test @inferred ndims(F) == 2
    @test @inferred size(F, 1) == 10
    @test @inferred length(F) == 100
    @test @inferred !issymmetric(F)
    @test @inferred !ishermitian(F)
    @test @inferred !ishermitian(FC)
    @test @inferred !isposdef(F)
    w = similar(v)
    mul!(w, F, v)
    @test w == F * v
    @test_throws MethodError F' * v
    @test_throws MethodError transpose(F) * v
    @test_throws MethodError mul!(w, adjoint(F), v)
    @test_throws MethodError mul!(w, transpose(F), v)

    # test composition of several maps with shared data #31
    global sizes = ( (5, 2), (3, 3), (3, 2), (2, 2), (9, 2), (7, 1) )
    N = length(sizes) - 1
    global Lf = []
    global Lt = []
    global Lc = []

    # build list of operators [LN, ..., L2, L1] for each mode
    for (fi, i) in [ (Symbol("f$i"), i) for i in 1:N]
        @eval begin
            function ($fi)(source)
                dest = ones(prod(sizes[$i + 1]))
                tmp = reshape(source, sizes[$i])
                return conj.($i * dest)
            end
            insert!(Lf, 1, LinearMap($fi, prod(sizes[$i + 1]), prod(sizes[$i])))
            insert!(Lt, 1, LinearMap(x -> x, $fi, prod(sizes[$i]), prod(sizes[$i + 1])))
            insert!(Lc, 1, LinearMap{ComplexF64}(x -> x, $fi, prod(sizes[$i]), prod(sizes[$i + 1])))
        end
    end
    @test size(prod(Lf[1:N])) == (prod(sizes[end]), prod(sizes[1]))
    @test isreal(prod(Lf[1:N]))
    # multiply as composition and as recursion
    v1 = ones(prod(sizes[1]))
    u1 = ones(prod(sizes[1]))
    w1 = im.*ones(ComplexF64, prod(sizes[1]))
    for i = N:-1:1
        v2 = prod(Lf[i:N]) * ones(prod(sizes[1]))
        u2 = transpose(LinearMap(prod(Lt[N:-1:i]))) * ones(prod(sizes[1]))
        w2 = adjoint(LinearMap(prod(Lc[N:-1:i]))) * ones(prod(sizes[1]))

        v1 = Lf[i] * v1
        u1 = transpose(Lt[i]) * u1
        w1 = adjoint(Lc[i]) * w1

        @test v1 == v2
        @test u1 == u2
        @test w1 == w2
    end
end


text:
# SteinSquaredFourierSumDensity
# This model is used in Section 5.1 of http://icml.cc/2012/papers/683.pdf
#
# The density has support [0,1] and is given by
#
# p(x) \propto (\sum_{i=1}^d a_i cos(2i\pi x) + b_i sin(2i\pi x))^2
#
# If d=0 then we reduce this to the uniform density

type SteinSquaredFourierSumDensity <: SteinDistribution
    a::Array{Float64, 1}
    b::Array{Float64, 1}
    c1::Float64
    c2::Float64
    c3::Float64
end

# Just the empty d=0 case
SteinSquaredFourierSumDensity() = SteinSquaredFourierSumDensity([], [])

# Scalar a,b
SteinSquaredFourierSumDensity(a::Float64, b::Float64) =
SteinSquaredFourierSumDensity([a], [b])

# Default (c1,c2,c3) = (1,1,1)
SteinSquaredFourierSumDensity(a::Array{Float64, 1}, b::Array{Float64, 1}) =
SteinSquaredFourierSumDensity(a, b, 1.0, 1.0, 1.0)

function supportlowerbound(d::SteinSquaredFourierSumDensity, j::Int64)
    0.0
end

function supportupperbound(d::SteinSquaredFourierSumDensity, j::Int64)
    1.0
end

function gradlogdensity(d::SteinSquaredFourierSumDensity, X::Array{Float64,2})
    p = length(d.a)
    if p == 0
        return zeros(size(X))
    end
    cosix = [cos(2*pi*i*x) for i=1:p, x=X]
    sinix = [sin(2*pi*i*x) for i=1:p, x=X]
    ia_i = d.a .* (1:p)
    ib_i = d.b .* (1:p)

    numer = cosix .* ib_i - sinix .* ia_i
    denom = cosix .* d.a  + sinix .* d.b

    gradlog = 4*pi*(sum(numer, 1) ./ sum(denom, 1))
    reshape(gradlog, size(X)...)
end

function numdimensions(d::SteinSquaredFourierSumDensity)
    length(d.a)
end


text:
using Test, ParametricAirfoils
include("testnaca.jl")
include("testcst.jl")