text:
import JLD

# The number of components in the color prior.
const NUM_COLOR_COMPONENTS = 8

# The number of types of celestial objects (here, stars and galaxies).
const NUM_SOURCE_TYPES = 2

mutable struct CatalogEntry
    pos::Vector{Float64}
    is_star::Bool
    star_fluxes::Vector{Float64}
    gal_fluxes::Vector{Float64}
    gal_frac_dev::Float64
    gal_ab::Float64
    gal_angle::Float64
    gal_scale::Float64
end

"""
Parameters of a single normal component of a galaxy.

Attributes:
  etaBar: The weight of the galaxy component
  nuBar: The scale of the galaxy component
"""
struct GalaxyComponent
    etaBar::Float64
    nuBar::Float64
end

const GalaxyPrototype = Vector{GalaxyComponent}

"""
Pre-defined shapes for galaxies.

Returns:
  dev_prototype: An array of GalaxyComponent for de Vaucouleurs galaxy types
  exp_prototype: An array of GalaxyComponent for exponenttial galaxy types
"""
function get_galaxy_prototypes()
    dev_amp = [
        4.26347652e-2, 2.40127183e-1, 6.85907632e-1, 1.51937350,
        2.83627243, 4.46467501, 5.72440830, 5.60989349]
    dev_amp /= sum(dev_amp)
    dev_var = [
        2.23759216e-4, 1.00220099e-3, 4.18731126e-3, 1.69432589e-2,
        6.84850479e-2, 2.87207080e-1, 1.33320254, 8.40215071]

    exp_amp = [
        2.34853813e-3, 3.07995260e-2, 2.23364214e-1,
        1.17949102, 4.33873750, 5.99820770]
    exp_amp /= sum(exp_amp)
    exp_var = [
        1.20078965e-3, 8.84526493e-3, 3.91463084e-2,
        1.39976817e-1, 4.60962500e-1, 1.50159566]

    # Adjustments to the effective radius hard-coded above.
    # (The effective radius is the distance from the center containing half
    # the light.)
    effective_radii = [1.078031, 0.928896]
    dev_var /= effective_radii[1]^2
    exp_var /= effective_radii[2]^2

    exp_prototype = [GalaxyComponent(exp_amp[j], exp_var[j]) for j in 1:6]
    dev_prototype = [GalaxyComponent(dev_amp[j], dev_var[j]) for j in 1:8]
    (dev_prototype, exp_prototype)
end

const galaxy_prototypes = get_galaxy_prototypes()

struct PriorParams
    is_star::Vector{Float64}
    flux_mean::Vector{Float64}
    flux_var::Vector{Float64}
    k::Matrix{Float64}
    color_mean::Array{Float64, 3}
    color_cov::Array{Float64, 4}
    gal_scale_mean::Float64
    gal_scale_var::Float64
end

function load_prior()
    # set a = [.99, .01] if stars are underrepresented
    # due to the greater flexibility of the galaxy model
    #is_star = [0.28, 0.72]
    is_star = [0.099, 0.001]
    k = Matrix{Float64}(NUM_COLOR_COMPONENTS, NUM_SOURCE_TYPES)
    color_mean = Array{Float64}(NUM_BANDS - 1, NUM_COLOR_COMPONENTS, NUM_SOURCE_TYPES)
    color_cov = Array{Float64}(NUM_BANDS - 1, NUM_BANDS - 1, NUM_COLOR_COMPONENTS, NUM_SOURCE_TYPES)

    prior_params = [JLD.load(joinpath(cfgdir, "star_prior.jld")),
                    JLD.load(joinpath(cfgdir, "gal_prior.jld"))]
    for i in 1:2
        k[:, i] = prior_params[i]["c_weights"]
        color_mean[:, :, i] = prior_params[i]["c_means"]
        color_cov[:, :, :, i] = prior_params[i]["c_covs"]
    end

    # log normal parameters for the r-band brightness prior.
    # these were fit by maximum likelihood to the output of primary
    # on one field.
    flux_mean = Float64[1.5035546, 1.07431]
    flux_var = Float64[1.9039063^2, 1.1177502^2]

    # If I remove this next statement, compile time for
    # benchmark_infer.jl jumps from 13 seconds to 300 seconds!
    # Really. It's crazy!
    log(42.)

    # Compile time is still over 20x as long if `log(42.)`
    # is replaced with `exp(2.)`. Insane! Julia 0.5.0 mac os x.
    # Update: happens on Linux with 0.5.0 binaries too.
    # exp(2.)

    # log normal prior parameters (location, scale) on galaxy scale.
    # determined by fitting a univariate log normal to primary's
    # output the region of stripe 82 we use for validation
    gal_scale_mean = 0.5015693
    gal_scale_var = 0.8590007^2

    PriorParams(is_star, flux_mean, flux_var, k, color_mean, color_cov, gal_scale_mean, gal_scale_var)
end

const prior = load_prior()


text:
using MatpowerCases, Base.Test

println("- extract_gencost check")
for c in casenames()
    if c != "case4gs" # case4gs has no gencost data
        mpc = loadcase(c; describe=false)
        @assert isa(extract_gencost(mpc), MatpowerCases.Gencost)
        @assert isa(extract_gencost(c), MatpowerCases.Gencost)
    end
end

ext: jl


text:
if nworkers()<=1
	addprocs(2);
end

using  jInv.Mesh
using  jInv.Utils
using  jInv.InverseSolve
using  EikonalInv
using  MAT

#############################################################################################################
modelDir = "../examples";

dataDir = pwd();
include("../drivers/prepareTravelTimeDataFiles.jl");
include("../drivers/readModelAndGenerateMeshMref.jl");
include("../drivers/setupTravelTimeTomography.jl");

########################################################################################################
######################################## for SEG 256X512 ###############################################
#######################################################################################################
dim     = 2;
pad     = 1;
jump    = 32;
offset  = 256;
(m,Minv,mref,boundsHigh,boundsLow) = readModelAndGenerateMeshMref(modelDir,"SEGmodel2Dsalt.dat",dim,pad,[0.0,13.5,0.0,4.2],[64,32]);
dataFilenamePrefix = string(dataDir,"/DATA_SEG",tuple((Minv.n+1)...));

prepareTravelTimeDataFiles(m,Minv,mref,boundsHigh,boundsLow,dataFilenamePrefix,pad,jump,offset);
resultsFilenamePrefix = "";
(Q,P,pMisRFs,SourcesSubInd,contDiv,Iact,sback,mref,boundsHigh,boundsLow,resultsFilename) = setupTravelTimeTomography(m,dataFilenamePrefix, resultsFilenamePrefix);

########################################################################################################
##### Set up Inversion #################################################################################
########################################################################################################
maxStep=0.1*maximum(boundsHigh);
a = minimum(boundsLow);
b = maximum(boundsHigh);
function modfun(m)
	bm,dbm = getBoundModel(m,a,b);
	bs,dbs = velocityToSlowSquared(bm);
	dsdm = dbm*dbs;
	return bs,dsdm;
end
mref = getBoundModelInv(mref,a,b);
boundsHigh = boundsHigh*10000000.0;
boundsLow = -boundsLow*100000000.0

cgit  = 2;
maxit = 1;
alpha = 1e+0;
pcgTol = 1e-1;

HesPrec=getSSORCGRegularizationPreconditioner(1.0,1e-5,2)

################################################# GIT VERSION OF JINV #################################################

regparams = [1.0,1.0,1.0,1e-2];
regfun(m, mref, M) = wdiffusionRegNodal(m, mref, M, Iact=Iact, C = regparams);

pInv = getInverseParam(Minv,modfun,regfun,alpha,mref[:],boundsLow,boundsHigh,
                     maxStep=maxStep,pcgMaxIter=cgit,pcgTol=pcgTol,
					 minUpdate=1e-3, maxIter = maxit,HesPrec=HesPrec);

mc,Dc,flag = projGNCG(mref[:],pInv,pMisRFs,indCredit = []);

clear!(pMisRFs);
#############################################################################################
rm("DATA_SEG(66, 33)_travelTime.dat");
rm("DATA_SEG(66, 33)_rcvMap.dat");
rm("DATA_SEG(66, 33)_srcMap.dat");
rm("DATA_SEG(66, 33)_PARAM.mat");
rm("jInv.out");


text:
# Fusion trees:
#==============================================================================#
"""
    struct FusionTree{G,N,M,L,T}

Represents a fusion tree of sectors of type `G<:Sector`, fusing (or splitting) `N` uncoupled
sectors to a coupled sector. This fusion tree has `M=max(0,N-2)` inner lines. Furthermore,
for `FusionStyle(G) isa DegenerateNonAbelian`, the `L=max(0,N-1)` corresponding vertices
carry a label of type `T`. If `FusionStyle(G) isa Union{Abelian,SimpleNonAbelian}`,
`T = Nothing`.
"""
struct FusionTree{G<:Sector,N,M,L,T}
    uncoupled::NTuple{N,G}
    coupled::G
    isdual::NTuple{N,Bool}
    innerlines::NTuple{M,G} # M = N-2
    vertices::NTuple{L,T} # L = N-1
    function FusionTree{G,N,M,L,T}(uncoupled::NTuple{N,G},
                                            coupled::G,
                                            isdual::NTuple{N,Bool},
                                            innerlines::NTuple{M,G},
                                            vertices::NTuple{L,T}) where
                                            {G<:Sector,N,M,L,T}
        new{G,N,M,L,T}(uncoupled, coupled, isdual, innerlines, vertices)
    end
end
FusionTree{G}(uncoupled::NTuple{N,Any},
                coupled,
                isdual::NTuple{N,Bool},
                innerlines,
                vertices = ntuple(n->nothing, StaticLength(N)-StaticLength(1))
                ) where {G<:Sector,N} =
    fusiontreetype(G, StaticLength(N))(map(s->convert(G,s),uncoupled),
        convert(G,coupled), isdual, map(s->convert(G,s), innerlines), vertices)
FusionTree(uncoupled::NTuple{N,G},
            coupled::G,
            isdual::NTuple{N,Bool},
            innerlines,
            vertices = ntuple(n->nothing, StaticLength(N)-StaticLength(1))
            ) where {G<:Sector,N} =
    fusiontreetype(G, StaticLength(N))(uncoupled, coupled, isdual, innerlines, vertices)

function FusionTree{G}(uncoupled::NTuple{N}, coupled = one(G),
                        isdual = ntuple(n->false, StaticLength(N))) where {G<:Sector, N}
    FusionStyle(G) isa Abelian ||
        error("fusion tree requires inner lines if `FusionStyle(G) <: NonAbelian`")
    FusionTree{G}(map(s->convert(G,s), uncoupled), convert(G, coupled), isdual,
                    _abelianinner(map(s->convert(G,s),(uncoupled..., dual(coupled)))))
end
function FusionTree(uncoupled::NTuple{N,G}, coupled::G = one(G),
                        isdual = ntuple(n->false, StaticLength(N))) where {G<:Sector, N}
    FusionStyle(G) isa Abelian ||
        error("fusion tree requires inner lines if `FusionStyle(G) <: NonAbelian`")
    FusionTree{G}(uncoupled, coupled, isdual, _abelianinner((uncoupled..., dual(coupled))))
end

# Properties
sectortype(::Type{<:FusionTree{G}}) where {G<:Sector} = G
FusionStyle(::Type{<:FusionTree{G}}) where {G<:Sector} = FusionStyle(G)
BraidingStyle(::Type{<:FusionTree{G}}) where {G<:Sector} = BraidingStyle(G)
Base.length(::Type{<:FusionTree{<:Sector,N}}) where {N} = N

sectortype(f::FusionTree) = sectortype(typeof(f))
FusionStyle(f::FusionTree) = FusionStyle(typeof(f))
BraidingStyle(f::FusionTree) = BraidingStyle(typeof(f))
Base.length(f::FusionTree) = length(typeof(f))

# Hashing, important for using fusion trees as key in a dictionary
function Base.hash(f::FusionTree{G}, h::UInt) where {G}
    h = hash(f.isdual, hash(f.coupled, hash(f.uncoupled, h)))
    if FusionStyle(G) isa SimpleNonAbelian
        h = hash(f.innerlines, h)
    end
    if FusionStyle(G) isa DegenerateNonAbelian
        h = hash(f.vertices, h)
    end
    return h
end
function Base.isequal(f1::FusionTree{G,N}, f2::FusionTree{G,N}) where {G<:Sector,N}
    f1.coupled == f2.coupled || return false
    @inbounds for i = 1:N
        f1.uncoupled[i] == f2.uncoupled[i] || return false
        f1.isdual[i] == f2.isdual[i] || return false
    end
    if FusionStyle(G) isa SimpleNonAbelian
        @inbounds for i=1:N-2
            f1.innerlines[i] == f2.innerlines[i] || return false
        end
    end
    if FusionStyle(G) isa DegenerateNonAbelian
        @inbounds for i=1:N-1
            f1.vertices[i] == f2.vertices[i] || return false
        end
    end
    return true
end
Base.isequal(f1::FusionTree, f2::FusionTree) = false

# Facilitate getting correct fusion tree types
Base.@pure fusiontreetype(::Type{G}, ::StaticLength{0}) where {G<:Sector} =
    FusionTree{G, 0, 0, 0, vertex_labeltype(G)}
Base.@pure fusiontreetype(::Type{G}, ::StaticLength{1}) where {G<:Sector} =
    FusionTree{G, 1, 0, 0, vertex_labeltype(G)}
Base.@pure fusiontreetype(::Type{G}, ::StaticLength{2}) where {G<:Sector} =
    FusionTree{G, 2, 0, 1, vertex_labeltype(G)}
Base.@pure fusiontreetype(::Type{G}, ::StaticLength{N}) where {G<:Sector, N} =
    _fusiontreetype(G, StaticLength(N),
        StaticLength(N) - StaticLength(2), StaticLength(N) - StaticLength(1))
Base.@pure _fusiontreetype(::Type{G}, ::StaticLength{N}, ::StaticLength{M},
                            ::StaticLength{L}) where {G<:Sector, N, M, L} =
    FusionTree{G,N,M,L,vertex_labeltype(G)}

# converting to actual array
function Base.convert(::Type{Array}, f::FusionTree{G, 0}) where {G}
    T = eltype(fusiontensor(one(G), one(G), one(G)))
    return fill(one(T), 1)
end
function Base.convert(::Type{Array}, f::FusionTree{G, 1}) where {G}
    c = f.coupled
    dc = dim(c)
    if f.isdual[1]
        Zcbartranspose = sqrt(dc)*reshape(fusiontensor(conj(c), c, one(c)), (dc, dc))
        return convert(Array, conj(Zcbartranspose))
    else
        convert(Array, reshape(fusiontensor(c, one(c), c), (dc, dc)))
    end
