text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["libgeos_c"], :libgeos),
    LibraryProduct(prefix, ["libgeos"], :libgeos_cpp),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaGeo/GEOSBuilder/releases/download/v3.7.2-0"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/GEOS.v3.7.2.aarch64-linux-gnu-gcc4.tar.gz", "5787a09eb5f1de0e8328b3fdb482e4cc798fd6b5a4689111f9a6befe15bdbb80"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/GEOS.v3.7.2.aarch64-linux-gnu-gcc7.tar.gz", "9b5ee2b297b10b296ea9c931a83bc3f299352e30a9ff696f0f8d402b3651f403"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/GEOS.v3.7.2.aarch64-linux-gnu-gcc8.tar.gz", "e48d44177d7b5a836860f04959878365a941ace156cfff06e718c37fd55b9c49"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/GEOS.v3.7.2.aarch64-linux-musl-gcc4.tar.gz", "90aec29158928a43a4e2bba4c3116a2d3fc5cc621c829146865d66a6c3cd888a"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/GEOS.v3.7.2.aarch64-linux-musl-gcc7.tar.gz", "7c57861224750ecd3d068a7d7d7e89b1171f4ba49cde4c6762c735b5520fc0cc"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/GEOS.v3.7.2.aarch64-linux-musl-gcc8.tar.gz", "0163dc3893c8ec3befae6227a564c0f26d4db15aab23783c17d70698aa78f29e"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/GEOS.v3.7.2.arm-linux-gnueabihf-gcc4.tar.gz", "99a6816b5503d173ef7c4b9ed1645abb1dabd91a8a7e051a7900e671ac7b62c6"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/GEOS.v3.7.2.arm-linux-gnueabihf-gcc7.tar.gz", "c0d893dd8ddc8fffed1979525eda0ad673f1a2db3cadbd1a05f14a20a8174ed6"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/GEOS.v3.7.2.arm-linux-gnueabihf-gcc8.tar.gz", "e752016544692fa174f157fe63069479888a871f88926a87345b8fb13c48400d"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/GEOS.v3.7.2.arm-linux-musleabihf-gcc4.tar.gz", "af0b3d62984c8bac1f0842124c6ac53a3476c963bdf1219c61215a326d60e7f1"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/GEOS.v3.7.2.arm-linux-musleabihf-gcc7.tar.gz", "6690373e7fc7963ea2c3edbe893001f0856dac27bfcf0550f087116ddaaf7668"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/GEOS.v3.7.2.arm-linux-musleabihf-gcc8.tar.gz", "fccda9b71fab947cefe9c900d383e15e7c4c08331b218a26fa5073b9b12f79ac"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/GEOS.v3.7.2.i686-linux-gnu-gcc4.tar.gz", "4e0b9390d289ca68be601e70a588949602025a459a95479dda35f63c3f29cc3d"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/GEOS.v3.7.2.i686-linux-gnu-gcc7.tar.gz", "82a71e47c60e6c25472519b1a09ef4aa6318a8bed9791f00d492507431e078f9"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/GEOS.v3.7.2.i686-linux-gnu-gcc8.tar.gz", "6d8a73c951dc7e30ebd2f9b8597eeb427f758a2e16854fe58fedcd8b3c2bb1e3"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/GEOS.v3.7.2.i686-linux-musl-gcc4.tar.gz", "aa03493edf4e4913da3d2b32b9237283f5fdc26e35b750deddfa2581f857e140"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/GEOS.v3.7.2.i686-linux-musl-gcc7.tar.gz", "d10d2f8ba7728d45291848c6d7377cf02110b65472580cc372230007f04a82e1"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/GEOS.v3.7.2.i686-linux-musl-gcc8.tar.gz", "ba1d6f8baabe5711f2e69c5076eeb6044595e6f6207524cd7542b23ca7ffa226"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/GEOS.v3.7.2.i686-w64-mingw32-gcc4.tar.gz", "dd98e9fc7fab8b8424f833731f06c40cacb43f1f3ba7036e1e3f79d8027a3453"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/GEOS.v3.7.2.i686-w64-mingw32-gcc7.tar.gz", "165770dbc66573d7a8e10030f7b8547107662cdd8cb7872b63db25d0e32de88c"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/GEOS.v3.7.2.i686-w64-mingw32-gcc8.tar.gz", "a01d8761434684c1f9123f4543467f0ad2d9920a9e659f80d6e5c9c6ddaffa79"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/GEOS.v3.7.2.powerpc64le-linux-gnu-gcc4.tar.gz", "d2121b5fde6163d8cd6b45fb00bd656196629b3dd1d7743634e66ed3b021b7a5"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/GEOS.v3.7.2.powerpc64le-linux-gnu-gcc7.tar.gz", "488df95ca419bbc527ef75a05a7d99a309c86304a42c0ed2313848d0ad2da843"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/GEOS.v3.7.2.powerpc64le-linux-gnu-gcc8.tar.gz", "7ad7f4efa44d0106f36b2cbe9be1fdd6a0ee7133700bc7bd3067d908649484ef"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/GEOS.v3.7.2.x86_64-apple-darwin14-gcc4.tar.gz", "46d2873ffef747e278af36c429bd14c46ee133aeb454e1d34d4a4915c321ab12"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/GEOS.v3.7.2.x86_64-apple-darwin14-gcc7.tar.gz", "e8aebc603654e84d4c3327b21273146b897bfd1694411316c31cd345ebf7ee4d"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/GEOS.v3.7.2.x86_64-apple-darwin14-gcc8.tar.gz", "8f154d5c5422a7d39e5a5da0dad2f35b7669e0b49683e7fe0604b75f4203bba3"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/GEOS.v3.7.2.x86_64-linux-gnu-gcc4.tar.gz", "be0fd5b0ca50092b27d64be31140871c68898ac635b59fffa3c192c78dab5034"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/GEOS.v3.7.2.x86_64-linux-gnu-gcc7.tar.gz", "d98b1f02d3d9591630207d75bf12242a80b98b28a2eb2dfe080c396fc8291040"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/GEOS.v3.7.2.x86_64-linux-gnu-gcc8.tar.gz", "19a4159ee0eb19dd71fc1ab2cd671083846cb6a33a355d9e850d1ebfad7be5f2"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/GEOS.v3.7.2.x86_64-linux-musl-gcc7.tar.gz", "1032d637fad87e201da9b8418f15dac2d6c4d905bc98ec041cebef718d52bc93"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/GEOS.v3.7.2.x86_64-linux-musl-gcc8.tar.gz", "bf98e8f95b95d531e6b0a4a6f2ead3d3c07192bdc9c4f6c5adc436b41f1d9177"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/GEOS.v3.7.2.x86_64-unknown-freebsd11.1-gcc4.tar.gz", "a1b35b13701604012b5f783d0187408d3dfe9825bb5a50815017da9556dc07c1"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/GEOS.v3.7.2.x86_64-unknown-freebsd11.1-gcc7.tar.gz", "efc48c1a34a38718dcc70ac860ce88ca0ae2ec71b5139d614f19eee602b4f827"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/GEOS.v3.7.2.x86_64-unknown-freebsd11.1-gcc8.tar.gz", "aaf02142353c250c531625f943f2a11cf52983c6f0a5ae24a2cd5251d8a13bed"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/GEOS.v3.7.2.x86_64-w64-mingw32-gcc7.tar.gz", "34ca72de29b79d9635fceed8764ec6428aa5c92f5ef4158fe70539b500db42c1"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/GEOS.v3.7.2.x86_64-w64-mingw32-gcc8.tar.gz", "e571095938420ec3422fc791040fc1970be5659cb54ab5b151bd81ef3d9c0077"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!")
