text:
export LLJITBuilder, LLJIT, ExecutionSession, JITDylib,
       ThreadSafeModule, ThreadSafeContext, OrcTargetAddress
export TargetMachineBuilder, targetmachinebuilder!, linkinglayercreator!
export mangle, lookup, intern
export ObjectLinkingLayer, register!

include("executionengine/utils.jl")

@checked struct TargetMachineBuilder
    ref::API.LLVMOrcJITTargetMachineBuilderRef
end
Base.unsafe_convert(::Type{API.LLVMOrcJITTargetMachineBuilderRef},
                    tmb::TargetMachineBuilder) = tmb.ref

function TargetMachineBuilder()
    ref = Ref{API.LLVMOrcJITTargetMachineBuilderRef}()
    @check API.LLVMOrcJITTargetMachineBuilderDetectHost(ref)
    TargetMachineBuilder(ref[])
end

function TargetMachineBuilder(tm::TargetMachine)
    tmb = API.LLVMOrcJITTargetMachineBuilderCreateFromTargetMachine(tm)
    TargetMachineBuilder(tmb)
end

function dispose(tmb::TargetMachineBuilder)
    API.LLVMOrcDisposeJITTargetMachineBuilder(tmb)
end

include("executionengine/lljit.jl")

@checked struct ExecutionSession
    ref::API.LLVMOrcExecutionSessionRef
end
Base.unsafe_convert(::Type{API.LLVMOrcExecutionSessionRef}, es::ExecutionSession) = es.ref

function ExecutionSession(lljit::LLJIT)
    es = API.LLVMOrcLLJITGetExecutionSession(lljit)
    ExecutionSession(es)
end

@checked struct ObjectLinkingLayer
    ref::API.LLVMOrcObjectLayerRef
end
Base.unsafe_convert(::Type{API.LLVMOrcObjectLayerRef}, oll::ObjectLinkingLayer) = oll.ref

function ObjectLinkingLayer(es::ExecutionSession)
    ref = LLVM.API.LLVMOrcCreateRTDyldObjectLinkingLayerWithSectionMemoryManager(es)
    ObjectLinkingLayer(ref)
end

function dispose(oll::ObjectLinkingLayer)
    API.LLVMOrcDisposeObjectLayer(oll)
end

function register!(oll::ObjectLinkingLayer, listener::JITEventListener)
    API.LLVMOrcRTDyldObjectLinkingLayerRegisterJITEventListener(oll, listener)
end

mutable struct ObjectLinkingLayerCreator
    cb
end

function ollc_callback(ctx::Ptr{Cvoid}, es::LLVM.API.LLVMOrcExecutionSessionRef, triple::Ptr{Cchar})
    es = ExecutionSession(es)
    triple = Base.unsafe_string(triple)

    ollc = Base.unsafe_pointer_to_objref(ctx)::ObjectLinkingLayerCreator
    oll = ollc.cb(es, triple)::ObjectLinkingLayer
    return oll.ref
end

"""
    linkinglayercreator!(builder::LLJITBuilder, creator::ObjectLinkingLayerCreator)

!!! warning
    The creator object needs to be rooted by the caller for the lifetime of the
    builder argument.
"""
function linkinglayercreator!(builder::LLJITBuilder, creator::ObjectLinkingLayerCreator)
    cb = @cfunction(ollc_callback,
                    LLVM.API.LLVMOrcObjectLayerRef,
                    (Ptr{Cvoid}, LLVM.API.LLVMOrcExecutionSessionRef, Ptr{Cchar}))
    linkinglayercreator!(builder, cb, Base.pointer_from_objref(creator))
end

include("executionengine/ts_module.jl")

@checked struct LLVMSymbol
    ref::API.LLVMOrcSymbolStringPoolEntryRef
end
Base.unsafe_convert(::Type{API.LLVMOrcSymbolStringPoolEntryRef}, sym::LLVMSymbol) = sym.ref
Base.convert(::Type{API.LLVMOrcSymbolStringPoolEntryRef}, sym::LLVMSymbol) = sym.ref

function Base.cconvert(::Type{Cstring}, sym::LLVMSymbol)
    return API.LLVMOrcSymbolStringPoolEntryStr(sym)
end

"""
    string(sym::LLVMSymbol)

Converts a `LLVMSymbol` to a Julia string without releasing the corresponding reference.
It is safe to continue using the symbol.
"""
function Base.string(sym::LLVMSymbol)
    cstr = API.LLVMOrcSymbolStringPoolEntryStr(sym)
    Base.unsafe_string(cstr)
end

"""
   String(sym::LLVMSymbol)

Converts a `LLVMSymbol` to a Julia string and releases the corresponding reference.
Likely you do not want to re-use the passed symbol.
"""
function Base.String(sym::LLVMSymbol)
    str = string(sym)
    release(sym)
    return str
end

function intern(es::ExecutionSession, string)
    entry = API.LLVMOrcExecutionSessionIntern(es, string)
    LLVMSymbol(entry)
end

function release(sym::LLVMSymbol)
    API.LLVMOrcReleaseSymbolStringPoolEntry(sym)
end

function retain(sym::LLVMSymbol)
    API.LLVMOrcRetainSymbolStringPoolEntry(sym)
end

# ORC always uses linker-mangled symbols internally (including for lookups, responsibility object maps, etc).
# IR uses non-linker-mangled names.
# If you're synthesizing IR from a requested-symbols map you'll need to demangle the name.
# Unfortunately we don't have a generic utility for that yet, but on MacOS it just means
# dropping the leading '_' if there is one, or prepending a \01 prefix (see https://llvm.org/docs/LangRef.html#identifiers)

function mangle(lljit::LLJIT, name)
    entry = API.LLVMOrcLLJITMangleAndIntern(lljit, name)
    return LLVMSymbol(entry)
end

@checked struct JITDylib
    ref::API.LLVMOrcJITDylibRef
end
Base.unsafe_convert(::Type{API.LLVMOrcJITDylibRef}, jd::JITDylib) = jd.ref

"""
    JITDylib(lljit::LLJIT)

Get the main JITDylib
"""
function JITDylib(lljit::LLJIT)
    ref = API.LLVMOrcLLJITGetMainJITDylib(lljit)
    JITDylib(ref)
end

"""
    JITDylib(es::ExecutionSession, name; bare=false)

Adds a new JITDylib to the ExecutionSession. The name must be unique and
the `bare=true` no standard platform symbols are made available.
"""
function JITDylib(es::ExecutionSession, name; bare=false)
    if bare
        ref = API.LLVMOrcExecutionSessionCreateBareJITDylib(es, name)
    else
        ref = Ref{API.LLVMOrcJITDylibRef}()
        @check API.LLVMOrcExecutionSessionCreateJITDylib(es, ref, name)
        ref = ref[]
    end
    JITDylib(ref)
end

@checked struct DefinitionGenerator
    ref::API.LLVMOrcDefinitionGeneratorRef
end
Base.unsafe_convert(::Type{API.LLVMOrcDefinitionGeneratorRef}, dg::DefinitionGenerator) = dg.ref

function add!(jd::JITDylib, dg::DefinitionGenerator)
    API.LLVMOrcJITDylibAddGenerator(jd, dg)
end

function CreateDynamicLibrarySearchGeneratorForProcess(prefix)
    ref = Ref{API.LLVMOrcDefinitionGeneratorRef}()
    @check API.LLVMOrcCreateDynamicLibrarySearchGeneratorForProcess(ref, prefix, C_NULL, C_NULL)
    DefinitionGenerator(ref[])
end

# LLVMOrcCreateCustomCAPIDefinitionGenerator(F, Ctx)

function lookup_dylib(es::ExecutionSession, name)
    ref = API.LLVMOrcExecutionSessionGetJITDylibByName(es, name)
    if ref == C_NULL
        return nothing
    end
    JITDylib(ref)
end

function add!(lljit::LLJIT, jd::JITDylib, obj::MemoryBuffer)
    @check API.LLVMOrcLLJITAddObjectFile(lljit, jd, obj)
    return nothing
end

# LLVMOrcLLJITAddObjectFileWithRT(J, RT, ObjBuffer)

function add!(lljit::LLJIT, jd::JITDylib, mod::ThreadSafeModule)
    @check API.LLVMOrcLLJITAddLLVMIRModule(lljit, jd, mod)
    return nothing
end

# LLVMOrcLLJITAddLLVMIRModuleWithRT(J, JD, TSM)

function Base.empty!(jd::JITDylib)
    API.LLVMOrcJITDylibClear(jd)
end

struct OrcTargetAddress
    ptr::API.LLVMOrcJITTargetAddress
end
Base.convert(::Type{API.LLVMOrcJITTargetAddress}, addr::OrcTargetAddress) = addr.ptr

Base.pointer(addr::OrcTargetAddress) = reinterpret(Ptr{Cvoid}, addr.ptr % UInt) # LLVMOrcTargetAddress is UInt64 even on 32-bit

OrcTargetAddress(ptr::Ptr{Cvoid}) = OrcTargetAddress(reinterpret(UInt, ptr))

"""
    lookup(lljit::LLJIT, name)

Takes an unmangled symbol names and searches for it in the LLJIT.
"""
function lookup(lljit::LLJIT, name)
    result = Ref{API.LLVMOrcJITTargetAddress}()
    @check API.LLVMOrcLLJITLookup(lljit, result, name)
    OrcTargetAddress(result[])
end

@checked struct IRTransformLayer
    ref::API.LLVMOrcIRTransformLayerRef
end
Base.unsafe_convert(::Type{API.LLVMOrcIRTransformLayerRef}, il::IRTransformLayer) = il.ref

function IRTransformLayer(lljit::LLJIT)
    ref = API.LLVMOrcLLJITGetIRTransformLayer(lljit)
    IRTransformLayer(ref)
end

function set_transform!(il::IRTransformLayer)
    API.LLVMOrcIRTransformLayerSetTransform(il)
end

@checked struct MaterializationResponsibility
    ref::API.LLVMOrcMaterializationResponsibilityRef
end
Base.unsafe_convert(::Type{API.LLVMOrcMaterializationResponsibilityRef}, mr::MaterializationResponsibility) = mr.ref

function emit(il::IRTransformLayer, mr::MaterializationResponsibility, tsm::ThreadSafeModule)
    API.LLVMOrcIRTransformLayerEmit(il, mr, tsm)
end

function get_requested_symbols(mr::MaterializationResponsibility)
    N = Ref{Csize_t}()
    ptr = API.LLVMOrcMaterializationResponsibilityGetRequestedSymbols(mr, N)
    syms = map(LLVMSymbol, Base.unsafe_wrap(Array, ptr, N[], own=false))
    API.LLVMOrcDisposeSymbols(ptr)
    return syms
end

abstract type AbstractMaterializationUnit end

function define(jd::JITDylib, mu::AbstractMaterializationUnit)
    API.LLVMOrcJITDylibDefine(jd, mu)
end

@checked struct MaterializationUnit <: AbstractMaterializationUnit
    ref::API.LLVMOrcMaterializationUnitRef
end
Base.unsafe_convert(::Type{API.LLVMOrcMaterializationUnitRef}, mu::MaterializationUnit) = mu.ref

mutable struct CustomMaterializationUnit <: AbstractMaterializationUnit
    materialize
    discard
    mu::MaterializationUnit
    function CustomMaterializationUnit(materialize, discard)
        new(materialize, discard)
    end
end
Base.cconvert(::Type{API.LLVMOrcMaterializationUnitRef}, mu::CustomMaterializationUnit) = mu.mu

const CUSTOM_MU_ROOTS = Base.IdSet{CustomMaterializationUnit}()

function __materialize(ctx::Ptr{Cvoid}, mr::API.LLVMOrcMaterializationResponsibilityRef)
    mu = Base.unsafe_pointer_to_objref(ctx)::CustomMaterializationUnit
    try
        mu.materialize(MaterializationResponsibility(mr))
    catch err
        bt = catch_backtrace()
        showerror(stderr, err, bt)
        API.LLVMOrcMaterializationResponsibilityFailMaterialization(mr)
    end
    nothing
end

function __discard(ctx::Ptr{Cvoid}, jd::API.LLVMOrcJITDylibRef, symbol::API.LLVMOrcSymbolStringPoolEntryRef)
    mu = Base.unsafe_pointer_to_objref(ctx)::CustomMaterializationUnit
    try
        mu.discard(JITDylib(jd), LLVMSymbol(symbol))
    catch err
        bt = catch_backtrace()
        showerror(stderr, err, bt)
    end
    nothing
end

function __destroy(ctx::Ptr{Cvoid})
    mu = Base.unsafe_pointer_to_objref(ctx)::CustomMaterializationUnit
    delete!(CUSTOM_MU_ROOTS, mu)
    nothing
end