end

function Base.convert(::Type{Array}, f::FusionTree{G,2}) where {G}
    a, b = f.uncoupled
    isduala, isdualb = f.isdual
    c = f.coupled
    da, db, dc = dim.((a,b,c))
    μ = f.vertices[1]
    X = reshape(fusiontensor(a, b, c, μ), da*db, dc)
    Za = convert(Array, FusionTree((a,), a, (isduala,), ()))
    Zb = convert(Array, FusionTree((b,), b, (isdualb,), ()))
    return convert(Array, reshape(kron(Zb, Za)*X, (da, db, dc)))
end

function Base.convert(::Type{Array}, f::FusionTree{G}) where {G}
    tailout = (f.innerlines[1], TupleTools.tail2(f.uncoupled)...)
    isdualout = (false, TupleTools.tail2(f.isdual)...)
    ftail = FusionTree(tailout, f.coupled, isdualout,
                        Base.tail(f.innerlines), Base.tail(f.vertices))
    Ctail = convert(Array, ftail)
    f1 = FusionTree((f.uncoupled[1], f.uncoupled[2]), f.innerlines[1],
                    (f.isdual[1], f.isdual[2]), (), (f.vertices[1],))
    C1 = convert(Array, f1)
    dtail = size(Ctail)
    d1 = size(C1)
    C = reshape(C1, d1[1]*d1[2], d1[3]) *
            reshape(Ctail, dtail[1], prod(Base.tail(dtail)))
    return reshape(C, (d1[1], d1[2], Base.tail(dtail)...))
end

# Show methods
function Base.show(io::IO, t::FusionTree{G,N,M,K,Nothing}) where {G<:Sector,N,M,K}
    print(IOContext(io, :typeinfo => G), "FusionTree{", G, "}(",
        t.uncoupled, ", ", t.coupled, ", ", t.isdual, ", ", t.innerlines, ")")
end
function Base.show(io::IO, t::FusionTree{G}) where {G<:Sector}
    print(IOContext(io, :typeinfo => G), "FusionTree{", G, "}(",
        t.uncoupled, ", ", t.coupled, ", ", t.isdual, ",",
        t.innerlines, ", ", t.vertices, ")")
end

# Manipulate fusion trees
include("manipulations.jl")

# Fusion tree iterators
include("iterator.jl")

# auxiliary routines
# _abelianinner: generate the inner indices for given outer indices in the abelian case
_abelianinner(outer::Tuple{}) = ()
_abelianinner(outer::Tuple{G}) where {G<:Sector} =
    outer[1] == one(G) ? () : throw(SectorMismatch())
_abelianinner(outer::Tuple{G,G}) where {G<:Sector} =
    outer[1] == dual(outer[2]) ? () : throw(SectorMismatch())
_abelianinner(outer::Tuple{G,G,G}) where {G<:Sector} =
    first(⊗(outer...)) == one(G) ? () : throw(SectorMismatch())
function _abelianinner(outer::NTuple{N,G}) where {G<:Sector,N}
    c = first(outer[1] ⊗ outer[2])
    return (c, _abelianinner((c, TupleTools.tail2(outer)...))...)
end


text:
"""
    Deflate(level=5)

Deflate/ZLIB lossless compression filter. `level` is an integer between 0 and 9,
inclusive, denoting the compression level, with 0 being no compression, 9 being the
highest compression (but slowest speed).

# External links
- $(h5doc("H5P_SET_DEFLATE"))
- [_Deflate_ on Wikipedia](https://en.wikipedia.org/wiki/Deflate)
"""
struct Deflate <: Filter
    level::Cuint
end
Deflate(;level=5) = Deflate(level)
Base.show(io::IO, deflate::Deflate) = print(io, Deflate, "(level=", Int(deflate.level), ")")

filterid(::Type{Deflate}) = API.H5Z_FILTER_DEFLATE
FILTERS[API.H5Z_FILTER_DEFLATE] = Deflate

function Base.push!(f::FilterPipeline, deflate::Deflate)
    API.h5p_set_deflate(f.plist, deflate.level)
    return f
end

"""
    Shuffle()

The shuffle filter de-interlaces a block of data by reordering the bytes. All the bytes
from one consistent byte position of each data element are placed together in one block;
all bytes from a second consistent byte position of each data element are placed together
a second block; etc. For example, given three data elements of a 4-byte datatype stored as
012301230123, shuffling will re-order data as 000111222333. This can be a valuable step in
an effective compression algorithm because the bytes in each byte position are often
closely related to each other and putting them together can increase the compression
ratio.

As implied above, the primary value of the shuffle filter lies in its coordinated use with
a compression filter; it does not provide data compression when used alone. When the
shuffle filter is applied to a dataset immediately prior to the use of a compression
filter, the compression ratio achieved is often superior to that achieved by the use of a
compression filter without the shuffle filter.

# External links
- $(h5doc("H5P_SET_SHUFFLE"))
"""
struct Shuffle <: Filter
end
filterid(::Type{Shuffle}) = API.H5Z_FILTER_SHUFFLE
FILTERS[API.H5Z_FILTER_SHUFFLE] = Shuffle

function Base.push!(f::FilterPipeline, ::Shuffle)
    API.h5p_set_shuffle(f.plist)
    return f
end

"""
    Fletcher32()

The Fletcher32 checksum filter. This doesn't perform compression, but instead checks the validity of the stored data.

This should be applied _after_ any lossy filters have been applied.

# External links
- $(h5doc("H5P_SET_FLETCHER32"))
- [_Fletcher's checksum_ on Wikipedia](https://en.wikipedia.org/wiki/Fletcher's_checksum)
"""
struct Fletcher32 <: Filter
end
filterid(::Type{Fletcher32}) = API.H5Z_FILTER_FLETCHER32
FILTERS[API.H5Z_FILTER_FLETCHER32] = Fletcher32
function Base.push!(f::FilterPipeline, ::Fletcher32)
    API.h5p_set_fletcher32(f.plist)
    return f
end

"""
    Szip(coding=:nn, pixels_per_block=8)

Szip compression lossless filter. Options:

- `coding`: the coding method: either `:ec` (entropy coding) or `:nn` (nearest neighbors,
  default)
- `pixels_per_block`: The number of pixels or data elements in each data block (typically
  8, 10, 16, or 32)

# External links
- $(h5doc("H5P_SET_SZIP"))
- [Szip Compression in HDF Products](https://support.hdfgroup.org/doc_resource/SZIP/)
"""
struct Szip <: Filter
    options_mask::Cuint
    pixels_per_block::Cuint
end
function Szip(;coding=:nn, pixels_per_block=8)
    options_mask = Cuint(0)
    if coding == :ec
        options_mask |= API.H5_SZIP_EC_OPTION_MASK
    elseif coding == :nn
        options_mask |= API.H5_SZIP_NN_OPTION_MASK
    else
        error("invalid coding option")
    end
    Szip(options_mask, pixels_per_block)
end
function Base.show(io::IO, szip::Szip)
    print(io, Szip, "(")
    if szip.options_mask & API.H5_SZIP_EC_OPTION_MASK != 0
        print(io, "coding=:ec,")
    elseif szip.options_mask & API.H5_SZIP_NN_OPTION_MASK != 0
        print(io, "coding=:nn,")
    end
    print(io, "pixels_per_block=", Int(szip.pixels_per_block))
end

filterid(::Type{Szip}) = API.H5Z_FILTER_SZIP
FILTERS[API.H5Z_FILTER_SZIP] = Szip
function Base.push!(f::FilterPipeline, szip::Szip)
    API.h5p_set_szip(f.plist, szip.options_mask, szip.pixels_per_block)
    return f
end

"""
    NBit()

The N-Bit filter.

# External links
- $(h5doc("H5P_SET_NBIT"))
"""
struct NBit <: Filter
end
filterid(::Type{NBit}) = API.H5Z_FILTER_NBIT
FILTERS[API.H5Z_FILTER_NBIT] = NBit
function Base.push!(f::FilterPipeline, ::NBit)
    API.h5p_set_nbit(f.plist)
    return f
end

"""
    ScaleOffset(scale_type::Integer, scale_offset::Integer)

The scale-offset filter.

# External links
- $(h5doc("H5P_SET_SCALEOFFSET"))
"""
struct ScaleOffset <: Filter
    scale_type::Cint
    scale_factor::Cint
end

filterid(::Type{ScaleOffset}) = API.H5Z_FILTER_SCALEOFFSET
FILTERS[API.H5Z_FILTER_SCALEOFFSET] = ScaleOffset
function Base.push!(f::FilterPipeline, scaleoffset::ScaleOffset)
    API.h5p_set_scaleoffset(f.plist, scaleoffset.scale_type, scaleoffset.scale_factor)
    return f
end


text:
push!(LOAD_PATH,"../src/")
using WebAssemblyText
using Documenter

makedocs(;
    modules=[WebAssemblyText],
    authors="Anders Gustafsson <andersgee@gmail.com> and contributors",
    repo="https://github.com/andersgee/WebAssemblyText.jl/blob/{commit}{path}#L{line}",
    sitename="WebAssemblyText.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://andersgee.github.io/WebAssemblyText.jl",
        assets=["assets/favicon.ico"],
    ),
    pages = [
        "Home"=>"index.md",
        "Internals" => "Internals.md",
    ]
)

deploydocs(;
    repo="github.com/Andersgee/WebAssemblyText.jl.git",
    target = "build",
    branch = "gh-pages",
)


text:
#
# This module is a RPC remote API. It contains Manager RPC
# functions identifiers. They are used for Remote Procedure Calls
# from any other processes. See Connection.Command type for details.
# It'salso important, that this module should contain only native
# Julia definition. It shouldn't contain other modules.
#
# @author DeadbraiN
#
module RpcApi
  export BLOCK_SIZE

  export RPC_GET_REGION
  export RPC_CREATE_ORGANISMS
  export RPC_CREATE_ORGANISM
  export RPC_SET_CONFIG
  export RPC_GET_CONFIG
  export RPC_SET_QUIET
  export RPC_MUTATE
  export RPC_GET_IPS
  export RPC_GET_ORGANISM
  export RPC_GET_AMOUNT
  export RPC_GET_ORGANISMS
  export RPC_SET_ENERGY
  export RPC_SET_ENERGY_RND
  export RPC_BACKUP
  export RPC_GET_STATISTICS
  export RPC_GET_BEST
  export RPC_SET_LEFT_WORLD
  export RPC_SET_RIGHT_WORLD
  export RPC_SET_UP_WORLD
  export RPC_SET_DOWN_WORLD
  export RPC_SET_WORLD_STREAMING
  export RPC_MARK_ORGANISM

  export RPC_ORG_STEP_LEFT
  export RPC_ORG_STEP_RIGHT
  export RPC_ORG_STEP_UP
  export RPC_ORG_STEP_DOWN

  export RPC_ORG_STEP_OK
  export RPC_ORG_STEP_FAIL

  export RPC_WORLD_CHANGE

  export Region
  export SimpleOrganism
  export Statistics
  export Org
  export Block
  #
  # Describes one organism's dot in a world
  #
  type Org
    offs::UInt8
    id::UInt32
    color::UInt16 # only 65536 colors are supported now
    age::UInt16
  end
  #
  # Describes one block of pixels of size 16x16
  #
  type Block
    energy::Array{UInt8, 1} # array of energy dots offsets within current 16x16 block
    orgs::Array{Org, 1}     # array of organisms within current 16x16 block
  end
  #
  # Describes pixels region in a world. This type will be transfered
  # to the client. Maximum supported sesolution 4096x4096
  #
  type Region
    xBlocks::UInt8
    yBlocks::UInt8
    width::UInt16
    height::UInt16
    blocks::Array{Block, 1}
    ips::Float64
  end
  #
  # This is an analog of Creature.Organism type. It's used for
  # transporting more simple version of organism data through
  # network.
  #
  type SimpleOrganism
    #
    # Unique organism id
    #
    id::UInt
    #
    # @inheritable
    # Code of organism. Array of atomic string blocks on Julia language.
    #
    code::Expr
    #
    # @inheritable
    # Mutations probability. Add, change, delete. e.g.: [1,10,2]
    # means, that "add" mutation will be 10 times rare then "change"
    # and 2 times rare then "delete" mutations.
    #
    mutationProbabilities::Array{Int, 1}
    #
    # @inheritable
    # Percent from code size of mutations, which will be applied to arganism after
    # clonning.
    #
    mutationsOnClonePercent::Float64
    #
    # @inheritable
    # Amount of iterations within organism's life loop, after that we
    # do mutations according to MUTATE_AMOUNT config amount. If 0, then
    # mutations will be disabled.
    #
    mutationPeriod::Int
    #
    # @inheritable
    # Value, which will be used like amount of mutations per
    # MUTATE_AFTER_TIMES iterations. 0 is a possible value if
    # we want to disable mutations. Percent means "percent from
    # code size".
    #
    mutationPercent::Float64
    #
    # Organism's energy. If it's zero, then organism is die.
    #
    energy::Int
    #
    # Amount of mutations starting from first ancestor
    #
    mutationsFromStart::Int
    #
    # @inheritable
    # Color of organism. Similar colors means relative organisms.
    # Example: 0xAABB
    #
    color::UInt16
    #
    # Organism's personal memory. Is used in any possible way.
    #
    mem::Dict{Float64, Float64}
    #
    # Organism's position in a 2D world. Starts from (1,1)
    # ends with (worldWidth, worldHeight) configurations.
    #
    pos::Array{Int, 1}
    #
    # Age of organism in iterations
    #
    age::Int
    #
    # Percent of energy, which is given to clonned organism.
    #
    cloneEnergyPercent::Float64
    #
    # Organism code size in lines
    #
    codeSize::Int
  end
  #
  # Structure of world statistics
  #
  type Statistics
    #
    # Amount of alive organisms
    #
    orgAmount::Int
    #
    # Iterations(cicles) per second
    #
    ips::Float64
    #
    # Total amount of organisms: alive+dead
    #
    orgTotalAmount::Int
    #
    # Width of the world
    #
    worldWidth::Int
    #
    # Height of the world
    #
    worldHeight::Int
    #
    # Global configuration in string representation
    #
    cfg::Array{String, 1}
    #
    # An organism with minimum energy
    #
    minOrg::SimpleOrganism
    #
    # An organism with maximum energy
    #
    maxOrg::SimpleOrganism
  end
  #
  # The size of one 16x16 block
  #
  const BLOCK_SIZE = 16
  #
  # RPC API unique identifiers. Only these functions may be called
  # remotely on the server.
  #
  const RPC_GET_REGION          = 1
  const RPC_CREATE_ORGANISMS    = 2
  const RPC_CREATE_ORGANISM     = 3
  const RPC_SET_CONFIG          = 4
  const RPC_GET_CONFIG          = 5
  const RPC_SET_QUIET           = 6
  const RPC_MUTATE              = 7
  const RPC_GET_IPS             = 8
  const RPC_GET_ORGANISM        = 9
  const RPC_GET_AMOUNT          = 10
  const RPC_GET_ORGANISMS       = 11
  const RPC_SET_ENERGY          = 12
  const RPC_SET_ENERGY_RND      = 13
  const RPC_BACKUP              = 14
  const RPC_GET_STATISTICS      = 15
  const RPC_GET_BEST            = 16
  const RPC_SET_LEFT_WORLD      = 17
  const RPC_SET_RIGHT_WORLD     = 18
  const RPC_SET_UP_WORLD        = 19
  const RPC_SET_DOWN_WORLD      = 20
  const RPC_SET_WORLD_STREAMING = 21
  const RPC_MARK_ORGANISM       = 22
  #
  # These IDs are for inter-instances communication. This is
  # something like internal protocol between Managers/servers.

  const RPC_ORG_STEP_LEFT       = 23
  const RPC_ORG_STEP_RIGHT      = 24
  const RPC_ORG_STEP_UP         = 25
  const RPC_ORG_STEP_DOWN       = 26
  #
  # Means that step from one instance to another has finished
  #
  const RPC_ORG_STEP_OK         = 27
  const RPC_ORG_STEP_FAIL       = 28
  #
  # These id's for request from server to client
  #
  const RPC_WORLD_CHANGE        = 29
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
using Restful
using Test
using HTTP
using JSON2
using Logging