end

# If we have a download, and we are unsatisfied (or the version we're
# trying to install is not itself installed) then load it up!
if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
    # Download and install binaries
    install(dl_info...; prefix=prefix, force=true, verbose=verbose)
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps_geos.jl"), products, verbose=verbose)


text:
# Copyright 2018 IBM Corp.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module OntologyJSON
export presentation_from_json, annotation_from_json

using Catlab
using ...Doctrine, ..Ontology

# Concepts
##########

""" Load Monocl concepts (as a presentation) from JSON documents.
"""
function presentation_from_json(docs)::Presentation
  # Be careful about load order because:
  # - To add a morphism, the domain/codomain objects must be added
  # - To add a subobject, the domain/codomain objects must already be added
  # - To add a submorphism, the domain/codomain morphisms must already be added
  presentation = Presentation{String}(Monocl)
  ob_docs = filter(doc -> doc["kind"] == "type", docs)
  hom_docs = filter(doc -> doc["kind"] == "function", docs)
  for doc in ob_docs
    add_ob_generator_from_json!(presentation, doc)
  end
  for doc in ob_docs
    add_subob_generators_from_json!(presentation, doc)
  end
  for doc in hom_docs
    add_hom_generator_from_json!(presentation, doc)
  end
  for doc in hom_docs
    add_subhom_generators_from_json!(presentation, doc)
  end
  presentation
end

""" Add object from JSON document to presentation.
"""
function add_ob_generator_from_json!(pres::Presentation, doc::AbstractDict)
  ob = Ob(Monocl, doc["id"])
  add_generator!(pres, ob)
end

""" Add subobjects from JSON document to presentation.
"""
function add_subob_generators_from_json!(pres::Presentation, doc::AbstractDict)
  ob = generator(pres, doc["id"])::Monocl.Ob
  names = get(doc, "is-a", [])
  names = isa(names, AbstractString) ? [ names ] : names
  for super_name in names
    super_ob = generator(pres, super_name)::Monocl.Ob
    add_generator!(pres, SubOb(ob, super_ob))
  end
end

""" Add morphism from JSON document to presentation.
"""
function add_hom_generator_from_json!(pres::Presentation, doc::AbstractDict)
  dom_ob = domain_ob_from_json(pres, doc["inputs"])
  codom_ob = domain_ob_from_json(pres, doc["outputs"])
  hom = Hom(doc["id"], dom_ob, codom_ob)
  add_generator!(pres, hom)
end

""" Add submorphisms from JSON document to presentation.
"""
function add_subhom_generators_from_json!(pres::Presentation, doc::AbstractDict)
  hom = generator(pres, doc["id"])::Monocl.Hom
  names = get(doc, "is-a", [])
  names = isa(names, AbstractString) ? [ names ] : names
  for super_name in names
    super_hom = generator(pres, super_name)::Monocl.Hom
    # FIXME: Do basic type inference to check these subobject relations hold.
    subob_dom = SubOb(dom(hom), dom(super_hom))
    subob_codom = SubOb(codom(hom), codom(super_hom))
    add_generator!(pres, SubHom(hom, super_hom, subob_dom, subob_codom))
  end
end

function domain_ob_from_json(pres::Presentation, docs)::Monocl.Ob
  if isempty(docs)
    munit(Monocl.Ob)
  else
    otimes([ generator(pres, doc["type"]) for doc in docs ])
  end
end

# Annotations
#############

const language_keys = [
  "class", "function", "method", "inputs", "outputs", "slots",
]

""" Load annotation from JSON document.
"""
function annotation_from_json(doc::AbstractDict, load_ref::Function)::Annotation
  parse_def = sexpr -> parse_json_sexpr(Monocl, sexpr;
    symbols=false, parse_head=parse_json_sexpr_term, parse_reference=load_ref)
  name = AnnotationID(doc["language"], doc["package"], doc["id"])
  lang = Dict{Symbol,Any}(
    Symbol(key) => doc[key] for key in language_keys if haskey(doc, key)
  )
  definition = parse_def(doc["definition"])
  if doc["kind"] == "type"
    slots = [ parse_def(slot["definition"]) for slot in get(doc, "slots", []) ]
    ObAnnotation(name, lang, definition, slots)
  elseif doc["kind"] == "function"
    HomAnnotation(name, lang, definition)
  else
    error("Invalid kind of annotation: $(doc["kind"])")
  end
end
function annotation_from_json(doc::AbstractDict, pres::Presentation)
  annotation_from_json(doc, name -> generator(pres, name))
end

""" Replace term names (S-exp head) in JSON S-expression.

Translates PLT terminology into category theory terminology.
"""
function parse_json_sexpr_term(x::String)
  get(json_sexpr_term_table, x, x)
end

const json_sexpr_term_table = Dict(
  "Type" => "Ob",
  "Function" => "Hom",
  "Subtype" => "SubOb",
  "Subfunction" => "SubHom",
  "product" => "otimes",
  "unit" => "munit",
  "swap" => "braid",
  "copy" => "mcopy",
  "apply" => "pair",
)

end

ext: jl

max_stars_repo_name: IBM/semanticflowgraph

max_issues_repo_name: epatters/semanticflowgraph

max_forks_repo_name: IBM/semanticflowgraph


text:
using PQPolygonSDK #test package -ycpan
using Dates
using DataFrames

# build a user model -
options = Dict{String,Any}()
options["email"] = "yp392@cornell.edu"
options["apikey"] = "abc1234" # do _not_ check in a real API key

# build the user model -
user_model = model(PQPolygonSDKUserModel, options)

# now that we have the user_model, let's build an endpoint model -
endpoint_options = Dict{String,Any}()
endpoint_options["ticker"] = "AAPL"
endpoint_options["ex_dividend_date"] = "N/A"
endpoint_options["record_date"] = "N/A"
endpoint_options["declaration_date"] = "N/A"
endpoint_options["pay_date"] = "N/A"
endpoint_options["frequency"] = "4"
endpoint_options["cash_amount"] = "N/A"
endpoint_options["dividend_type"] = "CD"
endpoint_options["order"] = "desc"
endpoint_options["limit"] = "N/A"
endpoint_options["sort"] = "N/A"
endpoint_model = model(PolygonDividendsEndpointModel, user_model, endpoint_options)

# build the url -
base = "https://api.polygon.io"
my_url_string = url(base, endpoint_model)

# make the call -
(header, data) = api(PolygonDividendsEndpointModel, my_url_string)


text:
module Tagset

export encode, decode

"""
    I: begin or internal, O: outside, E: end
"""
immutable IOE
    I::Int
    O::Int
    E::Int
end
IOE() = IOE(1,2,3)

function encode(tagset::IOE, ranges::Vector{UnitRange{Int}})
    tags = fill(tagset.O, last(ranges[end]))
    for r in ranges
        tags[r] = tagset.I
        tags[last(r)] = tagset.E
    end
    tags