function CustomMaterializationUnit(name, symbols, materialize, discard, init=C_NULL)
    this = CustomMaterializationUnit(materialize, discard)
    push!(CUSTOM_MU_ROOTS, this)

    ref = API.LLVMOrcCreateCustomMaterializationUnit(
        name,
        Base.pointer_from_objref(this), # escaping this, rooted in CUSTOM_MU_ROOTS
        symbols,
        length(symbols),
        init,
        @cfunction(__materialize, Cvoid, (Ptr{Cvoid}, API.LLVMOrcMaterializationResponsibilityRef)),
        @cfunction(__discard, Cvoid, (Ptr{Cvoid}, API.LLVMOrcJITDylibRef, API.LLVMOrcSymbolStringPoolEntryRef) ),
        @cfunction(__destroy, Cvoid, (Ptr{Cvoid},))
    )
    this.mu = MaterializationUnit(ref)
    return this
end

function absolute_symbols(symbols)
    ref = API.LLVMOrcAbsoluteSymbols(symbols, length(symbols))
    MaterializationUnit(ref)
end

@checked struct IndirectStubsManager
    ref::API.LLVMOrcIndirectStubsManagerRef
end
Base.unsafe_convert(::Type{API.LLVMOrcIndirectStubsManagerRef}, ism::IndirectStubsManager) = ism.ref

function LocalIndirectStubsManager(triple)
    ref = API.LLVMOrcCreateLocalIndirectStubsManager(triple)
    IndirectStubsManager(ref)
end

function dispose(ism::IndirectStubsManager)
    API.LLVMOrcDisposeIndirectStubsManager(ism)
end

@checked mutable struct LazyCallThroughManager
    ref::API.LLVMOrcLazyCallThroughManagerRef
end
Base.unsafe_convert(::Type{API.LLVMOrcLazyCallThroughManagerRef}, lcm::LazyCallThroughManager) = lcm.ref

function LocalLazyCallThroughManager(triple, es)
    ref = Ref{API.LLVMOrcLazyCallThroughManagerRef}()
    @check API.LLVMOrcCreateLocalLazyCallThroughManager(triple, es, C_NULL, ref)
    LazyCallThroughManager(ref[])
end

function dispose(lcm::LazyCallThroughManager)
    API.LLVMOrcDisposeLazyCallThroughManager(lcm)
end

function reexports(lctm::LazyCallThroughManager, ism::IndirectStubsManager, jd::JITDylib, symbols)
    ref = LLVM.API.LLVMOrcLazyReexports(lctm, ism, jd, symbols, length(symbols))
    MaterializationUnit(ref)
end


text:
## TODO Use Unitful to have currency weekly monthly annual counts as annotations

@with_kw mutable struct IncomeTaxSys
   earnings_rates :: RateBands = [19.0,20.0,21.0,41.0,46.0]
   earnings_bands :: RateBands = [2_049.0, 12_444.0, 30_930.0, 150_000.0]
   savings_rates  :: RateBands = [0.0, 20.0, 40.0, 45.0]
   savings_bands  :: RateBands = [5_000.0, 37_500.0, 150_000.0]
   dividend_rates :: RateBands = [0.0, 7.5,32.5,38.1]
   dividend_bands :: RateBands = [2_000.0, 37_500.0, 150_000.0]
   personal_allowance          = 12_500.00
   blind_persons_allowance     = 2_450.00
   married_couples_allowance   = 8_915.00
   mca_minimum                 = 3_450.00
   marriage_allowance          = 1_250.00
   personal_savings_allowance  = 1_000.00
end

function annualise!( it :: IncomeTaxSys )
   it.earnings_rates .*= 100.0
   it.earnings_bands .*= 100.0
   it.savings_rates .*= 100.0
   it.savings_bands .*= 100.0
   it.dividend_rates .*= 100.0
   it.earnings_bands .*= WEEKS_PER_YEAR
   it.savings_bands .*= WEEKS_PER_YEAR
   it.dividend_bands .*= WEEKS_PER_YEAR
   it.personal_allowance *= WEEKS_PER_YEAR
   it.blind_persons_allowance *= WEEKS_PER_YEAR
   it.married_couples_allowance *= WEEKS_PER_YEAR
   it.mca_minimum *= WEEKS_PER_YEAR
   it.marriage_allowance *= WEEKS_PER_YEAR
   it.personal_savings_allowance *= WEEKS_PER_YEAR
end

function weeklyise!( it :: IncomeTaxSys )

   it.earnings_rates ./= 100.0
   it.earnings_bands ./= 100.0
   it.savings_rates ./= 100.0
   it.savings_bands ./= 100.0
   it.dividend_rates ./= 100.0
   it.earnings_bands ./= WEEKS_PER_YEAR
   it.savings_bands ./= WEEKS_PER_YEAR
   it.dividend_bands ./= WEEKS_PER_YEAR
   it.personal_allowance /= WEEKS_PER_YEAR
   it.blind_persons_allowance /= WEEKS_PER_YEAR
   it.married_couples_allowance /= WEEKS_PER_YEAR
   it.mca_minimum /= WEEKS_PER_YEAR
   it.marriage_allowance /= WEEKS_PER_YEAR
   it.personal_savings_allowance /= WEEKS_PER_YEAR
end

"""
Map from
"""
function fromJSON( json :: Dict ) :: IncomeTaxSys
   it.earnings_rates = convert( RateBands, json["earnings_rates"] )
   it.earnings_bands  = convert( RateBands, json["earnings_bands"] )
   it.savings_rates = convert( RateBands, json["savings_rates"] )
   it.savings_bands = convert( RateBands, json["savings_bands"] )
   it.dividend_rates = convert( RateBands ,json["dividend_rates"] )
   it.earnings_bands = convert( RateBands ,json["earnings_bands"] )
   it.savings_bands = convert( RateBands ,json["savings_bands"] )
   it.dividend_bands = convert( RateBands ,json["dividend_bands"] )
   it.personal_allowance = json["personal_allowance"]
   it.blind_persons_allowance = json["blind_persons_allowance"]
   it.married_couples_allowance = json["married_couples_allowance"]
   it.mca_minimum = json["mca_minimum"]
   it.marriage_allowance = json["marriage_allowance"]
   it.personal_savings_allowance = json["personal_savings_allowance"]
   it
end

@with_kw mutable struct NationalInsuranceSys
   rates :: RateBands = [1.0]
end

@with_kw mutable struct TaxBenefitSystem
   name :: AbstractString = "Scotland 2919/20"
   it   :: IncomeTaxSys = IncomeTaxSys()
   ni   :: NationalInsuranceSys = NationalInsuranceSys()
end

function save( filename :: AbstractString, sys :: TaxBenefitSystem )
    JSON.print( filename, sys )
end

include( "../default_params/default2019_20.jl")
defsys = load()

ext: jl


text:
"""
Banking holidays for England and Wales.
"""
struct UKSettlement <: HolidayCalendar end
const UnitedKingdom = UKSettlement

# England and Wales Banking Holidays
function isholiday(::UKSettlement, dt::Dates.Date)

    yy = Dates.year(dt)
    mm = Dates.month(dt)
    dd = Dates.day(dt)

    # Bisection
    if mm >= 8
        # Fixed holidays
        if (
                # Late Summer Bank Holiday, August Bank Holiday
                adjustweekendholidayPost( findweekday(Dates.Monday, yy, 8, 1, false) ) == dt
                ||
                # Christmas
                adjustweekendholidayPost( Dates.Date(yy, 12, 25) ) == dt
                ||
                # Boxing
                adjustweekendholidayPost(adjustweekendholidayPost( Dates.Date(yy, 12, 25) ) + Dates.Day(1)) == dt
            )
            return true
        end

        # Fixed date holidays with mm >= 8
        if dt == Dates.Date(1999, 12, 31)
            return true
        end
    else
        # mm < 8
        # Fixed holidays
        if (
                # New Year's Day
                adjustweekendholidayPost( Dates.Date(yy, 01, 01) ) == dt
                ||
                # May Day, Early May Bank Holiday
                (adjustweekendholidayPost(findweekday(Dates.Monday, yy, 5, 1, true)) == dt && yy != 1995 && yy != 2020)
                ||
                # Spring Bank Holiday
                (adjustweekendholidayPost(findweekday(Dates.Monday, yy, 5, 1, false)) == dt && yy != 2012 && yy != 2002)
            )
            return true
        end

        # Easter occurs up to April, which is before August (mm < 8). See test/easter-min-max.jl .
        # Holidays based on easter date
        dt_rata::Int = Dates.days(dt)
        e_rata::Int = easter_rata(Dates.Year(yy))

        if (
                # Good Friday
                ( dt_rata == ( e_rata - 2 ) )
                ||
                # Easter Monday
                ( dt_rata == ( e_rata + 1 ) )
            )
            return true
        end

        # Fixed date holidays with mm < 8
        if (
                # Substitute date for Spring Bank Holiday
                (dt == Dates.Date(2012, 06, 04))
                ||
                # Diamond Jubilee of Queen Elizabeth II.
                (dt == Dates.Date(2012, 06, 05))
                ||
                # Golden Jubilee of Queen Elizabeth II.
                (dt == Dates.Date(2002, 06, 03))
                ||
                # Substitute date for Spring Bank Holiday
                (dt == Dates.Date(2002, 06, 04))
                ||
                # Wedding of Prince William and Catherine Middleton
                (dt == Dates.Date(2011, 04, 29))
                ||
                # Substitute date for Early May Bank Holiday in 1995
                (dt == Dates.Date(1995, 05, 08))
                ||
                # Substitute date for Early May Bank Holiday in 2020
                (dt == Dates.Date(2020, 05, 08))
            )
            return true
        end
    end

    return false
end


text:
mutable struct Config{S, T, U}
    population_size::Int
    num_generations::Int
    p_crossover::Float64
    p_mutation::Float64
    mutation_step::Float64
    min_depth::Int
    max_depth::Int
    tournament_size::Int
    pset::PrimitiveSet
    create_rand_tree::S
    create_population::T
    fitness::U

    p_mate::Float64
end

function Config(;population_size,
                 num_generations,
                 p_crossover,
                 p_mutation,
                 mutation_step,
                 min_depth,
                 max_depth,
                 tournament_size,
                 pset,
                 create_rand_tree,
                 create_population,
                 fitness,
                 p_mate = 0.25
                 )
    Config(population_size, num_generations, p_crossover, p_mutation, mutation_step, min_depth, max_depth, tournament_size, pset, create_rand_tree, create_population, fitness, p_mate)
end

max_stars_repo_name: kogad/JuliaGSGP

max_issues_repo_name: kogad/JuliaGSGP

max_forks_repo_name: kogad/JuliaGSGP


text:
using FromFile
@from "Core.jl" import Node, Options
@from "EquationUtils.jl" import countNodes

# Check if any binary operator are overly complex
function flagBinOperatorComplexity(tree::Node, ::Val{op}, options::Options)::Bool where {op}
    if tree.degree == 0
        return false
    elseif tree.degree == 1
        return flagBinOperatorComplexity(tree.l, Val(op), options)
    else
        if tree.op == op
            overly_complex::Bool = (
                    ((options.bin_constraints[op][1]::Int > -1) &&
                     (countNodes(tree.l) > options.bin_constraints[op][1]::Int))
                      ||
                    ((options.bin_constraints[op][2]::Int > -1) &&
                     (countNodes(tree.r) > options.bin_constraints[op][2]::Int))
                )
            if overly_complex
                return true
            end
        end
        return (flagBinOperatorComplexity(tree.l, Val(op), options) || flagBinOperatorComplexity(tree.r, Val(op), options))
    end
end

# Check if any unary operators are overly complex
function flagUnaOperatorComplexity(tree::Node, ::Val{op}, options::Options)::Bool where {op}
    if tree.degree == 0
        return false
    elseif tree.degree == 1
        if tree.op == op
            overly_complex::Bool = (
                      (options.una_constraints[op]::Int > -1) &&
                      (countNodes(tree.l) > options.una_constraints[op]::Int)
                )
            if overly_complex
                return true
            end
        end
        return flagUnaOperatorComplexity(tree.l, Val(op), options)
    else
        return (flagUnaOperatorComplexity(tree.l, Val(op), options) || flagUnaOperatorComplexity(tree.r, Val(op), options))
    end
end

"""Check if user-passed constraints are violated or not"""
function check_constraints(tree::Node, options::Options, maxsize::Int)::Bool
    if countNodes(tree) > maxsize
        return false
    end
    for i=1:options.nbin
        if options.bin_constraints[i] == (-1, -1)
            continue
        elseif flagBinOperatorComplexity(tree, Val(i), options)
            return false
        end
    end
    for i=1:options.nuna
        if options.una_constraints[i] == -1
            continue
        elseif flagUnaOperatorComplexity(tree, Val(i), options)
            return false
        end
    end
    return true
end

function check_constraints(tree::Node, options::Options)::Bool
    check_constraints(tree, options, options.maxsize)
end