import Restful: json

const db = Dict()
const app = Restful.app()

app.get("/") do req, res, route
    res.text("Hello World?")
end

app.get("/todo", json) do req, res, route
    res.json(keys(db) |> collect)
end

app.post("/todo", json) do req, res, route
    id = req.body |> hash |> string
    db[id] = req.body
    res.json(;id=id)
end

app.get("/todo/:id", json) do req, res, route
    route.id in keys(db) ? res.json(;content=db[route.id]) : res.code(404)
end

app.put("/todo/:id", json) do req, res, route
    db[route.id] = req.json.content
    res.code(200)
end

app.delete("/todo/:id", json) do req, res, route
    if route.id in keys(db)
        delete!(db, route.id)
        res.code(200)
    else
        res.code(404)
    end
end

@async with_logger(SimpleLogger(stderr, Logging.Warn)) do
    app.listen("127.0.0.1", 3001)
end

macro u_str(x)
    "http://127.0.0.1:3001$x"
end

@test String(HTTP.get(u"/").body) == "Hello World?"
@test HTTP.put(u"/todo/39", [], JSON2.write((content="task a",))).status == 200
@test String(HTTP.get(u"/todo/39").body) == "{\"content\":\"task a\"}"
@test String(HTTP.get(u"/todo").body) == "[\"39\"]"
@test HTTP.delete(u"/todo/39").status == 200
@test HTTP.get(u"/todo/39", status_exception=false).status == 404
@test String(HTTP.get(u"/todo").body) == "[]"
@test JSON2.read(String(HTTP.post(u"/todo", [], JSON2.write((content="task b",))).body)).id ==
      JSON2.read(String(HTTP.get(u"/todo").body))[1]
@test HTTP.put(u"/", status_exception=false).status == 405
# @test headers(put(url("/todo")))["Allow"] == "GET, POST"

app.close()


text:
using InvariantMeasures
using ValidatedNumerics

# Testing whether the "shooting method" beats taking successive preimages

# f = x -> 2*x + 0.1*InvariantMeasures.sinpi(2*x)
# preim1 = root(x ->f(x)-1., 0..4, 0)

f = x-> 2.5x*(1-x) + 0.01*x*x + 0.3*x*x*x + 0.003*x*x*x*x
preim1 = @interval(0.5)

mypreim(y) = root(x -> f(x)-y, hull(0..0,preim1), 0)

n = 15

y0 = @interval(0.21)

# method 1
ys = fill(y0, n+1)
ys[1] = y0
for i = 1:n
    ys[i+1] = mypreim(ys[i])
end

# method 2

zs = fill(hull(0..0,preim1), n)
fs = fill(f, n)

nthpreimage!(y0, fs, zs)

z = zs[1]

[zs[end:-1:1] ys[2:end]]
[diam.(zs[end:-1:1]) diam.(ys[2:end])]

diam.(zs[end:-1:1]) ./ diam.(ys[2:end])


text:
# Before running this, please make sure to activate and instantiate the environment
# corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Manifest.toml)
# so that you get an environment which matches the one used to generate the tutorials:
#
# ```julia
# cd("DataScienceTutorials") # cd to folder with the *.toml
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

using MLJ
using StableRNGs
import DataFrames

@load RidgeRegressor pkg=MultivariateStats

rng = StableRNG(6616) # for reproducibility
x1 = rand(rng, 300)
x2 = rand(rng, 300)
x3 = rand(rng, 300)
y = exp.(x1 - x2 -2x3 + 0.1*rand(rng, 300))
X = DataFrames.DataFrame(x1=x1, x2=x2, x3=x3)

test, train = partition(eachindex(y), 0.8);

Xs = source(X)
ys = source(y)

std_model = Standardizer()
stand = machine(std_model, Xs)
W = MLJ.transform(stand, Xs)

box_model = UnivariateBoxCoxTransformer()
box_mach = machine(box_model, ys)
z = MLJ.transform(box_mach, ys)

ridge_model = RidgeRegressor(lambda=0.1)
ridge = machine(ridge_model, W, z)
ẑ = predict(ridge, W)

ŷ = inverse_transform(box_mach, ẑ)

surrogate = Deterministic()
mach = machine(surrogate, Xs, ys; predict=ŷ)

fit!(mach)
predict(mach, X[test[1:5], :])

@from_network mach begin
    mutable struct CompositeModel
        regressor=ridge_model
    end
end

cm = machine(CompositeModel(), X, y)
res = evaluate!(cm, resampling=Holdout(fraction_train=0.8, rng=51),
                measure=rms)
round(res.measurement[1], sigdigits=3)

mutable struct CompositeModel2 <: DeterministicNetwork
    std_model::Standardizer
    box_model::UnivariateBoxCoxTransformer
    ridge_model::RidgeRegressor
end

function MLJ.fit(m::CompositeModel2, verbosity::Int, X, y)
    Xs = source(X)
    ys = source(y)
    W = MLJ.transform(machine(m.std_model, Xs), Xs)
    box = machine(m.box_model, ys)
    z = MLJ.transform(box, ys)
    ẑ = predict(machine(m.ridge_model, W, z), W)
    ŷ = inverse_transform(box, ẑ)
    mach = machine(Deterministic(), Xs, ys; predict=ŷ)
    fit!(mach, verbosity=verbosity - 1)
    return mach()
end

mdl = CompositeModel2(Standardizer(), UnivariateBoxCoxTransformer(),
                      RidgeRegressor(lambda=0.1))
cm = machine(mdl, X, y)
res = evaluate!(cm, resampling=Holdout(fraction_train=0.8), measure=rms)
round(res.measurement[1], sigdigits=3)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
"""
Base class for archive-specific component
of the `OptimizationResults`.
"""
abstract type ArchiveOutput end

"""
Base class for method-specific component
of the `OptimizationResults`.
"""
abstract type MethodOutput end

struct DummyMethodOutput <: MethodOutput end

# no method-specific output by default
MethodOutput(method::Optimizer) = DummyMethodOutput()

"""
The results of running optimization method.

Returned by `run!(oc::OptRunController)`.
Should be compatible (on the API level) with the `Optim` package.
See `make_opt_results()`.
"""
mutable struct OptimizationResults
    method::String           # FIXME symbol instead or flexible?
    stop_reason::String      # FIXME turn into type hierarchy of immutable reasons with their attached info
    iterations::Int
    start_time::Float64           # time (seconds) optimization started
    elapsed_time::Float64         # time (seconds) optimization finished
    parameters::Parameters        # all user-specified parameters to bboptimize()
    f_calls::Int                  # total number of fitness function evaluations
    fit_scheme::FitnessScheme     # fitness scheme used by the archive
    archive_output::ArchiveOutput # archive-specific output
    method_output::MethodOutput   # method-specific output

    function OptimizationResults(ctrl, oc)
        new(
            string(oc.parameters[:Method]),
            stop_reason(ctrl),
            num_steps(ctrl),
            start_time(ctrl), elapsed_time(ctrl),
            oc.parameters,
            num_func_evals(ctrl),
            fitness_scheme(evaluator(ctrl).archive),
            ArchiveOutput(evaluator(ctrl).archive),
            MethodOutput(ctrl.optimizer))
    end
end

stop_reason(or::OptimizationResults) = or.stop_reason
iterations(or::OptimizationResults) = or.iterations
start_time(or::OptimizationResults) = or.start_time
elapsed_time(or::OptimizationResults) = or.elapsed_time
parameters(or::OptimizationResults) = or.parameters
f_calls(or::OptimizationResults) = or.f_calls

fitness_scheme(or::OptimizationResults) = or.fit_scheme
best_candidate(or::OptimizationResults) = or.archive_output.best_candidate
best_fitness(or::OptimizationResults) = or.archive_output.best_fitness
# FIXME doesn't work if there's no best candidate
numdims(or::OptimizationResults) = length(best_candidate(or))

function general_stop_reason(or::OptimizationResults)
    detailed_reason = stop_reason(or)

    if occursin(r"Fitness .* within tolerance .* of optimum", detailed_reason)
        return "Within fitness tolerance of optimum"
    end

    if occursin(r"Delta fitness .* below tolerance .*", detailed_reason)
        return "Delta fitness below tolerance"
    end

    return detailed_reason
end

# Alternative nomenclature that mimics Optim.jl more closely.
# FIXME should be it be enabled only for MinimizingFitnessScheme?
Base.minimum(or::OptimizationResults) = best_candidate(or)
f_minimum(or::OptimizationResults) = best_fitness(or)
# FIXME lookup stop_reason
iteration_converged(or::OptimizationResults) = iterations(or) >= parameters(or)[:MaxSteps]
isinterrupted(or::OptimizationResults) = stop_reason(or) == "InterruptException()"

"""
`TopListArchive`-specific components of the optimization results.
"""
struct TopListArchiveOutput{F,C} <: ArchiveOutput
    best_fitness::F
    best_candidate::C

    TopListArchiveOutput(archive::TopListArchive{F}) where F =
        new{F,Individual}(best_fitness(archive), best_candidate(archive))
end

ArchiveOutput(archive::TopListArchive) = TopListArchiveOutput(archive)

"""
Wrapper for `FrontierIndividual` that allows easy access to the problem fitness.
"""
struct FrontierIndividualWrapper{F,FA} <: FitIndividual{F}
    inner::FrontierIndividual{FA}
    fitness::F

    FrontierIndividualWrapper{F,FA}(
            indi::FrontierIndividual{FA},
            fit_scheme::FitnessScheme{FA}) where {F, FA} =
        new{F, FA}(indi, convert(F, fitness(indi), fit_scheme))
end

params(indi::FrontierIndividualWrapper) = params(indi.inner)
archived_fitness(indi::FrontierIndividualWrapper) = fitness(indi.inner)

"""
`EpsBoxArchive`-specific components of the optimization results.
"""
struct EpsBoxArchiveOutput{N,F,FS<:EpsBoxDominanceFitnessScheme} <: ArchiveOutput
    best_fitness::NTuple{N,F}
    best_candidate::Individual
    frontier::Vector{FrontierIndividualWrapper{NTuple{N,F},IndexedTupleFitness{N,F}}} # inferred Pareto frontier
    fit_scheme::FS

    function EpsBoxArchiveOutput(archive::EpsBoxArchive{N,F}) where {N,F}
        fit_scheme = fitness_scheme(archive)
        FIW = FrontierIndividualWrapper{NTuple{N,F}, IndexedTupleFitness{N,F}}
        new{N,F,typeof(fit_scheme)}(convert(NTuple{N,F}, best_fitness(archive), fit_scheme),
                                    best_candidate(archive),
                                    FIW[FIW(frontel, fit_scheme) for frontel in pareto_frontier(archive)],
                                    fit_scheme)
    end
end

ArchiveOutput(archive::EpsBoxArchive) = EpsBoxArchiveOutput(archive)

pareto_frontier(or::OptimizationResults) = or.archive_output.frontier

"""
`PopulationOptimizer`-specific components of the `OptimizationResults`.
Stores the final population.
"""
struct PopulationOptimizerOutput{P} <: MethodOutput
    population::P

    # FIXME PO is only required so that julia doesn't think it's an old inner ctor syntax
    PopulationOptimizerOutput(method::PO) where {PO<:PopulationOptimizer} =
        new{typeof(population(method))}(population(method))
end

MethodOutput(optimizer::PopulationOptimizer) = PopulationOptimizerOutput(optimizer)

population(or::OptimizationResults) = or.method_output.population


text:
module OnlineStatsBase

using Statistics, Dates, OrderedCollections, LinearAlgebra

import LearnBase: nobs, value, fit!
import StatsBase: StatsBase

export
    OnlineStat, Weight,
    nobs, value, fit!, eachrow, eachcol,
    EqualWeight, ExponentialWeight, LearningRate, LearningRate2, HarmonicWeight, McclainWeight,
    # Stats
    Counter, CountMap, CovMatrix, Extrema, FTSeries, Group, GroupBy, Mean, Moments, Series, Sum, Variance

#-----------------------------------------------------------------------# OnlineStat
abstract type OnlineStat{T} end

# Stats that hold a collection of other stats
abstract type StatCollection{T} <: OnlineStat{T} end
function Base.show(io::IO, o::StatCollection)
    print(io, name(o, false, false))
    print_stat_tree(io, o.stats)
end
function print_stat_tree(io::IO, stats)
    for (i, stat) in enumerate(stats)
        char = i == length(stats) ? '└' : '├'
        print(io, "\n  $(char)── $stat")
    end
end

nobs(o::OnlineStat) = o.n

"""
    value(stat::OnlineStat)

Calculate the value of `stat` from its "sufficient statistics".
"""
@generated function value(o::OnlineStat)
    r = first(fieldnames(o))
    return :(o.$r)