end

function decode(tagset::IOE, tags::Vector{Int})
    bpos = 0
    ranges = UnitRange{Int}[]
    for i = 1:length(tags)
        t = tags[i]
        t != tagset.O && bpos == 0 && (bpos = i)
        if t == tagset.E
            push!(ranges, bpos:i)
            bpos = 0
        end
    end
    ranges
end

end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
## Procedures for Controlling DReal Environment
## +===========================================

is_dreal_initialised = false

function init_dreal!()
  global is_dreal_initialised
  is_dreal_initialised == true && return
  # dlopen_shared_deps!()
  opensmt_init()
  is_dreal_initialised = true
end

# @doc "Set Verboisty of DReal output" ->
set_verbosity_level!(ctx::Context, level::Int) = opensmt_set_verbosity(ctx.ctx, @compat Int32(level))
set_verbosity_level!(level::Int) = set_verbosity_level!(global_context(), level)

# @doc "Set precision of a DReal" ->
set_precision!(ctx::Context, p::Float64) = opensmt_set_precision(ctx.ctx, p)


text:
using GLM, Distributions, DataFrames

data1 = parse.(Float64, readlines("../data/machine1.csv"))
data2 = parse.(Float64, readlines("../data/machine2.csv"))
data3 = parse.(Float64, readlines("../data/machine3.csv"))

function manualANOVA(allData)
    nArray = length.(allData)
    d = length(nArray)

    xBarTotal = mean(vcat(allData...))
    xBarArray = mean.(allData)

    ssBetween = sum( [nArray[i]*(xBarArray[i] - xBarTotal)^2 for i in 1:d] )
    ssWithin = sum([sum([(ob - xBarArray[i])^2 for ob in allData[i]])
				for i in 1:d])
    dfBetween = d-1
    dfError = sum(nArray)-d

    msBetween = ssBetween/dfBetween
    msError = ssWithin/dfError
    fStat = msBetween/msError
    pval = ccdf(FDist(dfBetween,dfError),fStat)
    return (fStat,pval)
end

function glmANOVA(allData)
    nArray = length.(allData)
    d = length(nArray)

    treatment = vcat([fill(k,nArray[k]) for k in 1:d]...)
    response = vcat(allData...)
    dataFrame = DataFrame(Response=response, Treatment=categorical(treatment))
    modelH0  = lm(@formula(Response ~ 1), dataFrame)
    modelH1a = lm(@formula(Response ~ 1 + Treatment), dataFrame)
    res = ftest(modelH1a.model, modelH0.model)
    (res.fstat[1],res.pval[1])
end

println("Manual ANOVA: ", manualANOVA([data1, data2, data3]))
println("GLM ANOVA: ", glmANOVA([data1, data2, data3]))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using Base.BitUnsigned_types

for U in BitUnsigned_types
    @eval begin
        fullshift(::Type{$U}) = $(unsigned(8*sizeof(U)))
        halfshift(::Type{$U}) = $(unsigned(4*sizeof(U)))
        halfmax(::Type{$U}) = $(typemax(U) >> 4*sizeof(U))
    end
end

fullshift(a::BitUnsigned) = fullshift(typeof(a))
halfshift(a::BitUnsigned) = halfshift(typeof(a))
halfmax(a::BitUnsigned) = halfmax(typeof(a))
_split(a::BitUnsigned) = a & halfmax(a), a >> halfshift(a)

_add(a::U, b::U) where {U<:BitUnsigned} = a+b, U(a > ~b)
_sub(a::U, b::U) where {U<:BitUnsigned} = a-b, U(b > a)
function _mul(a::U, b::U) where {U<:BitUnsigned}
    a0, a1 = _split(a)
    b0, b1 = _split(b)

    c00 = a0 * b0
    c11 = a1 * b1

    c01 = a0 * b1
    c10 = a1 * b0

    d1 = c11 + (c01 >> halfshift(U)) + (c10 >> halfshift(U))
    c01 <<= halfshift(U)
    c10 <<= halfshift(U)
    d1 += (c00 > ~c01)
    d0 = c00+c01
    d1 += (d0 > ~c10)
    d0 += c10

    return d0, d1
end

_divrem(n::U, d::U) where {U<:Unsigned} = divrem(n, d)
# divrem for a two word integer n1n0, provided that q is one word, i.e. n1 < d
function _divrem(n0::U, n1::U, d::U) where {U<:Unsigned}
    n1 > d && throw(OverflowError()) # TODO: remove when checked that we never call it otherwise
    iszero(n1) && return divrem(n0,d)

    # by splitting, the problem becomes: m3m2m1m0 / d1d0 -> q1q0 and r1r0
    if leading_zeros(d) >= halfshift(U)
        # simple case d1 = 0, also m3 = 0 given the promise n1 = m3m2 < d = d1d0
        # reduces to short division m2m1m0 / d0 -> q1q0 and r0
        d0 = d
        m2 = n1
        m0, m1 = _split(n0)
        q1, r1 = divrem( (m2 << halfshift(U)) + m1, d0)
        q0, r0 = divrem( (r1 << halfshift(U)) + m0, d0)
        return (q1 << halfshift(U)) + q0, r0
    else
        # first normalize
        s = leading_zeros(d)
        d = d << s
        n1 = (n1 << s) + (n0 >> (fullshift(U)-s)) # no overflow since n1 < d
        n0 = (n0 << s)

        # now split
        d0, d1 = _split(d) # d1 will be properly normalized
        # m3, m2 = _split(n1)
        m0, m1 = _split(n0)
        b = one(U) << halfshift(U)

        q̂, r̂ = divrem(n1, d1)
        # q̂ - q <= 2 (because of normalization) and q <= b-1
        # therefore q̂ <= b+1 and we can compute q̂*d0 without overflow
        if q̂*d0 > b*r̂+m1
            q̂ = q̂-one(q̂)
            r̂ += d1
            # repeat test (at most once since q̂ - q <= 2)
            if r̂ < b && q̂*d0 > b*r̂+m1
                q̂ = q̂-one(q̂)
                r̂ += d1
            end
        end
        q1 = q̂ # should now be exact

        # compute m3m2m1 - q*d1d0
        # = m3m2*b + m1 - q*d1*b - q*d0
        # = r*b + m1 - q*d0 # r*b might overflow, therefore subtract from r the high part of q*d0
        t0, t1 = _split(q1*d0)
        n1 = (r̂-t1) << halfshift(U) + (m1-t0)

        # now we want n1*b+m0 / (d1*b+d0) -> q0, r0
        q̂, r̂ = divrem(n1, d1)
        if q̂*d0 > b*r̂+m0
            q̂ = q̂-one(q̂)
            r̂ += d1
            # repeat test (at most once since q̂ - q <= 2)
            if r̂ < b && q̂*d0 > b*r̂+m0
                q̂ = q̂-one(q̂)
                r̂ += d1
            end
        end
        q0 = q̂ # should now be exact

        t0, t1 = _split(q0*d0)
        r = (r̂-t1) << halfshift(U) + (m0-t0)
        q = q1 << halfshift(U) + q0
        return q, (r >> s)
    end
end