text:
# Create an image
x = LinRange(-3, 3, 150)
y = LinRange(-2, 2, 100)
data = (sin.(exp.(x'.^2 .+ y.^2)) .+ sin.(x') .+ cos.(y) .+ 3) ./ 6
# Draw the image as a color scale
imshow(data)


text:
using RLInterface
using POMDPModels
using Test
using Random

include("zmq.jl")

function sim(env, nsteps=100, rng=MersenneTwister(0))
    o = reset!(env)
    step = 1
    done = false
    r_tot = 0.0
    na = length(actions(env))
    dims = obs_dimensions(env)
    while !done && step <= nsteps
        action = sample_action(env)
        obs, rew, done, info = step!(env, action)
        r_tot += rew
        step += 1
    end
    return r_tot
end

@testset "Sim" begin
    envs = [MDPEnvironment(SimpleGridWorld()),
            MDPEnvironment(InvertedPendulum()),
            MDPEnvironment(MountainCar()),
            POMDPEnvironment(TMaze()),
            POMDPEnvironment(BabyPOMDP()),
            POMDPEnvironment(TigerPOMDP()),
            KMarkovEnvironment(TMaze(), k=4),
            KMarkovEnvironment(BabyPOMDP(), k=4),
            KMarkovEnvironment(TigerPOMDP(), k=4)]

    for env in envs
        r = sim(env)
        process(env)
    end
end

if VERSION >= v"1.1"
    @testset "type stability" begin
        env = MDPEnvironment(SimpleGridWorld())
        @inferred reset!(env)
        @inferred step!(env, :up)
        @inferred sim(env)

        env = POMDPEnvironment(TigerPOMDP())
        @inferred reset!(env)
        @inferred step!(env, 0)
        @inferred sim(env)

        env = KMarkovEnvironment(TigerPOMDP(), k=4)
        @inferred reset!(env)
        @inferred step!(env, 0)
        @inferred sim(env)
    end
end


text:
module TestCBF

import VecMathOptInterface
using Test

const MOI = VecMathOptInterface
const MOIU = MOI.Utilities
const CBF = MOI.FileFormats.CBF
const CBF_TEST_FILE = "test.cbf"
const MODELS_DIR = joinpath(@__DIR__, "models")

function _set_var_and_con_names(model::MOI.ModelLike)
    variable_names = String[]
    for j in MOI.get(model, MOI.ListOfVariableIndices())
        var_name_j = "v" * string(j.value)
        push!(variable_names, var_name_j)
        MOI.set(model, MOI.VariableName(), j, var_name_j)
    end

    idx = 0
    constraint_names = String[]
    single_variable_constraints = Tuple[]
    for i in MOI.get(
        model,
        MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.Integer}(),
    )
        idx += 1
        x = MOI.get(model, MOI.VariableName(), MOI.VariableIndex(i.value))
        push!(single_variable_constraints, (x, MOI.Integer()))
    end
    for S in [
            MOI.Reals,
            MOI.Zeros,
            MOI.Nonnegatives,
            MOI.Nonpositives,
            MOI.SecondOrderCone,
            MOI.RotatedSecondOrderCone,
            MOI.PositiveSemidefiniteConeTriangle,
            MOI.ExponentialCone,
            MOI.DualExponentialCone,
            MOI.PowerCone{Float64},
            MOI.DualPowerCone{Float64},
        ],
        F in [MOI.VectorOfVariables, MOI.VectorAffineFunction{Float64}]

        for i in MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
            idx += 1
            con_name_i = "c" * string(idx)
            push!(constraint_names, con_name_i)
            MOI.set(model, MOI.ConstraintName(), i, con_name_i)
        end
    end

    return variable_names, constraint_names, single_variable_constraints
end

function Base.isapprox(
    f::MOI.VectorOfVariables,
    g::MOI.VectorAffineFunction{Float64};
    kwargs...,
)
    return isapprox(MOI.VectorAffineFunction{Float64}(f), g)
end

function _test_write_then_read(model_string::String)
    model1 = CBF.Model()
    MOIU.loadfromstring!(model1, model_string)
    args = _set_var_and_con_names(model1)

    MOI.write_to_file(model1, CBF_TEST_FILE)
    model2 = CBF.Model()
    MOI.read_from_file(model2, CBF_TEST_FILE)
    _set_var_and_con_names(model2)

    return MOI.Test.util_test_models_equal(model1, model2, args...)
end

function _test_read(filename::String, model_string::String)
    model1 = CBF.Model()
    MOIU.loadfromstring!(model1, model_string)
    args = _set_var_and_con_names(model1)

    model2 = CBF.Model()
    MOI.read_from_file(model2, filename)
    _set_var_and_con_names(model2)

    return MOI.Test.util_test_models_equal(model1, model2, args...)
end

function test_show()
    @test sprint(show, CBF.Model()) == "A Conic Benchmark Format (CBF) model"
end

function test_support_errors()
    for set in [
        MOI.EqualTo(1.0),
        MOI.LessThan(1.0),
        MOI.GreaterThan(1.0),
        MOI.Interval(1.0, 2.0),
        MOI.Semiinteger(1.0, 2.0),
        MOI.Semicontinuous(1.0, 2.0),
        MOI.ZeroOne(),
    ]
        model_string = """
        variables: x
        minobjective: x
        x in $set
        """
        model = CBF.Model()
        err = MOI.UnsupportedConstraint{MOI.VariableIndex,typeof(set)}
        @test_throws err MOIU.loadfromstring!(model, model_string)
    end
end

function test_read_nonempty()
    model = CBF.Model()
    MOI.add_variable(model)
    @test_throws(
        ErrorException("Cannot read in file because model is not empty."),
        MOI.read_from_file(model, joinpath(MODELS_DIR, "example_A.cbf")),
    )
end

function test_read_incompatible()
    model = CBF.Model()
    @test_throws Exception MOI.read_from_file(
        model,
        joinpath(MODELS_DIR, "incompatible_version.cbf"),
    )
end

function test_read_badcones()
    for filename in [
        "bad_cone_string_A.cbf",
        "bad_cone_string_B.cbf",
        "bad_cone_string_C.cbf",
        "bad_cone_string_D.cbf",
    ]
        model = CBF.Model()
        @test_throws Exception MOI.read_from_file(
            model,
            joinpath(MODELS_DIR, filename),
        )
    end
end

function test_read_badpowerdim()
    for filename in ["bad_power_dim_A.cbf", "bad_power_dim_B.cbf"]
        model = CBF.Model()
        @test_throws Exception MOI.read_from_file(
            model,
            joinpath(MODELS_DIR, filename),
        )
    end
end

function test_read_corrupt()
    for filename in [
        "corrupt_line_A.cbf",
        "corrupt_line_B.cbf",
        "corrupt_line_C.cbf",
        "corrupt_line_D.cbf",
        "corrupt_line_E.cbf",
    ]
        model = CBF.Model()
        @test_throws Exception MOI.read_from_file(
            model,
            joinpath(MODELS_DIR, filename),
        )
    end
end

const _WRITE_READ_MODELS = [
    (
        "min VariableIndex",
        """
    variables: x
    minobjective: x
""",
    ),
    (
        "min ScalarAffine",
        """
    variables: x, y
    minobjective: 1.2x + -1y + 1
""",
    ),
    (
        "max VariableIndex",
        """
    variables: x
    maxobjective: x
""",
    ),
    (
        "max ScalarAffine",
        """
    variables: x, y
    maxobjective: 1.2x + -1y + 1
""",
    ),
    (
        "VariableIndex in Integer",
        """
    variables: x, y
    minobjective: 1.2x
    y in Integer()
""",
    ),
    (
        "VectorOfVariables in Zeros",
        """
    variables: x, y
    minobjective: x
    c1: [x, y] in Zeros(2)
""",
    ),
    (
        "VectorOfVariables in Reals",
        """
    variables: x, y
    minobjective: x
    c1: [x, y] in Reals(2)
""",
    ),
    (
        "VectorOfVariables in Nonnegatives",
        """
    variables: x, y
    minobjective: x
    c1: [x, y] in Nonnegatives(2)
""",
    ),
    (
        "VectorOfVariables in Nonpositives",
        """
    variables: x, y
    minobjective: x
    c1: [y, x] in Nonpositives(2)
""",
    ),
    (
        "VectorAffineFunction in Zeros",
        """
    variables: x, y
    minobjective: 1.2x
    c1: [x + 2y + -1.1, 0] in Zeros(2)
""",
    ),
    (
        "VectorAffineFunction in Reals",
        """
    variables: x, y
    minobjective: 1.2x
    c1: [1x, 2y] in Reals(2)
""",
    ),
    (
        "VectorAffineFunction in Nonnegatives",
        """
    variables: x, y
    minobjective: 1.2x
    c1: [1.1 * x, y + 1] in Nonnegatives(2)
""",
    ),
    (
        "VectorAffineFunction in Nonpositives",
        """
    variables: x
    minobjective: 1.2x
    c1: [-1.1 * x + 1] in Nonpositives(1)
""",
    ),
    (
        "VectorOfVariables in SecondOrderCone",
        """
    variables: x, y, z
    minobjective: x
    c1: [x, y, z] in SecondOrderCone(3)
""",
    ),
    (
        "VectorOfVariables in RotatedSecondOrderCone",
        """
    variables: x, y, z
    minobjective: x
    c1: [x, y, z] in RotatedSecondOrderCone(3)
""",
    ),
    (
        "VectorOfVariables in ExponentialCone",
        """
    variables: x, y, z
    minobjective: x
    c1: [x, y, z] in ExponentialCone()
""",
    ),
    (
        "VectorOfVariables in DualExponentialCone",
        """
    variables: x, y, z
    minobjective: x
    c1: [x, y, z] in DualExponentialCone()
""",
    ),
    (
        "VectorOfVariables in PowerCone",
        """
    variables: x, y, z
    minobjective: x
    c1: [x, y, z] in PowerCone(2.0)
""",
    ),
    (
        "VectorOfVariables in DualPowerCone",
        """
    variables: x, y, z
    minobjective: x
    c1: [x, y, z] in DualPowerCone(2.0)
""",
    ),
    (
        "VectorOfVariables in PositiveSemidefiniteConeTriangle",
        """
    variables: x, y, z
    minobjective: x
    c1: [x, y, z] in PositiveSemidefiniteConeTriangle(2)
""",
    ),
    (
        "VectorAffineFunction in SecondOrderCone",
        """
    variables: x, y, z
    minobjective: 1.2x
    c1: [1.1x, y + 1, 2x + z] in SecondOrderCone(3)
""",
    ),
    (
        "VectorAffineFunction in RotatedSecondOrderCone",
        """
    variables: x, y, z
    minobjective: 1.2x
    c1: [1.1x, y + 1, 2x + z] in RotatedSecondOrderCone(3)
""",
    ),
    (
        "VectorAffineFunction in ExponentialCone",
        """
    variables: x, y, z
    minobjective: 1.2x
    c1: [1.1x, y + 1, 2x + z] in ExponentialCone()
""",
    ),
    (
        "VectorAffineFunction in DualExponentialCone",
        """
    variables: x, y, z
    minobjective: 1.2x
    c1: [1.1x, y + 1, 2x + z] in DualExponentialCone()
""",
    ),
    (
        "VectorAffineFunction in PowerCone",
        """
    variables: x, y, z
    minobjective: 1.2x
    c1: [1.1x, y + 1, 2x + z] in PowerCone(2.0)
""",
    ),
    (
        "VectorAffineFunction in DualPowerCone",
        """
    variables: x, y, z
    minobjective: 1.2x
    c1: [1.1x, y + 1, 2x + z] in DualPowerCone(2.0)
""",
    ),
    (
        "VectorAffineFunction in PositiveSemidefiniteConeTriangle",
        """
    variables: x, y, z
    minobjective: 1.2x
    c1: [1.1x, y + 1, 2x + z] in PositiveSemidefiniteConeTriangle(2)
""",
    ),
]

function test_write_read_models()
    for (model_name, model_string) in _WRITE_READ_MODELS
        _test_write_then_read(model_string)
    end
end

const _EXAMPLE_MODELS = [
    (
        "example_A.cbf",
        """
    variables: U, V, W, X, Y, Z, x, y, z
    minobjective: y + 2U + 2V + 2W + 2Y + 2Z
    c1: [U, V, W, X, Y, Z] in PositiveSemidefiniteConeTriangle(3)
    c2: [y + U + W + Z + -1, x + z + U + 2V + W + 2X + 2Y + Z + -0.5] in Zeros(2)
    c3: [y, x, z] in SecondOrderCone(3)
""",
    ),
    (
        "example_B.cbf",
        """
    variables: X, Y, Z, x, y
    minobjective: 1 + x + y + X + Z
    c1: [X, Y, Z] in PositiveSemidefiniteConeTriangle(2)
    c2: [2Y + -1x + -1y] in Nonnegatives(1)
    c3: [3y + -1, x + y, 3x + -1] in PositiveSemidefiniteConeTriangle(2)
""",
    ),
    (
        "example_C.cbf",
        """
    variables: a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p
    maxobjective: a + b + c + d + e + f + g + h + i + j + k + l + m + n + o + p + -1
    c1: [b] in Zeros(1)
    c2: [c] in Nonnegatives(1)
    c3: [d] in Nonpositives(1)
    c4: [e, f, g] in SecondOrderCone(3)
    c5: [h, i, j] in RotatedSecondOrderCone(3)
    c6: [m, l, k] in ExponentialCone()
    c7: [p, o, n] in DualExponentialCone()
""",
    ),
    (
        "example_D.cbf",
        """
    variables: u, v, w, x, y, z
    maxobjective: w + z
    c1: [u, v, w] in PowerCone(0.5)
    c2: [x, y, z] in DualPowerCone(0.1)
    c3: [1, u, u + v] in PowerCone(0.8)
    c4: [1, y, x + y] in DualPowerCone(0.75)
    u in Integer()
    w in Integer()
    y in Integer()
""",
    ),
]