end

#-----------------------------------------------------------------------# Base
Base.:(==)(o::OnlineStat, o2::OnlineStat) = false
function Base.:(==)(o1::T, o2::T) where {T<:OnlineStat}
    nms = fieldnames(typeof(o1))
    all(getfield.(Ref(o1), nms) .== getfield.(Ref(o2), nms))
end

Base.copy(o::OnlineStat) = deepcopy(o)

"""
    merge!(stat1, stat2)

Merge `stat1` into `stat2` (supported by most `OnlineStat` types).

# Example

    a = fit!(Mean(), 1:10)
    b = fit!(Mean(), 11:20)
    merge!(a, b)
"""
function Base.merge!(o::OnlineStat, o2::OnlineStat)
    (nobs(o) > 0 || nobs(o2) > 0) && _merge!(o, o2)
    o
end
_merge!(o, o2) = @warn("Merging $(name(o2)) into $(name(o)) is not well-defined.  No merging occurred.")
Base.merge(o::OnlineStat, o2::OnlineStat) = merge!(copy(o), o2)

#-----------------------------------------------------------------------# Show
function Base.show(io::IO, o::OnlineStat)
    print(io, name(o, false, false), ": ")
    print(io, "n=", nobs(o))
    print(io, " | value=")
    show(IOContext(io, :compact => true), value(o))
end
function name(T::Type, withmodule = false, withparams = true)
    s = string(T)
    if !withmodule
        s = replace(s, r"([a-zA-Z]*\.)" => "")  # remove text that ends in period
    end
    if !withparams
        s = replace(s, r"\{(.*)" => "")  # remove "{" to the end of the string
    end
    s
end
name(o, args...) = name(typeof(o), args...)

#-----------------------------------------------------------------------# fit!
"""
    fit!(stat::OnlineStat, data)

Update the "sufficient statistics" of a `stat` with more data.   If `typeof(data)` is not
the type of a single observation for the provided `stat`, `fit!` will attempt to iterate
through and `fit!` each item in `data`.  Therefore, `fit!(Mean(), 1:10)` translates
roughly to:

```
o = Mean()
for x in 1:10
    fit!(o, x)
end
```
"""
fit!(o::OnlineStat{T}, yi::T) where {T} = (_fit!(o, yi); return o)

function fit!(o::OnlineStat{I}, y::T) where {I, T}
    T == eltype(y) && error("The input for $(name(o,false,false)) is a $I.  Found $T.")
    for yi in y
        fit!(o, yi)
    end
    o
end

#-----------------------------------------------------------------------# utils
function _fit!(o::OnlineStat{T}, arg) where {T}
    error("A $(typeof(arg)) is not a single observation for $((name(o, false, true)))")
end

"""
    smooth(a, b, γ)

Weighted average of `a` and `b` with weight `γ`.
"""
smooth(a, b, γ) = a + γ * (b - a)

"""
    smooth!(a, b, γ)

Update `a` in place by applying the [`smooth`](@ref) function elementwise with `b`.
"""
function smooth!(a, b, γ)
    for (i, bi) in zip(eachindex(a), b)
        a[i] = smooth(a[i], bi, γ)
    end
end

"""
    smooth_syr!(A::AbstractMatrix, x, γ::Number)

Weighted average of symmetric rank-1 update.  Updates the upper triangle of:

`A = (1 - γ) * A + γ * x * x'`
"""
function smooth_syr!(A::AbstractMatrix, x, γ::Number)
    for j in 1:size(A, 2), i in 1:j
        A[i, j] = smooth(A[i,j], x[i] * conj(x[j]), γ)
    end
end

# bessel correction
bessel(o) = nobs(o) / (nobs(o) - 1)

Statistics.std(o::OnlineStat; kw...) = sqrt.(var(o; kw...))

input(o::OnlineStat{T}) where {T} = T

const TwoThings{T,S} = Union{Tuple{T,S}, Pair{T,S}, NamedTuple{names, Tuple{T,S}}} where names

#-----------------------------------------------------------------------# Compat
@static if VERSION < v"1.1.0"
    export eachrow, eachcol
    eachrow(A::Union{AbstractVector, AbstractMatrix}) = (view(A, i, :) for i in axes(A, 1))
    eachcol(A::Union{AbstractVector, AbstractMatrix}) = (view(A, :, i) for i in axes(A, 2))
else
    import Base: eachrow, eachcol
end

@deprecate eachrow(x::AbstractMatrix, y::AbstractVector) zip(eachrow(x), y)
@deprecate eachcol(x::AbstractMatrix, y::AbstractVector) zip(eachcol(x), y)

include("weight.jl")
include("stats.jl")
end


text:
# test will live in Caesar.jl
using Base.Test
using IncrementalInference, CloudGraphs
#
#
# IncrementalInference.setDataLayerAPI( DataLayerAPI(
# CloudGraphs.add_vertex!,
# CloudGraphs.get_vertex,
# CloudGraphs.make_edge,
# CloudGraphs.add_edge!,
# CloudGraphs.out_neighbors,
# CloudGraphs.update_vertex!,
# CloudGraphs.update_edge!
# ) )

# Example using CloudGraphs DB functionality
# Do fourdoorstest

# initCloudGrphMagicSomething()....

#CloudGraph.registerPackedType!(cloudGraph, IncrementalInference.PackedFunctionNodeData{PackedOdo},
#                               coolencode, cooldecode)

@test false


text:
###########################
# An n-dimensional simplex
###########################

"Supertype of an N-dimensional simplex."
abstract type Simplex{T,C} <: Domain{T} end

isclosedset(::Simplex{T,:closed}) where {T} = true
isclosedset(::Simplex{T,:open}) where {T} = false

isopenset(d::Simplex) = !isclosedset(d)

"The unit simplex is a polytope with the origin and all unit vectors as vertices."
abstract type UnitSimplex{T,C} <: Simplex{T,C} end

const ClosedUnitSimplex{T} = UnitSimplex{T,:closed}
const OpenUnitSimplex{T} = UnitSimplex{T,:open}

UnitSimplex(n::Int) = DynamicUnitSimplex(n)
UnitSimplex(::Val{N}) where {N} = EuclideanUnitSimplex{N}()

UnitSimplex{T}(n::Int) where {T <: StaticTypes} = StaticUnitSimplex{T}(n)
UnitSimplex{T}(::Val{N}) where {N,T} = StaticUnitSimplex{T}(Val(N))
UnitSimplex{T}() where {T <: StaticTypes} = StaticUnitSimplex{T}()
UnitSimplex{T}(n::Int) where {T} = DynamicUnitSimplex{T}(n)

UnitSimplex{T,C}(n::Int) where {T <: StaticTypes,C} = StaticUnitSimplex{T,C}(n)
UnitSimplex{T,C}(::Val{N}) where {N,T,C} = StaticUnitSimplex{T,C}(Val(N))
UnitSimplex{T,C}() where {T <: StaticTypes,C} = StaticUnitSimplex{T,C}()
UnitSimplex{T,C}(n::Int) where {T,C} = DynamicUnitSimplex{T,C}(n)

insimplex_closed(x) = mapreduce( t-> t >= 0, &, x) && norm(x,1) <= 1
insimplex_open(x) = mapreduce( t-> t > 0, &, x) && norm(x,1) < 1
insimplex_closed(x, tol) = mapreduce( t-> t >= -tol, &, x) && norm(x,1) <= 1+tol
insimplex_open(x, tol) = mapreduce( t-> t > -tol, &, x) && norm(x,1) < 1+tol

indomain(x, d::ClosedUnitSimplex) = length(x)==dimension(d) && insimplex_closed(x)
indomain(x, d::OpenUnitSimplex) = length(x)==dimension(d) && insimplex_open(x)
approx_indomain(x, d::ClosedUnitSimplex, tolerance) = length(x)==dimension(d) && insimplex_closed(x, tolerance)
approx_indomain(x, d::OpenUnitSimplex, tolerance) = length(x)==dimension(d) && insimplex_open(x, tolerance)

isempty(::UnitSimplex) = false

==(d1::UnitSimplex, d2::UnitSimplex) =
    isclosedset(d1)==isclosedset(d2) && dimension(d1)==dimension(d2)
hash(d::UnitSimplex, h::UInt) = hashrec("UnitSimplex", isclosedset(d), dimension(d), h)

boundingbox(d::UnitSimplex{T}) where {T} = UnitCube{T}(dimension(d))

boundary(d::UnitSimplex{T,:open}) where {T} = EmptySpace{T}()

distance_to(d::UnitSimplex, x) = x ∈ d ? zero(prectype(d)) : minimum(distance_to(el, x) for el in components(boundary(d)))

function normal(d::UnitSimplex, x)
    z = similar(x)
    fill!(z, 0)
    if sum(x) ≈ 1
        fill!(z, one(eltype(z))/sqrt(length(z)))
    else
        index = findmin(x)[2]
        z[index] = -1
    end
    return convert(eltype(d), z/norm(z))
end

"A unit simplex whose dimension is determined by its element type."
struct StaticUnitSimplex{T,C} <: UnitSimplex{T,C}
end

StaticUnitSimplex(::Val{N}) where {N} = StaticUnitSimplex{SVector{N,Float64}}()

StaticUnitSimplex{T}() where {T} = StaticUnitSimplex{T,:closed}()

StaticUnitSimplex{T}(n::Int) where {T} =
    (@assert n == euclideandimension(T); StaticUnitSimplex{T}())
StaticUnitSimplex{T}(::Val{N}) where {N,T} =
    (@assert N == euclideandimension(T); StaticUnitSimplex{T}())

StaticUnitSimplex{T,C}(n::Int) where {T,C} =
    (@assert n == euclideandimension(T); StaticUnitSimplex{T,C}())
StaticUnitSimplex{T,C}(::Val{N}) where {N,T,C} =
    (@assert N == euclideandimension(T); StaticUnitSimplex{T,C}())

const EuclideanUnitSimplex{N,T,C} = StaticUnitSimplex{SVector{N,T},C}

EuclideanUnitSimplex{N}() where {N} = EuclideanUnitSimplex{N,Float64}()

## A StaticUnitSimplex{<:Number} equals the interval [0,1]  (open or closed)
convert(::Type{Interval}, d::StaticUnitSimplex{T,:closed}) where {T <: Number} =
    UnitInterval{T}()
convert(::Type{Interval}, d::StaticUnitSimplex{T,:open}) where {T <: Number} =
    OpenInterval{T}(0, 1)

canonicaldomain(::Equal, d::StaticUnitSimplex{T}) where {T<:Number} = convert(Interval, d)

boundary(d::StaticUnitSimplex{T}) where {T<:Number} = boundary(convert(Interval, d))

"A unit simplex with vector elements with variable dimension determined by a field."
struct DynamicUnitSimplex{T,C} <: UnitSimplex{T,C}
    dimension   ::  Int

    DynamicUnitSimplex{T,C}(n::Int) where {T,C} = new(n)
    DynamicUnitSimplex{T,C}(n::Int) where {T<:StaticTypes,C} =
        (@assert n == euclideandimension(T); new(n))
end

DynamicUnitSimplex(n::Int) = DynamicUnitSimplex{Vector{Float64}}(n)
DynamicUnitSimplex{T}(n::Int) where {T} = DynamicUnitSimplex{T,:closed}(n)

dimension(d::DynamicUnitSimplex) = d.dimension

const VectorUnitSimplex{T,C} = DynamicUnitSimplex{Vector{T},C}

VectorUnitSimplex(dimension) = VectorUnitSimplex{Float64}(dimension)

show(io::IO, d::EuclideanUnitSimplex{N,Float64,:closed}) where {N} = print(io, "UnitSimplex(Val($(N)))")
show(io::IO, d::VectorUnitSimplex{Float64,:closed}) = print(io, "UnitSimplex($(dimension(d)))")

center(d::EuclideanUnitSimplex{N,T}) where {N,T} = ones(SVector{N,T})/N
center(d::VectorUnitSimplex{T}) where {T} = ones(T, dimension(d))/dimension(d)

corners(d::UnitSimplex) = vcat([origin(d)], [ unitvector(d, i) for i in 1:dimension(d)])

interior(d::EuclideanUnitSimplex{N,T}) where {N,T} = EuclideanUnitSimplex{N,T,:open}()
closure(d::EuclideanUnitSimplex{N,T}) where {N,T} = EuclideanUnitSimplex{N,T,:closed}()
interior(d::VectorUnitSimplex{T}) where {T} = VectorUnitSimplex{T,:open}(dimension(d))
closure(d::VectorUnitSimplex{T}) where {T} = VectorUnitSimplex{T,:closed}(dimension(d))

# We pick the center point, because it belongs to the domain regardless of
# whether it is open or closed.
point_in_domain(d::UnitSimplex) = center(d)

similardomain(d::StaticUnitSimplex{S,C}, ::Type{T}) where {S,T,C} =
    StaticUnitSimplex{T,C}()
similardomain(d::DynamicUnitSimplex{S,C}, ::Type{T}) where {S,T,C} =
    DynamicUnitSimplex{T,C}(d.dimension)

simplex_face_map(a::Number, b::Number, c::SVector{2}, d::SVector{2}) =
	AffineMap((d-c)/(b-a), c - (d-c)/(b-a)*a)
function simplex_face_map(a::Number, b::Number, c::Vector, d::Vector)
	@assert length(c) == length(d) == 2
	AffineMap((d-c)/(b-a), c - (d-c)/(b-a)*a)
end

function boundary(d::StaticUnitSimplex{SVector{2,T},:closed}) where {T}
    d0 = UnitInterval{T}()
	T0 = zero(T)
	T1 = one(T)
	maps = [
		simplex_face_map(T0, T1, SVector(T0,T0), SVector(T1,T0)),
		simplex_face_map(T0, T1, SVector(T1,T0), SVector(T0,T1)),
		simplex_face_map(T0, T1, SVector(T0,T1), SVector(T0,T0))
	]
	faces = map(m -> ParametricDomain(m, d0), maps)
	UnionDomain(faces)
end

# function boundary(d::StaticUnitSimplex{SVector{N,T},:closed}) where {N,T}
# 	left2 = infimum(d)
# 	right2 = supremum(d)
# 	d0 = UnitSimplex{SVector{N-1,T},:closed}()
# 	T0 = zero(T)
# 	T1 = one(T)
#
# 	map1 = cube_face_map(left1, right1, left2, right2, 1, left2[1])
# 	MAP = typeof(map1)
# 	maps = MAP[]
# 	for dim in 1:N
# 		push!(maps, cube_face_map(left1, right1, left2, right2, dim, left2[dim]))
# 		push!(maps, cube_face_map(left1, right1, left2, right2, dim, right2[dim]))
# 	end
# 	faces = map(m -> ParametricDomain(m, d_unit), maps)
# 	UnionDomain(faces)
# end