# alternative implementation using widening: works up to UInt64
function _mul2(a::U, b::U) where {U<:Union{UInt8,UInt16,UInt32,UInt64}}
    t = widen(a)*widen(b)
    c = convert(U, t & typemax(U))
    d = convert(U, (t >> fullshift(U)))
    return c, d
end
function _divrem2(n0::U, n1::U, d::U) where {U<:Union{UInt8,UInt16,UInt32,UInt64}}
    q, r = divrem((widen(n1)<<fullshift(U))+widen(n0), widen(d))
    return convert(U, q), convert(U, r)
end

# estimate single word divisor q from first tree words of dividend and first two
# words of divisor
function _estimateq(u′′, u′, u, v′, v)
    if v == u
        q̂ = typemax(u)
        r̂′, r̂ = _add(u′, v)
    else
        q̂, r̂′ = _divrem(u′, u, v)
        r̂ = zero(u)
    end
    p′, p = _mul(q̂, v′)
    if iszero(r̂) && (p > r̂′ || (p == r̂′ && p′ > u′′))
        q̂ -= one(q̂)
        r̂′, carry = _add(r̂′, v)
        r̂ += carry
        p′, borrow = _sub(p′, v′)
        p -= borrow
    end
    if iszero(r̂) && (p > r̂′ || (p == r̂′ && p′ > u′′))
        q̂ -= one(q̂)
    end
    return q̂
end


text:
abstract type ConfigStruct end

#TODO: this type is probably not needed anymore
BSum = Union{Symbol, Tuple{Int,Int}}

@enum ChiFillType zero_χ_fill lambda_χ_fill χ_fill

"""
    ModelParameters <: ConfigStruct

Contains model parameters for the Hubbard model.
This is typically generated from a config.toml file using
the [`readConfig`](@ref readConfig) function.
"""
struct ModelParameters <: ConfigStruct
    U::Float64              # Hubbard U
    μ::Float64              # chemical potential
    β::Float64              # inverse temperature
    n::Float64              # number density
end
# TODO: save kgrid info here: kGrid::String           # String encoding information about the grid

"""
    SimulationParameters <: ConfigStruct

Contains simulation parameters for the ladder DGA computations.
This is typically generated from a config.toml file using
the [`readConfig`](@ref readConfig) function.

Fields
-------------
- **`n_iω`**    : Number of positive bosonic frequencies (full number will be `2*n_iω+1`
- **`n_iν`**    : Number of positive fermionic frequencies (full number will be `2*n_iν`
- **`shift`**   : Flag specifying if `-n_iν:n_iν-1` is shifted by `-ωₙ/2` at each `ωₙ` slice (centering the main features)
- **`tc_type_f`** : Symbol specifying the type of sum extrapolation for fermionic sums. Implemented for `:nothing` (naive sum) and `:richardson`. See `SeriesAcceleration` package for more details.
- **`tc_type_b`** : Symbol specifying the type of sum extrapolation for bosonic sums. Implemented for `:nothing` (naive sum), `:coeffs` (subtracted known tail coefficients) and `:richardson`. See `SeriesAcceleration` package for more details.
- **`λc_type`** : Symbol specifying the type of λ-correction. `:nothing` for no correction, `:sp` for correction in spin channel, using the filling of the model as reference for the fit, `:sp_ch` for λ-corrections in spin an charge channel, `:sp_ch_q` for k-vector dependent λ correction in both channels.
- **`ωsum_type`**  :
- **`λ_rhs`**   :
- **`fullChi`** :
- **`χFillType`**  :
- **`bosonic_tail_coeffs`**   :
- **`fermionic_tail_coeffs`** :
- **`usable_prct_reduction`** :
"""
struct SimulationParameters <: ConfigStruct
    n_iω::Int64             # number of bosonic frequencies
    n_iν::Int64             # number of fermionic frequencies
    shift::Bool            # shift of center for interval of bosonic frequencies
    tc_type_f::Symbol  # use correction for finite ν sums.
    tc_type_b::Symbol  # use correction for finite ω sums.
    λc_type::Symbol  # which type of lambda correction to use (currecntly: nothing, sp, sp_ch, TOOD: sp_ch_q
    ωsum_type::BSum
    λ_rhs::Symbol
    fullChi::Bool
    χFillType::ChiFillType # values to be set outside the usable interval
    bosonic_tail_coeffs::Array{Int,1}   # tail
    fermionic_tail_coeffs::Array{Int,1}
    usable_prct_reduction::Float64      # safety cutoff for usable ranges
    ω_smoothing::Symbol                 # nothing, range, full
    dbg_full_eom_omega::Bool
end

"""
    EnvironmentVars <: ConfigStruct

Contains various settings, controlling the I/O behaviour
of this module.
This is typically generated from a config.toml file using
the [`readConfig`](@ref readConfig) function.
"""
struct EnvironmentVars <: ConfigStruct
    inputDataType::String
    writeFortran::Bool
    loadAsymptotics::Bool
    inputDir::String
    freqFile::String
    inputVars::String
    asymptVars::String
    cast_to_real::Bool
    loglevel::String      # disabled, error, warn, info, debug
    logfile::String       # STDOUT, STDERR, filename
    progressbar::Bool     # true/false enable or disable progress bar
end


text:
using AsymptoticPosteriors
using Base.Test

# write your own tests here
@test 1 == 2

include(joinpath(dirname(pathof(AsymptoticPosteriors)), "../examples/NGSCov.jl"));
# include("/home/chris/Documents/progwork/julia/NGStan.jl");
@time ap = AsymptoticPosterior(n_big_4, init8, Val(8));
@time apc = AsymptoticPosterior(n_big_4c, initCorr10, Val(10));
summary(ap)
summary(apc)

using BenchmarkTools
# @benchmark AsymptoticPosterior($n_big_4, $init8, Val(8))
# @benchmark four_num_sum.(($ap,), 1:8)
@benchmark AsymptoticPosterior($n_big_2c, $initCorr8, Val(8))
@benchmark four_num_sum.(($apc,), 1:10)

low_sensitivity = CorrErrors((0.1,0.2,0.3,0.4), (0.8, 0.8), (0.98,0.98));
lowS1 = NoGoldDataCorr(low_sensitivity,100,200,16);
ap_lows1 = AsymptoticPosterior(lowS1, initCorr10, Val(10));
summary(ap_lows1)
rstan(lowS1, 10_000)
using RCall
R"options(width=160)"
summary(ap_lows1)
R"summary(res)$summary"

low_specificity = CorrErrors((0.1,0.2,0.3,0.4), (0.98, 0.98), (0.8,0.8));
lowC1 = NoGoldDataCorr(low_specificity,100,200,16);
ap_lowc1 = AsymptoticPosterior(lowC1, initCorr10, Val(10));
summary(ap_lowc1)
rstan(lowC1, 10_000)
summary(ap_lowc1)
R"summary(res)$summary"

n =  [13   12; 17   20; 11   11; 59   57; 279  308; 721  692; 11   15; 19   15];
test_data = NoGoldDataCorr{8}(n =n,
αS1 = 16.0, βS1 = 4.0,
αS2 = 16.0, βS2 = 4.0,
αC1 = 16.0, βC1 = 4.0,
αC2 = 16.0, βC2 = 4.0);
@time apc = AsymptoticPosterior(test_data, initCorr8, Val(8));
summary(apc)