function test_example_models()
    for (model_name, model_string) in _EXAMPLE_MODELS
        _test_read(joinpath(MODELS_DIR, model_name), model_string)
        _test_write_then_read(model_string)
    end
end

function runtests()
    for name in names(@__MODULE__, all = true)
        if startswith("$(name)", "test_")
            @testset "$name" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    sleep(1.0)  # Allow time for unlink to happen.
    rm(CBF_TEST_FILE, force = true)
    return
end

end

TestCBF.runtests()


text:
using Test

using MemoryBasedCF
using SparseArrays

R = sparse([
    2.0  1.0  0.0
    1.0  0.0  2.0
    0.0  1.0  2.0
    0.0  0.0  2.0
])

@testset "default" begin
    m = memorize(R)

    @test m.nu == 4
    @test m.ni == 3
    @test m.bu == [1.5, 1.5, 1.5, 2.0]
    @test m.bi == [1.5, 1.0, 2.0]
    @test m.Dui == sparse([
         0.5  -0.5   0.0
        -0.5   0.0   0.5
         0.0  -0.5   0.5
         0.0   0.0   0.0
    ])
    @test m.Diu == sparse([
         0.5  -0.5   0.0   0.0
         0.0   0.0   0.0   0.0
         0.0   0.0   0.0   0.0
    ])
    @test isapprox(m.Sii, sparse([
         1.0  -0.5  -0.5
        -0.5   1.0  -0.5
        -0.5  -0.5   1.0
    ]), atol = 1e-2)
    @test isapprox(m.Suu, sparse([
         1.0  -1.0   0.0   0.0
        -1.0   1.0   0.0   0.0
         0.0   0.0   0.0   0.0
         0.0   0.0   0.0   0.0
    ]), atol = 1e-2)

    @testset "itembased" begin
        expected_scores = [
            2.0  1.0  1.5
            1.0  1.5  2.0
            1.5  1.0  2.0
            2.0  2.0  2.0
        ]
        expected_ranked_items = [
            1  3  2
            3  2  1
            3  1  2
            1  2  3
        ]
        expected_ranked_scores = [
            2.0  1.5  1.0
            2.0  1.5  1.0
            2.0  1.5  1.0
            2.0  2.0  2.0
        ]

        @test isapprox(itembased_scores(m, [1, 2, 3, 4]), expected_scores, atol = 1e-2)
        @test isapprox(itembased_scores(m, [1, 2], [1, 3]), expected_scores[[1, 2], [1, 3]], atol = 1e-2)

        ranked_items, ranked_scores = itembased_rankings(m, 3, [1, 2, 3, 4])
        @test ranked_items == expected_ranked_items
        @test isapprox(ranked_scores, expected_ranked_scores, atol = 1e-2)

        ranked_items, ranked_scores = itembased_rankings(m, 2, [1, 2, 3, 4])
        @test ranked_items == expected_ranked_items[:, 1:2]
        @test isapprox(ranked_scores, expected_ranked_scores[:, 1:2], atol = 1e-2)

        ranked_items, ranked_scores = itembased_rankings(m, 2, [1, 2], [1, 3])
        @test ranked_items == [1 3; 3 1]
        @test isapprox(ranked_scores, [2.0 1.5; 2.0 1.0], atol = 1e-2)
    end

    @testset "userbased" begin
        expected_scores = [
            2.0  1.0  2.0
            1.0  1.0  2.0
            1.5  1.0  2.0
            1.5  1.0  2.0
        ]
        expected_ranked_items = [
            1  3  2
            3  1  2
            3  1  2
            3  1  2
        ]
        expected_ranked_scores = [
            2.0  2.0  1.0
            2.0  1.0  1.0
            2.0  1.5  1.0
            2.0  1.5  1.0
        ]

        @test isapprox(userbased_scores(m, [1, 2, 3, 4]), expected_scores, atol = 1e-2)
        @test isapprox(userbased_scores(m, [1, 2], [1, 3]), expected_scores[[1, 2], [1, 3]], atol = 1e-2)

        ranked_items, ranked_scores = userbased_rankings(m, 3, [1, 2, 3, 4])
        @test ranked_items == expected_ranked_items
        @test isapprox(ranked_scores, expected_ranked_scores, atol = 1e-2)

        ranked_items, ranked_scores = userbased_rankings(m, 2, [1, 2, 3, 4])
        @test ranked_items == expected_ranked_items[:, 1:2]
        @test isapprox(ranked_scores, expected_ranked_scores[:, 1:2], atol = 1e-2)

        ranked_items, ranked_scores = userbased_rankings(m, 2, [1, 2], [1, 3])
        @test ranked_items == [1 3; 3 1]
        @test isapprox(ranked_scores, [2.0  2.0; 2.0  1.0], atol = 1e-2)
    end
end

@testset "centering=false" begin
    c = Config(centering = false)
    m = memorize(R, c)

    @test m.nu == 4
    @test m.ni == 3
    @test m.bu == [0.0, 0.0, 0.0, 0.0]
    @test m.bi == [0.0, 0.0, 0.0]
    @test m.Dui == sparse([
        2.0  1.0  0.0
        1.0  0.0  2.0
        0.0  1.0  2.0
        0.0  0.0  2.0
    ])
    @test m.Diu == sparse([
        2.0  1.0  0.0  0.0
        1.0  0.0  1.0  0.0
        0.0  2.0  2.0  2.0
    ])

    r2, r5, r12 = sqrt(2), sqrt(5), sqrt(12)
    @test isapprox(m.Sii, sparse([
               1.0   2/(r5*r2)  2/(r5*r12)
         2/(r2*r5)         1.0  2/(r2*r12)
        2/(r12*r5)  2/(r12*r2)         1.0
    ]), atol = 1e-2)
    @test isapprox(m.Suu, sparse([
              1.0  2/(r5*r5)  1/(r5*r5)       0.0
        2/(r5*r5)        1.0  4/(r5*r5)  4/(r5*2)
        1/(r5*r5)  4/(r5*r5)        1.0  4/(r5*2)
              0.0   4/(r5*2)   4/(r5*2)       1.0
    ]), atol = 1e-2)
end


text:
# Autogenerated wrapper script for Xorg_libXt_jll for powerpc64le-linux-gnu
export libXt

using Xorg_libSM_jll
using Xorg_libX11_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"