text:
using LabReports, Test
using LabReports: replace_powers, to_origin, comment_value

@testset "replace_power" begin
    @test replace_powers("300 mV s^-1") == "300 mV/s"
    @test replace_powers("10 m^2") == "10 m\\+(2)"
    @test replace_powers("1 μA^-2") == "1 μA\\+(-2)"
    @test replace_powers("1 μA^2") == "1 μA\\+(2)"
    @test replace_powers("15 mA cm^-2") == "15 mA/cm\\+(2)"
end


text:
"""
    substitute(expr, dict; fold=true)

substitute any subexpression that matches a key in `dict` with
the corresponding value. If `fold=false`,
expressions which can be evaluated won't be evaluated.

```julia
julia> substitute(1+sqrt(y), Dict(y => 2), fold=true)
2.414213562373095
julia> substitute(1+sqrt(y), Dict(y => 2), fold=false)
1 + sqrt(2)
```
"""
function substitute(expr, dict; fold=true)
    haskey(dict, expr) && return dict[expr]

    if istree(expr)
        op = substitute(operation(expr), dict; fold=fold)
        if fold
            canfold = !(op isa Symbolic)
            args = map(unsorted_arguments(expr)) do x
                x′ = substitute(x, dict; fold=fold)
                canfold = canfold && !(x′ isa Symbolic)
                x′
            end
            canfold && return op(args...)
            args
        else
            args = map(x->substitute(x, dict, fold=fold), unsorted_arguments(expr))
        end

        similarterm(expr,
                    op,
                    args,
                    symtype(expr);
                    metadata=metadata(expr))
    else
        expr
    end
end

"""
    occursin(needle::Symbolic, haystack::Symbolic)

Determine whether the second argument contains the first argument. Note that
this function doesn't handle associativity, commutativity, or distributivity.
"""
Base.occursin(needle::Symbolic, haystack::Symbolic) = _occursin(needle, haystack)
Base.occursin(needle, haystack::Symbolic) = _occursin(needle, haystack)
Base.occursin(needle::Symbolic, haystack) = _occursin(needle, haystack)
function _occursin(needle, haystack)
    isequal(needle, haystack) && return true

    if istree(haystack)
        args = arguments(haystack)
        for arg in args
            occursin(needle, arg) && return true
        end
    end
    return false
end


text:
# ----------------------------------------------------
# --- Import modules
# ----------------------------------------------------
using DigitalComm
using Test
# ----------------------------------------------------
# --- Tests
# ----------------------------------------------------
println("Tests for Hard constellation");

for alphabet = [4;16;64;256]
	@testset  "$alphabet-QAM" begin
		# Minimal Tx
		mcs	  = alphabet;		  # Constellation size
		n	  = Int(log2(mcs));		  # Bit per symbol
		N	  = 4096;	  # Numbver of test symbols
		bitSeq	= genBitSequence(N*n);
		qamSeq	= bitMappingQAM(mcs,bitSeq);
		# Minimal Rx
		qamDec	= hardConstellation(mcs,qamSeq);
		# Checking size of decoding vector
		@test length(qamDec) == N;
		# Checking that bang and classic methods are the same
		buff  = zeros(Complex{Float64},N);
		hardConstellation!(buff,mcs,qamSeq);
		@test all(qamDec .== buff);
		# Checking that decoded method is initial method
		@test all(qamDec .≈ qamSeq)
		# Checking that introduces a small biais does not change voronoi region
		biais	= 1/sqrt(2/3*(mcs-1)) * 0.1 * (1+1im)
		@test all( hardConstellation(mcs, qamSeq .+ biais) .≈ qamSeq)
	end
end


text:
module PMTK4BookDemos

    export unpairedTestHandednessDemo, amazonSellersDelta, amazonSellersTheta
    include("BayesianStatistics/amazonSellerDemo.jl")
    export drawBayesChangeOfVar
    include("BayesianStatistics/bayesChangeOfVar.jl")
    export drawBetaHPD
    include("BayesianStatistics/betaHPD.jl")
    export betaCredibleInt
    include("BayesianStatistics/betaCredibleInt.jl")
    export drawBimodalDemo
    include("BayesianStatistics/bimodalDemo.jl")
    export coinsModelSelDemo
    include("BayesianStatistics/coinsModelSelDemo.jl")
    export drawGammaDistribution
    include("BayesianStatistics/gammaPlotDemo.jl")
    export lossDemo, hingeLossDemo
    include("BayesianStatistics/lossFunctionDemo.jl")
    export mcBetaQuantileDemo, mcNormalQuantileDemo
    include("BayesianStatistics/mcQuantileDemo.jl")
    export mixBetaDemo, mixBetaBiased
    include("BayesianStatistics/mixBetaDemo.jl")
    export PRhandDemo, ROChandDemo
    include("BayesianStatistics/PRhand.jl")
    export robustPriorDemo
    include("BayesianStatistics/robustPriorDemo.jl")
    export shrinkageDemoBaseball, shrinkageMSEBaseball
    include("BayesianStatistics/shrinkageDemoBaseball.jl")

end


text:
function deconv_boot_test(X::Vector, Y::Vector, σ_X::Vector, σ_Y::Vector, t::Function, method::DeconvolutionMethod; niter=100)
    X_and_Y = cat(X, Y, dims=1)
    σ_X_and_Y = cat(σ_X, σ_Y, dims=1)

    ϵ_X = Normal.(0.0, σ_X)
    ϵ_Y = Normal.(0.0, σ_Y)
    ϵ_X_and_Y = Normal.(0.0, σ_X_and_Y)

    t_obs = t(X, Y, σ_X, σ_Y)
    n_X = length(X)
    n_Y = length(Y)
    n_XY = n_X + n_Y
    null_decon = decon(method, X_and_Y, ϵ_X_and_Y)
    nabove = 0
    for i in 1:niter
        Xboot = rand(null_decon, n_X)
        Yboot = rand(null_decon, n_Y)
        Xtilde = Xboot .+ rand.(ϵ_X)
        Ytilde = Yboot .+ rand.(ϵ_Y)
        t_boot = t(Xtilde, Ytilde, σ_X, σ_Y)
        if t_boot > t_obs
            nabove += 1
        end
    end
    return nabove / niter
end

function deconv_boot_stat(X::Vector, Y::Vector, σ_X::Vector, σ_Y::Vector, t::Function, method::DeconvolutionMethod; niter=100)
    X_and_Y = cat(X, Y, dims=1)
    σ_X_and_Y = cat(σ_X, σ_Y, dims=1)

    ϵ_X = Normal.(0.0, σ_X)
    ϵ_Y = Normal.(0.0, σ_Y)
    ϵ_X_and_Y = Normal.(0.0, σ_X_and_Y)

    t_obs = t(X, Y, σ_X, σ_Y)
    n_X = length(X)
    n_Y = length(Y)
    n_XY = n_X + n_Y
    t_record = Vector{typeof(t_obs)}(niter)
    null_decon = decon(method, X_and_Y, ϵ_X_and_Y)
    for i in 1:niter
        Xboot = rand(null_decon, n_X)
        Yboot = rand(null_decon, n_Y)
        Xtilde = Xboot .+ rand.(ϵ_X)
        Ytilde = Yboot .+ rand.(ϵ_Y)
        t_boot = t(Xtilde, Ytilde, σ_X, σ_Y)
        t_record[i] = t_boot
    end
    return t_record
end


text:
using DataFrames
using Gadfly
using Mamba
using Distributions
using StatsBase

dat = readtable("test.txt.gz", separator = '\t', header=false, names=[:Length, :DeltaG])

hex=plot(dat, x=:Length, y=:DeltaG, Geom.hexbin(xbincount=35),
         Coord.cartesian(xmin=0, ymax=15),
         Guide.xlabel("Duplex Length"),
         Guide.ylabel("Gibbs Free Energy Change (ΔG)"))

dens=plot(dat[dat[:Length] .>= 5,:], color=:Length, x=:DeltaG, Geom.density, Guide.colorkey("Duplex Length"),
          Guide.xlabel("Gibbs Free Energy Change (ΔG)"))

dd = DuplexDist( dat )
model=plotcdf( dd, dat[:Length], dat[:DeltaG] )

draw(PDF("deltaG.pdf", 8inch, 8inch), vstack(hstack(hex,dens), model))

type DuplexDist
   length::Distributions.Poisson
   deltag::Vector{Distributions.Normal}
end

# G[L] ~ Normal(mu[L], s[L])
# L ~ Poisson(lambda)
function DuplexDist( df::DataFrame )
   # Fit Model
   dgdat = Vector{DataArray}()
   for i in 1:maximum(df[:Length])
      push!( dgdat, dat[dat[:,1] .== i, 2] )
   end
   lendat = df[:Length]

   len = fit(Poisson, lendat)
   dg  = Vector{Distributions.Normal}(length(dgdat))
   for i in 1:N
      dg[i] = fit(Normal, dgdat[i])
   end

   DuplexDist( len, dg )
end

function cdf( dd::DuplexDist, len::Int, dg::Float64 )
   prob = 0.0
   for i in len:length(dd.deltag)
      prob += pdf( dd.length, i ) * cdf( dd.deltag[i], dg )
   end
   prob += ccdf( dd.length, length(dd.deltag)+1 )
   prob
end

# P(G,L) = P(G|L)P(L)
pdf( dd::DuplexDist, len::Int, dg::Float64 ) = pdf( dd.length, len ) * pdf( dd.deltag[len], dg )

function plotcdf( dd::DuplexDist, lenvec, dgvec )
   l = Vector{Int}()
   d = Vector{Float64}()
   c = Vector{Float64}()
   p = Vector{Float64}()
   irange = 1:maximum(lenvec)
   jrange = minimum(dgvec):1.0:maximum(dgvec)
   for i in irange
      for j in jrange
         push!( l, i )
         push!( d, j )
         push!( c, cdf( dd, i, j ) )
         push!( p, pdf( dd, i, j ) )
      end
   end
   z = reshape( p, (length(jrange), length(irange)))

   pval = plot(layer(x=l, y=d, color=c, Geom.rectbin),
         Coord.cartesian(xmin=0, ymax=15),
         Guide.xlabel("Duplex Length"),
         Guide.ylabel("Gibbs Free Energy Change (ΔG)"),
         Guide.colorkey("p-value"), Scale.color_log10)

   prob =plot(layer(x=l, y=d, color=p, Geom.rectbin),
         Coord.cartesian(xmin=0, ymax=15),
         Guide.xlabel("Duplex Length"),
         Guide.ylabel("Gibbs Free Energy Change (ΔG)"),
         Guide.colorkey("Pdf(x)"))

   hstack(prob, pval)
end

#= MAMBA

fitdat = Dict{Symbol, Any}(
  :len => len,
  :dg => dg
)
N = length(dg)

model = Model(
   dg = Stochastic(1, (len, mu, s2) -> UnivariateDistribution[ Normal(mu[i], sqrt(s2[i])) for i in 1:N ], false),

   mu = Stochastic(1, () -> Normal(0, 1)),
   s2 = Stochastic(1, () -> InverseGamma(0.01, 0.01)),

   len = Stochastic(alpha -> begin
                               lambda = max(0, 1/alpha)
                               Poisson(lambda)
                             end, false),
   alpha = Stochastic(() -> Exponential(1/10)),
)

inits = [
  Dict{Symbol, Any}(
     :alpha => rand(Exponential(1/10)),
     :len => mean(fitdat[:len]),
     :s2 => [rand(InverseGamma(0.01, 0.01)) for i in 1:N],
     :mu => [mean(dg[i]) for i in 1:N],
     :dg => [mean(dg[i]) for i in 1:N]
  )
for j in 1:5
]

scheme = [Slice(:alpha, 1.0, Univariate),
          Slice(:len, 1.0, Univariate),
          NUTS([:mu, :s2])]

setinputs!(model, fitdat)
setinits!(model, inits[1])
setsamplers!(model, scheme)

sim = mcmc(model, fitdat, inits, 10000, burnin=250, thin=2, chains=1)

=#


text:
function move(model::CRRJITContractPremiumLatticeModel;
    from::PQContractPremiumLatticePoint, to::PQContractPremiumLatticePoint)::Array{Float64,1}

    # initialize -
    premium_array = Array{Float64,1}()
    p_from = premium(model,from)
    p_to = premium(model,to)

    push!(premium_array, p_from)
    push!(premium_array, p_to)

    # return -
    return premium_array
end

function roll(model::CRRJITContractPremiumLatticeModel, origin::PQContractPremiumLatticePoint,
    move::Pair{PQContractPremiumLatticePoint, PQContractPremiumLatticePoint})::DataFrame

    # initialize -
    premium_array = Array{Float64,2}(undef,3,5)

    # how much was the original positon?
    premium_array[1,1] = premium(model, origin)
    premium_array[2,1] = premium(model, move.first)
    premium_array[3,1] = premium(model, move.second)

    # compute delta -
    premium_array[1,2] = δ(model, origin)
    premium_array[2,2] = δ(model, move.first)
    premium_array[3,2] = δ(model, move.second)

    # compute θ -
    premium_array[1,3] = θ(model, origin)
    premium_array[2,3] = θ(model, move.first)
    premium_array[3,3] = θ(model, move.second)

    # compute vega -
    premium_array[1,4] = vega(model, origin)
    premium_array[2,4] = vega(model, move.first)
    premium_array[3,4] = vega(model, move.second)

    # compute γ -
    premium_array[1,5] = γ(model, origin)
    premium_array[2,5] = γ(model, move.first)
    premium_array[3,5] = γ(model, move.second)

    # package -
    df = DataFrame(P = premium_array[:,1] , δ = premium_array[:,2], θ = premium_array[:,3],
        vega = premium_array[:,4], γ = premium_array[:,5])

    # return -
    return df
end

function roll(model::CRRJITContractPremiumLatticeModel,
    move::Pair{PQContractPremiumLatticePoint, PQContractPremiumLatticePoint})::DataFrame

    # initialize -
    premium_array = Array{Float64,2}(undef,2,5)

    # how much was the original positon?
    premium_array[1,1] = premium(model, move.first)
    premium_array[2,1] = premium(model, move.second)

    # compute delta -
    premium_array[1,2] = δ(model, move.first)
    premium_array[2,2] = δ(model, move.second)

    # compute delta -
    premium_array[1,3] = θ(model, move.first)
    premium_array[2,3] = θ(model, move.second)

    # compute vega -
    premium_array[1,4] = vega(model, move.first)
    premium_array[2,4] = vega(model, move.second)

    # compute γ -
    premium_array[1,5] = γ(model, move.first)
    premium_array[2,5] = γ(model, move.second)

    # package -
    df = DataFrame(P = premium_array[:,1] , δ = premium_array[:,2], θ = premium_array[:,3],
        vega = premium_array[:,4], γ = premium_array[:,5])

    # return -
    return df