using AsymptoticPosteriors
include("/home/chris/Documents/progwork/julia/NGS.jl")
f3(x) = NGS(cdb4, x, Val{8}())
@time ap = AsymptoticPosterior(f3, initial_x, Val(8));
quantile(ap, 0.025)

interval(ap, ind) = (quantile(ap, 0.025, ind), quantile(ap, 0.975, ind))
@benchmark interval($ap, 1)
@benchmark interval($ap, 2)
@benchmark interval($ap, 3)
@benchmark interval($ap, 4)
@benchmark interval($ap, 5)
@benchmark interval($ap, 6)
@benchmark interval($ap, 7)
@benchmark interval($ap, 8)

@benchmark AsymptoticPosterior(f3, $initial_x, Val(8))
@benchmark AsymptoticPosteriors.fit!($ap.pl.map, $initial_x)

using Revise, AsymptoticPosteriors
using Traceur
@code_warntype AsymptoticPosteriors.AsymptoticConfig(f2, max_density_params, Val{8}())
config = @trace AsymptoticPosteriors.AsymptoticConfig(f2, max_density_params, Val{8}());
# config = AsymptoticPosteriors.AsymptoticConfig(f2, max_density_params, Val{8}());

@code_warntype AsymptoticPosteriors.AsymptoticMAP(max_density_params, config, Val{8}())
# map = @trace AsymptoticPosteriors.AsymptoticMAP(max_density_params, config, Val{8}());
map = AsymptoticPosteriors.AsymptoticMAP(max_density_params, config, Val{8}());

@code_warntype AsymptoticPosterior(map)
# ap = AsymptoticPosterior(map)

# ap = AsymptoticPosterior(f2, max_density_params, Val{8}())

function odgen(map::M) where {T, N, M <: AsymptoticPosteriors.AsymptoticMAP{N,T}}
    F = x -> begin
        map::M
        AsymptoticPosteriors.arrays_equal(x, map.last_x) ? map.lastval[] : AsymptoticPosteriors.update_grad!(map, x)
    end

    G = (out, x) -> begin
        AsymptoticPosteriors.update_grad!(out, map::M, x)
        out
    end

    FG = (out, x) -> begin
        map::M
        AsymptoticPosteriors.update_grad!(out, map, x)
        map.lastval[]
    end

    initial_x = Vector{T}(undef, N+1)

    @inbounds begin
        for i ∈ 1:N
            initial_x[i] = map.θhat[i]
        end
        initial_x[end] = zero(T)
    end
    OnceDifferentiable(F, G, FG, initial_x, zero(T), Val{true}())
end

odmax = odgen(map);

backtrack = BFGS(LineSearches.InitialStatic(), LineSearches.BackTracking(), x -> Matrix{eltype(x)}(I, length(x), length(x)), Optim.Flat())

options = Optim.InternalUseOptions(backtrack)

# state = Optim.initial_state(backtrack, options, odmax, initial_x)

fieldnames(typeof(odmax))

xtest = vcat(max_density_params, [1]);
theta = @view(xtest[1:8]);
out = similar(xtest);

AsymptoticPosteriors.update_hessian!(map, theta)
AsymptoticPosteriors.qb(map, theta )

@code_warntype odmax.f(xtest)

@code_warntype AsymptoticPosteriors.update_grad!(map, xtest)

# Purterbation confusion!?!
AsymptoticPosteriors.update_grad!(map, max_density_params)
# ERROR: Cannot determine ordering of Dual tags Void and ForwardDiff.Tag{#f2,Float64}
# Where does the Tag{#f2,Float64} dual tag come from?
# julia> typeof(map.configuration.theta)
# SubArray{ForwardDiff.Dual{ForwardDiff.Tag{#f2,Float64},Float64,8},1,Array{ForwardDiff.Dual{ForwardDiff.Tag{#f2,Float64},Float64,8},1},Tuple{UnitRange{Int64}},true}

# Run through WTF theta is supposed to be in each use case! Seems wrong.

ap = AsymptoticPosterior(f2, max_density_params, Val{8}());

ap.map.θhat |> detransform_res

AsymptoticPosteriors.bound(ap, 1, 0.025)

precise_options =  Optim.Options(g_tol = 1e-30,x_tol=0.0,f_tol=0.0)
opt_precise = optimize(f2, BigFloat.(max_density_params), hage, precise_options)
prex_max_dens = Optim.minimizer(opt_precise)

xbf = prex_max_dens[1]
mdp_bf = prex_max_dens[2:end]

opt_1 = optimize(y -> f2(vcat(xbf-eps()/2,y)), mdp_bf, hage, precise_options)
min1 = Optim.minimizer(opt_1)
minimum_1 = Optim.minimum(opt_1)
opt_2 = optimize(y -> f2(vcat(xbf+eps()/2,y)), mdp_bf, hage, precise_options)
min2 = Optim.minimizer(opt_2)
minimum_2 = Optim.minimum(opt_2)

grad = @. (min2 - min1)/eps()
dlp = (minimum_2 - minimum_1) / eps()

precise_hessian = ForwardDiff.hessian(f2, Float64.(prex_max_dens))
estimated_grad = -inv(BigFloat.(precise_hessian[2:end,2:end]))*BigFloat.(precise_hessian[2:end,1])
grad .- estimated_grad

precise_grad = ForwardDiff.gradient(f2, Float64.(prex_max_dens))
estimated_dlp = precise_grad' * vcat([1.0],estimated_grad)
dlp - estimated_dlp

xbf2 = xbf - 1.96 * inv(precise_hessian)[1]
opt_precise2 = optimize(y -> f2(vcat(xbf2,y)), mdp_bf, hage, precise_options)
prex_max_dens2 = Optim.minimizer(opt_precise2)
prex_max_dens2_ = vcat(xbf2, prex_max_dens2)

opt_3 = optimize(y -> f2(vcat(xbf2-eps()/2,y)), mdp_bf, hage, precise_options)
min3 = Optim.minimizer(opt_3)
minimum_3 = Optim.minimum(opt_3)
opt_4 = optimize(y -> f2(vcat(xbf2+eps()/2,y)), mdp_bf, hage, precise_options)
min4 = Optim.minimizer(opt_4)
minimum_4 = Optim.minimum(opt_4)

grad2 = @. (min4 - min3)/eps()
dlp2 = (minimum_4 - minimum_3) / eps()

precise_hessian2 = ForwardDiff.hessian(f2, Float64.(prex_max_dens2_))
estimated_grad2 = -inv(BigFloat.(precise_hessian2[2:end,2:end]))*BigFloat.(precise_hessian2[2:end,1])
grad2 .- estimated_grad2

precise_grad2 = ForwardDiff.gradient(f2, Float64.(prex_max_dens2_))
estimated_dlp2 = precise_grad2' * vcat([1.0],estimated_grad2)
dlp2 - estimated_dlp2

phi = inv(precise_hessian2);
estimated_dlp3 = precise_grad2[1] / (phi[1] * precise_hessian2[1] )


text:
#######################################################################
# Debugging IJulia

# in the Jupyter front-end, enable verbose output via IJulia.set_verbose()
verbose = false
"""
    set_verbose(v=true)

This function enables (or disables, for `set_verbose(false)`) verbose
output from the IJulia kernel, when called within a running notebook.
This consists of log messages printed to the terminal window where
`jupyter` was launched, displaying information about every message sent
or received by the kernel.   Used for debugging IJulia.
"""
function set_verbose(v::Bool=true)
    global verbose = v