# Relative path to `libXt`
const libXt_splitpath = ["lib", "libXt.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libXt_path = ""

# libXt-specific global declaration
# This will be filled out by __init__()
libXt_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libXt = "libXt.so.6"

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"Xorg_libXt")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # We first need to add to LIBPATH_list the libraries provided by Julia
    append!(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)])
    # From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
    # then append them to our own.
    foreach(p -> append!(PATH_list, p), (Xorg_libSM_jll.PATH_list, Xorg_libX11_jll.PATH_list,))
    foreach(p -> append!(LIBPATH_list, p), (Xorg_libSM_jll.LIBPATH_list, Xorg_libX11_jll.LIBPATH_list,))

    global libXt_path = normpath(joinpath(artifact_dir, libXt_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libXt_handle = dlopen(libXt_path)
    push!(LIBPATH_list, dirname(libXt_path))

    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(LIBPATH_list, ':')

    # Add each element of LIBPATH to our DL_LOAD_PATH (necessary on platforms
    # that don't honor our "already opened" trick)
    #for lp in LIBPATH_list
    #    push!(DL_LOAD_PATH, lp)
    #end
end  # __init__()


text:
using Documenter, ForwardDiff

makedocs(modules=[ForwardDiff],
         doctest = false,
         sitename = "ForwardDiff",
         pages = ["Introduction" => "index.md",
                  "User Documentation" => [
                    "Limitations of ForwardDiff" => "user/limitations.md",
                    "Differentiation API" => "user/api.md",
                    "Advanced Usage Guide" => "user/advanced.md",
                    "Upgrading from Older Versions" => "user/upgrade.md"],
                  "Developer Documentation" => [
                    "How ForwardDiff Works" => "dev/how_it_works.md",
                    "How to Contribute" => "dev/contributing.md"]])

deploydocs(
    repo = "github.com/JuliaDiff/ForwardDiff.jl.git"
)


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ f81b2c54-c541-41db-bdd2-fc9d7caa354c
using Pkg; Pkg.develop(path=".."); Pkg.activate("..")

# ╔═╡ 7c28cf4c-cd0c-11eb-0507-b10f2a883e8c
begin
    using RDatasets
    using LinearAlgebra
    using Plots
    using StatsPlots
    plotly(size=(680, 500))
end

# ╔═╡ 93477f7b-aef4-44e5-ba3b-8285a6261500
md"
# K-means

JingYu Ning
"

# ╔═╡ 64f0d05d-9c4e-4935-bf33-a5ee5835f99a
md"
## IRIS dataset
"

# ╔═╡ 7f59bdb6-c6de-4a25-9ff0-8a81457a4242
iris = dataset("datasets", "iris")

# ╔═╡ 63073482-7add-4234-ab34-6c8068b2f5a6
begin
    cols = names(iris)

    function gen_pots()
        plots = Matrix{Plots.Plot}(undef, 4, 4)
        for i in 1:4, j in 1:4
            if i == j
                plots[i, j] = @df(iris, histogram(cols(i)))
            elseif i < j
                plots[i, j] = @df(iris, histogram2d(
                    cols(i), cols(j),
                    nbins=20, color=:coolwarm, colorbar=false
                ))
            else
                plots[i, j] =  @df(iris, scatter(cols(i), cols(j)))
            end

            (i==1) && (plots[i, j] = plot!(ylabel=cols[j])) # add x label
            (j==4) && (plots[i, j] = plot!(xlabel=cols[i])) # add y label
        end

        return plots
    end
end

# ╔═╡ fde3d98a-be02-47fe-a8c1-293092298659
plot(gen_pots()..., layout=(4, 4), legend=false)

# ╔═╡ 30f8c663-e43b-44e2-9c91-b3e8dd760e5f
md"
## Initial center points
"

# ╔═╡ b1b771a2-3e94-4666-a668-01dee57f9832
md"
**Data size and numbers of classes**
"

# ╔═╡ 52e9ada2-cff4-4fe4-abb0-ec0ebc2d3fa1
begin
    n = size(iris, 1)
    n_classes = 3
end;

# ╔═╡ e6f38859-fd04-4dd6-84f3-991fdc72317d
md"
**Square of distance between two points**
"

# ╔═╡ 9ce95067-b54f-4a9d-8158-8e25b4a2480f
begin
    function d²(data::DataFrame, p₁_i::Integer, p₂_i::Integer)
        p₁₂ = Vector(data[p₁_i, 1:4]) - Vector(data[p₂_i, 1:4])

        return p₁₂' * p₁₂
    end

    function d²(data::DataFrame, p₁_i::Integer, p₂::Vector)
        p₁₂ = Vector(data[p₁_i, 1:4]) - p₂

        return p₁₂' * p₁₂
    end
end;

# ╔═╡ cbb9a475-36d9-4c9f-a12a-d2cc5c339be7
md"
**Find three center points**

```julia
# pick one point randomly
c₁_i = rand(1:n)

# calculate probabilities for candidate points
d²s = [d²(iris, i, c₁_i) for i in 1:n]
ps = d²s ./ sum(d²s)

# second point is the farest one from c₁
c₂_i = argmax(ps)

# calculate probabilities for candidate points
d²s = [min(d²(iris, i, c₁_i), d²(iris, i, c₂_i)) for i in 1:n]
ps = d²s ./ sum(d²s)

# third point is the farest one from c₁ and from c₂
c₃_i = argmax(ps)
```
"

# ╔═╡ ad7bfec2-3fe8-428c-ac92-ad56183cbcdd
begin
    # pick one point randomly
    c₁_i = rand(1:n)
    # calculate probabilities for candidate points
    d²s = [d²(iris, i, c₁_i) for i in 1:n]
    ps = d²s ./ sum(d²s)

    # second point is the farest one from c₁
    c₂_i = argmax(ps)
    # calculate probabilities for candidate points
    d²s = [min(d²(iris, i, c₁_i), d²(iris, i, c₂_i)) for i in 1:n]
    ps = d²s ./ sum(d²s)

    # third point is the farest one from c₁ and from c₂
    c₃_i = argmax(ps)

    cs_i = [c₁_i, c₂_i, c₃_i]
end

# ╔═╡ 2c6511ca-e070-47ac-8f2e-c1308f491f98
let
    plots = gen_pots()
    for i in 2:4, j in 1:(i-1)
        plots[i, j] = @df iris[cs_i, :] scatter!(plots[i, j], cols(i), cols(j))
    end
    plot(plots..., layout=(4, 4), legend=false)
end

# ╔═╡ bb476b47-ba34-41d7-a32f-fc204f4a328c
md"
## Training
"

# ╔═╡ 9cf597b7-9ea2-4c5d-b810-83b5d027d758
begin
    signal = "⚡"

    iris[!, :Predict] = zeros(n)
    cs = Matrix(iris[cs_i, 1:4])
    for _ in 1:10
        iris[!, :Predict] = map(
            i->argmin([d²(iris, i, cs[ci, :]) for ci in 1:n_classes]),
            1:n
        )

        for c in 1:n_classes
            class_data = iris[iris.Predict.==c, :]
            cs[c, :] = sum(Matrix(class_data[:, 1:4]), dims=1) ./ size(class_data, 1)
        end
    end
end

# ╔═╡ 7b5ea310-f5c2-4f9d-91b0-2ebae7b8639e
md"
**Prediction**
"

# ╔═╡ e0b73f20-18ba-4dff-9b75-0fe82d072167
let
    signal
    plots = gen_pots()
    for i in 2:4, j in 1:(i-1)
        plots[i, j] = @df iris scatter!(plots[i, j], cols(i), cols(j), group=:Predict)
        plots[i, j] = @df iris[cs_i, :] scatter!(plots[i, j], cols(i), cols(j))
    end
    plot(plots..., layout=(4, 4), legend=false)
end

# ╔═╡ 413aa844-5ce1-4349-a715-b62c8db9a1f8
md"
**Ground truth**
"

# ╔═╡ 1ae9328e-72c1-46e4-ba1e-742205c4bd88
let
    plots = gen_pots()
    for i in 2:4, j in 1:(i-1)
        plots[i, j] = @df iris scatter!(plots[i, j], cols(i), cols(j), group=:Species)
    end
    plot(plots..., layout=(4, 4), legend=false)
end

# ╔═╡ Cell order:
# ╟─93477f7b-aef4-44e5-ba3b-8285a6261500
# ╟─f81b2c54-c541-41db-bdd2-fc9d7caa354c
# ╠═7c28cf4c-cd0c-11eb-0507-b10f2a883e8c
# ╟─64f0d05d-9c4e-4935-bf33-a5ee5835f99a
# ╠═7f59bdb6-c6de-4a25-9ff0-8a81457a4242
# ╟─63073482-7add-4234-ab34-6c8068b2f5a6
# ╟─fde3d98a-be02-47fe-a8c1-293092298659
# ╟─30f8c663-e43b-44e2-9c91-b3e8dd760e5f
# ╟─b1b771a2-3e94-4666-a668-01dee57f9832
# ╠═52e9ada2-cff4-4fe4-abb0-ec0ebc2d3fa1
# ╟─e6f38859-fd04-4dd6-84f3-991fdc72317d
# ╠═9ce95067-b54f-4a9d-8158-8e25b4a2480f
# ╟─cbb9a475-36d9-4c9f-a12a-d2cc5c339be7
# ╠═ad7bfec2-3fe8-428c-ac92-ad56183cbcdd
# ╟─2c6511ca-e070-47ac-8f2e-c1308f491f98
# ╟─bb476b47-ba34-41d7-a32f-fc204f4a328c
# ╠═9cf597b7-9ea2-4c5d-b810-83b5d027d758
# ╟─7b5ea310-f5c2-4f9d-91b0-2ebae7b8639e
# ╟─e0b73f20-18ba-4dff-9b75-0fe82d072167
# ╟─413aa844-5ce1-4349-a715-b62c8db9a1f8
# ╟─1ae9328e-72c1-46e4-ba1e-742205c4bd88


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.DLR_Visualization (Modia3D/renderer/DLR_Visualization/_module.jl)
#

function Composition.initializeVisualization(renderer::Modia3D.AbstractDLR_VisualizationRenderer, velements::Vector{Composition.Object3D})::Nothing
    simVis = renderer.simVis
    @assert(length(velements) > 0)
    if simVis.isInitialized
        Composition.closeVisualization(renderer)
    end

    # Initialize SimVis
    SimVis_init(simVis)

    # Determine visualize functions and ids of visual elements
    for i in eachindex(velements)
        push!(simVis.ids, SimVis_getObjectID(simVis, 0))
    end
    simVis.velements = copy(velements)
    simVis.isInitialized = true
    return nothing
end

function Composition.visualize!(renderer::Modia3D.AbstractDLR_VisualizationRenderer, time::Float64)
    simVis = renderer.simVis
    @assert(length(simVis.velements) > 0)
    if simVis.isInitialized
        velements = simVis.velements
        ids       = simVis.ids
        for i in eachindex(velements)
            visualizeObject(velements[i], ids[i], simVis)
        end
        SimVis_setTime(simVis, time)
    else
        error("visualize! called without initializing visualization first.")
    end
end

function Composition.closeVisualization(renderer::Modia3D.AbstractDLR_VisualizationRenderer)
    simVis = renderer.simVis
    if simVis.isInitialized
        sleep(0.01)
        for id in simVis.ids
        SimVis_freeObjectID(simVis, id)
        end
        SimVis_shutdown(simVis)
    end
    empty!(simVis.ids)
    empty!(simVis.velements)
    simVis.isInitialized = false
end

function Composition.isVisible(feature::Shapes.Solid, renderer::Modia3D.AbstractRenderer)
    return !isnothing(feature.visualMaterial) && !isnothing(feature.shape)
end

function Composition.isVisible(feature::Shapes.Visual, renderer::Modia3D.AbstractRenderer)
    if isnothing(feature.shape)
        return false
    elseif typeof(feature.shape) == Shapes.TextShape
        return typeof(renderer) == ProfessionalEdition
    else
        return !isnothing(feature.visualMaterial)
    end
end


text:
export attach_gate!, contract, circuit2tn

_tmat(::Type{T}, blk) where T = transpose(Matrix(mat(T, blk)))

function attach_gate!(gb::GraphBuilder, tree::ChainBlock{N}) where N
    for b in subblocks(tree)
        attach_gate!(gb, b)
    end
end

function attach_gate!(gb::GraphBuilder{T}, tree::PutBlock{N}) where {T,N}
    tn_attach!(gb, tree.locs=>_tmat(T, content(tree)))
end

function attach_gate!(gb::GraphBuilder, tree::PutBlock{N,2,<:SWAPGate}) where N
    swaplines!(gb, tree.locs...)
end

function attach_gate!(gb::GraphBuilder{T}, tree::PutBlock{N,C,<:RotationGate{C,TR,<:KronBlock{C,C2,<:NTuple{C2,AbstractBlock{1}}}}}) where {N,C,C2,TR,T}
    tensors = Any[]
    if seperatecontrol()
        rt = content(tree)
        kr = rt.block
        for (i,loc_) in enumerate(kr.locs)
            loc = tree.locs[loc_]
            mi = IMatrix{2}()
            mk = _tmat(T, kr[loc_])
            if i==1
                mi = cos(rt.theta/2)*mi
                mk = -im*sin(rt.theta/2)*mk
            end
            if i==1 || i==length(kr.locs)
                tensor = zeros(T,2,2,2)
                tensor[:,:,1] = mi
                tensor[:,:,2] = mk
            else
                tensor = zeros(T,2,2,2,2)
                tensor[:,:,1,1] = mi
                tensor[:,:,2,2] = mk
            end
            push!(tensors, (loc...,)=>tensor)
        end
        tn_attach!(gb, tensors...)
    else
        tn_attach!(gb, tree.locs=>_tmat(T, content(tree)))
    end
end

function attach_gate!(gb::GraphBuilder{T}, tree::ControlBlock{N,GT,C,M}) where {N,GT,C,M,T}
    if seperatecontrol()
        tensors = Any[]
        for i in 1:C
            cbit = tree.ctrl_locs[i]
            cval = tree.ctrl_config[i]
            tensor = i == 1 ? control3(T;inverse=cval==0) : control4(T;inverse=cval==0)
            push!(tensors, (cbit,)=>tensor)
        end
        tensor = zeros(T, 2^M, 2^M, 2)
        tensor[1:2^(2*M)] .= vec(IMatrix{1<<M}())
        tensor[end-2^(2*M)+1:end] .= vec(_tmat(T, content(tree)))
        tn_attach!(gb, tensors..., tree.locs=>tensor)
    else
        locs = (tree.ctrl_locs..., tree.locs...)
        address = zeros(Int, N)
        address[[locs...]] .= 1:C+M
        c = map_address(tree, AddressInfo(C+M, address))
        tn_attach!(gb, locs=>Matrix(c))
    end
end

function circuit2tn(::Type{T}, c::AbstractBlock; initial_config=zeros(Int,nqubits(c)), final_config=zeros(Int, nqubits(c))) where T
    gb = GraphBuilder(T, initial_config)
    attach_gate!(gb, c)
    finish!(gb, final_config)
    return generate_eingraph(gb)
end

circuit2tn(c::AbstractBlock; initial_config=zeros(Int,nqubits(c)), final_config=zeros(Int, nqubits(c))) = circuit2tn(ComplexF64, c; initial_config=initial_config, final_config=final_config)

function opennetwork(c::AbstractBlock{N}) where N
    gb = GraphBuilder(zeros(Int, N))
    gb.tensors = gb.tensors[N+1:end]
    gb.labels = gb.tensors[N+1:end]
    attach_gate!(gb, c)
    return generate_eingraph(gb)
end


text:
using Random

struct RandomAgent{T}
    actions::T
end
start!(agent::RandomAgent, env_s_tp1, rng::AbstractRNG=Random.GLOBAL_RNG) = rand(rng, agent.actions)
step!(agent::RandomAgent, env_s_tp1, r, terminal, rng::AbstractRNG=Random.GLOBAL_RNG) = rand(rng, agent.actions)


text:
module RocketStartWithOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: start_with()" begin

    println("Testing: operator start_with()")

    run_proxyshowcheck("StartWith", start_with(0))

    run_testset([
        (
            source      = from(1:5) |> start_with(0),
            values      = @ts([ 0, 1, 2, 3, 4, 5, c ]),
            source_type = Int
        ),
        (
            source      = from(1:5) |> start_with(0.0),
            values      = @ts([ 0.0, 1, 2, 3, 4, 5, c ]),
            source_type = Union{Float64, Int}
        ),
        (
            source      = completed(Int) |> start_with("Hello"),
            values      = @ts([ "Hello", c]),
            source_type = Union{String, Int}
        ),
        (
            source = faulted(Int, "e") |> start_with(1),
            values = @ts([ 1, e("e") ]),
            source_type = Int
        ),
        (
            source = never(Int) |> start_with(0),
            values = @ts([ 0 ])
        )
    ])

end

end


text:
export iauAper13
"""
In the star-independent astrometry parameters, update only the
Earth rotation angle.  The caller provides UT1, (n.b. not UTC).

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   ut11    double      UT1 as a 2-part...
   ut12    double      ...Julian Date (Note 1)
   astrom  iauASTROM*  star-independent astrometry parameters:
   pmt    double       not used
   eb     double[3]    not used
   eh     double[3]    not used
   em     double       not used
   v      double[3]    not used
   bm1    double       not used
   bpn    double[3][3] not used
   along  double       longitude + s' (radians)
   xpl    double       not used
   ypl    double       not used
   sphi   double       not used
   cphi   double       not used
   diurab double       not used
   eral   double       not used
   refa   double       not used
   refb   double       not used

Returned:
   astrom  iauASTROM*  star-independent astrometry parameters:
   pmt    double       unchanged
   eb     double[3]    unchanged
   eh     double[3]    unchanged
   em     double       unchanged
   v      double[3]    unchanged
   bm1    double       unchanged
   bpn    double[3][3] unchanged
   along  double       unchanged
   xpl    double       unchanged
   ypl    double       unchanged
   sphi   double       unchanged
   cphi   double       unchanged
   diurab double       unchanged
   eral   double       "local" Earth rotation angle (radians)
   refa   double       unchanged
   refb   double       unchanged

Notes:

1) The UT1 date (n.b. not UTC) ut11+ut12 is a Julian Date,
   apportioned in any convenient way between the arguments ut11 and
   ut12.  For example, JD(UT1)=2450123.7 could be expressed in any
   of these ways, among others:

         ut11           ut12

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in cases
   where the loss of several decimal digits of resolution is
   acceptable.  The J2000 and MJD methods are good compromises
   between resolution and convenience.  The date & time method is
   best matched to the algorithm used:  maximum precision is
   delivered when the ut11 argument is for 0hrs UT1 on the day in
   question and the ut12 argument lies in the range 0 to 1, or vice
   versa.

2) If the caller wishes to provide the Earth rotation angle itself,
   the function iauAper can be used instead.  One use of this
   technique is to substitute Greenwich apparent sidereal time and
   thereby to support equinox based transformations directly.

3) This is one of several functions that inserts into the astrom
   structure star-independent parameters needed for the chain of
   astrometric transformations ICRS <-> GCRS <-> CIRS <-> observed.

   The various functions support different classes of observer and
   portions of the transformation chain:

         functions         observer        transformation

      iauApcg iauApcg13    geocentric      ICRS <-> GCRS
      iauApci iauApci13    terrestrial     ICRS <-> CIRS
      iauApco iauApco13    terrestrial     ICRS <-> observed
      iauApcs iauApcs13    space           ICRS <-> GCRS
      iauAper iauAper13    terrestrial     update Earth rotation
      iauApio iauApio13    terrestrial     CIRS <-> observed

   Those with names ending in "13" use contemporary SOFA models to
   compute the various ephemerides.  The others accept ephemerides
   supplied by the caller.

   The transformation from ICRS to GCRS covers space motion,
   parallax, light deflection, and aberration.  From GCRS to CIRS
   comprises frame bias and precession-nutation.  From CIRS to
   observed takes account of Earth rotation, polar motion, diurnal
   aberration and parallax (unless subsumed into the ICRS <-> GCRS
   transformation), and atmospheric refraction.

Called:
   iauAper      astrometry parameters: update ERA
   iauEra00     Earth rotation angle, IAU 2000

This revision:   2013 September 25

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
function iauAper13(ut11::Real, ut12::Real, astrom::iauASTROM)
   ref_astrom = Ref{iauASTROM}(astrom)
   ccall((:iauAper13, libsofa_c), Cvoid,
            (Cdouble, Cdouble, Ref{iauASTROM}),
            convert(Float64, ut11),
            convert(Float64, ut12),
            ref_astrom)

   return ref_astrom[]
end


text:
# instantiate project
import Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

const OUTPUT_FOLDER = "output"
const RESULTS_FILE = "runtimes.csv"

LONG = length(ARGS) == 1 && ARGS[1] ∈ ("LONG", "long")

function main()

    # --------------------------
    # Initialization
    # --------------------------
    if !isdir(OUTPUT_FOLDER)
        mkdir(OUTPUT_FOLDER)
        mkdir(joinpath(OUTPUT_FOLDER, "SDOF"))
        mkdir(joinpath(OUTPUT_FOLDER, "Clamped"))
        mkdir(joinpath(OUTPUT_FOLDER, "Heat1D"))
        mkdir(joinpath(OUTPUT_FOLDER, "Heat3D"))
    end
    global io = open(joinpath(OUTPUT_FOLDER, RESULTS_FILE), "w")
    global TARGET_FOLDER
    println("Running examples...")

    # --------------------------
    # Single-degree-of-freedom
    # --------------------------
    println("\n (1/5) Running single-degree-of-freedom (SDOF) model... \n")
    TARGET_FOLDER = joinpath(OUTPUT_FOLDER, "SDOF")
    include(joinpath("examples", "SDOF", "SDOF_Preliminaries.jl"))
    include(joinpath("examples", "SDOF", "SDOF_Time.jl"))
    include(joinpath("examples", "SDOF", "SDOF_PE_AD.jl"))

    # --------------------------
    # Clamped bar
    # --------------------------
    TARGET_FOLDER = joinpath(OUTPUT_FOLDER, "Clamped")
    println("\n (2/5) Running Clamped bar model\n")
    include(joinpath("examples", "Clamped", "Clamped.jl"))

    # --------------------------
    # Lambs Wave
    # --------------------------
    if LONG
        TARGET_FOLDER = joinpath(OUTPUT_FOLDER, "LambsWave")
        println("\n (3/5) Running Lambs wave propagation problem...\n")
        include(joinpath("examples", "LambsWave", "lambs_results.jl"))
    end

    # --------------------------
    # Heat 1D
    # --------------------------
    TARGET_FOLDER = joinpath(OUTPUT_FOLDER, "Heat1D")
    println("\n (4/5) Running one-dimensional heat transfer model...\n")
    include(joinpath("examples", "Heat1D", "Heat1D_TempProfile.jl"))
    include(joinpath("examples", "Heat1D", "Heat1D_TempGradient.jl"))

    # --------------------------
    # Heat 3D
    # --------------------------
    TARGET_FOLDER = joinpath(OUTPUT_FOLDER, "Heat3D")
    println("\n (5/5) Running three-dimensional concrete casting model...\n")
    include(joinpath("examples", "Heat3D", "Heat3D.jl"))

    print(io, "\n")
    println("Finished running benchmarks.")
    close(io)
    nothing
end

main()

max_stars_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_issues_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_forks_repo_name: JuliaReach/SetPropagation-FEM-Examples


text:
module NeuroCore

using ImageCore, ImageAxes, Unitful, ImageMetadata, Markdown
using StaticArrays, Rotations, CoordinateTransformations

export CoordinateList,
       NeuroMetadata,
       NeuroMetaArray,
       BIDSMetadata,
       # methods
       calmax,
       calmax!,
       calmin,
       calmin!,
       description,
       description!,
       freqdim,
       freqdim!,
       is_anatomical,
       is_functional,
       is_electrophysiology,
       neurohelp,
       phasedim,
       phasedim!,
       stream_offset,
       slicedim,
       slicedim!,
       slice_end,
       slice_start,
       spatial_offset,
       spatial_units,
       time_units

const OneF64Sec = 1.0u"s"

"F64Sec - A `Float64` type with seconds units."
const F64Sec = typeof(OneF64Sec)

const OneF64Tesla = 1.0u"T"

"F64Tesla - A `Float64` type with tesla units."
const F64Tesla = typeof(OneF64Tesla)

"F64Hertz - A `Float64` type with hertz units."
const OneF64Hz = 1.0u"Hz"
const F64Hz = typeof(OneF64Hz)

"IntHertz - A `Float64` type with degree units."
const OneIntDeg = 1.0u"°"
const IntDeg = typeof(OneIntDeg)

const OnekOhm = 1.0u"kΩ"
const F64kOhm = typeof(OnekOhm)

const CoordinateList = Dict{Symbol,NTuple{3,Float64}}

const NeuroAffine{R} = AffineMap{R,SArray{Tuple{3},Float64,1,3}}

include("enums.jl")
include("properties.jl")
include("coordinates.jl")
include("bids.jl")
include("getproperty.jl")
include("metadata.jl")
include("array.jl")
include("traits.jl")

"""
    neurohelp(func[;])

"""
neurohelp(func) = neurohelp(stdout, func)
neurohelp(io::IO, input::Symbol) = neurohelp(io, getproperty(NeuroCore, input))
function neurohelp(io::IO, input)
    buffer = IOBuffer()
    println(buffer, Base.Docs.doc(input))
    Markdown.parse(String(take!(buffer)))
end

end


text:
#5) Finds observations connected for every cluster

function check_clustering(clustering_var)
    NT = length(clustering_var)
    counter = ones(size(clustering_var,1));

    #Indexing obs number per worked/id
    gcs = Int.(@transform(groupby(DataFrame(counter = counter, clustering_var = clustering_var), :clustering_var), gcs = cumsum(:counter)).gcs)
    maxD = maximum(gcs)

    index = collect(1:NT)

    #This will be the output, and we will append observations to it in the loop
    list_final=DataFrame(row = Int64[],col = Int64[], id_cluster = Int64[])

    for t=1:maxD
        rowsel =  findall(x->x==true,gcs.==t)
        list_base = DataFrame( id_cluster= [clustering_var[x] for x in rowsel], row =
        [index[x] for x in rowsel]  )

            for tt=t:maxD
                colsel =  findall(x->x==true,gcs.==tt)
                list_sel =  DataFrame( id_cluster= [clustering_var[x] for x in colsel], col =
                [index[x] for x in colsel] )

                merge = outerjoin(list_base, list_sel, on = :id_cluster)
                merge = dropmissing(merge)
                merge = merge[:,[:row, :col, :id_cluster]]
                append!(list_final, merge)

            end
        end

        sort!(list_final, (:row))
        list_final = Matrix(list_final)

        return (list_final = list_final , nnz_2 = size(list_final,1) )

end

# Assumes that obs_leaveoutset is an array of zeros
function whole_prunning!(obs_leaveoutset,y,first_id,second_id;settings)
    obs,  y  , first_id , second_id  = find_connected_set(y,first_id,second_id,settings)
    obs,  y  , first_id , second_id  = prunning_connected_set(y,first_id,second_id,obs,settin)
    obs,  y  , first_id , second_id  = drop_single_obs(y,first_id,second_id,obs)
    for x in obs
        obs_leaveoutset[x] = 1
    end
    return nothing
end

#7) Computes exact Pii to do inference (may be deprecated in the future)
function do_Pii(X, clustering_var)
    n=size(X,1)

    #If no clustering, Lambda_P is just diagonal matrix.
    if clustering_var == nothing
        clustering_var = collect(1:n)
    end

    #Set matrices for parallel environment.
    xx=X'*X
    #P = aspreconditioner(ruge_stuben(xx))
    ldli, la = computeLDLinv(xx)
    buffs = zeros(size(la)[1],Threads.nthreads())
    compute_sol = []
    for i in 1:Threads.nthreads()
        P = approxcholOperator(ldli,buffs[:,i])
        push!(compute_sol,approxcholSolver(P,la))
    end

    #Return the structure of the indexes associated with the clustering variable
    elist, nnz_2 = check_clustering(clustering_var)
    M=size(elist,1)

    #Set elist
    elist_1 = elist[:,1]
    elist_2 = elist[:,2]
    Pii=zeros(M)

    Threads.@threads for i=1:M
        #zexact = zeros(size(X,2))
        col = elist_2[i]
        row = elist_1[i]
        #cg!(zexact, xx, X[col,:], Pl = P , log=true, maxiter=300)
        zexact = compute_sol[Threads.threadid()]([X[col,:]...],verbose=false)

        Pii[i]= (SparseMatrixCSC{Float64,Int64}(X[col,:])'*zexact)[1]
    end

    Lambda_P = sparse(elist[:,1],elist[:,2],Pii,n,n)

    return Lambda_P
end

#8) Performs Statistical Inference on Results
function lincom_KSS(y,X,Z,Transform,clustering_var,Lambda_P; joint_test =false, labels=nothing, restrict=nothing, nsim = 10000, settings = settings)
    #SET DIMENSIONS
    n=size(X,1)
    K=size(X,2)
    #Add Constant
    Z=hcat(ones(size(Z,1)), Z)

    # PART 1: ESTIMATE HIGH DIMENSIONAL MODEL
    xx=X'*X
    xy=X'*y
    compute_sol = approxcholSolver(xx;verbose = settings.print_level > 0)
    beta = compute_sol([xy...];verbose=false)
    eta=y-X*beta

    # PART 1B: VERIFY LEAVE OUT COMPUTATION
    if Lambda_P == nothing
        Lambda_P=do_Pii(X,clustering_var)
    end

    if Lambda_P != nothing && clustering_var !=nothing
        nnz_1=nnz(Lambda_P)
        nnz_2=check_clustering(clustering_var).nnz_2

        if nnz_1 == nnz_2
            println("The structure of the specified Lambda_P is consistent with the level of clustering required by the user.")
        elseif nnz_1 != nnz_2
            error("The user wants cluster robust inference but the Lambda_P provided by the user is not consistent with the level of clustering asked by the user. Try to omit input Lambda_P when running lincom_KSS")
        end
    end
    I_Lambda_P = I-Lambda_P
    eta_h = I_Lambda_P\eta

    #PART 2: SET UP MATRIX FOR SANDWICH FORMULA
    rows,columns, V = findnz(Lambda_P)

    aux= 0.5*(y[rows].*eta_h[columns] + y[columns].*eta_h[rows])
    sigma_i=sparse(rows,columns,aux,n,n)

    aux= 0.5*(eta[rows].*eta[columns] + eta[columns].*eta[rows])
    sigma_i_res=sparse(rows,columns,aux,n,n)

    r=size(Z,2);
    wy=Transform*beta
    zz=Z'*Z

    numerator=Z\wy
    chet=wy-Z*numerator
    aux= 0.5*(chet[rows].*chet[columns] + chet[columns].*chet[rows])
    sigma_i_chet=sparse(rows,columns,aux,n,n)

    #PART 3: COMPUTE
    denominator=zeros(r,1)
    denominator_RES=zeros(r,1)

    for q=1:r
        v=sparse([q],[1.0],[1.0],r,1)
        v=zz\[v...]
        v=Z*v
        v=Transform'*v

        right = compute_sol(v;verbose=false)

        left=right'

        denominator[q]=left*(X'*sigma_i*X)*right
        denominator_RES[q]=left*(X'*sigma_i_res*X)*right
    end

    test_statistic=numerator./(sqrt.(denominator))
    #zz_inv=zz^(-1)
    SE_linear_combination_NAI=zz\(Z'*sigma_i_chet*Z)/zz

    #PART 4: REPORT
    println("Inference on Linear Combinations:")
    if labels == nothing
        for q=2:r
            if q <= r
                println("Linear Combination - Column Number ", q-1," of Z: ", numerator[q] )
                println("Standard Error of the Linear Combination - Column Number ", q-1," of Z: ", sqrt(denominator[q]) )
                println("T-statistic - Column Number ", q-1, " of Z: ", test_statistic[q])
            end
        end
    else
        for q=2:r
            tell_me = labels[q-1]
            println("Linear Combination associated with ", tell_me,": ", numerator[q] )
            println("Standard Error  associated with ", tell_me,": ", sqrt(denominator[q]) )
            println("T-statistic  associated with ", tell_me,": ", test_statistic[q])
        end
    end

    # PART 5: Joint-test. Quadratic form beta'*A*beta
    if joint_test == true

        if restrict  == nothing
            restrict=sparse(collect(1:r-1),collect(2:r),1.0,r-1,r)
        end

        v=restrict*(zz\(Z'*Transform))
        v=v'
        #v=sparse(v) #ldiv doesn't work for sparse RHS
        r=size(v,2)

        #Auxiliary
        aux=xx\v[:,:]
        opt_weight=v'*aux
        opt_weight=opt_weight^(-1)
        opt_weight=(1/r)*(opt_weight+opt_weight')/2

        #Eigenvalues, eigenvectors, and relevant components
        lambda , Qtilde = eigs( v*opt_weight*v', xx; nev=r,ritzvec=true)
        lambda = Real.(lambda)
        Qtilde = Real.(Qtilde)
        #lambdaS, QtildeS = eigs(v*opt_weight*v', xx; nev=1,which=:SM,ritzvec=true)
        #lambdaS = [lambdaL; lambdaS]
        #Qtilde = hcat(QtildeL,QtildeS)

        W=X*Qtilde
        V_b=W'*sigma_i*W
        V_b = (1/2)*(V_b + V_b')

        #Now focus on obtaining matrix Lambda_B with the A test associated with a joint hypothesis testing.
        Bii=opt_weight^(0.5)*aux';
        Bii=Bii*X'
        Bii=Bii'
        Bii = 0.5*(Bii[rows,:].*Bii[columns,:] + Bii[columns,:].*Bii[rows,:])
        Bii = sum(Bii,dims=2)[:]
        Lambda_B=sparse(rows,columns,Bii,n,n)

        #Leave Out Joint-Statistic
        stat=(v'*beta)'*opt_weight*(v'*beta)-y'*Lambda_B*eta_h

        #Now simulate critical values under the null.
        mu=zeros(r)
        sigma = V_b
        b_sim = MvNormal(mu,sigma)
        b_sim = rand(b_sim, nsim)

        #theta_star_sim=sum(lambda'.*(b_sim.^2 - diag(V_b)'),2)
        theta_star_sim = sum(lambda.*b_sim.^2 .- lambda.* diag(V_b),dims=1)
        pvalue=mean(theta_star_sim.>stat)

        #Report
        println("Joint-Test Statistic: ", stat)
        println("p-value: ", pvalue)

    end

    test_statistic=test_statistic[2:end]
    linear_combination=numerator[2:end]
    SE_linear_combination_KSS=sqrt.(denominator[2:end])
    SE_linear_combination_RES=sqrt.(denominator_RES[2:end])
    SE_linear_combination_NAI=diag(SE_linear_combination_NAI)
    SE_linear_combination_NAI=sqrt.(SE_linear_combination_NAI[2:end])

    return nothing
end


text:
module BayesTests

using Distributions, Gadfly,
      KernelDensity, Cubature, Mamba

export BinomTest, CorrelationTest, TTest,
       plot, update, posterior, ppv, hpd

include("common.jl")
include("TTest.jl")
include("BinomTest.jl")
include("CorrelationTest.jl")

end # module


text:
#==============================================================================#
# sign.jl
#
# AWS Request Signing.
#
# Copyright OC Technology Pty Ltd 2014 - All rights reserved
#==============================================================================#

using MbedTLS

function sign!(r::AWSRequest, t = now(Dates.UTC))

    if r[:service] in ["sdb", "importexport"]
        sign_aws2!(r, t)
    else
        sign_aws4!(r, t)
    end
end

# Create AWS Signature Version 2 Authentication query parameters.
# http://docs.aws.amazon.com/general/latest/gr/signature-version-2.html

function sign_aws2!(r::AWSRequest, t)

    query = Dict{AbstractString,AbstractString}()
    for elem in split(r[:content], '&', keepempty=false)
        (n, v) = split(elem, "=")
        query[n] = HTTP.unescapeuri(v)
    end

    r[:headers]["Content-Type"] =
        "application/x-www-form-urlencoded; charset=utf-8"

    creds = check_credentials(r[:creds])
    query["AWSAccessKeyId"] = creds.access_key_id
    query["Expires"] = Dates.format(t + Dates.Second(120),
                                   dateformat"yyyy-mm-dd\THH:MM:SS\Z")
    query["SignatureVersion"] = "2"
    query["SignatureMethod"] = "HmacSHA256"
    if creds.token != ""
        query["SecurityToken"] = creds.token
    end

    query = Pair[k => query[k] for k in sort(collect(keys(query)))]

    u = HTTP.URI(r[:url])
    to_sign = "POST\n$(u.host)\n$(u.path)\n$(HTTP.escapeuri(query))"

    secret = creds.secret_key
    push!(query, "Signature" =>
                  digest(MD_SHA256, to_sign, secret) |> base64encode |> strip)

    r[:content] = HTTP.escapeuri(query)
end

# Create AWS Signature Version 4 Authentication Headers.
# http://docs.aws.amazon.com/general/latest/gr/signature-version-4.html

function sign_aws4!(r::AWSRequest, t)

    # ISO8601 date/time strings for time of request...
    date = Dates.format(t,"yyyymmdd")
    datetime = Dates.format(t, dateformat"yyyymmdd\THHMMSS\Z")

    # Authentication scope...
    scope = [date, r[:region], r[:service], "aws4_request"]

    creds = check_credentials(r[:creds])
    # Signing key generated from today's scope string...
    signing_key = string("AWS4", creds.secret_key)
    for element in scope
        signing_key = digest(MD_SHA256, element, signing_key)
    end

    # Authentication scope string...
    scope = join(scope, "/")

    # SHA256 hash of content...
    content_hash = bytes2hex(digest(MD_SHA256, r[:content]))

    # HTTP headers...
    delete!(r[:headers], "Authorization")
    merge!(r[:headers], Dict(
        "x-amz-content-sha256" => content_hash,
        "x-amz-date"           => datetime,
        "Content-MD5"          => base64encode(digest(MD_MD5, r[:content]))
    ))
    if creds.token != ""
        r[:headers]["x-amz-security-token"] = creds.token
    end

    # Sort and lowercase() Headers to produce canonical form...
    canonical_headers = ["$(lowercase(k)):$(strip(v))" for (k,v) in r[:headers]]
    signed_headers = join(sort([lowercase(k) for k in keys(r[:headers])]), ";")

    # Sort Query String...
    uri = HTTP.URI(r[:url])
    query = HTTP.URIs.queryparams(uri.query)
    query = Pair[k => query[k] for k in sort(collect(keys(query)))]

    # Create hash of canonical request...
    canonical_form = string(r[:verb], "\n",
                            r[:service] == "s3" ? uri.path
                                                : HTTP.escapepath(uri.path), "\n",
                            HTTP.escapeuri(query), "\n",
                            join(sort(canonical_headers), "\n"), "\n\n",
                            signed_headers, "\n",
                            content_hash)
    if debug_level > 2
        println("canonical_form:")
        println(canonical_form)
        println("")
    end

    canonical_hash = bytes2hex(digest(MD_SHA256, canonical_form))

    # Create and sign "String to Sign"...
    string_to_sign = "AWS4-HMAC-SHA256\n$datetime\n$scope\n$canonical_hash"
    signature = bytes2hex(digest(MD_SHA256, string_to_sign, signing_key))

    if debug_level > 2
        println("string_to_sign:")
        println(string_to_sign)
        println("")
        println("signature:")
        println(signature)
        println("")
    end

    # Append Authorization header...
    r[:headers]["Authorization"] = string(
        "AWS4-HMAC-SHA256 ",
        "Credential=$(creds.access_key_id)/$scope, ",
        "SignedHeaders=$signed_headers, ",
        "Signature=$signature"
    )
end

#==============================================================================#
# End of file.
#==============================================================================#


text:
""" 2-category of categories, functors, and natural transformations.

Categories in mathematics appear in the large, often as categories of sets with
extra structure, and in the small, as algebraic structures that generalize
groups, monoids, preorders, and graphs. This division manifests in Catlab as
well. Large categories (in spirit, if not in the [technical
sense](https://ncatlab.org/nlab/show/large+category)) occur throughout Catlab as
`@instance`s of the theory of categories. For computational reasons, small
categories are usually presented by generators and relations.

This module provides a minimal interface to accomodate both situations. Category
instances are supported through the wrapper type [`TypeCat`](@ref). Finitely
presented categories are provided by another module, [`FinCats`](@ref).
"""
module Categories
export Cat, TypeCat, Ob, Functor, Transformation,
  dom, codom, compose, id, ob, hom, is_hom_equal,
  ob_map, hom_map, dom_ob, codom_ob, component

using AutoHashEquals

using ...GAT, ..Sets
import ...Theories: Category2, Ob, ob, hom, dom, codom, compose, ⋅, ∘, id,
  composeH, *

# Categories
############

""" Abstract base type for a category.

The objects and morphisms in the category have Julia types `Ob` and `Hom`,
respectively. Note that these types do *not* necessarily form an `@instance` of
the theory of categories, as they may not meaningfully form a category outside
the context of this object. For example, a finite category regarded as a
reflexive graph with a composition operation might have type `Cat{Int,Int}`,
where the objects and morphisms are numerical identifiers for vertices and edges
in the graph.

The basic operations available in any category are: [`dom`](@ref),
[`codom`](@ref), [`id`](@ref), [`compose`](@ref).
"""
abstract type Cat{Ob,Hom} end

""" Coerce or look up object in category.
"""
ob(::Cat, x) = x

""" Coerce or look up morphism in category.
"""
hom(::Cat, f) = f

""" Domain of morphism in category.
"""
dom(C::Cat, f) = dom(f)

""" Codomain of morphism in category.
"""
codom(C::Cat, f) = codom(f)

""" Identity morphism on object in category.
"""
id(C::Cat, x) = id(x)

""" Compose morphisms in a category.
"""
compose(C::Cat, fs...) = compose(fs...)

""" Are two morphisms in a category equal?

By default, just checks for equality of Julia objects using ``==``. In some
categories, checking equality of morphisms may involve nontrivial reasoning.
"""
is_hom_equal(C::Cat, f, g) = is_hom_equal(f, g)
is_hom_equal(f, g) = f == g

# Instances
#----------

""" Pair of Julia types regarded as a category.

The Julia types should form an `@instance` of the theory of categories
(`Theories.Category`).
"""
struct TypeCat{Ob,Hom} <: Cat{Ob,Hom} end

TypeCat(Ob::Type, Hom::Type) = TypeCat{Ob,Hom}()

Ob(::TypeCat{T}) where T = TypeSet{T}()

# FIXME: This isn't practical because types are often too tight.
#ob(::TypeCat{Ob,Hom}, x) where {Ob,Hom} = convert(Ob, x)
#hom(::TypeCat{Ob,Hom}, f) where {Ob,Hom} = convert(Hom, f)

# Functors
##########

""" Abstract base type for a functor between categories.

A functor has a domain and a codomain ([`dom`](@ref) and [`codom`](@ref)), which
are categories, and object and morphism maps, which can be evaluated using
[`ob_map`](@ref) and [`hom_map`](@ref). The functor object can also be called
directly when the objects and morphisms have distinct Julia types. This is
sometimes but not always the case (see [`Cat`](@ref)), so when writing generic
code one should prefer the `ob_map` and `hom_map` functions.
"""
abstract type Functor{Dom<:Cat,Codom<:Cat} end

""" Evaluate functor on object.
"""
function ob_map end

""" Evaluate functor on morphism.
"""
function hom_map end

""" Forgetful functor Ob: Cat → Set.

Sends a category to its set of objects and a functor to its object map.
"""
function Ob end

@auto_hash_equals struct IdentityFunctor{Dom<:Cat} <: Functor{Dom,Dom}
  dom::Dom
end

codom(F::IdentityFunctor) = F.dom

ob_map(F::IdentityFunctor, x) = ob(F.dom, x)
hom_map(F::IdentityFunctor, f) = hom(F.dom, f)

# Instances
#----------

(F::Functor{TypeCat{Ob,Hom}})(x::Ob) where {Ob,Hom} = ob_map(F, x)
(F::Functor{TypeCat{Ob,Hom}})(f::Hom) where {Ob,Hom} = hom_map(F, f)

# Natural transformations
#########################

""" Abstract base type for a natural transformation between functors.

A natural transformation ``α: F ⇒ G`` has a domain ``F`` and codomain ``G``
([`dom`](@ref) and [`codom`](@ref)), which are functors ``F,G: C → D`` having
the same domain ``C`` and codomain ``D``. The transformation consists of a
component ``αₓ: Fx → Gx`` in ``D`` for each object ``x ∈ C``, accessible using
[`component`](@ref) or indexing notation (`Base.getindex`).
"""
abstract type Transformation{C<:Cat,D<:Cat,Dom<:Functor{C,D},Codom<:Functor{C,D}} end

""" Component of natural transformation.
"""
function component end

@inline Base.getindex(α::Transformation, c) = component(α, c)

""" Domain object of natural transformation.

Given ``α: F ⇒ G: C → D``, this function returns ``C``.
"""
dom_ob(α::Transformation) = dom(dom(α)) # == dom(codom(α))

""" Codomain object of natural transformation.

Given ``α: F ⇒ G: C → D``, this function returns ``D``.
"""
codom_ob(α::Transformation) = codom(dom(α)) # == codom(codom(α))

@auto_hash_equals struct IdentityTransformation{C<:Cat,D<:Cat,Dom<:Functor{C,D}} <:
     Transformation{C,D,Dom,Dom}
  dom::Dom
end

codom(α::IdentityTransformation) = α.dom

const IdIdTransformation{C<:Cat} = IdentityTransformation{C,C,IdentityFunctor{C}}

# 2-category of categories
##########################

@instance Category2{Cat,Functor,Transformation} begin
  dom(F::Functor) = F.dom
  codom(F::Functor) = F.codom
  id(C::Cat) = IdentityFunctor(C)

  function compose(F::Functor, G::Functor)
    codom(F) == dom(G) || error("Domain mismatch in composition $F ⋅ $G")
    compose_id(F, G)
  end

  dom(α::Transformation) = α.dom
  codom(α::Transformation) = α.codom
  id(F::Functor) = IdentityTransformation(F)

  function compose(α::Transformation, β::Transformation)
    codom(α) == dom(β) || error("Domain mismatch in composition $α ⋅ $β")
    compose_id(α, β)
  end
  function composeH(α::Transformation, β::Transformation)
    codom_ob(α) == dom_ob(β) || error("Domain mismatch in composition $α * $β")
    composeH_id(α, β)
  end
  function composeH(α::Transformation, H::Functor)
    codom_ob(α) == dom(H) || error("Domain mismatch in whiskering $α * $H")
    composeH_id(α, H)
  end
  function composeH(F::Functor, β::Transformation)
    codom(F) == dom_ob(β) || error("Domain mismatch in whiskering $F * $β")
    composeH_id(F, β)
  end
end

# XXX: Is this normalization of identities using multiple dispatch a good idea?
# In contrast to `Sets`, it requires a lot of boilerplate here.

@inline compose_id(F::Functor, G::Functor) = do_compose(F, G)
@inline compose_id(F::Functor, ::IdentityFunctor) = F
@inline compose_id(::IdentityFunctor, G::Functor) = G
@inline compose_id(F::IdentityFunctor, ::IdentityFunctor) = F

@inline compose_id(α::Transformation, β::Transformation) = do_compose(α, β)
@inline compose_id(α::Transformation, ::IdentityTransformation) = α
@inline compose_id(::IdentityTransformation, β::Transformation) = β
@inline compose_id(α::IdentityTransformation, ::IdentityTransformation) = α

@inline composeH_id(α::Transformation, β::Transformation) = do_composeH(α, β)
@inline composeH_id(α::Transformation, ::IdIdTransformation) = α
@inline composeH_id(::IdIdTransformation, β::Transformation) = β
@inline composeH_id(α::IdIdTransformation, ::IdIdTransformation) = α

@inline composeH_id(α::Transformation, H::Functor) = do_composeH(α, H)
@inline composeH_id(α::Transformation, ::IdentityFunctor) = α
@inline composeH_id(α::IdentityTransformation, H::Functor) =
  id(compose_id(dom(α), H))
@inline composeH_id(α::IdentityTransformation, ::IdentityFunctor) = α

@inline composeH_id(F::Functor, β::Transformation) = do_composeH(F, β)
@inline composeH_id(::IdentityFunctor, β::Transformation) = β
@inline composeH_id(F::Functor, β::IdentityTransformation) =
  id(compose_id(F, dom(β)))
@inline composeH_id(::IdentityFunctor, β::IdentityTransformation) = β

function do_compose end

@inline function do_composeH(α::Transformation, β::Transformation)
  do_composeH(α, β, Val{:covariant})
end
function do_composeH(α::Transformation, β::Transformation,
                     ::Type{Val{:covariant}})
  G, H = codom(α), dom(β)
  compose_id(composeH_id(α, H), composeH_id(G, β))
end
function do_composeH(α::Transformation, β::Transformation,
                     ::Type{Val{:contravariant}})
  F, K = dom(α), codom(β)
  compose_id(composeH_id(F, β), composeH_id(α, K))
end

end


text:
module NSGAIII

using ProgressMeter
using Requires

include("indivs.jl")
include("functions.jl")
include("crossover.jl")
include("mutation.jl")
include("realcoding.jl")

export nsga, RealCoding, decode, encode, DennisDas

function nsga(popSize::Integer, nbGen::Integer, init::Function, z::Function, H::Int=15; args...)
    nbobj = length(z(init()))
    nsga(popSize, nbGen, init, z, DennisDas(nbobj, H) ; args...)
end

function nsga(popSize::Integer, nbGen::Integer, init::Function, z::Function, references::Vector{Vector{Float64}} ;
    fCV=(x)->0. , pmut=0.05, fmut=default_mutation!, fcross = default_crossover, seed=typeof(init())[], fplot = (x)->nothing)

    popSize = max(popSize, length(references))
    if popSize % 4 != 0
        popSize += 4 - (popSize % 4)
    end

    X = typeof(init())
    P = [indiv(init(), z, fCV) for _=1:popSize-length(seed)]
    append!(P, indiv.(convert.(X, seed),z, fCV))
    fast_non_dominated_sort!(P)
    associate_references!(P, references)
    Q = similar(P)

    @showprogress 0.1 for gen = 1:nbGen

        ind_Q = 1
        for pass = 1:2
            shuffle!(P)
            for i = 1:4:popSize
                pa = niching_based_selection(P[i], P[i+1], references)
                pb = niching_based_selection(P[i+2], P[i+3], references)
                ca,cb = crossover(pa, pb, fcross)

                rand() < pmut && mutate!(ca, fmut)
                rand() < pmut && mutate!(cb, fmut)

                eval!(ca, z, fCV)
                eval!(cb, z, fCV)

                Q[ind_Q] = ca
                Q[ind_Q+1] = cb

                ind_Q += 2
            end
        end

        @assert length(Q) == length(P)

        F = fast_non_dominated_sort!(vcat(P, Q))
        i = 1
        empty!(P)
        while length(P) + length(F[i]) <= popSize
            append!(P, F[i])
            i += 1
        end

        if length(P) != popSize
            T = vcat(P, F[i])
            normalize_pop!(T)
            associate_references!(T, references)
            n_count = niche_count(P, references)
            while length(P) < popSize
                niching!(P, F[i], n_count)
            end
        else
            associate_references!(P, references)
        end

        fplot(P)
    end
    P
end

@require vOptGeneric begin
    function nsga(popSize, nbGen, H, m ; kwargs...)

        vd = @eval Main getvOptData(m)
        @assert all(isfinite, m.colLower) "All variables must be bounded"
        @assert all(isfinite, m.colUpper) "All variables must be bounded"

        @assert !(:Int in m.colCat) "Only continuous and binary variables are supported"
        ϵ = map(x -> x==:Cont ? 5 : 0, m.colCat)
        d = RealCoding(ϵ, m.colLower, m.colUpper)

        init = () -> bitrand(d.nbbitstotal)

        function evaluate(obj, x)
            dot(obj.aff.coeffs, x[map(v-> getfield(v, :col), obj.aff.vars)]) + obj.aff.constant
        end

        function z(bits)
            x = decode(bits, d)
            ((evaluate(obj, x) for obj in vd.objs)...)
        end

        function CV(bits)
            x = decode(bits, d)
            res = 0.
            for CSTR in m.linconstr

                if CSTR.lb != -Inf && CSTR.lb != typemin(Float64)
                    if CSTR.lb == 0
                        g = dot(CSTR.terms.coeffs, x[map(v-> getfield(v, :col), CSTR.terms.vars)])
                        res += max(0, -g)
                    elseif CSTR.lb > 0
                        g = dot(CSTR.terms.coeffs, x[map(v-> getfield(v, :col), CSTR.terms.vars)]) / CSTR.lb - 1
                        res += max(0, -g)
                    else
                        g = dot(CSTR.terms.coeffs, x[map(v-> getfield(v, :col), CSTR.terms.vars)]) / CSTR.lb - 1
                        res += max(0, g)
                    end

                end

                if CSTR.ub != Inf && CSTR.lb != typemax(Float64)
                     if CSTR.ub == 0
                        g = dot(CSTR.terms.coeffs, x[map(v-> getfield(v, :col), CSTR.terms.vars)])
                        res += max(0, g)
                    elseif CSTR.ub > 0
                        g = dot(CSTR.terms.coeffs, x[map(v-> getfield(v, :col), CSTR.terms.vars)]) / CSTR.ub - 1
                        res += max(0, g)
                    else
                        g = dot(CSTR.terms.coeffs, x[map(v-> getfield(v, :col), CSTR.terms.vars)]) / CSTR.ub - 1
                        res += max(0, -g)
                    end
                end
            end
            res
        end

        for i = 1:length(vd.objs)
            if vd.objSenses[i] == :Max
                vd.objs[i] = vd.objs[i] * -1
            end
        end

        res = nsga(popSize, nbGen, init, z, H ; fCV = CV, kwargs...)

        for i = 1:length(vd.objs)
            if vd.objSenses[i] == :Max
                vd.objs[i] = vd.objs[i] * -1
            end
        end

        signs = Tuple(s == :Min ? 1 : -1 for s in vd.objSenses)

        @show signs

        for indiv in res
            indiv.y = indiv.y .* signs
        end

        [(decode(ind.x, d), ind.y, ind.CV) for ind in res]

    end
end

end # module


text:
module Factors

import Base.Iterators.product
const FactorTable = Dict{NamedTuple,Float64}

struct Variable
    name::Symbol
    m::Int # number of possible values
end

struct Factor
    vars::Vector{Variable}
    table::FactorTable
end

variablenames(ϕ::Factor) = [var.name for var in ϕ.vars]

function namedtuple(n::Array{Symbol,1}, v::Tuple)
    return NamedTuple{Tuple(n)}(v)
end

function assignments(vars::AbstractVector{Variable})
    n = [var.name for var in vars]
    #n = Tuple([var.name for var in vars])
    #return [NamedTuple{n}(v) for v in product((1:v.m for v in vars)...)]
    return [namedtuple(n, v) for v in product((1:v.m for v in vars)...)]
end

X = Variable(:x, 2)
Y = Variable(:y, 2)
Z = Variable(:z, 2)

ϕ = Factor([X, Y, Z], FactorTable(
    (x=1, y=1, z=1) => 0.08,
    (x=1, y=1, z=2) => 0.31,
    (x=1, y=2, z=1) => 0.09,
    (x=1, y=2, z=2) => 0.37,
    (x=2, y=1, z=1) => 0.01,
    (x=2, y=1, z=2) => 0.05,
    (x=2, y=2, z=1) => 0.02,
    (x=2, y=2, z=2) => 0.07,
))

display(ϕ.table)

display(ϕ.vars)

display(variablenames(ϕ))

"""
    normalize!(ϕ::Factor)
Normalize a factor ϕ, which divides all the entries in the factor by the
same scalar so they sum to 1.
"""
function normalize!(ϕ::Factor)

    z = sum(p for (a,p) in ϕ.table)
    for (a,p) in ϕ.table
        ϕ.table[a] = p/z
    end
    return ϕ
end

k = assignments(ϕ.vars)
@show k
@show k[2,1,1]
@show k[2,1,1].x
@show k[2,1,1].y
@show k[2,1,1].z

println(k)
println(typeof(k))
println(size(k))

sum([2,3,2])

end


text:
module Capstan
    println("There's nothing here!")
end # module