end


text:
struct ExactNonlinearSolver <: NonlinearCombinatorialAlgorithm
    solver::Any # MINLP solver with support for the nonlinear function.
end

function solve(nli::NonlinearCombinatorialInstance, algo::ExactNonlinearSolver)
    return solve(nli, algo, Val(nli.nonlinear_function))
end

function solve(
    nli::NonlinearCombinatorialInstance,
    algo::ExactNonlinearSolver,
    ::Val{SquareRoot},
)
    # This formulation is valid only for maximising (concave objective
    # function). Otherwise, think about general MINLP solvers...
    @assert nli.combinatorial_structure.objective == Maximise()

    # Base formulation for the combinatorial problem.
    m, x = formulation(nli.combinatorial_structure, nli.formulation)

    # Formulation for the nonlinear term.
    dot_product = @variable(m, lower_bound = 0.0) # b^T x
    nonlinear_term = @variable(m, lower_bound = 0.0) # t

    set_name(dot_product, "dot_product")
    set_name(nonlinear_term, "nonlinear_term")

    @constraint(
        m,
        eq,
        dot_product == sum(
            nli.nonlinear_coefficients[i] * x[i] for
            i in eachindex(nli.nonlinear_coefficients)
        )
    )
    @constraint(
        m,
        cone,
        [1 + dot_product, 1 - dot_product, 2 * nonlinear_term] in
        SecondOrderCone()
    )

    # New objective.
    @objective(
        m,
        Max,
        sum(
            nli.linear_coefficients[i] * x[i] for
            i in eachindex(nli.nonlinear_coefficients)
        ) + nonlinear_term
    )

    # Solve the new formulation.
    set_silent(m)
    t0 = time_ns()
    optimize!(m)
    t1 = time_ns()

    if termination_status(m) != MOI.OPTIMAL
        error("The model was not solved correctly.")
    end

    # Retrieve the solution.
    dict_sol =
        Dict(k => JuMP.value(x[k]) for k in eachindex(nli.linear_coefficients))
    return make_solution(nli.combinatorial_structure, dict_sol)
end

function solve(
    nli::NonlinearCombinatorialInstance,
    algo::ExactNonlinearSolver,
    ::Val{Square},
)
    # This formulation is valid only for minimising (convex objective
    # function). Otherwise, think about general MINLP solvers...
    @assert nli.combinatorial_structure.objective == Minimise()

    # Base formulation for the combinatorial problem.
    m, x = formulation(nli.combinatorial_structure, solver=algo.solver)

    # Formulation for the nonlinear term. Quite easy, JuMP supports quadratic
    # objective functions.
    dot_product = @variable(m, lower_bound = 0.0) # b^T x
    @constraint(
        m,
        eq,
        dot_product == sum(
            nli.nonlinear_coefficients[i] * x[i] for
            i in eachindex(nli.nonlinear_coefficients)
        )
    )

    # New objective.
    @objective(
        m,
        Max,
        sum(
            nli.linear_coefficients[i] * x[i] for
            i in eachindex(nli.nonlinear_coefficients)
        ) + nonlinear_term^2
    )

    # Solve the new formulation.
    set_silent(m)
    t0 = time_ns()
    optimize!(m)
    t1 = time_ns()

    if termination_status(m) != MOI.OPTIMAL
        error("The model was not solved correctly.")
    end

    # Retrieve the solution.
    dict_sol =
        Dict(k => JuMP.value(x[k]) for k in eachindex(nli.linear_coefficients))
    return make_solution(nli.combinatorial_structure, dict_sol)
end


text:
# WARNING: This file was generated from seq/fastq/parser.rl using ragel. Do not edit!
const fastqparser_start  = 18
const fastqparser_first_final  = 18
const fastqparser_error  = 0
const fastqparser_en_main  = 18
const _fastqparser_nfa_targs = Int8[0, 0 , ]
const _fastqparser_nfa_offsets = Int8[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , ]
const _fastqparser_nfa_push_actions = Int8[0, 0 , ]
const _fastqparser_nfa_pop_trans = Int8[0, 0 , ]
# FIXME: output type may be too loose
Ragel.@generate_read!_function(
"fastqparser",
FASTQReader,
SeqRecord,
begin
	if p == pe
		@goto _test_eof

	end
	if cs  == 18
		@goto st_case_18
	elseif cs  == 0
		@goto st_case_0
	elseif cs  == 1
		@goto st_case_1
	elseif cs  == 2
		@goto st_case_2
	elseif cs  == 3
		@goto st_case_3
	elseif cs  == 4
		@goto st_case_4
	elseif cs  == 5
		@goto st_case_5
	elseif cs  == 6
		@goto st_case_6
	elseif cs  == 7
		@goto st_case_7
	elseif cs  == 8
		@goto st_case_8
	elseif cs  == 9
		@goto st_case_9
	elseif cs  == 10
		@goto st_case_10
	elseif cs  == 11
		@goto st_case_11
	elseif cs  == 12
		@goto st_case_12
	elseif cs  == 19
		@goto st_case_19
	elseif cs  == 13
		@goto st_case_13
	elseif cs  == 14
		@goto st_case_14
	elseif cs  == 15
		@goto st_case_15
	elseif cs  == 16
		@goto st_case_16
	elseif cs  == 17
		@goto st_case_17
	end
	@goto st_out
	@label ctr39
	state.linenum += 1
	@goto st18
	@label st18
	p+= 1
	if p == pe
		@goto _test_eof18

	end
	@label st_case_18
	if (data[1+(p )]) == 10
		@goto ctr39
	elseif (data[1+(p )]) == 32
		@goto st18
	elseif (data[1+(p )]) == 64
		ck = 0;
		if (length(input.qualbuf) == length(input.seqbuf)
			)
			ck += 1

		end
		if 1 <= ck
			@goto st1

		end
		@goto st0
	end
	if 9 <= (data[1+(p )])&& (data[1+(p )])<= 13
		@goto st18

	end
	@goto st0
	@label st_case_0
	@label st0
	cs = 0
	@goto _out
	@label ctr41
	if length(input.seqbuf) != length(input.qualbuf)
		error("Error parsing FASTQ: sequence and quality scores must be of equal length")
	end

	# if a name and s description are repeated (i.e. old-fashioned fastq)
		# make sure they match.
		if (!isempty(input.name2buf) && input.name2buf != output.name) ||
			(!isempty(input.desc2buf) && input.desc2buf != output.metadata.description)
			error("Error parsing FASTQ: sequence and quality scores have non-matching identifiers")
		end

		# sequence
		resize!(output.seq, input.seqbuf.position - 1)
		if !isnull(input.fill_ambiguous)
			byte = convert(UInt8, convert(Char, get(input.fill_ambiguous)))
			for i in 1:input.seqbuf.position-1
			b = input.seqbuf.buffer[i]
			if isambiguous(convert(DNANucleotide, convert(Char, b)))
				input.seqbuf.buffer[i] = byte
			end
		end
	end
	encode_copy!(output.seq, 1, input.seqbuf.buffer, 1, input.seqbuf.position - 1)

	# quality
	check_quality_string(input.quality_encoding, input.qualbuf.buffer, 1, input.qualbuf.position - 1)
	decode_quality_string!(input.quality_encoding, input.qualbuf.buffer,
	output.metadata.quality, 1,
	input.qualbuf.position - 1)

	# reset temporaries for the next run
	empty!(input.qualbuf)
	empty!(input.seqbuf)
	empty!(input.name2buf)
	empty!(input.desc2buf)

	Ragel.@yield 1
	@goto st1
	@label st1
	p+= 1
	if p == pe
		@goto _test_eof1

	end
	@label st_case_1
	if (data[1+(p )])== 32
		@goto st0

	end
	if 9 <= (data[1+(p )])&& (data[1+(p )])<= 13
		@goto st0

	end
	@goto ctr0
	@label ctr0
	Ragel.@anchor!; @goto st2
	@label st2
	p+= 1
	if p == pe
		@goto _test_eof2

	end
	@label st_case_2
	if (data[1+(p )]) == 10
		@goto ctr4
	elseif (data[1+(p )]) == 12
		@goto st0
	elseif (data[1+(p )]) == 13
		@goto ctr5
	elseif (data[1+(p )]) == 32
		@goto ctr3
	end
	if 9 <= (data[1+(p )])&& (data[1+(p )])<= 11
		@goto ctr3

	end
	@goto st2
	@label ctr3
	Ragel.@copy_from_anchor!(output.name); @goto st3
	@label st3
	p+= 1
	if p == pe
		@goto _test_eof3

	end
	@label st_case_3
	if (data[1+(p )]) == 10
		@goto st0
	elseif (data[1+(p )]) == 32
		@goto st3
	end
	if (data[1+(p )])> 11
		if (data[1+(p )])<= 13
			@goto st0

		end

	elseif (data[1+(p )])>= 9
		@goto st3

	end
	@goto ctr6
	@label ctr6
	Ragel.@anchor!; @goto st4
	@label st4
	p+= 1
	if p == pe
		@goto _test_eof4

	end
	@label st_case_4
	if (data[1+(p )]) == 10
		@goto ctr9
	elseif (data[1+(p )]) == 13
		@goto ctr10
	end
	@goto st4
	@label ctr4
	Ragel.@copy_from_anchor!(output.name); state.linenum += 1
	@goto st5
	@label ctr9
	Ragel.@copy_from_anchor!(output.metadata.description); state.linenum += 1
	@goto st5
	@label ctr37
	state.linenum += 1
	@goto st5
	@label st5
	p+= 1
	if p == pe
		@goto _test_eof5

	end
	@label st_case_5
	if (data[1+(p )]) == 10
		@goto ctr11
	elseif (data[1+(p )]) == 13
		@goto st7
	end
	if 65 <= (data[1+(p )])&& (data[1+(p )])<= 122
		@goto ctr13

	end
	@goto st0
	@label ctr11
	state.linenum += 1
	@goto st6
	@label ctr34
	Ragel.@append_from_anchor!(input.seqbuf); state.linenum += 1
	@goto st6
	@label st6
	p+= 1
	if p == pe
		@goto _test_eof6

	end
	@label st_case_6
	if (data[1+(p )]) == 10
		@goto ctr11
	elseif (data[1+(p )]) == 13
		@goto st7
	elseif (data[1+(p )]) == 43
		@goto st8
	end
	if 65 <= (data[1+(p )])&& (data[1+(p )])<= 122
		@goto ctr13

	end
	@goto st0
	@label ctr35
	Ragel.@append_from_anchor!(input.seqbuf); @goto st7
	@label st7
	p+= 1
	if p == pe
		@goto _test_eof7

	end
	@label st_case_7
	if (data[1+(p )])== 10
		@goto ctr11

	end
	@goto st0
	@label st8
	p+= 1
	if p == pe
		@goto _test_eof8

	end
	@label st_case_8
	if (data[1+(p )]) == 10
		@goto ctr16
	elseif (data[1+(p )]) == 13
		@goto st15
	elseif (data[1+(p )]) == 32
		@goto st0
	end
	if 9 <= (data[1+(p )])&& (data[1+(p )])<= 12
		@goto st0

	end
	@goto ctr15
	@label ctr15
	Ragel.@anchor!; @goto st9
	@label st9
	p+= 1
	if p == pe
		@goto _test_eof9

	end
	@label st_case_9
	if (data[1+(p )]) == 10
		@goto ctr20
	elseif (data[1+(p )]) == 12
		@goto st0
	elseif (data[1+(p )]) == 13
		@goto ctr21
	elseif (data[1+(p )]) == 32
		@goto ctr19
	end
	if 9 <= (data[1+(p )])&& (data[1+(p )])<= 11
		@goto ctr19

	end
	@goto st9
	@label ctr19
	Ragel.@copy_from_anchor!(input.name2buf); @goto st10
	@label st10
	p+= 1
	if p == pe
		@goto _test_eof10

	end
	@label st_case_10
	if (data[1+(p )]) == 10
		@goto st0
	elseif (data[1+(p )]) == 32
		@goto st10
	end
	if (data[1+(p )])> 11
		if (data[1+(p )])<= 13
			@goto st0

		end

	elseif (data[1+(p )])>= 9
		@goto st10

	end
	@goto ctr22
	@label ctr22
	Ragel.@anchor!; @goto st11
	@label st11
	p+= 1
	if p == pe
		@goto _test_eof11

	end
	@label st_case_11
	if (data[1+(p )]) == 10
		@goto ctr25
	elseif (data[1+(p )]) == 13
		@goto ctr26
	end
	@goto st11
	@label ctr16
	state.linenum += 1
	@goto st12
	@label ctr20
	Ragel.@copy_from_anchor!(input.name2buf); state.linenum += 1
	@goto st12
	@label ctr25
	Ragel.@copy_from_anchor!(input.desc2buf); state.linenum += 1
	@goto st12
	@label st12
	p+= 1
	if p == pe
		@goto _test_eof12

	end
	@label st_case_12
	if (data[1+(p )]) == 10
		@goto ctr27
	elseif (data[1+(p )]) == 13
		@goto st13
	end
	if 33 <= (data[1+(p )])&& (data[1+(p )])<= 126
		ck = 0;
		if (length(input.qualbuf) + input.qualcount < length(input.seqbuf)
			)
			ck += 1

		end
		if 1 <= ck
			@goto ctr29

		end
		@goto st0

	end
	@goto st0
	@label ctr27
	state.linenum += 1
	@goto st19
	@label ctr31
	Ragel.@append_from_anchor!(input.qualbuf)
	input.qualcount = 0
	state.linenum += 1
	@goto st19
	@label st19
	p+= 1
	if p == pe
		@goto _test_eof19

	end
	@label st_case_19
	if (data[1+(p )]) == 10
		@goto ctr27
	elseif (data[1+(p )]) == 13
		@goto st13
	elseif (data[1+(p )]) == 64
		ck = 0;
		if (length(input.qualbuf) + input.qualcount < length(input.seqbuf)
			)
			ck += 1

		end
		if (length(input.qualbuf) == length(input.seqbuf)
			)
			ck += 2

		end
		if ck < 2
			if 1 <= ck
				@goto ctr29

			end

		elseif ck > 2
			@goto ctr29

		else
			@goto ctr41

		end
		@goto st0
	end
	if 33 <= (data[1+(p )])&& (data[1+(p )])<= 126
		ck = 0;
		if (length(input.qualbuf) + input.qualcount < length(input.seqbuf)
			)
			ck += 1

		end
		if 1 <= ck
			@goto ctr29

		end
		@goto st0

	end
	@goto st0
	@label ctr32
	Ragel.@append_from_anchor!(input.qualbuf)
	input.qualcount = 0
	@goto st13
	@label st13
	p+= 1
	if p == pe
		@goto _test_eof13

	end
	@label st_case_13
	if (data[1+(p )])== 10
		@goto ctr27

	end
	@goto st0
	@label ctr29
	Ragel.@anchor!; input.qualcount += 1
	@goto st14
	@label ctr33
	input.qualcount += 1
	@goto st14
	@label st14
	p+= 1
	if p == pe
		@goto _test_eof14

	end
	@label st_case_14
	if (data[1+(p )]) == 10
		@goto ctr31
	elseif (data[1+(p )]) == 13
		@goto ctr32
	end
	if 33 <= (data[1+(p )])&& (data[1+(p )])<= 126
		ck = 0;
		if (length(input.qualbuf) + input.qualcount < length(input.seqbuf)
			)
			ck += 1

		end
		if 1 <= ck
			@goto ctr33

		end
		@goto st0

	end
	@goto st0
	@label ctr21
	Ragel.@copy_from_anchor!(input.name2buf); @goto st15
	@label ctr26
	Ragel.@copy_from_anchor!(input.desc2buf); @goto st15
	@label st15
	p+= 1
	if p == pe
		@goto _test_eof15

	end
	@label st_case_15
	if (data[1+(p )])== 10
		@goto ctr16

	end
	@goto st0
	@label ctr13
	Ragel.@anchor!; @goto st16
	@label st16
	p+= 1
	if p == pe
		@goto _test_eof16

	end
	@label st_case_16
	if (data[1+(p )]) == 10
		@goto ctr34
	elseif (data[1+(p )]) == 13
		@goto ctr35
	end
	if 65 <= (data[1+(p )])&& (data[1+(p )])<= 122
		@goto st16

	end
	@goto st0
	@label ctr5
	Ragel.@copy_from_anchor!(output.name); @goto st17
	@label ctr10
	Ragel.@copy_from_anchor!(output.metadata.description); @goto st17
	@label st17
	p+= 1
	if p == pe
		@goto _test_eof17

	end
	@label st_case_17
	if (data[1+(p )])== 10
		@goto ctr37

	end
	@goto st0
	@label st_out
	@label _test_eof18
	cs = 18
	@goto _test_eof
	@label _test_eof1
	cs = 1
	@goto _test_eof
	@label _test_eof2
	cs = 2
	@goto _test_eof
	@label _test_eof3
	cs = 3
	@goto _test_eof
	@label _test_eof4
	cs = 4
	@goto _test_eof
	@label _test_eof5
	cs = 5
	@goto _test_eof
	@label _test_eof6
	cs = 6
	@goto _test_eof
	@label _test_eof7
	cs = 7
	@goto _test_eof
	@label _test_eof8
	cs = 8
	@goto _test_eof
	@label _test_eof9
	cs = 9
	@goto _test_eof
	@label _test_eof10
	cs = 10
	@goto _test_eof
	@label _test_eof11
	cs = 11
	@goto _test_eof
	@label _test_eof12
	cs = 12
	@goto _test_eof
	@label _test_eof19
	cs = 19
	@goto _test_eof
	@label _test_eof13
	cs = 13
	@goto _test_eof
	@label _test_eof14
	cs = 14
	@goto _test_eof
	@label _test_eof15
	cs = 15
	@goto _test_eof
	@label _test_eof16
	cs = 16
	@goto _test_eof
	@label _test_eof17
	cs = 17
	@goto _test_eof
	@label _test_eof
	if p == eof
		if cs  == 19
			if length(input.seqbuf) != length(input.qualbuf)
				error("Error parsing FASTQ: sequence and quality scores must be of equal length")
			end

			# if a name and s description are repeated (i.e. old-fashioned fastq)
				# make sure they match.
				if (!isempty(input.name2buf) && input.name2buf != output.name) ||
					(!isempty(input.desc2buf) && input.desc2buf != output.metadata.description)
					error("Error parsing FASTQ: sequence and quality scores have non-matching identifiers")
				end

				# sequence
				resize!(output.seq, input.seqbuf.position - 1)
				if !isnull(input.fill_ambiguous)
					byte = convert(UInt8, convert(Char, get(input.fill_ambiguous)))
					for i in 1:input.seqbuf.position-1
					b = input.seqbuf.buffer[i]
					if isambiguous(convert(DNANucleotide, convert(Char, b)))
						input.seqbuf.buffer[i] = byte
					end
				end
			end
			encode_copy!(output.seq, 1, input.seqbuf.buffer, 1, input.seqbuf.position - 1)

			# quality
			check_quality_string(input.quality_encoding, input.qualbuf.buffer, 1, input.qualbuf.position - 1)
			decode_quality_string!(input.quality_encoding, input.qualbuf.buffer,
			output.metadata.quality, 1,
			input.qualbuf.position - 1)

			# reset temporaries for the next run
			empty!(input.qualbuf)
			empty!(input.seqbuf)
			empty!(input.name2buf)
			empty!(input.desc2buf)

			Ragel.@yield 0
		end

	end
	@label _out
end)