end

# set this to false for debugging, to disable stderr redirection
"""
The IJulia kernel captures all [stdout and stderr](https://en.wikipedia.org/wiki/Standard_streams)
output and redirects it to the notebook.   When debugging IJulia problems,
however, it can be more convenient to *not* capture stdout and stderr output
(since the notebook may not be functioning). This can be done by editing
`IJulia.jl` to set `capture_stderr` and/or `capture_stdout` to `false`.
"""
const capture_stdout = true
const capture_stderr = true

include("init.jl")
include("stdio.jl")
include("display.jl")
include("execute_request.jl")
include("inline.jl")


text:
using ClineSimulation
using Test

# Tests for the initialfreqs function
@testset "initial frequencies function" begin
    @test_throws "Total number of demes must be even and positive" initialfreqs(-1, 0, 1)
    @test_throws "Total number of demes must be even and positive" initialfreqs(0, 0, 1)
    @test_throws "Total number of demes must be even and positive" initialfreqs(1, 0, 1)
    @test_throws "Allelic frequencies must be different" initialfreqs(2, 0, 0)
    @test_throws "Allelic frequencies must be different" initialfreqs(2, 0.51, 0.51)
    @test_throws "Allelic frequencies must be different" initialfreqs(2, 1, 1)
    @test_throws "Allelic frequencies must be between 0 and 1" initialfreqs(2, 0, 1.1)
    @test_throws "Allelic frequencies must be between 0 and 1" initialfreqs(2, -0.1, 1)
    @test initialfreqs(4, 0, 1) == [0 0 1 1]
    @test initialfreqs(2, 0.51, 0.52) == [0.51 0.52]
    @test length(initialfreqs(100, 0, 1)) == 100
end


text:
using ParameterizedFunctions, OrdinaryDiffEq, DiffEqParamEstim
using BlackBoxOptim, NLopt, Plots, RecursiveArrayTools, QuadDIRECT
gr(fmt=:png)

loc_bounds = Tuple{Float64, Float64}[(0, 5), (0, 5), (0, 5), (0, 5)]
glo_bounds = Tuple{Float64, Float64}[(0, 10), (0, 10), (0, 10), (0, 10)]
loc_init = [1,0.5,3.5,1.5]
glo_init = [5,5,5,5]

f = @ode_def LotkaVolterraTest begin
    dx = a*x - b*x*y
    dy = -c*y + d*x*y
end a b c d

u0 = [1.0,1.0]                          #initial values
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1,0]                   #parameters used, these need to be estimated from the data
tspan = (0.0, 30.0)                     # sample of 3000 observations over the (0,30) timespan
prob = ODEProblem(f, u0, tspan,p)
tspan2 = (0.0, 3.0)                     # sample of 3000 observations over the (0,30) timespan
prob_short = ODEProblem(f, u0, tspan2,p)

dt = 30.0/3000
tf = 30.0
tinterval = 0:dt:tf
t  = collect(tinterval)

h = 0.01
M = 300
tstart = 0.0
tstop = tstart + M * h
tinterval_short = 0:h:tstop
t_short = collect(tinterval_short)

#Generate Data
data_sol_short = solve(prob_short,Tsit5(),saveat=t_short,reltol=1e-9,abstol=1e-9)
data_short = convert(Array, data_sol_short)
data_sol = solve(prob,Tsit5(),saveat=t,reltol=1e-9,abstol=1e-9)
data = convert(Array, data_sol)

p1 = plot(data_sol_short)

p2 = plot(data_sol)

obj_short = build_loss_objective(prob_short,Tsit5(),L2Loss(t_short,data_short),tstops=t_short)
res1 = bboptimize(obj_short;SearchRange = glo_bounds, MaxSteps = 7e3)
# Lower tolerance could lead to smaller fitness (more accuracy)

obj_short = build_loss_objective(prob_short,Tsit5(),L2Loss(t_short,data_short),tstops=t_short,reltol=1e-9)
res1 = bboptimize(obj_short;SearchRange = glo_bounds, MaxSteps = 7e3)
# Change in tolerance makes it worse

obj_short = build_loss_objective(prob_short,Vern9(),L2Loss(t_short,data_short),tstops=t_short,reltol=1e-9,abstol=1e-9)
res1 = bboptimize(obj_short;SearchRange = glo_bounds, MaxSteps = 7e3)
# using the moe accurate Vern9() reduces the fitness marginally and leads to some increase in time taken

obj_short = build_loss_objective(prob_short,Vern9(),L2Loss(t_short,data_short),tstops=t_short,reltol=1e-9,abstol=1e-9)

opt = Opt(:GN_ORIG_DIRECT_L, 4)
lower_bounds!(opt,[0.0,0.0,0.0,0.0])
upper_bounds!(opt,[10.0,10.0,10.0,10.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,glo_init)

opt = Opt(:GN_CRS2_LM, 4)
lower_bounds!(opt,[0.0,0.0,0.0,0.0])
upper_bounds!(opt,[10.0,10.0,10.0,10.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,glo_init)

opt = Opt(:GN_ISRES, 4)
lower_bounds!(opt,[0.0,0.0,0.0,0.0])
upper_bounds!(opt,[10.0,10.0,10.0,10.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,glo_init)

opt = Opt(:GN_ESCH, 4)
lower_bounds!(opt,[0.0,0.0,0.0,0.0])
upper_bounds!(opt,[10.0,10.0,10.0,10.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,glo_init)

opt = Opt(:LN_BOBYQA, 4)
lower_bounds!(opt,[0.0,0.0,0.0,0.0])
upper_bounds!(opt,[5.0,5.0,5.0,5.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,loc_init)

opt = Opt(:LN_NELDERMEAD, 4)
lower_bounds!(opt,[0.0,0.0,0.0,0.0])
upper_bounds!(opt,[5.0,5.0,5.0,5.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,loc_init)

opt = Opt(:LD_SLSQP, 4)
lower_bounds!(opt,[0.0,0.0,0.0,0.0])
upper_bounds!(opt,[5.0,5.0,5.0,5.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,loc_init)

opt = Opt(:LN_COBYLA, 4)
lower_bounds!(opt,[0.0,0.0,0.0,0.0])
upper_bounds!(opt,[5.0,5.0,5.0,5.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,loc_init)

opt = Opt(:LN_NEWUOA_BOUND, 4)
lower_bounds!(opt,[0.0,0.0,0.0,0.0])
upper_bounds!(opt,[5.0,5.0,5.0,5.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,loc_init)

opt = Opt(:LN_PRAXIS, 4)
lower_bounds!(opt,[0.0,0.0,0.0,0.0])
upper_bounds!(opt,[5.0,5.0,5.0,5.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,loc_init)

opt = Opt(:LN_SBPLX, 4)
lower_bounds!(opt,[0.0,0.0,0.0,0.0])
upper_bounds!(opt,[5.0,5.0,5.0,5.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,loc_init)

opt = Opt(:LD_MMA, 4)
lower_bounds!(opt,[0.0,0.0,0.0,0.0])
upper_bounds!(opt,[5.0,5.0,5.0,5.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,loc_init)

opt = Opt(:LD_TNEWTON_PRECOND_RESTART, 4)
lower_bounds!(opt,[0.0,0.0,0.0,0.0])
upper_bounds!(opt,[5.0,5.0,5.0,5.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,loc_init)

obj = build_loss_objective(prob,Vern9(),L2Loss(t,data),tstops=t,reltol=1e-9,abstol=1e-9)
res1 = bboptimize(obj;SearchRange = glo_bounds, MaxSteps = 4e3)

opt = Opt(:GN_ORIG_DIRECT_L, 4)
lower_bounds!(opt,[0.0,0.0,0.0,0.0])
upper_bounds!(opt,[10.0,10.0,10.0,10.0])
min_objective!(opt, obj.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,glo_init)

opt = Opt(:GN_CRS2_LM, 4)
lower_bounds!(opt,[0.0,0.0,0.0,0.0])
upper_bounds!(opt,[10.0,10.0,10.0,10.0])
min_objective!(opt, obj.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 20000)
@time (minf,minx,ret) = NLopt.optimize(opt,glo_init)

opt = Opt(:GN_ISRES, 4)
lower_bounds!(opt,[0.0,0.0,0.0,0.0])
upper_bounds!(opt,[10.0,10.0,10.0,10.0])
min_objective!(opt, obj.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 50000)
@time (minf,minx,ret) = NLopt.optimize(opt,glo_init)

opt = Opt(:GN_ESCH, 4)
lower_bounds!(opt,[0.0,0.0,0.0,0.0])
upper_bounds!(opt,[10.0,10.0,10.0,10.0])
min_objective!(opt, obj.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 20000)
@time (minf,minx,ret) = NLopt.optimize(opt,glo_init)

opt = Opt(:LN_BOBYQA, 4)
lower_bounds!(opt,[0.0,0.0,0.0,0.0])
upper_bounds!(opt,[5.0,5.0,5.0,5.0])
min_objective!(opt, obj.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,loc_init)

opt = Opt(:LN_NELDERMEAD, 4)
lower_bounds!(opt,[0.0,0.0,0.0,0.0])
upper_bounds!(opt,[5.0,5.0,5.0,5.0])
min_objective!(opt, obj.cost_function2)
xtol_rel!(opt,1e-9)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,loc_init)

opt = Opt(:LD_SLSQP, 4)
lower_bounds!(opt,[0.0,0.0,0.0,0.0])
upper_bounds!(opt,[5.0,5.0,5.0,5.0])
min_objective!(opt, obj.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,loc_init)

obj_short = build_loss_objective(prob_short,Tsit5(),L2Loss(t_short,data_short),tstops=t_short)
lower = [0.0,0.0,0.0,0.0]
upper = [5.0,5.0,5.0,5.0]
splits = ([0.0,1.0,3.0],[0.0,1.0,3.0],[0.0,1.0,3.0],[0.0,1.0,3.0])
root, x0 = analyze(obj_short,splits,lower,upper)

minimum(root)

obj = build_loss_objective(prob,Vern9(),L2Loss(t,data),tstops=t,reltol=1e-9,abstol=1e-9)
lower = [0.0,0.0,0.0,0.0]
upper = [10.0,10.0,10.0,10.0]
splits = ([0.0,3.0,6.0],[0.0,3.0,6.0],[0.0,3.0,6.0],[0.0,3.0,6.0])
root, x0 = analyze(obj,splits,lower,upper)

minimum(root)


text:
@inline fractional_x_index(x, ::Type{Cell}, grid::RegularCartesianGrid) = @inbounds (x - grid.xC[1]) / grid.Δx
@inline fractional_x_index(x, ::Type{Face}, grid::RegularCartesianGrid) = @inbounds (x - grid.xF[1]) / grid.Δx

@inline fractional_y_index(y, ::Type{Cell}, grid::RegularCartesianGrid) = @inbounds (y - grid.yC[1]) / grid.Δy
@inline fractional_y_index(y, ::Type{Face}, grid::RegularCartesianGrid) = @inbounds (y - grid.yF[1]) / grid.Δy

@inline fractional_z_index(z, ::Type{Cell}, grid::RegularCartesianGrid) = @inbounds (z - grid.zC[1]) / grid.Δz
@inline fractional_z_index(z, ::Type{Face}, grid::RegularCartesianGrid) = @inbounds (z - grid.zF[1]) / grid.Δz

"""
    fractional_indices(x, y, z, loc, grid::RegularCartesianGrid)

Convert the coordinates `(x, y, z)` to _fractional_ indices on a regular Cartesian grid located at `loc`
where `loc` is a 3-tuple of `Cell` and `Face`. Fractional indices are floats indicating a location between
grid points.
"""
@inline function fractional_indices(x, y, z, loc, grid::RegularCartesianGrid)
    i = fractional_x_index(x, loc[1], grid)
    j = fractional_y_index(y, loc[2], grid)
    k = fractional_z_index(z, loc[3], grid)
    return i, j, k
end

# Trilinear Lagrange polynomials
@inline ϕ₁(ξ, η, ζ) = (1 - ξ) * (1 - η) * (1 - ζ)
@inline ϕ₂(ξ, η, ζ) = (1 - ξ) * (1 - η) * ζ
@inline ϕ₃(ξ, η, ζ) = (1 - ξ) * η * (1 - ζ)
@inline ϕ₄(ξ, η, ζ) = (1 - ξ) * η * ζ
@inline ϕ₅(ξ, η, ζ) = ξ * (1 - η) * (1 - ζ)
@inline ϕ₆(ξ, η, ζ) = ξ * (1 - η) * ζ
@inline ϕ₇(ξ, η, ζ) = ξ * η * (1 - ζ)
@inline ϕ₈(ξ, η, ζ) = ξ * η * ζ

@inline _interpolate(field, ξ, η, ζ, i, j, k) =
    @inbounds (  ϕ₁(ξ, η, ζ) * field[i,   j,   k  ]
               + ϕ₂(ξ, η, ζ) * field[i,   j,   k+1]
               + ϕ₃(ξ, η, ζ) * field[i,   j+1, k  ]
               + ϕ₄(ξ, η, ζ) * field[i,   j+1, k+1]
               + ϕ₅(ξ, η, ζ) * field[i+1, j,   k  ]
               + ϕ₆(ξ, η, ζ) * field[i+1, j,   k+1]
               + ϕ₇(ξ, η, ζ) * field[i+1, j+1, k  ]
               + ϕ₈(ξ, η, ζ) * field[i+1, j+1, k+1])

"""
    interpolate(field, x, y, z)

Interpolate `field` to the physical point `(x, y, z)` using trilinear interpolation.
"""
@inline function interpolate(field, x, y, z)
    i, j, k = fractional_indices(x, y, z, location(field), field.grid)

    # Convert fractional indices to unit cell coordinates 0 <= (ξ, η, ζ) <=1
    # and integer indices (with 0-based indexing).
    ξ, i = modf(i)
    η, j = modf(j)
    ζ, k = modf(k)

    # Convert indices to proper integers and shift to 1-based indexing.
    return _interpolate(field.data, ξ, η, ζ, Int(i+1), Int(j+1), Int(k+1))
end


text:
## Reproducibility instructions

# 1) Open Julia Pro/Atom
# 2) Juno: Work in folder (right click on folder)
# 3) Juno: Activate environment in folder (right click on folder)

## Create environments

# 1) set in Projects folder
# 2) julia -e 'Pkg.generate("Project_Name")'
# 3) set in project folder
# 4) julia --project=. / Reproducibility instructions

## Saving package versions in Project.toml and Manifest.toml
Pkg.add("Distributed", preserve=PRESERVE_DIRECT)
Pkg.add(["DataFrames","CSV",
         "DataFramesMeta", "LinearAlgebra",
         "Parameters", "JuMP", "Ipopt", "SharedArrays",
         "Statistics","BlackBoxOptim", "NLopt", "Distributions",
         "StatsBase"], preserve=PRESERVE_DIRECT)

max_stars_repo_name: hans-mtz/SCR

max_issues_repo_name: hans-mtz/SCR

max_forks_repo_name: hans-mtz/SCR


text:
# -----------------------------------------------------------------------
# Decimal coded genetic algorithm
# -----------------------------------------------------------------------
# Reference: JENKINS, W. M. A decimal-coded evolutionary algorithm for
# constrained optimization. Computers and Structures, v. 80,
# p. 471–480, 2002.
# -----------------------------------------------------------------------
# Required packages
using DataFrames, Statistics
include("random_init.jl")

# Algorithm implemented according to the following paper:
# "A decimal-coded evolutionary algorithm for constrained optimization", W.M. Jenkins, 2002
# The algorithm is decimal coded, mutation based. There is no crossover.

# Mutation

"""
    αmax(gen::Int; α0::Real, maxgen::Int)

## Description
Function to compute the maximum shift for the generation
## Arguments
- gen: current generation
- α0: shift for the first generation
- maxgen: total number of generations
"""
αmax(gen::Int; α0::Real, maxgen::Int) = α0*(1-gen/maxgen)^2 + 1

"""
αind(gen::Int, ind::Int; α0::Real, maxgen::Int, ucb::Int, popsize::Int, rounded=true)

## Description
Function to compute the shift for each individual acconding to the current generation
## Arguments
- gen: current generation
- ind: individual id
- α0: shift for the first generation
- maxgen: total number of generations
- ucb: upper class boundary (elitism)
- popsize: population size
- digits: decimal places
"""
function αind(gen::Int, ind::Int; α0::Real, maxgen::Int, ucb::Int, popsize::Int, digits=0)

    shift_gen = αmax(gen, α0=α0, maxgen=maxgen)

    α = (ind - ucb)*(shift_gen - 1)/(popsize - ucb) + 1
    α = round(α, digits=digits)

    return α
end

# Elitism

"""
    elite!(pop::DataFrame; lcb::Int, cols= :fitness)

## Description
Order population from the most well adapted induvidual to the least adapted.
Also, remove individuals from the lower boundary and replace them by the most well adapted
## Arguments
- pop: population dataframe with coordinates and fitness
- lcb: lower class boundary
- cols: idicates the finess column
"""
function elite!(pop::DataFrame; lcb::Int, cols= :fitness)

    sort!(pop, cols)
    delete!(pop, 1:lcb)
    sort!(pop, cols, rev=true)
    append!(pop, pop[1:lcb,:])

    return pop
end

# Boundaries

"""
    isindomain(coord::Real, lims)

## Description
Identify if a coordinate is inside a domain
## Arguments
- coord: coordinate
- lims: boundaries (range or 1D array)
## Example
```jldoctest
julia> isindomain(10, 1:4)
false

julia> isindomain(10, 1:11)
true
```
"""
isindomain(coord::Real, lims) = coord >= lims[1] && coord <= lims[end]

"""
    findnearest(val::Real, lims)

## Description
Find the nearest neighbor in the case the coordinate is outside the domain
## Arguments
- coord: coordinate
- lims: noundaries (range or 1D array)
"""
function findnearest(coord::Real, lims)

    if coord < lims[1]
        return lims[1]
    elseif coord > lims[end]
        return lims[end]
    else
        error("The coordinate should be outside the domain!")
    end
end

"""
    indomain_coord!(coord::Real, boundaries)

## Description
Make all the coordinates stay inside the boundary
## Arguments
- coord: coordinate
- boundaries: domain limits for the direction
"""
function indomain_coord(coord::Real, boundaries)

    if !isindomain(coord, boundaries)
        coord_new = findnearest(coord, boundaries)
    else
        coord_new = coord
    end

    return coord_new
end

"""
    gadecimal_step(pop::DataFrame; Tm::AbstractVector, α0::AbstractVector, gen::int, maxgen::Int, ucb::Int, lcb::Int, digits=0)

## Description
One step of decimal-coded GA algorithm
## Arguments
- pop: input population
- Tm: mutation rate for each direction
- α0: shift for the first generation
- gen: current generation
- maxgen: total number of generations
- ucb: upper class boundary (elitism)
- lcb: lower class boundary
- digits: decimal places
"""
function gadecimal_step(pop::DataFrame; Tm::AbstractVector, α0::AbstractVector, gen::Int, maxgen::Int,
                        ucb::Int, lcb::Int, digits=0, boundaries::AbstractVector)

    n = nrow(pop)

    pop_new = deepcopy(pop)
    elite!(pop_new, lcb=lcb)
    new_individuals = []

    for i in (ucb+1):n

        ncoord = length(α0)

        α = [αind(gen, i, α0=a, maxgen=maxgen, ucb=ucb, popsize=n, digits=digits) for a in α0]
        mutate = rand(ncoord) .> Tm

        if sum(mutate) != 0
            push!(new_individuals, i)
        end

        for j in 1:ncoord
            if !mutate[j]
                nothing
            else
                op = rand([+,-])
                pop_new[i,j] = op(pop_new[i,j], α[j])
                pop_new[i,j] = indomain_coord(pop_new[i,j], boundaries[j])
            end
        end
    end

    return pop_new, new_individuals
end

"""
    gadecimal_criteria(df::AbstractVector; ucb::Int, mean_tol::AbstractVector, std_tol::AbstractVector)

## Description
Criteria for the decimal-coded GA algorithm.
The criteria is based on the mean stability of the upper class boundary (3 last generations are considered) and on its standard deviation.
## Arguments
- df: vector with all the dataframes obtained during the process (check experiments.jl)
- ucb: upper class boundary
- mean_tol: vector with the tolerance for the mean values
- std_tol: vector with the tolerance for the std values
"""
function gadecimal_criteria(df::AbstractVector; ucb::Int, mean_tol::AbstractVector, std_tol::AbstractVector)

    n = length(df)
    i = n - 1
    k = n - 2

    if k <= 0
        return false
    else
        ncoord = ncol(df[end]) - 1
        mean_end = [mean(df[end][1:ucb,q]) for q in 1:ncoord]
        std_end = [std(df[end][1:ucb,q]) for q in 1:ncoord]
        mean_i = [mean(df[i][1:ucb,q]) for q in 1:ncoord]
        mean_k = [mean(df[k][1:ucb,q]) for q in 1:ncoord]

        Δmean_i = abs.(mean_end - mean_i)
        Δmean_k = abs.(mean_end - mean_k)

        if all(Δmean_i .< mean_tol) && all(Δmean_k .< mean_tol) && all(std_end .< std_tol)
            return true
        else
            return false
        end
    end
end