text:
using Aqua
using MetaDataGraphs
using Documenter
using Graphs
using JET
using LinearAlgebra
using Test

DocMeta.setdocmeta!(MetaDataGraphs, :DocTestSetup, :(using MetaDataGraphs); recursive=true)

@testset verbose = true "MetaDataGraphs.jl" begin
    @testset verbose = true "Doctests" begin
        doctest(MetaDataGraphs)
    end
    @testset verbose = true "Code quality" begin
        Aqua.test_all(MetaDataGraphs)
    end
    @testset verbose = true "Structure" begin
        include("structure.jl")
    end
end


text:
############################################################
# FactCheck.jl
# A testing framework for Julia
# http://github.com/JuliaLang/FactCheck.jl
# MIT Licensed
############################################################

module TestFactCheck

using FactCheck
using Base.Test

############################################################
# Before we excerse the other various parts of FactCheck,
# check we actually catch and report errors correctly. This
# also allows us to test printing code for the Failure and
# Error cases, which wouldn't be tested otherwise.
print_with_color(:blue,"Testing Result counting and printing, not actual errors!\n")
facts("Test error pathways") do
    a_success = @fact 1 => 1 "I will never be seen"
    println(a_success)
    a_failure = @fact 1 => 2 "one doesn't equal two!"
    a_error   = @fact 2^-1 => 0.5 "domains are tricky"
    a_pending = @pending not_really_pending() "sorta pending"
    println(a_pending)
end
stats = getstats()
FactCheck.clear_results()
@test stats["nSuccesses"] == 1
@test stats["nFailures"] == 1
@test stats["nErrors"] == 1
@test stats["nPending"] == 1
@test stats["nNonSuccessful"] == 2
print_with_color(:blue,"Done, begin actual FactCheck tests\n")

############################################################
# Begin actual tests
type Foo a end
type Bar a end
type Baz end
type Bazz a end
==(x::Foo, y::Foo) = x.a == y.a

facts("Testing core functionality") do
    @fact 1 => 1
    @fact 2*2 => 4
    @fact uppercase("foo") => "FOO"
    @fact_throws 2^-1
    @fact 2*[1,2,3] => [2,4,6]
    @fact Foo(1) => Foo(1)
end

facts("FactCheck assertion helper functions") do

    context("`not` works for values and functions") do
        notone = not(1)
        @fact notone(2) => true
        @fact notone(1) => false
        @fact 2 => not(1)

        noteven = not(iseven)
        @fact noteven(3) => true
        @fact noteven(2) => false
        @fact not(iseven)(2) => false
        @fact 3 => not(iseven)
    end

    context("`truthy` is anything other than nothing or false (which is 0)") do
        @fact truthy(-1) => true
        @fact truthy("") => true
        @fact truthy([]) => true
        @fact truthy(Dict())  => true
        @fact truthy(nothing) => false
        @fact truthy(false)   => false
        @fact truthy(0)       => false
    end

    context("`anything` is always true") do
        @fact anything(false)   => true
    end

    context("`exactly` can be used to check object equality") do
        @fact exactly(exactly)(exactly) => true

        x() = ()
        @fact exactly(x)(x) => true

        # types with no fields return a singleton object when instantiated
        @fact exactly(Baz())(Baz()) => true

        @fact exactly(Bazz(1))(Bazz(1)) => false
    end

    context("`roughly` compares numbers... roughly") do
        @fact 2.4999999999999 => roughly(2.5)
        @fact 9.5 => roughly(10; atol=1.0)
        @fact 10.5 => roughly(10; atol=1.0)
    end

    context("`roughly` compares matrixes... roughly") do
        X = [1.1 1.2; 2.1 2.2]
        Y = X + [0 0.000001; -0.00000349 0.00001]
        Z = [1 1; 2 2]
        @fact X => roughly(Y)
        @fact X => roughly(Z; atol=0.2)
    end

end

exitstatus()

end # module


text:
struct HeightModel <: Function
    formula::Function
    b
end

##named equations
#2 parameter equations, mainly from the LMFOR R package
#sorted alphabetically
curtis=(x,b)->4.5+b[1]*((x/1+x)^b[2])
michailoff=(x,b)->4.5+b[1]*exp(-b[2]*x^(-1))
meyer=(x,b)->4.5+b[1]*(1-exp(-b[2]*x))
micment=(x,b)->4.5+b[1]*x/(b[2]+x)
micment2=(x,b)-> 4.5+x/(b[1]+b[2]*x)
naslund=(x,b)->4.5+(x^2/(b[1]+b[2]*x)^2)
naslund2=(x,b)-> 4.5+x^2/(b[1]+exp(b[2])*x)^2
naslund3=(x,b)-> 4.5+x^2/(exp(b[1]) + b[2]*x)^2
naslund4=(x,b)-> 4.5+x^2/(exp(b[1]) +exp(b[2])*x)^2
power=(x,b)-> 4.5+b[1]*x^b[2]
wyckoff=(x,b)->4.5+exp(b[1]+(b[2]/(x+1)))

#3 parameter equations, mainly from LMFOR R package
#sorted alphabetically
chapman=(x,b)->4.5+b[1]*(1-exp(-b[2]*x))^b[3]
gompertz=(x,b)->4.5+b[1]*exp(-b[2]*exp(-b[3]*x))
hossfeldIV=(x,b)->4.5+b[1]/(1+1/(b[2]*x^b[3]))
korf=(x,b)->4.5+b[1]*exp(-b[2]*x^-b[3])
logistic=(x,b)->4.5+b[1]/(1+b[2]*exp(-b[3]*x))
monserud=(x,b)->4.5+exp(b[1]+(b[2]*x*b[3]))
prodan=(x,b)->4.5+x^2/(b[1]+b[2]*x+b[3]*x^2)
ratkowsky=(x,b)->4.5+b[1]*exp(-b[2]/(x+b[3]))
sibbesen=(x,b)->4.5+b[1]*x^(b[2]*x^(-b[3]))
weibull=(x,b)->4.5+b[1]*(1-exp(-b[2]*x^b[3]))

function estimate_height(t::Tree,hd::HeightModel)
    hd.formula(t.diameter, hd.b[t.species])
end

function estimate_height!(t::Tree,hd::HeightModel)
    t.height = ismissing(t.height) ? hd.formula(t.diameter, hd.b[t.species]) : t.height
end


text:
using Combinatorics

function solve(n::Vector{<:AbstractString}, pred::Vector{<:Function})
    rst = Vector{typeof(n)}(0)
    for candidate in permutations(n)
        if all(p(candidate) for p in predicates)
            push!(rst, candidate)
        end
    end
    return rst
end

Names = ["Baker", "Cooper", "Fletcher", "Miller", "Smith"]
predicates = [
    (s) -> last(s) != "Baker",
    (s) -> first(s) != "Cooper",
    (s) -> first(s) != "Fletcher" && last(s) != "Fletcher",
    (s) -> findfirst(s, "Miller") > findfirst(s, "Cooper"),
    (s) -> abs(findfirst(s, "Smith") - findfirst(s, "Fletcher")) != 1,
    (s) -> abs(findfirst(s, "Cooper") - findfirst(s, "Fletcher")) != 1]

solutions = solve(Names, predicates)
foreach(x -> println(join(x, ", ")), solutions)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
## Example 1 from paper
# "Lasso Regularization Paths for NARMAX Models via Coordinate Descent"
# Antonio H. Ribeiro and Luis A. Aguirre
using NarmaxLasso
using Plots
plotlyjs()

## Signal generation
# Defining nonlinear system
function armax_mdl(y, u, v)
    return (0.5*y[1] - 0.5*u[1] + 0.5*v[1] )
end
yterms0 = [1]
uterms0 = [1]
vterms0 = [1]

# Generate signals
u = generate_random_input(2000, 5)
y = simulate(u, armax_mdl, yterms0, uterms0, vterms0;
             σv=0.3, seed=1)

## Model Estimation
# Generate all candidate terms
ny = 10
nu = 10
nv = 10
mdl = generate_all(NarmaxLasso.NarmaxRegressors,
                                NarmaxLasso.Monomial,
                                ny, nu, nv, 1)

# Lasso regression
result = narmax_lasso(y, u, mdl; λ_min_ratio=1e-6)

# Plot result
plot(result, xscale=:log10)

## Validation
# Generate signals
uv = generate_random_input(1000, 5, seed=14)
yv = simulate(uv, armax_mdl, yterms0, uterms0, vterms0;
                σv=0.3, seed=11)

# Compute Mean Absolute Erro
n = length(result.λ)
mae = Vector{Float64}(n)
for i = 1:n
    ysv = simulate(uv, result.mdl, result.β[:, i])
    mae[i] = mean(abs(ysv-yv))
end
vline!([result.λ[indmin(mae)]], color=:black, line=:dash)


text:
# TODO add Foldable class to refer to the reduce, foldr, foldl?

module TypeClasses
export neutral, combine, ⊕, reduce_monoid, foldr_monoid, foldl_monoid,
  orelse, ⊘,
  foreach, @syntax_foreach,
  map, @syntax_map,
  pure, ap, curry, mapn, @mapn, tupled, neutral_applicative, combine_applicative, orelse_applicative,
  flatmap, flatten, ↠, @syntax_flatmap,
  flip_types

using Compat
using Reexport
using Requires

using Monadic
export @pure  # a user only needs @pure from Monadi

@reexport using DataTypesBasic

include("Utils/Utils.jl")
using .Utils

include("DataTypes/DataTypes.jl")
@reexport using .DataTypes

# TypeClasses
# ===========
# (they may already include default implementations in terms of other typeclasses)

# we decided for a flat module without submodules to simplify overloading the functions
include("TypeClasses/MonoidAlternative.jl")
include("TypeClasses/FunctorApplicativeMonad.jl")
include("TypeClasses/FlipTypes.jl")  # depends on both Monoid and Applicative
# TODO add Arrow typeclass (composition)

# Instances
# =========

include("TypeInstances/AbstractVector.jl")
include("TypeInstances/Callable.jl")
include("TypeInstances/Dict.jl")
include("TypeInstances/Future.jl")
include("TypeInstances/Iterable.jl")  # only supplies default functions, no actual dispatch is done (Reason: there had been too many conflicts with Dict already)
include("TypeInstances/Monoid.jl")
include("TypeInstances/Pair.jl")
include("TypeInstances/State.jl")
include("TypeInstances/String.jl")
include("TypeInstances/Task.jl")
include("TypeInstances/Tuple.jl")
include("TypeInstances/Writer.jl")

include("TypeInstances_DataTypesBasic/Const.jl")
include("TypeInstances_DataTypesBasic/Identity.jl")
include("TypeInstances_DataTypesBasic/Either.jl")
include("TypeInstances_DataTypesBasic/MultipleExceptions.jl")
include("TypeInstances_DataTypesBasic/ContextManager.jl")

# extra interactions between ContainerTypes
include("convert.jl")

function __init__()
  # there are known issues with require. See this issue for updates https://github.com/JuliaLang/Pkg.jl/issues/1285
  @require Dictionaries="85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4" include("TypeInstances/AbstractDictionary.jl")
end

end # module


text:
module GISSMOReader

import Unicode
import JLD
import InfoZIP
import DelimitedFiles

import HTTP
import JSON3
import JSON

include("../src/database/download_edata.jl")
include("../src/database/helpers.jl")
include("../src/database/BMRB_rest.jl")
include("../src/database/GISSMO_entries.jl")
include("../src/database/metadata.jl")

include("../src/assemble.jl")
include("../src/utils.jl")

end


text:
function sumbyterm(dtm :: AbstractMatrix)
  sum(dtm, 2)[:]
end

function sumbydocument(dtm :: AbstractMatrix)
  sum(dtm, 1)[:]
end


text:
# All the types and methods specific to our data storage hierarchy are here. The root of the image
# data is assumed to have subdirectories in the form "SS_/visitV/tT", where SS is a zero-padded
# subject number, V is a visit number, and T is a trial number.

struct Subject
	dataroot::String
	number::Integer
end

struct Visit
	subject::Subject
	number::Integer
end
Visit(dataroot,s::Integer,number::Integer) = Visit(Subject(dataroot,s),number)

struct Trial
	visit::Visit
	number::Integer
	resultdir::String
end
Trial(v::Visit,number::Integer) = Trial(v,number,".")
Trial(dataroot,s::Integer,v::Integer,number::Integer,resultdir::String=".") = Trial(Visit(Subject(dataroot,s),v),number,resultdir)

#
# Subject methods
#

dataroot(s::Subject) = s.dataroot
number(S::Subject) = S.number
String(s::Subject) = s.number < 10 ? "0$(s.number)" : "$(s.number)"
dirname(s::Subject,full=false) = full ? joinpath(s.dataroot,String(s)*"_") : String(s)*"_"
fullname(s::Subject) = dirname(s,true)
function visits(s::Subject)
	vis = Visit[]
	screen = t->startswith(t,"visit")
	contents = filter(screen,readdir(fullname(s)))
	for v in 1:2
		newvis = Visit(s,v)
		if String(newvis) in contents
			push!(vis,newvis)
		end
	end
	return vis
end
numvisits(s::Subject) = length(visits(s))
trials(s::Subject) = reduce(vcat,trials(v) for v in visits(s))
numtrials(s::Subject) = sum( numtrials(v) for v in visits(s) )
length(s::Subject) = sum( numframes(t) for t in trials(s) )

#
# Visit methods
#

subject(v::Visit) = v.subject
dataroot(v::Visit) = dataroot(v.subject)
number(V::Visit) = [number(V.subject),V.number]
String(v::Visit) = "visit$(v.number)"
dirname(v::Visit,full=false) = joinpath(dirname(v.subject,full),String(v))
fullname(v::Visit) = dirname(v,true)
function trials(v::Visit)
	tri = Trial[]
	screen = t->startswith(t,"t")
	contents = filter(screen,readdir(fullname(v)))
	for t in 1:10
		newtri = Trial(v,t)
		if String(newtri) in contents
			push!(tri,newtri)
		end
	end
	return tri
end
numtrials(v::Visit) = length(trials(v))
length(v::Visit) = sum( numframes(t) for t in trials(v) )

#
# Trial methods
#

visit(t::Trial) = t.visit
subject(t::Trial) = subject(visit(t))
dataroot(t::Trial) = dataroot(t.visit)
number(T::Trial) = append!(number(T.visit),T.number)
String(t::Trial) = "t$(t.number)"
dirname(t::Trial,full=false) = joinpath(dirname(t.visit,full),String(t))
fullname(t::Trial) = dirname(t,true)
function shortname(t::Trial)
	v = visit(t)
	s = subject(v)
	"S$(String(s))_V$(v.number)_T$(t.number)"
end
numframes(t::Trial) = count(isimg,readdir(fullname(t)))
filenames(t::Trial;join=false) = filter(isimg,readdir(fullname(t),join=join))
length(t::Trial) = numframes(t)

#
# Versions of key methods that work with trials
#

function resize(T::Trial,destroot::String,sz=(2824÷2,4240÷2))
	summary = resize(filenames(T,join=true),joinpath(destroot,dirname(T)),sz)
	T2 = Trial(destroot,number(T)...)
	summary = (fname=filenames(T2),summary...)
	CSV.write(joinpath(fullname(T2),"summary.csv"),summary)
	tmp = joinpath(tempname(),"summary.jld2")
	save(tmp,shortname(T2),summary)
	cp(tmp,joinpath(fullname(T2),"summary.jld2"),force=false)
	return T
end

summarize(root,subj,vis,tri) = summarize(Trial(root,subj,vis,tri))
function summarize(T::Trial;dosave=true)
	summary = summarize(filenames(T,join=true))
	summary = (fname=filenames(T),summary...)
	if dosave
		outfile = shortname(T)
		CSV.write(joinpath(dataroot(T),"summaries","$(outfile).csv"),summary)
		tmp = joinpath(tempname(),"$(outfile).jld2")
		save(tmp,shortname(T),summary)
		cp(tmp,joinpath(dataroot(T),"summaries","$(outfile).jld2"),force=false)
	end
	return summary
end

function detect(T::Trial;dosave=true,size=[])
	result = detect(filenames(T,join=true),size)
	result = (fname=filenames(T),result...)
	if dosave
		outfile = shortname(T)
		save("$(outfile).jld2","result",result)
		CSV.write("$(outfile).csv",result)
	end
	return result
end

#
# Accessing results/summary data for a Trial
#

results(root,subject,visit,trial) = results(Trial(root,subject,visit,trial))
results(subject,visit,trial) = results(Trial(".",subject,visit,trial))
#results(t::Trial,dir=t.resultdir) = load(joinpath(dir,shortname(t)*".jld2"))["result"] |> DataFrame
results(t::Trial,dir=t.resultdir) = load(joinpath(dir,shortname(t)*".csv")) |> DataFrame

summary(root,subject,visit,trial) = summary(Trial(root,subject,visit,trial))
function summary(T::Trial)
	sn = shortname(T)
	load(joinpath(dataroot(T),"summaries","$(sn).jld2"))[sn] |> DataFrame
end

function makemovie(T::Trial,sz=(470,706);numframes=Inf)
	folder = filenames(T,join=true)
	imgstack = makemovie(folder,results(T),sz,numframes=numframes)
	encodevideo(shortname(T)*".mp4",imgstack,framerate=10)
	return shortname(T)*".mp4"
end

#
# Utilities
#

intensity(T::Trial) = intensity(summary(T))
darkframes(T::Trial) = darkframes(summary(T))
goodframes(T::Trial) = goodframes(summary(T))

#
# ImageFolder: iterator of images for a given trial
#
struct ImageFolder
	t::Trial
	preprocessor
	file::AbstractVector
end

ImageFolder(t::Trial,prep=missing) = ImageFolder(t,prep,filenames(t,join=true))
images(t::Trial,prep=missing) = ImageFolder(t,prep)
ImageFolder(dataroot,s::Integer,v::Integer,t::Integer) = ImageFolder(Trial(dataroot,s,v,t))

function getindex(img::ImageFolder,n::Integer)
	X = load(img.file[n])
	if !ismissing(img.preprocessor)
		X = img.preprocessor(X)
	end
	return X
end
getindex(f::ImageFolder,inds) = ImageFolder(f.t,f.preprocessor,f.file[inds])

iterate(f::ImageFolder) = isempty(f.file) ? nothing : (getindex(f,1),1)
iterate(f::ImageFolder,state) = state==length(f.file) ? nothing : (getindex(f,state+1),state+1)
length(f::ImageFolder) = length(f.file)
isempty(f::ImageFolder) = isempty(f.file)
IteratorEltype(ImageFolder) = EltypeUnknown()

max_stars_repo_name: tobydriscoll/cornea-detection

max_issues_repo_name: tobydriscoll/cornea-detection

max_forks_repo_name: tobydriscoll/cornea-detection


text:
# Copyright 2020 Google LLC

# 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

#     https://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.

function part1(filename="day15input.txt", nth=2020)
    starting_numbers = map(s->parse(Int, s), split(String(read(filename)), ","))
    last_said = Dict{Int, Int}()
    prev = nothing
    saying = nothing
    for i ∈ 1:nth
        if i <= length(starting_numbers)
            saying = starting_numbers[i]
        elseif prev ∈ keys(last_said)
            saying = i - last_said[prev] - 1
        else
            saying = 0
        end
        if prev !== nothing
            last_said[prev] = i-1
        end
        prev = saying
    end
    return saying
end

println("Part 1: ", part1())
part2() = part1("day15input.txt", 30000000)
println("Part 2: ", part2())


text:
using Test
using Mongoc

include("../../src/CholeskyExperiment.jl")
using ..CholeskyExperiment:ExperimentManager,
                           gendecomposition,
                           gendecomposition_all,
                           gendecomposition_all_json,
                           gendecomposition_mult_json,
                           gendecomposition_all_jsons,
                           genoriginal,
                           genoriginal_all,
                           DecompositionDB.isdecomposition,
                           DecompositionDB.getdecomposition,
                           load_instance_by_paths
include("init_manager.jl")

testdirectory = dirname(@__FILE__)
cd(testdirectory)
files = readdir()
for file in files
    if split(file, ".")[end] == "jl" && file != "runtest.jl" && file != "test_generate.jl"
        include(joinpath(testdirectory, file))
    end
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
function _dot(x::AbstractVector, A::AbstractBandedMatrix, y::AbstractVector)
    require_one_based_indexing(x, y)
    # This should include a message but matches LinearAlgebra.dot
    (axes(x)..., axes(y)...) == axes(A) || throw(DimensionMismatch())
    l,u = bandwidths(A)
    M = size(A,1)
    T = typeof(dot(first(x), first(A), first(y)))
    s = zero(T)
    i₁ = first(eachindex(x))
    x₁ = first(x)
    @inbounds for j in colsupport(y)
        yj = y[j]
        if !iszero(yj)
            temp = zero(adjoint(A[i₁,j]) * x₁)
            @simd for i in max(1,j-u):min(M,j+l)
                temp += adjoint(A[i,j]) * x[i]
            end
            s += dot(temp, yj)
        end
    end
    s
end

dot(x::AbstractVector, A::AbstractBandedMatrix, y::AbstractVector) =
    _dot(x, A, y)

function dot(x::AbstractVector, A::BandedMatrix{<:Any,<:AbstractFill}, y::AbstractVector)
    require_one_based_indexing(x, y)
    l,u = bandwidths(A)
    l == -u || return _dot(x, A, y)

    # This should include a message but matches LinearAlgebra.dot
    (axes(x)..., axes(y)...) == axes(A) || throw(DimensionMismatch())
    T = typeof(dot(first(x), first(A), first(y)))
    Av = unique_value(A.data)
    iszero(Av) && return zero(T)

    M,N = size(A)
    @inbounds begin
        xv = view(x, max(1,l+1):min(M,N+l))
        yv = view(y, max(1,u+1):min(M+u,N))

        Av*dot(xv, yv)
    end
end


text:
"""
        greedy_star2_coloring

Find a coloring of a given input graph such that
no two vertices connected by an edge have the same
color using greedy approach. The number of colors
used may be equal or greater than the chromatic
number `χ(G)` of the graph.

A star coloring is a special type of distance - 1  coloring,
For a coloring to be called a star coloring, it must satisfy
two conditions:

1. every pair of adjacent vertices receives distinct  colors
(a distance-1 coloring)

2. For any vertex v, any color that leads to a two-colored path
involving v and three other vertices  is  impermissible  for  v.
In other words, every path on four vertices uses at least three
colors.

Reference: Gebremedhin AH, Manne F, Pothen A. **What color is your Jacobian? Graph coloring for computing derivatives.** SIAM review. 2005;47(4):629-705.

TODO: add text explaining the difference between star1 and star2
"""
function color_graph(g::Graphs.AbstractGraph, ::GreedyStar2Color)
    v = nv(g)
    colorvec = zeros(Int, v)

    forbidden_colors = zeros(Int, v+1)

    for vertex_i = vertices(g)
        for w in inneighbors(g, vertex_i)
            if colorvec[w] != 0
                forbidden_colors[colorvec[w]] = vertex_i
            end

            for x in inneighbors(g, w)
                if colorvec[x] != 0
                    if colorvec[w] == 0
                        forbidden_colors[colorvec[x]] = vertex_i
                    else
                        if colorvec[x] < colorvec[w]
                            forbidden_colors[colorvec[x]] = vertex_i
                        end
                    end
                end
            end
        end
        colorvec[vertex_i] = find_min_color(forbidden_colors, vertex_i)
    end
    return colorvec
end