text:
"""
A tool for composing and sequencing music using Julia.
"""
module SirenSeq

include("midi/Midi.jl")
include("scales/Scales.jl")

using SirenSeq.Midi
using SirenSeq.Scales

include("constants/GenFluidSF144.jl")
include("constants/GenMuseScoreSF144.jl")
include("constants/GmDrums.jl")

## use Julia commandline ? option for exported function descritions
export

## from basics.jl
IR, Exp, Atom, Duratom,
atomOrExp, offsets,
C, S, R,
ftoi14, ftoi7, itof14, itof7,

## from shortcuts.jl
D, F, A, T, V, Cha, Sca,
setDefaultScale, setDefaultOctave,
N, B,
Bank, Prog, Vol,
Wh, Aft, Cn7, Cn14,
Sus0, Sus1,
SoundsOff, CntrlOff, NotesOff,
Mark, Sig, Tpo,
Inst, Da,
Con, Kl, Spl,

## from exportmidi.jl
setDefaultBpm, makeMidi,

## submodules
Play, Render, Midi, Scales,
GenFluidSF144, GenMuseScoreSF144, GmDrums

include("basics.jl")
include("atoms.jl")
include("shortcuts.jl")
include("display.jl")
include("exportmidi.jl")
include("playmidi.jl")
include("midirender.jl")
include("csound.jl")

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## client.jl - frontend handling command line options, environment setup,
##             and REPL

have_color = false
default_color_warn = :yellow
default_color_error = :light_red
default_color_info = :cyan
default_color_debug = :blue
default_color_input = :normal
default_color_answer = :normal
color_normal = text_colors[:normal]

function repl_color(key, default)
    env_str = get(ENV, key, "")
    c = tryparse(Int, env_str)
    c_conv = something(c, Symbol(env_str))
    haskey(text_colors, c_conv) ? c_conv : default
end

error_color() = repl_color("JULIA_ERROR_COLOR", default_color_error)
warn_color()  = repl_color("JULIA_WARN_COLOR" , default_color_warn)
info_color()  = repl_color("JULIA_INFO_COLOR" , default_color_info)
debug_color()  = repl_color("JULIA_DEBUG_COLOR" , default_color_debug)

input_color()  = text_colors[repl_color("JULIA_INPUT_COLOR", default_color_input)]
answer_color() = text_colors[repl_color("JULIA_ANSWER_COLOR", default_color_answer)]

stackframe_lineinfo_color() = repl_color("JULIA_STACKFRAME_LINEINFO_COLOR", :bold)
stackframe_function_color() = repl_color("JULIA_STACKFRAME_FUNCTION_COLOR", :bold)

function repl_cmd(cmd, out)
    shell = shell_split(get(ENV, "JULIA_SHELL", get(ENV, "SHELL", "/bin/sh")))
    shell_name = Base.basename(shell[1])

    # Immediately expand all arguments, so that typing e.g. ~/bin/foo works.
    cmd.exec .= expanduser.(cmd.exec)

    if isempty(cmd.exec)
        throw(ArgumentError("no cmd to execute"))
    elseif cmd.exec[1] == "cd"
        new_oldpwd = pwd()
        if length(cmd.exec) > 2
            throw(ArgumentError("cd method only takes one argument"))
        elseif length(cmd.exec) == 2
            dir = cmd.exec[2]
            if dir == "-"
                if !haskey(ENV, "OLDPWD")
                    error("cd: OLDPWD not set")
                end
                dir = ENV["OLDPWD"]
            end
            cd(dir)
        else
            cd()
        end
        ENV["OLDPWD"] = new_oldpwd
        println(out, pwd())
    else
        @static if !Sys.iswindows()
            if shell_name == "fish"
                shell_escape_cmd = "begin; $(shell_escape_posixly(cmd)); and true; end"
            else
                shell_escape_cmd = "($(shell_escape_posixly(cmd))) && true"
            end
            cmd = `$shell -c $shell_escape_cmd`
        end
        run(ignorestatus(cmd))
    end
    nothing
end

# deprecated function--preserved for DocTests.jl
function ip_matches_func(ip, func::Symbol)
    for fr in StackTraces.lookup(ip)
        if fr === StackTraces.UNKNOWN || fr.from_c
            return false
        end
        fr.func === func && return true
    end
    return false
end

function scrub_repl_backtrace(bt)
    if bt !== nothing && !(bt isa Vector{Any}) # ignore our sentinel value types
        bt = stacktrace(bt)
        # remove REPL-related frames from interactive printing
        eval_ind = findlast(frame -> !frame.from_c && frame.func === :eval, bt)
        eval_ind === nothing || deleteat!(bt, eval_ind:length(bt))
    end
    return bt
end

function display_error(io::IO, er, bt)
    printstyled(io, "ERROR: "; bold=true, color=Base.error_color())
    bt = scrub_repl_backtrace(bt)
    showerror(IOContext(io, :limit => true), er, bt, backtrace = bt!==nothing)
    println(io)
end
function display_error(io::IO, stack::Vector)
    printstyled(io, "ERROR: "; bold=true, color=Base.error_color())
    bt = Any[ (x[1], scrub_repl_backtrace(x[2])) for x in stack ]
    show_exception_stack(IOContext(io, :limit => true), bt)
end
display_error(stack::Vector) = display_error(stderr, stack)
display_error(er, bt=nothing) = display_error(stderr, er, bt)

function eval_user_input(errio, @nospecialize(ast), show_value::Bool)
    errcount = 0
    lasterr = nothing
    while true
        try
            if have_color
                print(color_normal)
            end
            if lasterr !== nothing
                invokelatest(display_error, errio, lasterr)
                errcount = 0
                lasterr = nothing
            else
                ast = Meta.lower(Main, ast)
                value = Core.eval(Main, ast)
                ccall(:jl_set_global, Cvoid, (Any, Any, Any), Main, :ans, value)
                if !(value === nothing) && show_value
                    if have_color
                        print(answer_color())
                    end
                    try
                        invokelatest(display, value)
                    catch
                        @error "Evaluation succeeded, but an error occurred while displaying the value" typeof(value)
                        rethrow()
                    end
                    println()
                end
            end
            break
        catch
            if errcount > 0
                @error "SYSTEM: display_error(errio, lasterr) caused an error"
            end
            errcount += 1
            lasterr = catch_stack()
            if errcount > 2
                @error "It is likely that something important is broken, and Julia will not be able to continue normally" errcount
                break
            end
        end
    end
    isa(stdin, TTY) && println()
    nothing
end

function _parse_input_line_core(s::String, filename::String)
    ex = ccall(:jl_parse_all, Any, (Ptr{UInt8}, Csize_t, Ptr{UInt8}, Csize_t),
               s, sizeof(s), filename, sizeof(filename))
    if ex isa Expr && ex.head === :toplevel
        if isempty(ex.args)
            return nothing
        end
        last = ex.args[end]
        if last isa Expr && (last.head === :error || last.head === :incomplete)
            # if a parse error happens in the middle of a multi-line input
            # return only the error, so that none of the input is evaluated.
            return last
        end
    end
    return ex
end

function parse_input_line(s::String; filename::String="none", depwarn=true)
    # For now, assume all parser warnings are depwarns
    ex = if depwarn
        _parse_input_line_core(s, filename)
    else
        with_logger(NullLogger()) do
            _parse_input_line_core(s, filename)
        end
    end
    return ex
end
parse_input_line(s::AbstractString) = parse_input_line(String(s))

function parse_input_line(io::IO)
    s = ""
    while !eof(io)
        s *= readline(io, keep=true)
        e = parse_input_line(s)
        if !(isa(e,Expr) && e.head === :incomplete)
            return e
        end
    end
end

# detect the reason which caused an :incomplete expression
# from the error message
# NOTE: the error messages are defined in src/julia-parser.scm
incomplete_tag(ex) = :none
function incomplete_tag(ex::Expr)
    Meta.isexpr(ex, :incomplete) || return :none
    msg = ex.args[1]
    occursin("string", msg) && return :string
    occursin("comment", msg) && return :comment
    occursin("requires end", msg) && return :block
    occursin("\"`\"", msg) && return :cmd
    occursin("character", msg) && return :char
    return :other
end

# call include() on a file, ignoring if not found
include_ifexists(mod::Module, path::AbstractString) = isfile(path) && include(mod, path)

function exec_options(opts)
    if !isempty(ARGS)
        idxs = findall(x -> x == "--", ARGS)
        length(idxs) > 0 && deleteat!(ARGS, idxs[1])
    end
    quiet                 = (opts.quiet != 0)
    startup               = (opts.startupfile != 2)
    history_file          = (opts.historyfile != 0)
    color_set             = (opts.color != 0) # --color!=auto
    global have_color     = (opts.color == 1) # --color=on
    global is_interactive = (opts.isinteractive != 0)

    # pre-process command line argument list
    arg_is_program = !isempty(ARGS)
    repl = !arg_is_program
    cmds = unsafe_load_commands(opts.commands)
    for (cmd, arg) in cmds
        if cmd == 'e'
            arg_is_program = false
            repl = false
        elseif cmd == 'E'
            arg_is_program = false
            repl = false
        elseif cmd == 'L'
            # nothing
        else
            @warn "Unexpected command -$cmd'$arg'"
        end
    end

    # remove filename from ARGS
    global PROGRAM_FILE = arg_is_program ? popfirst!(ARGS) : ""

    # Load Distributed module only if any of the Distributed options have been specified.
    distributed_mode = (opts.worker == 1) || (opts.nprocs > 0) || (opts.machine_file != C_NULL)
    if distributed_mode
        let Distributed = require(PkgId(UUID((0x8ba89e20_285c_5b6f, 0x9357_94700520ee1b)), "Distributed"))
            Core.eval(Main, :(const Distributed = $Distributed))
            Core.eval(Main, :(using .Distributed))
        end

        invokelatest(Main.Distributed.process_opts, opts)
    end

    # load ~/.julia/config/startup.jl file
    startup && load_julia_startup()

    # process cmds list
    for (cmd, arg) in cmds
        if cmd == 'e'
            Core.eval(Main, parse_input_line(arg))
        elseif cmd == 'E'
            invokelatest(show, Core.eval(Main, parse_input_line(arg)))
            println()
        elseif cmd == 'L'
            # load file immediately on all processors
            if !distributed_mode
                include(Main, arg)
            else
                # TODO: Move this logic to Distributed and use a callback
                @sync for p in invokelatest(Main.procs)
                    @async invokelatest(Main.remotecall_wait, include, p, Main, arg)
                end
            end
        end
    end

    # load file
    if arg_is_program
        # program
        if !is_interactive
            ccall(:jl_exit_on_sigint, Cvoid, (Cint,), 1)
        end
        try
            include(Main, PROGRAM_FILE)
        catch
            invokelatest(display_error, catch_stack())
            if !is_interactive
                exit(1)
            end
        end
    end
    repl |= is_interactive
    if repl
        interactiveinput = isa(stdin, TTY)
        if interactiveinput
            global is_interactive = true
            banner = (opts.banner != 0) # --banner!=no
        else
            banner = (opts.banner == 1) # --banner=yes
        end
        run_main_repl(interactiveinput, quiet, banner, history_file, color_set)
    end
    nothing
end

function load_julia_startup()
    # If the user built us with a specific Base.SYSCONFDIR, check that location first for a startup.jl file
    #   If it is not found, then continue on to the relative path based on Sys.BINDIR
    BINDIR = Sys.BINDIR::String
    SYSCONFDIR = Base.SYSCONFDIR::String
    if !isempty(SYSCONFDIR) && isfile(joinpath(BINDIR, SYSCONFDIR, "julia", "startup.jl"))
        include(Main, abspath(BINDIR, SYSCONFDIR, "julia", "startup.jl"))
    else
        include_ifexists(Main, abspath(BINDIR, "..", "etc", "julia", "startup.jl"))
    end
    !isempty(DEPOT_PATH) && include_ifexists(Main, abspath(DEPOT_PATH[1], "config", "startup.jl"))
    return nothing
end

const repl_hooks = []

"""
    atreplinit(f)

Register a one-argument function to be called before the REPL interface is initialized in
interactive sessions; this is useful to customize the interface. The argument of `f` is the
REPL object. This function should be called from within the `.julia/config/startup.jl`
initialization file.
"""
atreplinit(f::Function) = (pushfirst!(repl_hooks, f); nothing)

function __atreplinit(repl)
    for f in repl_hooks
        try
            f(repl)
        catch err
            showerror(stderr, err)
            println(stderr)
        end
    end
end
_atreplinit(repl) = invokelatest(__atreplinit, repl)

# The REPL stdlib hooks into Base using this Ref
const REPL_MODULE_REF = Ref{Module}()

# run the requested sort of evaluation loop on stdio
function run_main_repl(interactive::Bool, quiet::Bool, banner::Bool, history_file::Bool, color_set::Bool)
    global active_repl
    # load interactive-only libraries
    if !isdefined(Main, :InteractiveUtils)
        try
            let InteractiveUtils = require(PkgId(UUID(0xb77e0a4c_d291_57a0_90e8_8db25a27a240), "InteractiveUtils"))
                Core.eval(Main, :(const InteractiveUtils = $InteractiveUtils))
                Core.eval(Main, :(using .InteractiveUtils))
            end
        catch ex
            @warn "Failed to import InteractiveUtils into module Main" exception=(ex, catch_backtrace())
        end
    end

    if interactive && isassigned(REPL_MODULE_REF)
        invokelatest(REPL_MODULE_REF[]) do REPL
            term_env = get(ENV, "TERM", @static Sys.iswindows() ? "" : "dumb")
            term = REPL.Terminals.TTYTerminal(term_env, stdin, stdout, stderr)
            color_set || (global have_color = REPL.Terminals.hascolor(term))
            banner && Base.banner(term)
            if term.term_type == "dumb"
                active_repl = REPL.BasicREPL(term)
                quiet || @warn "Terminal not fully functional"
            else
                active_repl = REPL.LineEditREPL(term, have_color, true)
                active_repl.history_file = history_file
            end
            # Make sure any displays pushed in .julia/config/startup.jl ends up above the
            # REPLDisplay
            pushdisplay(REPL.REPLDisplay(active_repl))
            _atreplinit(active_repl)
            REPL.run_repl(active_repl, backend->(global active_repl_backend = backend))
        end
    else
        # otherwise provide a simple fallback
        if interactive && !quiet
            @warn "REPL provider not available: using basic fallback"
        end
        banner && Base.banner()
        let input = stdin
            if isa(input, File) || isa(input, IOStream)
                # for files, we can slurp in the whole thing at once
                ex = parse_input_line(read(input, String))
                if Meta.isexpr(ex, :toplevel)
                    # if we get back a list of statements, eval them sequentially
                    # as if we had parsed them sequentially
                    for stmt in ex.args
                        eval_user_input(stderr, stmt, true)
                    end
                    body = ex.args
                else
                    eval_user_input(stderr, ex, true)
                end
            else
                while isopen(input) || !eof(input)
                    if interactive
                        print("julia> ")
                        flush(stdout)
                    end
                    try
                        eval_user_input(stderr, parse_input_line(input), true)
                    catch err
                        isa(err, InterruptException) ? print("\n\n") : rethrow()
                    end
                end
            end
        end
    end
    nothing
end

# MainInclude exists to hide Main.include and eval from `names(Main)`.
baremodule MainInclude
using ..Base
# We inline the definition of include from loading.jl/include_relative to get one-frame stacktraces.
# include(fname::AbstractString) = Main.Base.include(Main, fname)
function include(fname::AbstractString)
    mod = Main
    isa(fname, String) || (fname = Base.convert(String, fname))
    path, prev = Base._include_dependency(mod, fname)
    for callback in Base.include_callbacks # to preserve order, must come before Core.include
        Base.invokelatest(callback, mod, path)
    end
    tls = Base.task_local_storage()
    tls[:SOURCE_PATH] = path
    local result
    try
        result = ccall(:jl_load_, Any, (Any, Any), mod, path)
    finally
        if prev === nothing
            Base.delete!(tls, :SOURCE_PATH)
        else
            tls[:SOURCE_PATH] = prev
        end
    end
    return result
end
eval(x) = Core.eval(Main, x)
end

"""
    eval(expr)

Evaluate an expression in the global scope of the containing module.
Every `Module` (except those defined with `baremodule`) has its own 1-argument
definition of `eval`, which evaluates expressions in that module.
"""
MainInclude.eval

"""
    include(path::AbstractString)

Evaluate the contents of the input source file in the global scope of the containing module.
Every module (except those defined with `baremodule`) has its own 1-argument
definition of `include`, which evaluates the file in that module.
Returns the result of the last evaluated expression of the input file. During including,
a task-local include path is set to the directory containing the file. Nested calls to
`include` will search relative to that path. This function is typically used to load source
interactively, or to combine files in packages that are broken into multiple source files.

Use [`Base.include`](@ref) to evaluate a file into another module.
"""
MainInclude.include

function _start()
    empty!(ARGS)
    append!(ARGS, Core.ARGS)
    if ccall(:jl_generating_output, Cint, ()) != 0 && JLOptions().incremental == 0
        # clear old invalid pointers
        PCRE.__init__()
    end
    try
        exec_options(JLOptions())
    catch
        invokelatest(display_error, catch_stack())
        exit(1)
    end
    if is_interactive && have_color
        print(color_normal)
    end
end

ext: jl


text: Dict{String,Any}("longpi" => Dict{String,Any}("value" => "3.141592653589793","type" => "float"),"neglongpi" => Dict{String,Any}("value" => "-3.141592653589793","type" => "float"))

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
"""
This is an empty module.
Run the script MarketBigPictureWatch_run.jl instead.
"""
module MarketBigPictureWatch

end


text:
@static if Sys.iswindows()
    using WinRPM
    WinRPM.update()
end


text:
Name = "coname"
DividendYieldPercent = "yie"
LongTermDebtToEquity = "qto"
MarketCapitalizationInMillion = "mkt"
NetProfitMarginPercent = "qpm"
OneDayPriceChangePercent = "prl"
PriceEarningsRatio = "pee"
PriceToBookValue = "pri"
PriceToFreeCashFlow = "prf"
ReturnOnEquityPercent = "ttm"


text:
#Some of these tests require building the full system to have a valid PM object
@testset "AC Power Flow Monitored Line Flow Constraints and bounds" begin
    system = build_system("c_sys5_ml")
    line = PSY.get_component(Line, system, "1")
    PSY.convert_component!(MonitoredLine, line, system)
    devices = Dict{Symbol, DeviceModel}(
        :Generators => DeviceModel(ThermalStandard, ThermalDispatch),
        :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
    )
    branches = Dict{Symbol, DeviceModel}(
        :ML => DeviceModel(MonitoredLine, StaticLineBounds),
        :L => DeviceModel(Line, StaticLineBounds),
    )
    template = OperationsProblemTemplate(StandardPTDFModel, devices, branches, services)
    limits = PSY.get_flowlimits(PSY.get_component(MonitoredLine, system, "1"))
    op_problem_m = OperationsProblem(
        PSI.GenericOpProblem,
        template,
        system;
        optimizer = OSQP_optimizer,
        PTDF = build_PTDF5(),
    )
    for b in PSI.get_variable(op_problem_m.psi_container, :Fp__Line)
        @test JuMP.has_lower_bound(b)
        @test JuMP.has_upper_bound(b)
    end
    for b in PSI.get_variable(op_problem_m.psi_container, :Fp__MonitoredLine)
        @test JuMP.has_lower_bound(b)
        @test JuMP.has_upper_bound(b)
    end
    monitored = solve!(op_problem_m)
    flow = monitored.variable_values[:Fp__MonitoredLine][1, 1]
    @test isapprox(flow, limits.from_to, atol = 1e-3)
end

@testset "AC Power Flow Monitored Line Flow Constraints" begin
    system = build_system("c_sys5_ml")
    devices = Dict{Symbol, DeviceModel}(
        :Generators => DeviceModel(ThermalStandard, ThermalDispatch),
        :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
    )
    branches = Dict{Symbol, DeviceModel}(
        :ML => DeviceModel(MonitoredLine, FlowMonitoredLine),
        :L => DeviceModel(Line, StaticLineBounds),
    )
    template = OperationsProblemTemplate(ACPPowerModel, devices, branches, services)
    line = PSY.get_component(Line, system, "1")
    PSY.convert_component!(MonitoredLine, line, system)
    line = PSY.get_component(MonitoredLine, system, "1")
    limits = PSY.get_flowlimits(line)
    op_problem_m =
        OperationsProblem(TestOpProblem, template, system; optimizer = ipopt_optimizer)
    monitored = solve!(op_problem_m)
    fq = monitored.variable_values[:FqFT__MonitoredLine][1, 1]
    fp = monitored.variable_values[:FpFT__MonitoredLine][1, 1]
    flow = sqrt((fp[1])^2 + (fq[1])^2)
    @test isapprox(flow, limits.from_to, atol = 1e-3)
end

@testset "DC PowerFlow Monitored Line Branch Flow constraints" begin
    devices = Dict{Symbol, DeviceModel}(
        :Generators => DeviceModel(ThermalStandard, ThermalDispatch),
        :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
    )
    branches = Dict{Symbol, DeviceModel}(
        :ML => DeviceModel(MonitoredLine, FlowMonitoredLine),
        :L => DeviceModel(Line, StaticLineBounds),
    )
    template = OperationsProblemTemplate(DCPPowerModel, devices, branches, services)
    system = build_system("c_sys5_ml")
    line = PSY.get_component(Line, system, "1")
    PSY.convert_component!(MonitoredLine, line, system)
    line = PSY.get_component(MonitoredLine, system, "1")
    limits = PSY.get_flowlimits(PSY.get_component(MonitoredLine, system, "1"))
    rate = PSY.get_rate(PSY.get_component(MonitoredLine, system, "1"))
    op_problem_m =
        OperationsProblem(TestOpProblem, template, system; optimizer = ipopt_optimizer)
    monitored = solve!(op_problem_m)
    fp = monitored.variable_values[:Fp__MonitoredLine][1, 1]
    @test isapprox(fp, limits.from_to, atol = 1e-3)
    @test isapprox(fp, rate, atol = 1e-3)
end


text:
# %% imports
using Cambrian
using EvolutionaryStrategies
using JLD
using PyCall
using Statistics
gym = pyimport("gym")
include("../Modules/Networks/net_arch.jl")
include("../Modules/Poisson/PoissonStateSpike.jl")
include("cartpole_fitness.jl")

# %% read the results file
#checker results
elites = load("Julia/Actors/cartpoleElites/sNESelites_cartpole_1.jld", "elites")
# @show typeof(elites)

# %% evaluate functions

function create_net(elites::Array{AbstractESIndividual,1}, indiv::Int64)
    params = YAML.load(open("Julia/Actors/cartpole_cfglif.yml"))
    arch = params["arch"]

    genes = elites[indiv].genes

    # Input linear combination
    genes_alpha = map(x->borne(x,-1.0, 1.0),genes[1:arch[1]*4])
    matalpha = reshape(genes_alpha, (arch[1],4))

    # Weights of the network
    genes_w = genes[arch[1]*4+1:end]

    # Initialize network
    n = Network(arch, params)
    exc, inh = ei_idx(arch)
    set_weight!(n,genes_w, exc, inh)

    return n, matalpha
end

# run one episode with no seeds
function sim(elites::Array{AbstractESIndividual,1}, indiv::Int64)
    net, matalpha = create_net(elites, indiv)
    reward = play_episode(net, matalpha, rng_seed = false)
    return reward
end

# %% evaluate the individuals

# test robustness

function test_random(elites::Array{AbstractESIndividual,1})
    mean_reward = Float64[]
    for i in 1:length(elites)
        tmp_reward = []
        for ep in 1:30
            reward = sim(elites,i)
            push!(tmp_reward, reward)
        end
        push!(mean_reward, mean(tmp_reward))
    end
    return mean_reward
end

function test_random(elites::Array{AbstractESIndividual,1}, indiv::Int64)
    tmp_reward = []
    for ep in 1:30
        reward = sim(res,i)
        push!(tmp_reward, reward)
    end
    mean_r = mean(tmp_reward)
    return mean_r
end

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
#####
##### `show` methods
#####

function Base.show(io::IO, w::TimeSpan)
    start_string = format_duration(first(w))
    stop_string = format_duration(last(w))
    return print(io, "TimeSpan(", start_string, ", ", stop_string, ')')
end

function Base.show(io::IO, samples::Samples)
    return if get(io, :compact, false)
        print(io, "Samples(", summary(samples.data), ')')
    else
        duration_in_seconds = size(samples.data, 2) / samples.signal.sample_rate
        duration_in_nanoseconds = round(Int, duration_in_seconds * 1_000_000_000)
        println(io, "Samples (", format_duration(duration_in_nanoseconds), "):")
        println(io, "  signal.channel_names: ",
                channel_names_string(samples.signal.channel_names))
        println(io, "  signal.start_nanosecond: ", samples.signal.start_nanosecond, " (",
                format_duration(samples.signal.start_nanosecond), ")")
        println(io, "  signal.stop_nanosecond: ", samples.signal.stop_nanosecond, " (",
                format_duration(samples.signal.stop_nanosecond), ")")
        println(io, "  signal.sample_unit: ", repr(samples.signal.sample_unit))
        println(io, "  signal.sample_resolution_in_unit: ",
                samples.signal.sample_resolution_in_unit)
        println(io, "  signal.sample_offset_in_unit: ",
                samples.signal.sample_offset_in_unit)
        println(io, "  signal.sample_type: ", samples.signal.sample_type)
        println(io, "  signal.sample_rate: ", samples.signal.sample_rate, " Hz")
        println(io, "  signal.file_extension: ", repr(samples.signal.file_extension))
        println(io, "  signal.file_options: ", repr(samples.signal.file_options))
        println(io, "  encoded: ", samples.encoded)
        println(io, "  data:")
        show(io, "text/plain", samples.data)
    end
end

function Base.show(io::IO, signal::Signal)
    return if get(io, :compact, false)
        print(io, "Signal(", channel_names_string(signal.channel_names), ")")
    else
        println(io, "Signal:")
        println(io, "  channel_names: ", channel_names_string(signal.channel_names))
        println(io, "  start_nanosecond: ", signal.start_nanosecond, " (",
                format_duration(signal.start_nanosecond), ")")
        println(io, "  stop_nanosecond: ", signal.stop_nanosecond, " (",
                format_duration(signal.stop_nanosecond), ")")
        println(io, "  sample_unit: ", repr(signal.sample_unit))
        println(io, "  sample_resolution_in_unit: ", signal.sample_resolution_in_unit)
        println(io, "  sample_offset_in_unit: ", signal.sample_offset_in_unit)
        println(io, "  sample_type: ", signal.sample_type)
        println(io, "  sample_rate: ", signal.sample_rate, " Hz")
        println(io, "  file_extension: ", repr(signal.file_extension))
        print(io, "  file_options: ", repr(signal.file_options))
    end
end

function Base.show(io::IO, recording::Recording)
    return if get(io, :compact, false)
        duration_string = isempty(recording.signals) ? "<no signals>" :
                          format_duration(duration(recording))
        print(io, "Recording(", duration_string, ')')
    else
        if isempty(recording.signals)
            duration_string = "<no signals>"
        else
            duration_string = string(format_duration(duration(recording)), "; ",
                                     duration(recording).value / 1_000_000_000, " seconds")
        end
        println(io, "Recording (", duration_string, ')')
        println(io, "  signals:")
        compact_io = IOContext(io, :compact => true)
        for (name, signal) in recording.signals
            println(compact_io, "    :", name, " => ", signal)
        end
        print(io, "  annotations: (", length(recording.annotations), " total)")
    end
end

function Base.show(io::IO, dataset::Dataset)
    return print(io, "Dataset(", dataset.path, ", ", length(dataset.recordings),
                 " recordings)")
end

#####
##### utilities
#####

function channel_names_string(channel_names)
    return string('[', join(map(repr, channel_names), ", "), ']')
end

function nanosecond_to_periods(ns::Integer)
    μs, ns = divrem(ns, 1000)
    ms, μs = divrem(μs, 1000)
    s, ms = divrem(ms, 1000)
    m, s = divrem(s, 60)
    hr, m = divrem(m, 60)
    return (hr, m, s, ms, μs, ns)
end

format_duration(t::Period) = format_duration(convert(Nanosecond, t).value)

function format_duration(ns::Integer)
    hr, m, s, ms, μs, ns = nanosecond_to_periods(ns)
    hr = lpad(hr, 2, '0')
    m = lpad(m, 2, '0')
    s = lpad(s, 2, '0')
    ms = lpad(ms, 3, '0')
    μs = lpad(μs, 3, '0')
    ns = lpad(ns, 3, '0')
    return string(hr, ':', m, ':', s, '.', ms, μs, ns)
end


text:
module Matching
import Base: >=, <=, >, <, copy, show

export
  Agent, PreferenceRanking, Preference, Agent, RankedPreferences, prefers, add_rank!,
  gale_shapley, ttc_with_counter

include("types.jl")
include("deferred_acceptance.jl")
include("ttc.jl")

end


text:
"""
```
Module      : FinanceLib.FixedIncomes.Bonds
Description : Implement Fixed Incomes modules for the FinanceLib library
Copyright   : (c) 2022 Kishaloy Neogi
License     : MIT
Maintainer  : Kishaloy Neogi
Email       : nkishaloy@yahoo.com
```
The module describes the base modules of Bonds.

You may see the github repository at <https://github.com/n-kishaloy/FinanceLib.jl>
"""
module Bonds

import FinanceLib as Fl
import Roots

include("Rates.jl")

import FinanceLib.FixedIncomes.Bonds.Rates as Rt

"""
`CouponBond : struct defining a Coupon bond`
* c     = Coupon per period
* freq  = Frequency of coupon payment per period
* T     = Life of the Bond
"""
struct CouponBond
  c     :: Float64
  freq  :: Int64
  T     :: Float64
end

"""
`priceCouponBond(r, bd) = Price of coupon bonds`

* r   = rate as Float64, where rate is constant
* bd  = CouponBond
"""
priceCouponBond(r::Float64, bd) = (bd.c/bd.freq) * sum( 1.0 ./ (1+r/bd.freq).^(1:(bd.freq*bd.T)) ) + 1/(1+r/bd.freq)^(bd.T * bd.freq)

"""
`priceCouponBond(r, bd) = Price of coupon bonds`

* r   = RateCurve of spot rates
* bd  = CouponBond
"""
priceCouponBond(r::Fl.RateCurve{Fl.NomRate}, bd) = (bd.c/bd.freq) * sum( (i -> 1/(1+ Fl.estimR(r,i/bd.freq)/bd.freq)^i ).(1:(bd.freq*bd.T)) ) + 1/(1 + Fl.estimR(r, bd.T)/bd.freq)^(bd.T * bd.freq)

"""
`ytmCoupleBonds(bd :: CouponBond, P) = YTM of Coupon Bonds`

* bd = CouponBond
* P = Price of CouponBond
"""
function ytmCouponBond(bd::CouponBond, P)
  cf = (i -> (i, bd.c/bd.freq)).(1:(bd.freq*bd.T))
  cf[end] = (cf[end][1], cf[end][2] + 1.0); push!(cf, (0.0, -P))
  Fl.irr(cf)*bd.freq
end

end


text:
# This file is part of Kpax3. License is MIT.

function test_mcmc_logmarglikmerge()
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  settings = Kpax3.KSettings(ifile, ofile)

  data = UInt8[0 0 0 0 0 1;
               1 1 1 1 1 0;
               0 0 1 0 1 1;
               1 1 0 1 0 0;
               1 1 0 0 0 0;
               0 0 0 1 1 0;
               1 1 1 0 0 0;
               0 0 0 1 1 1;
               0 0 1 0 0 0;
               1 0 0 1 0 1;
               0 1 0 0 1 0;
               0 0 0 0 0 1;
               1 1 1 0 0 0;
               0 0 0 1 1 0;
               1 1 0 0 1 1;
               0 0 1 1 0 0;
               1 1 0 1 0 0;
               0 0 1 0 1 1]

  (m, n) = size(data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

  # [1; 1; 1; 2; 3; 3] => [1; 1; 1; 3; 3; 3]
  state1 = Kpax3.AminoAcidState(data, [1; 1; 1; 2; 3; 3], priorR, priorC, settings)
  support1 = Kpax3.MCMCSupport(state1, priorC)

  ni = vec(sum(float(data[:, [4; 5; 6]]), dims=2))
  vi = 3

  Kpax3.updatelogmargliki!(ni, vi, priorC, support1)

  Kpax3.logmarglikmerge!(state1.cl, state1.k, 3, 2, priorC, support1)

  state2 = Kpax3.AminoAcidState(data, [1; 1; 1; 2; 2; 2], priorR, priorC, settings)
  support2 = Kpax3.MCMCSupport(state2, priorC)

  @test isapprox(support1.logmlikcandidate, support2.logmlik, atol=ε)

  nothing
end

test_mcmc_logmarglikmerge()

function test_mcmc_logmargliksplit()
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  settings = Kpax3.KSettings(ifile, ofile)

  data = UInt8[0 0 0 0 0 1;
               1 1 1 1 1 0;
               0 0 1 0 1 1;
               1 1 0 1 0 0;
               1 1 0 0 0 0;
               0 0 0 1 1 0;
               1 1 1 0 0 0;
               0 0 0 1 1 1;
               0 0 1 0 0 0;
               1 0 0 1 0 1;
               0 1 0 0 1 0;
               0 0 0 0 0 1;
               1 1 1 0 0 0;
               0 0 0 1 1 0;
               1 1 0 0 1 1;
               0 0 1 1 0 0;
               1 1 0 1 0 0;
               0 0 1 0 1 1]

  (m, n) = size(data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

  # [1; 1; 1; 2; 2; 2] => [1; 1; 1; 2; 3; 3]
  state1 = Kpax3.AminoAcidState(data, [1; 1; 1; 2; 2; 2], priorR, priorC, settings)
  support1 = Kpax3.MCMCSupport(state1, priorC)

  support1.ni = vec(sum(float(data[:, 4]), dims=2))
  support1.vi = 1

  support1.nj = vec(sum(float(data[:, [5; 6]]), dims=2))
  support1.vj = 2

  Kpax3.updatelogmargliki!(priorC, support1)
  Kpax3.updatelogmarglikj!(priorC, support1)

  Kpax3.logmargliksplit!(state1.cl, state1.k, 2, priorC, support1)

  state2 = Kpax3.AminoAcidState(data, [1; 1; 1; 2; 3; 3], priorR, priorC, settings)
  support2 = Kpax3.MCMCSupport(state2, priorC)

  @test isapprox(support1.logmlikcandidate, support2.logmlik, atol=ε)

  nothing
end

test_mcmc_logmargliksplit()


text:
module RBD

export
SolverBody, SolverConstraint3DofLinear, SolverConstraint1Dof, SolverConstraint1DofAngular, SolverConstraintRigid,
solve, projectToManifoldVelocity, projectToR3Velocity!,
setupPointToPointConstraint, setup1DofConstraint, setup1DofAngularConstraint, setupRigidConstraint,
skew, translation, create_transform, transform, inertia_tensor, rotation,
Matrix3x3, Vector3, Matrix4x4, Transform, step_world, RigidBody, to_solver_body,
RigidBodyConstraint, PointToPointConstraint, Linear1DofConstraint, ContactConstraint, Angular1DofConstraint, RigidConstraint,
ndof, bodyA, bodyB, body_ids,
to_solver_constraint, linearJa, linearJb, angularJa, angularJb, rhs

using LinearAlgebra, StaticArrays, ReferenceFrameRotations, SparseArrays

Scalar = Float64
Vector3 = SVector{3, Scalar}
Vector6 = SVector{6, Scalar}
Matrix3x3 = SMatrix{3,3, Scalar, 9}
Matrix4x4 = SMatrix{4,4, Scalar, 16}
Matrix6x6 = SMatrix{6,6, Scalar, 36}
Matrix3x6 = SMatrix{3,6, Scalar, 18}
Matrix6x3 = SMatrix{6,3, Scalar, 18}
Transform = Matrix4x4

transform(A::AT, v::VT) where{T, S, AT<:SMatrix{4,4,T,16}, VT<:SVector{3,S}} = SVector{3,S}((A*vcat(v, 1))[1:3])
transform(A::AT, v::VT) where{T, S, AT<:Matrix{T}, VT<:SVector{3,S}} = SVector{3,S}((A*vcat(v, 1))[1:3])
rotate(A::AT, v) where{T, AT<:SMatrix{4,4,T,16}} = SVector{3,T}(A[1:3,1:3]*v)
rotate(A::AT, v) where{T, AT<:Matrix{T}} = SVector{3,T}(A[1:3,1:3]*v)
translation(wTb) = Vector3(wTb[4,1:3])
rotation(A::Transform) = Matrix3x3(A[1:3,1:3])

create_transform(R, t) = vcat(hcat(R, zeros(3)), hcat(t', 1)) |> Transform

Base.:>>>(q::T, v) where {T<:Quaternion} = vect(q*v*inv(q))

skew(ω::Vector3) = [0.0 -ω[3] ω[2];
         ω[3] 0.0  -ω[1];
         -ω[2] ω[1] 0.0]

function inertia_tensor(box::Vector3, m::Scalar)
    x,y,z = box
    (m/12.0) * Vector3(y^2+z^2, x^2+z^2, x^2+y^2)
end

function generate_frame(z)
    y = Vector3(0,1,0)
    if abs(y⋅z) < 0.9
        x = normalize(y×z)
        y = z×x
    else
        x = Vector3(1,0,0)
        y = normalize(z×x)
        x = y×z
    end
    (x,y,z)
end

import LinearAlgebra.normalize

normalize(q::Quaternion) = q/norm(q)

##=========== Constraint Solver ==============##

#== Solver types ==#

mutable struct SolverBody
    v::Vector3
    ω::Vector3

    #Iinv::Vector3   # hmmm... if we want in localspace need transform stored
    #priciple_axis_transform::Quaternion   # transform into space where inertia tensor I is diagonal
    wTb::Transform

    Iinv::Matrix3x3   # in worldspace
    massInv::Scalar
end

abstract type SolverConstraint end
abstract type SolverConstraintClamped <: SolverConstraint end

linearJa(c::T) where T<:SolverConstraint = zeros(3)
linearJb(c::T) where T<:SolverConstraint = zeros(3)
angularJa(c::T) where T<:SolverConstraint = zeros(3)
angularJb(c::T) where T<:SolverConstraint = zeros(3)
ndof(c::T) where T<:SolverConstraint = 1
bodyA(c::T) where T<:SolverConstraint = c.bodyA
bodyB(c::T) where T<:SolverConstraint = c.bodyB
rhs(c::T) where T<:SolverConstraint = c.manifoldVelocityToCorrect + projectToManifoldVelocity(c)

mutable struct SolverConstraint1Dof <: SolverConstraint

    manifoldVelocityToCorrect::Scalar
    solutionImpulse::Scalar

    bodyA::SolverBody
    bodyB::SolverBody

    Ja::Vector3  # Jacobian for body A
    Jb::Vector3  # Jacobian for body B

    axis::Vector3  # world space axis of constraint

    Dinv::Scalar # effective mass inverse

end

SolverConstraint1Dof(bodyA, bodyB) = SolverConstraint1Dof(0, 0, bodyA, bodyB, Vector3(zeros(3)), Vector3(zeros(3)), Vector3(zeros(3)), 0)

linearJa(c::T) where T<:SolverConstraint1Dof = c.axis
linearJb(c::T) where T<:SolverConstraint1Dof = -c.axis
angularJa(c::T) where T<:SolverConstraint1Dof = c.Ja
angularJb(c::T) where T<:SolverConstraint1Dof = c.Jb

mutable struct SolverConstraint1DofAngular <: SolverConstraint

    manifoldVelocityToCorrect::Scalar
    solutionImpulse::Scalar

    bodyA::SolverBody
    bodyB::SolverBody

    Ja::Vector3  # Jacobian for body A
    Jb::Vector3  # Jacobian for body B

    Dinv::Scalar # effective mass inverse

end

SolverConstraint1DofAngular(bodyA, bodyB) = SolverConstraint1DofAngular(0, 0, bodyA, bodyB, Vector3(zeros(3)), Vector3(zeros(3)), 0)

linearJa(c::T) where T<:SolverConstraint1DofAngular = zeros(3)
linearJb(c::T) where T<:SolverConstraint1DofAngular = zeros(3)
angularJa(c::T) where T<:SolverConstraint1DofAngular = c.Ja
angularJb(c::T) where T<:SolverConstraint1DofAngular = c.Jb

mutable struct SolverConstraintContact <: SolverConstraintClamped
    normal_constraint::SolverConstraint1Dof
    friction_constraint_id::Int64
end

SolverConstraintContact(bodyA::T, bodyB::T) where T<:SolverBody = SolverConstraintContact(
    SolverConstraint1Dof(0, 0, bodyA, bodyB, Vector3(zeros(3)), Vector3(zeros(3)), Vector3(zeros(3)), 0), 0)

unclamped_constraint(c:: SolverConstraintContact) = c.normal_constraint

mutable struct SolverConstraint3DofLinear <: SolverConstraint

    manifoldVelocityToCorrect::Vector3
    solutionImpulse::Vector3

    bodyA::SolverBody
    bodyB::SolverBody

    angularJa::Matrix3x3  # Jacobian for body A
    angularJb::Matrix3x3  # Jacobian for body B

    Dinv::Matrix3x3 # effective mass inverse

end

SolverConstraint3DofLinear(bodyA, bodyB) = SolverConstraint3DofLinear(
    Vector3(zeros(3)), Vector3(zeros(3)),
    bodyA, bodyB,
    Matrix3x3(zeros(3,3)), Matrix3x3(zeros(3,3)), Matrix3x3(zeros(3,3)))

ndof(c::SolverConstraint3DofLinear) = 3
linearJa(c::T) where T<:SolverConstraint3DofLinear = Matrix3x3(diagm(ones(3)))
linearJb(c::T) where T<:SolverConstraint3DofLinear = -Matrix3x3(diagm(ones(3)))
angularJa(c::T) where T<:SolverConstraint3DofLinear = c.angularJa
angularJb(c::T) where T<:SolverConstraint3DofLinear = c.angularJb

mutable struct SolverConstraintRigid <: SolverConstraint

    manifoldVelocityToCorrect::Vector6
    solutionImpulse::Vector6

    bodyA::SolverBody
    bodyB::SolverBody

    angularJa::Matrix3x3  # Jacobian for body A
    angularJb::Matrix3x3  # Jacobian for body B

    frame::Matrix3x3

    Dinv::Matrix6x6 # effective mass inverse

end

SolverConstraintRigid(bodyA, bodyB) = SolverConstraintRigid(
    Vector6(zeros(6)), Vector6(zeros(6)),
    bodyA, bodyB,
    zero(Matrix3x3), zero(Matrix3x3), zero(Matrix3x3), zero(Matrix6x6))

ndof(c::SolverConstraintRigid) = 6
linearJa(c::T) where T<:SolverConstraintRigid = Matrix6x3(vcat(diagm(ones(3)), zeros(3,3)))
linearJb(c::T) where T<:SolverConstraintRigid = Matrix6x3(vcat(-diagm(ones(3)), zeros(3,3)))
angularJa(c::T) where T<:SolverConstraintRigid = Matrix6x3(vcat(c.angularJa, c.frame'))
angularJb(c::T) where T<:SolverConstraintRigid = Matrix6x3(vcat(c.angularJb, -c.frame'))

#== end Solver types ==#

function setup1DofConstraint( dt_inv, wTa::Transform, wTb::Transform, pointOnAInA::Vector3, pointOnBInB::Vector3, worldAxis::Vector3, c::SolverConstraint1Dof )
    r1 = rotate(wTa, pointOnAInA);
    r2 = rotate(wTb, pointOnBInB);

    c.Ja = r1×worldAxis  # non-permissable angular axis
    c.Jb = -r2×worldAxis

    c.Dinv = inv((c.bodyA.Iinv*c.Ja)⋅c.Ja + (c.bodyB.Iinv*c.Jb)⋅c.Jb + c.bodyA.massInv + c.bodyB.massInv)
    c.axis = worldAxis

    x1 = translation(wTa) + r1;
    x2 = translation(wTb) + r2;

    # Compute the positional constraint error (scaled by the Baumgarte coefficient 'beta')
    positionError = (x1 - x2);
    scaleError = dt_inv*0.2;
    positionError = positionError*scaleError;

    c.manifoldVelocityToCorrect = positionError⋅c.axis;

    return c
end

function projectToManifoldVelocity( c::SolverConstraint1Dof )
    p = c.bodyA.v ⋅ c.axis;
    p -= c.bodyB.v ⋅ c.axis;

    p += c.Ja ⋅ c.bodyA.ω;
    p += c.Jb ⋅ c.bodyB.ω;

    return p;
end

function projectToR3Velocity!(  c::SolverConstraint1Dof, dp::Scalar )
    dv = dp * c.bodyA.massInv;
    c.bodyA.v += dv * c.axis;
    dv = dp * c.bodyB.massInv;
    c.bodyB.v -= dv * c.axis;
    dω = c.Ja*dp;
    dω = c.bodyA.Iinv*dω;
    c.bodyA.ω += dω;
    dω = c.Jb*dp;
    dω = c.bodyB.Iinv*dω;
    c.bodyB.ω += dω;
    c
end

"""
    provide two mutually perpendicular axis which are also perpendicular to the axis of the prohibited rotation.  One on body A and one on body B.
    this relation should be met when the constraint is at rest a⟂×b⟂ = constrained
"""
function setup1DofAngularConstraint( dt_inv, perp_axisa::Vector3, perp_axisb::Vector3, c::SolverConstraint1DofAngular )

    axis = normalize(perp_axisa×perp_axisb)
    c.Ja = axis
    c.Jb = -c.Ja

    c.Dinv = inv((c.bodyA.Iinv*c.Ja)⋅c.Ja + (c.bodyB.Iinv*c.Jb)⋅c.Jb)

    # Compute the constraint error (scaled by the Baumgarte coefficient 'beta')
    angularError = perp_axisa⋅perp_axisb;
    scaleError = dt_inv*0.2;
    angularError = angularError*scaleError;

    c.manifoldVelocityToCorrect = angularError;

    return c
end

function projectToManifoldVelocity( c::SolverConstraint1DofAngular )
    p = c.Ja ⋅ c.bodyA.ω;
    p += c.Jb ⋅ c.bodyB.ω;

    return p;
end

function projectToR3Velocity!(  c::SolverConstraint1DofAngular, dp::Scalar )
    dω = c.Ja*dp;
    dω = c.bodyA.Iinv*dω;
    c.bodyA.ω += dω;
    dω = c.Jb*dp;
    dω = c.bodyB.Iinv*dω;
    c.bodyB.ω += dω;
    c
end

#== Contact ==#

function setupContactConstraint( dt_inv, wTa::Transform, wTb::Transform, pointOnAInA::Vector3, pointOnBInB::Vector3, worldAxis::Vector3, c::SolverConstraintContact )
    c.normal_constraint = setup1DofConstraint( dt_inv, wTa, wTb, pointOnAInA, pointOnBInB, worldAxis, c.normal_constraint )
    c
end

projectToManifoldVelocity( c::SolverConstraintContact ) = projectToManifoldVelocity(c.normal_constraint)

function projectToR3Velocity!( c::SolverConstraintContact, dp::Scalar )
    c.normal_constraint = projectToR3Velocity!(c.normal_constraint, dp)
    c
end

#== Point to Point ===#

function projectToManifoldVelocity( c::SolverConstraint3DofLinear )
    p = c.bodyA.v;
    p -= c.bodyB.v;

    ppart = c.bodyA.ω;
    p += c.angularJa*ppart;

    ppart = c.bodyB.ω;
    ppart = c.angularJb*ppart;
    p += ppart;

    return p;
end

function projectToR3Velocity!( c::SolverConstraint3DofLinear, dp::Vector3 )
    dv = dp * c.bodyA.massInv;
    c.bodyA.v += dv;
    dv = dp * c.bodyB.massInv;
    c.bodyB.v -= dv;
    dω = dp;
    dω = c.angularJa'*dω;
    dω = c.bodyA.Iinv*dω;
    c.bodyA.ω += dω;
    dω = dp;
    dω = c.angularJb'*dω;
    dω = c.bodyB.Iinv*dω;
    c.bodyB.ω += dω;
    c
end

projectToR3Velocity!( c::SolverConstraint3DofLinear, dp::V ) where V<:AbstractVector = projectToR3Velocity!(c, Vector3(dp))

# output Dinv and solutionImpulse->zero
# to be more Julian should change to method returning Dinv
function init( c::SolverConstraint3DofLinear )
    Dinv = Matrix3x3(diagm(repeat([c.bodyA.massInv+c.bodyB.massInv], 3)))

    Jt = c.angularJa'
    dpartM = c.angularJa*c.bodyA.Iinv;
    dpartM *= Jt;
    Dinv += dpartM;

    Jt = c.angularJb'
    dpartM = c.angularJb*c.bodyB.Iinv;
    dpartM *= Jt;
    Dinv += dpartM;

    c.Dinv = Matrix3x3(inv(Dinv));

    c.solutionImpulse = zero(c.solutionImpulse);

    return c;
end

function setupPointToPointConstraint( dt_inv, wTa::Transform, wTb::Transform, pointOnAInA::Vector3, pointOnBInB::Vector3, c )
    r1 = rotate(wTa, pointOnAInA);
    r2 = rotate(wTb, pointOnBInB);

    x1 = translation(wTa) + r1;
    x2 = translation(wTb) + r2;

    r1 = -r1;
    c.angularJa = skew(r1)
    c.angularJb = skew(r2)

    # Compute the positional constraint error (scaled by the Baumgarte coefficient 'beta')
    positionError = (x1 - x2);
    scaleError = dt_inv*0.2;
    positionError = positionError*scaleError;

    c.manifoldVelocityToCorrect = positionError;

    return init(c)
end

#== 6DOF all (rigid) fixed ===#

function projectToManifoldVelocity( c::SolverConstraintRigid )
    p = c.bodyA.v;
    p -= c.bodyB.v;

    p += c.angularJa*c.bodyA.ω;
    p += c.angularJb*c.bodyB.ω;

    pω = c.frame'*c.bodyA.ω;
    pω -= c.frame'*c.bodyB.ω;

    return Vector6(p..., pω...);
end

function projectToR3Velocity!( c::SolverConstraintRigid, dj::Vector6 )

    # linear dof part
    dp = Vector3(dj[1:3])
    dpω = Vector3(dj[4:6])
    dv = dp * c.bodyA.massInv;
    c.bodyA.v += dv;
    dv = dp * c.bodyB.massInv;
    c.bodyB.v -= dv;
    dω = c.angularJa'*dp;
    dω = c.bodyA.Iinv*dω;
    c.bodyA.ω += dω;
    dω = c.angularJb'*dp;
    dω = c.bodyB.Iinv*dω;
    c.bodyB.ω += dω;

    # angular dof part
    dω = c.frame*dpω;
    dω = c.bodyA.Iinv*dω;
    c.bodyA.ω += dω;
    dω = -c.frame*dpω;
    dω = c.bodyB.Iinv*dω;
    c.bodyB.ω += dω;

    c
end

projectToR3Velocity!( c::SolverConstraintRigid, dp::V ) where V<:AbstractVector = projectToR3Velocity!(c, Vector6(dp))

# output Dinv and solutionImpulse->zero
# to be more Julian should change to method returning Dinv
function init( c::SolverConstraintRigid )
    D11 = Matrix3x3(diagm(repeat([c.bodyA.massInv+c.bodyB.massInv], 3)))

    F = c.frame;

    Jat = c.angularJa'
    JaI = c.angularJa*c.bodyA.Iinv;
    D11 += JaI*Jat;

    Jbt = c.angularJb'
    JbI = c.angularJb*c.bodyB.Iinv;
    D11 += JbI*Jbt;

    D12 = JaI*F-JbI*F;

    D21 = F'*JaI'-F'*JbI'  # inertia tensor is symmetric so (J*I)' = I*J' which goes in D21

    D22 = F'*c.bodyA.Iinv*F + F'*c.bodyB.Iinv*F;

    D = [hcat(D11, D12) ;
         hcat(D21, D22)]
    if det(D) != 0
        c.Dinv = Matrix6x6(inv(D));
    end

    c.solutionImpulse = zero(c.solutionImpulse);

    return c;
end

function setupRigidConstraint( dt_inv, wTa::Transform, wTb::Transform, pointOnAInA::Vector3, pointOnBInB::Vector3, Fa::Matrix3x3, Fb::Matrix3x3, c )
    r1 = rotate(wTa, pointOnAInA);
    r2 = rotate(wTb, pointOnBInB);

    x1 = translation(wTa) + r1;
    x2 = translation(wTb) + r2;

    r1 = -r1;
    c.angularJa = skew(r1)
    c.angularJb = skew(r2)

    # Compute the positional constraint error (scaled by the Baumgarte coefficient 'beta')
    positionError = (x1 - x2);
    scaleError = dt_inv*0.2;
    positionError = positionError*scaleError;

    # angular part
    xa, ya, za = eachcol(Fa)
    xb, yb, zb = eachcol(Fb)
    angularError = Vector3( ya⋅zb, za⋅xb, xa⋅yb )
    scaleError = dt_inv*0.2;
    angularError = angularError*scaleError;

    c.manifoldVelocityToCorrect = Vector6(positionError..., angularError...);

    xw = normalize(ya×zb)
    yw = normalize(za×xw)
    c.frame = Matrix3x3(hcat(xw, yw, xw×yw))

    return init(c)
end

#== end Constraint type methods ==#

function solve( c::SolverConstraint )
    # project velocities onto constraint manifold. i.e. this gives us the magnitude of the velocity violating the constraint conditions
    manV = projectToManifoldVelocity(c);

    manV += c.manifoldVelocityToCorrect;

    Δimpulse = -c.Dinv*manV;  # momentum now

    projectToR3Velocity!( c, Δimpulse );

    c.solutionImpulse += Δimpulse;
end

# clamps accumulated solution impulse to be > 0
function solve( cc::SolverConstraintClamped )
    c = unclamped_constraint(cc)
    # project velocities onto constraint manifold. i.e. this gives us the magnitude of the velocity violating the constraint conditions
    manV = projectToManifoldVelocity(c);

    manV += c.manifoldVelocityToCorrect;

    Δimpulse = -c.Dinv*manV
    solutionImpulse′ = max.(zero(Scalar), c.solutionImpulse + Δimpulse);
    Δimpulse = solutionImpulse′ - c.solutionImpulse

    projectToR3Velocity!( c, Δimpulse );

    c.solutionImpulse += Δimpulse;
end

##=========== Stepper ==============##

mutable struct RigidBody
    x::Vector3
    m::Scalar
    I::Vector3   # bodyspace
    v::Vector3
    ω::Vector3   # world space
    q::Quaternion{Scalar}
    id::Int64
end

RigidBody(x, m=1.0,inertia=ones(3)) = RigidBody(x, m, inertia,
                                                   zeros(3), zeros(3), Quaternion(1.0,0,0,0), -1)

apply_inertia(q::Quaternion, I, v) = q>>>(I.*inv(q)>>>v)

function ω_implicit(ω₀::Vector3, It, q::Quaternion, dt)
    # one step of newtons method to solve for new angular velocity = f(ω′) = I(ω′-ω)+ω′xIω′*dt = 0
    # df(ω′)/ω′ = I + (1xIω′+ω′xI)*dt
    # df(ω) = I + (ωxI - Iωx1)*dt
    ω = inv(q)>>>ω₀
    Iω = It*ω
    f = ω×Iω*dt
    df = It + (skew(ω)*It - skew(Iω))*dt

    ω′ = ω - df\f
    q>>>ω′
end

function integrate_implicit!(rb::RigidBody, dt::Scalar)
    rb.x += rb.v*dt
    rb.ω = ω_implicit(rb.ω, Diagonal(rb.I), rb.q, dt)
    q2 = rb.q + 0.5*dt*rb.ω*rb.q
    rb.q = q2/norm(q2)

    return rb
end

function step_world(world, dt::Scalar)
    for rb in world
        integrate_implicit!(rb, dt)
    end
end

function step_world(bodies, constraints, dt::Scalar; external_force = rb->zero(Vector3), external_torque = rb->zero(Vector3), constraint_iterations=8, use_global_solver=true)
    # Symplectic Euler.  update velocity before constraints and position after constraints
    # integrate velocity
    for rb in bodies
        rb.v += external_force(rb)*dt
        rb.ω += external_torque(rb)*dt
        rb.ω = ω_implicit(rb.ω, Diagonal(rb.I), rb.q, dt)
    end

    if use_global_solver
        solver_bodies = solve_global(constraints, bodies, dt) |> first
        for (i, bodᵢ) in enumerate(solver_bodies)
            bodies[i].v = bodᵢ.v
            bodies[i].ω = bodᵢ.ω
        end
    else
        solver_bodies = map(to_solver_body, bodies)
        c = map(cᵢ->to_solver_constraint(cᵢ, solver_bodies, dt), constraints)

        # solve constraints for updated velocity
        for i in 1:constraint_iterations
            for cᵢ in c
                solve(cᵢ)
            end
        end
    end

    for (i, bodᵢ) in enumerate(solver_bodies)
        bodies[i].v = bodᵢ.v
        bodies[i].ω = bodᵢ.ω
    end

    # integrate positions with updated velocities
    for rb in bodies
        rb.x += rb.v*dt
        rb.q = normalize(rb.q + 0.5*dt*rb.ω*rb.q)
    end
end

#===========  Entities ==============#

abstract type RigidBodyConstraint end

struct PointToPointConstraint <: RigidBodyConstraint
    bodyA::Int64
    bodyB::Int64
    ra::Vector3  # offset from com to pivot in space of body A
    rb::Vector3
end

struct Linear1DofConstraint <: RigidBodyConstraint
    bodyA::Int64
    bodyB::Int64
    ra::Vector3  # offset from com to pivot in space of body A
    rb::Vector3
    axisb::Vector3   # axis of constraint in body B space
end

struct Angular1DofConstraint <: RigidBodyConstraint
    bodyA::Int64
    bodyB::Int64
    perp_axisa::Vector3  # offset from com to pivot in space of body A
    perp_axisb::Vector3
end

Angular1DofConstraint(bodyA, bodyB, axis::Vector3, qa::Quaternion, qb::Quaternion) = Angular1DofConstraint(bodyA, bodyB, ([inv(qa), inv(qb)].>>>generate_frame(axis)[1:2])...)

struct ContactConstraint <: RigidBodyConstraint
    bodyA::Int64
    bodyB::Int64
    ra::Vector3  # offset from com to pivot in space of body A
    rb::Vector3
    normal::Vector3   # contact normal in body B space
end

"""
    rigidly constrain all degrees of freedom
"""
struct RigidConstraint <: RigidBodyConstraint
    bodyA::Int64
    bodyB::Int64
    ra::Vector3  # offset from com to pivot in space of body A
    rb::Vector3
    aRb::Matrix3x3  # relative frame from b to a
end

body_ids(c::T, rb::Vector{RigidBody}) where T<:RigidBodyConstraint = (rb[c.bodyA].id, rb[c.bodyB].id)
ndof(c::T) where T<:RigidBodyConstraint = 1
ndof(c::T) where T<:PointToPointConstraint = 3
ndof(c::T) where T<:RigidConstraint = 6

function to_solver_body(rb::RigidBody)
    R = quat_to_dcm(rb.q)'
    Iinv_world = R*diagm(inv.(rb.I))*R'

    wTb = create_transform(Array(quat_to_dcm(rb.q)'), rb.x)
    SolverBody(rb.v, rb.ω, wTb, Iinv_world, inv(rb.m))
end

function to_solver_constraint(p2p::PointToPointConstraint, solver_bodies, dt)
    rba = solver_bodies[p2p.bodyA]
    rbb = solver_bodies[p2p.bodyB]
    setupPointToPointConstraint( 1.0/dt, rba.wTb, rbb.wTb, p2p.ra, p2p.rb, SolverConstraint3DofLinear(rba, rbb) )
end

function to_solver_constraint(c::Linear1DofConstraint, solver_bodies, dt)
    rba = solver_bodies[c.bodyA]
    rbb = solver_bodies[c.bodyB]
    setup1DofConstraint( 1.0/dt, rba.wTb, rbb.wTb, c.ra, c.rb, rotate(rbb.wTb, c.axisb), SolverConstraint1Dof(rba, rbb) )
end

function to_solver_constraint(c::Angular1DofConstraint, solver_bodies, dt)
    rba = solver_bodies[c.bodyA]
    rbb = solver_bodies[c.bodyB]
    setup1DofAngularConstraint( 1.0/dt, rotate(rba.wTb, c.perp_axisa), rotate(rbb.wTb, c.perp_axisb), SolverConstraint1DofAngular(rba, rbb) )
end

function to_solver_constraint(c::ContactConstraint, solver_bodies, dt)
    rba = solver_bodies[c.bodyA]
    rbb = solver_bodies[c.bodyB]
    setupContactConstraint( 1.0/dt, rba.wTb, rbb.wTb, c.ra, c.rb, rotate(rbb.wTb, c.normal), SolverConstraintContact(rba, rbb) )
end

function to_solver_constraint(c::RigidConstraint, solver_bodies, dt)
    rba = solver_bodies[c.bodyA]
    rbb = solver_bodies[c.bodyB]

    # the b in wTb is for "body" not body "b". everywhere else here it's for body "b"
    wRb_a = rotation(rba.wTb)*c.aRb
    wRa_b = rotation(rbb.wTb)*c.aRb'
    setupRigidConstraint( 1.0/dt, rba.wTb, rbb.wTb, c.ra, c.rb, wRb_a, wRa_b, SolverConstraintRigid(rba, rbb) )
end

function forcepsd(A, diagfactor=1.0e-8)
    for i in 1:size(A)[1]
        A[i, i+1:end] = A[i+1:end, i]
    end
    A = A + Diagonal(A) * diagfactor
    A
end

#==
end # module

module BlockPSDSolver

using RBD, LinearAlgebra, SparseArrays, StaticArrays, ReferenceFrameRotations
==#

"""
    inefficient full matrix from scratch solve.  for reference/debugging purposes.
"""
function solve_global(constraints::Vector{T}, bodies::Vector{RigidBody}, dt) where T<:RigidBodyConstraint

    rbids = (sort ∘ collect)(mapreduce(x->body_ids(x, bodies), ∪, constraints; init = Set{Int64}()))

    id2active = zeros(Int64, last(rbids))
    id2active[rbids] = 1:length(rbids)

    i2original = Vector{Int64}(undef, length(rbids));
    for i in 1:length(rbids)
        j = 1
        while bodies[i].id != rbids[j]
            j += 1
        end
        i2original[j] = i
    end

    solver_bodies = map(to_solver_body, bodies)
    c = map(cᵢ->to_solver_constraint(cᵢ, solver_bodies, dt), constraints)

    I = zeros(Int64, sum(map(ndof, constraints))*12)
    row = 1
    for cᵢ in constraints
        for dofi in 1:ndof(cᵢ)
            I[((row-1)*12+1):(row*12)] .= row
            row += 1
        end
    end

    J = map(constraints) do cᵢ
            a,b = id2active[[body_ids(cᵢ, bodies)...]]
            a = (a-1)*6
            b = (b-1)*6
            repeat([(+).(a, 1:6)..., (+).(b, 1:6)...], ndof(cᵢ))
        end

    # either 1 or two vectors
    V = map(c) do cᵢ
            nd = ndof(cᵢ)
            Jax = linearJa(cᵢ)
            Jbx = linearJb(cᵢ)
            Jaω = angularJa(cᵢ)
            Jbω = angularJb(cᵢ)
            if nd == 1
                [Jax, Jaω, Jbx, Jbω]
            else
                map(1:nd) do i
                    [Jax[i,:], Jaω[i,:], Jbx[i,:], Jbω[i,:]]
                end
            end

        end

    ∂C∂x = sparse(I, reduce(vcat, J), reduce(vcat, reduce(vcat, reduce(vcat, V))))

    M⁻¹ = sparse(diagm(zeros(6*length(rbids))))
    if length(rbids) != length(i2original)
        throw("failed: length(rbids) == length(i2original)")
    end

    for (i, b) in enumerate(solver_bodies[i2original])
        mi = (i-1)*6+1
        M⁻¹[mi:(mi+2),mi:(mi+2)] = diagm(repeat([b.massInv], 3))
        Ii = mi+3
        M⁻¹[Ii:(Ii+2), Ii:(Ii+2)] = b.Iinv
    end

    A = ∂C∂x*M⁻¹*(∂C∂x') |> factorize

    b = -mapreduce(rhs, vcat, c; init=Vector{Float64}([]))

    # solve for constraint impulses
    λ = A\b

    i = 1
    for cᵢ in c
        nd = ndof(cᵢ)
        if nd == 1
            projectToR3Velocity!(cᵢ, λ[i])
            cᵢ.solutionImpulse += λ[i];
        else
            projectToR3Velocity!(cᵢ, λ[i:(i+nd-1)])
            cᵢ.solutionImpulse += λ[i:(i+nd-1)];
        end
        i += nd
    end

    (solver_bodies, map(cᵢ->cᵢ.solutionImpulse, c))
end

end

#==
using .RBD
using LinearAlgebra, StaticArrays, ReferenceFrameRotations

drawit = true
using MeshCat, GeometryBasics, CoordinateTransformations

dt = 0.0167
box_radius = Vector3(0.05,0.2,0.4)
offset = Vector3(0,0.0,0)

if drawit
    vis = Visualizer(); anim = Animation()
end

function add_box( offset, box_radius, id)
    x = offset
    rb = RigidBody(x)
    rb.q = Quaternion(1.0,0,0,0); rb.ω = [0.0,0,0]; rb.x = x; rb.v = [0.0,0.0,0.0]
    rb.m = 10
    rb.I = inertia_tensor(Vector3(2*box_radius), rb.m)
    rb.id = id
    symname = Symbol("box", id)
    global drawit
    if drawit
        setobject!(vis[symname], Rect(-Vec(box_radius), Vec(box_radius)*2));
        atframe(anim, 0) do
            settransform!(vis[symname], compose(Translation(rb.x...), LinearMap(quat_to_dcm(rb.q)')))
        end
    end
    (rb, symname)
end

function constrain_boxes( rb, rb2, id, id2, p )
    #p2p = PointToPointConstraint(id, id2, p-rb.x, p-rb2.x)
    #c1 = Angular1DofConstraint(id, id2, Vector3(0,0,1), rb.q, rb2.q)
    #c2 = Angular1DofConstraint(id, id2, Vector3(1,0,0), rb.q, rb2.q)
    #c3 = Angular1DofConstraint(id, id2, Vector3(0,1,0), rb.q, rb2.q)
    #[p2p, c1, c2, c3]

    cr = RigidConstraint(id, id2, p-rb.x, p-rb2.x, Matrix3x3(diagm(ones(3))))
    [cr]
end

nbox = 50
bodies = Vector{RigidBody}([])
drawbod = []
constraints = Vector{RigidBodyConstraint}([])
for i in 1:nbox
    rb, symname = add_box(Vector3(0,2*(i-1)*box_radius[2]-nbox*box_radius[2],0), box_radius, i)
    rb.ω = Vector3(0,0,1000)
    global bodies
    global drawbod
    global constraints
    bodies = vcat(bodies, [rb])
    drawbod = vcat(drawbod, (symname, box_radius, rb))

    if i > 1
        constraints = vcat(constraints, constrain_boxes(bodies[i-1], rb, i-1, i, rb.x+Vector3(0,box_radius[2],0)))
    end
end

for i in 1:250
    step_world(bodies, constraints, dt; use_global_solver=(length(constraints)>1), constraint_iterations=32 )
    if drawit
        atframe(anim, i) do
            for i in 1:nbox
                rb = drawbod[i][3]
                settransform!(vis[drawbod[i][1]], compose(Translation(rb.x...), LinearMap(quat_to_dcm(rb.q)')))
            end
        end
    end
end

# DRAW
if drawit
    setanimation!(vis, anim)
end

sleep(3)  # seems to need a bit of time to connect to renderer
==#

#==

==#


text:
#=

  "A puzzle" in ConstraintSolver.jl

  From "God plays dice"
  "A puzzle"
  http://gottwurfelt.wordpress.com/2012/02/22/a-puzzle/

  And the sequel "Answer to a puzzle"
  http://gottwurfelt.wordpress.com/2012/02/24/an-answer-to-a-puzzle/

  This problem instance was taken from the latter blog post.

  """
  8809 = 6
  7111 = 0
  2172 = 0
  6666 = 4
  1111 = 0
  3213 = 0
  7662 = 2
  9312 = 1
  0000 = 4
  2222 = 0
  3333 = 0
  5555 = 0
  8193 = 3
  8096 = 5
  7777 = 0
  9999 = 4
  7756 = 1
  6855 = 3
  9881 = 5
  5531 = 0

  2581 = ?
  """

  Note:
  This model yields 10 solutions, since x[5] is not
  restricted in the constraints.
  All solutions has x assigned to the correct result.

  The problem stated in "A puzzle"
  http://gottwurfelt.wordpress.com/2012/02/22/a-puzzle/
  is
  """
  8809 = 6
  7662 = 2
  9312 = 1
  8193 = 3
  8096 = 5
  7756 = 1
  6855 = 3
  9881 = 5

  2581 = ?
  """

  This problem instance - using the same principle - yields
  two different solutions of x, one is the same (correct) as
  for the above problem instance, and one is not.
  This is because here both x[5] and x[2] are underdefined.

  Note:
  This problem has another non-algebraic and - let's say - topological
  approach which yield the same solution as the first problem and one
  of the two solutions of the second problem.

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function a_puzzle(problem,print_solutions=true,all_solutions=true,timeout=6)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>timeout,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))
    m = problem[:m]
    p = problem[:p] # the unknown

    n = 10
    @variable(model, 0 <= x[1:n] <= 9, Int)
    @variable(model, 0 <= z <= 9, Int)

    for i in 1:length(m)
        num, a = m[i]
        ds = split(num,"").|>d->parse(Int,d)
        @constraint(model, a == sum([x[ds[j]+1] for j in 1:4 ]))
    end
    # Find the unknown
    @constraint(model, z == sum([x[p[j]+1] for j in 1:4 ]))

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    num_sols = 0
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                # println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                z_val = convert.(Integer,JuMP.value.(z; result=sol))
                println("x:$x_val  z:$z_val")

            end
        end
    else
        println("status:$status")
    end

    return status, num_sols
end

problems = Dict(
    :1 => Dict(
        :m =>  [["8809",6],
                ["7111", 0],
                ["2172", 0],
                ["6666", 4],
                ["1111", 0],
                ["3213", 0],
                ["7662", 2],
                ["9312", 1],
                ["0000", 4],
                ["2222", 0],
                ["3333", 0],
                ["5555", 0],
                ["8193", 3],
                ["8096", 5],
                ["7777", 0],
                ["9999", 4],
                ["7756", 1],
                ["6855", 3],
                ["9881", 5],
                ["5531", 0]],
        :p => [2,5,8,1], # The unknown

    ),

    :2 => Dict(
        :m => [
            ["8809",6],
            ["7662",2],
            ["9312",1],
            ["8193",3],
            ["8096",5],
            ["7756",1],
            ["6855",3],
            ["9881",5],
        ],
        :p => [2,5,8,1]

    )

)

println("Problem I:")
@time a_puzzle(problems[:1],true,true)

println("\nProblem II:")
@time a_puzzle(problems[:2],true,true)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using Distributions

using Random
Random.seed!(12); # Set seed for reproducibility

using StatsPlots

p_true = 0.5;

N = 100;

data = rand(Bernoulli(p_true), N);

data[1:5]

prior_belief = Beta(1, 1);

function updated_belief(prior_belief::Beta, data::AbstractArray{Bool})
    # Count the number of heads and tails.
    heads = sum(data)
    tails = length(data) - heads

    # Update our prior belief in closed form (this is possible because we use a conjugate prior).
    return Beta(prior_belief.α + heads, prior_belief.β + tails)
end

# Show updated belief for increasing number of observations
@gif for n in 0:N
    plot(
        updated_belief(prior_belief, data[1:n]);
        size=(500, 250),
        title="Updated belief after $n observations",
        xlabel="probability of heads",
        ylabel="",
        legend=nothing,
        xlim=(0, 1),
        fill=0,
        α=0.3,
        w=3,
    )
    vline!([p_true])
end

using Turing

using MCMCChains

# Unconditioned coinflip model with `N` observations.
@model function coinflip(; N::Int)
    # Our prior belief about the probability of heads in a coin toss.
    p ~ Beta(1, 1)

    # Heads or tails of a coin are drawn from `N` independent and identically
    # distributed Bernoulli distributions with success rate `p`.
    y ~ filldist(Bernoulli(p), N)

    return y
end

# Convenience method for constructing a coinflip model
# that is conditioned on observations `y`.
coinflip(y::AbstractVector{<:Real}) = coinflip(; N=length(y)) | (; y);

model = coinflip(data);

sampler = HMC(0.05, 10);

chain = sample(model, sampler, 1_000; progress=false);

histogram(chain)

@assert isapprox(mean(chain, :p), 0.5; atol=0.1)

# Visualize a blue density plot of the approximate posterior distribution using HMC (see Chain 1 in the legend).
density(chain; xlim=(0, 1), legend=:best, w=2, c=:blue)

# Visualize a green density plot of the posterior distribution in closed-form.
plot!(
    0:0.01:1,
    pdf.(updated_belief(prior_belief, data), 0:0.01:1);
    xlabel="probability of heads",
    ylabel="",
    title="",
    xlim=(0, 1),
    label="Closed-form",
    fill=0,
    α=0.3,
    w=3,
    c=:lightgreen,
)

# Visualize the true probability of heads in red.
vline!([p_true]; label="True probability", c=:red)

max_stars_repo_name: TuringLang/TuringTutorialsOutput

max_issues_repo_name: TuringLang/TuringTutorialsOutput

max_forks_repo_name: TuringLang/TuringTutorialsOutput


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module AccessorsTest

import Dates
using ExtendedDates
using Test

@testset "year/semester/quarter/month/week/day" begin
    # year, month, and day return the indivial components
    # of yearmonthday, avoiding additional calculations when possible
    @test year(-1) == 0
    @test semester(-1) == 2
    @test quarter(-1) == 4
    @test month(-1) == 12
    #    @test week(-1) == 52 # 0 and negative days not working !!
    @test day(-1) == 30

    @test year(0) == 0
    @test semester(0) == 2
    @test quarter(0) == 4
    @test month(0) == 12
    #    @test week(0) == 52 # 0 and negative days not working !!
    @test day(0) == 31

    @test year(1) == 1
    @test semester(1) == 1
    @test quarter(1) == 1
    @test month(1) == 1
    @test week(1) == 1
    @test day(1) == 1

    @test year(730120) == 2000
    @test month(730120) == 1
    @test day(730120) == 1
end

@testset "year, month, day, hour, minute, second over many years" begin
    let from = 1, to = 2, y = 0, s = 0, q = 0, m = 0, w = 1, d = 0
        dd = Dates.dayofweek(Dates.Date(from, 1, 1))
        for y = from:to
            maxweek = ExtendedDates.weeksinyear(y)
            for m = 1:12
                for d = 1:Dates.daysinmonth(y, m)
                    dt = DayDate(y, m, d)
                    s = div(rem(m - 1, 12), 6) + 1
                    q = div(rem(m - 1, 12), 3) + 1
                    @test y == year(dt)
                    @test m == month(dt)
                    @test d == day(dt)
                    @test s == semester(dt)
                    @test q == quarter(dt)
                    @test w == week(dt)
                    if dd > 6
                        w = (w + 1 > maxweek) ? 1 : w + 1
                        dd = 1
                    else
                        dd += 1
                    end
                end
            end
        end
    end
end
@testset "week" begin
    # Tests from https://en.wikipedia.org/wiki/ISO_week_date
    @test week(DayDate(2005, 1, 1)) == 53
    @test week(DayDate(2005, 1, 2)) == 53
    @test week(DayDate(2005, 12, 31)) == 52
    @test week(DayDate(2007, 1, 1)) == 1
    @test week(DayDate(2007, 12, 30)) == 52
    @test week(DayDate(2007, 12, 31)) == 1
    @test week(DayDate(2008, 1, 1)) == 1
    @test week(DayDate(2008, 12, 28)) == 52
    @test week(DayDate(2008, 12, 29)) == 1
    @test week(DayDate(2008, 12, 30)) == 1
    @test week(DayDate(2008, 12, 31)) == 1
    @test week(DayDate(2009, 1, 1)) == 1
    @test week(DayDate(2009, 12, 31)) == 53
    @test week(DayDate(2010, 1, 1)) == 53
    @test week(DayDate(2010, 1, 2)) == 53
    @test week(DayDate(2010, 1, 2)) == 53
    # Tests from http://www.epochconverter.com/date-and-time/weeknumbers-by-year.php?year=1999
    dt = DayDate(1999, 12, 27)
    check = (52, 52, 52, 52, 52, 52, 52, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2)
    for i = 1:21
        @test week(dt) == check[i]
        dt = dt + Dates.Day(1)
    end
    # Tests from http://www.epochconverter.com/date-and-time/weeknumbers-by-year.php?year=2000
    dt = DayDate(2000, 12, 25)
    for i = 1:21
        @test week(dt) == check[i]
        dt = dt + Dates.Day(1)
    end
    # Test from http://www.epochconverter.com/date-and-time/weeknumbers-by-year.php?year=2030
    dt = DayDate(2030, 12, 23)
    for i = 1:21
        @test week(dt) == check[i]
        dt = dt + Dates.Day(1)
    end
    # Tests from http://www.epochconverter.com/date-and-time/weeknumbers-by-year.php?year=2004
    dt1 = DayDate(2004, 12, 20)
    check = (52, 52, 52, 52, 52, 52, 52, 53, 53, 53, 53, 53, 53, 53, 1, 1, 1, 1, 1, 1, 1)
    for i = 1:21
        @test week(dt1) == check[i]
        dt1 = dt1 + Dates.Day(1)
    end
end
@testset "Vectorized accessors" begin
    a = DayDate(2014, 1, 1)
    dr = [a, a, a, a, a, a, a, a, a, a]
    @test year.(dr) == repeat([2014], 10)
    @test month.(dr) == repeat([1], 10)
    @test day.(dr) == repeat([1], 10)
end

end


text:
println("Dv) Test of the PairAnnihilation1Photon module with ASF from an internally generated initial- and final-state multiplet.")

@warn("\n\n !!! This example does not work properly at present !!! \n\n")
#
wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=JAC.Radial.Grid(true), nuclearModel=Nuclear.Model(36.) )

wb = perform(wa)


text:
using SparseSensors
using LinearAlgebra
using Gadfly
using DataFrames

#Setup the experiment
r = 11; # Number of basis modes
n = 1000;
x = collect(0.0:1/n:1.0);
vde_basis = VandermondeBasis(x,r);
Ψ = vde_basis.Ψ; #Get the vandermonde basis matrix from the Basis struct

#Make design matrix
X = copy(transpose(Ψ));
n_samples,n_features = size(X);

#Setup QR pivot sensor placement algorithm
qr_pivot = QRPivot(X);
fit(qr_pivot);
pivots = qr_pivot.pivots;

#Select the top 15 sensor locations
f = abs.(x.^2 .- 0.5);
selected_sensors = get_sensors(pivots,15);
x_sensed = x[selected_sensors];
y_sensed = f[selected_sensors];

#Ground truth
df_true = DataFrame();
df_true[!,"x_true"] = x
df_true[!,"y_true"] = f

#Sensed
df = DataFrame()
df[!,"x_sensed"] = x_sensed;
df[!,"y_sensed"] = y_sensed;

#Plot the results
p1 = plot(df,
    layer(x=:x_sensed,y=:y_sensed,color=["Optimized Sensors"]),
    layer(df_true,x=:x_true,y=:y_true,Geom.line,Geom.point,color=["True Function"]),
    Guide.xlabel("x"),Guide.ylabel("y"),
    Theme(background_color=colorant"white"))

#Save image
draw(SVG("example.svg",7inch,3.326inch),p1)


text:
module LatticeVars

using ..Points
import Distributions
using Distributions: Distribution

export LatticeVal, AbstractLatticeVar, LatticeVar, LatticeVarParam, init!,
    get_num_sites_visited

###
### LatticeVal
###

"""
    LatticeVal{N, T, D}

A container for values of type `T` attached to sites
on an `N`-dimensional hyper-cubic lattice.
Values are stored in `dict::D`, typically a `Dict`.
"""
struct LatticeVal{N, T, D}
    dict::D
end

"""
    LatticeVal{N=1, T=Float64}()

Construct a `LatticeVal` object. `dict` is
of type `Dict{Tuple(Int, Int, ...), T}`.
"""
function LatticeVal{N, T}() where {N, T}
    dict = Dict{Tuple{(Int for i in 1:N)...}, T}()
    return LatticeVal{N, T, typeof(dict)}(dict)
end

LatticeVal{N}() where N = LatticeVal{N, Float64}()
LatticeVal() = LatticeVal{1}()
Base.eltype(v::LatticeVal{<:Any, T}) where T = T

"""
    init!(vals::LatticeVal)

Remove all values stored in `vals`.
"""
init!(vals::LatticeVal) = empty!(vals)

import Base: empty!, length

for f in (:empty!, :length)
    @eval Base.$(f)(vals::LatticeVal) = ($f)(vals.dict)
end

get_num_sites_visited(vals::LatticeVal) = length(vals)

Base.setindex!(vals::LatticeVal, val, inds::Int...) = vals.dict[inds] = val
Base.getindex(vals::LatticeVal, inds::Int...) = vals.dict[inds]
Base.setindex!(vals::LatticeVal, val, p::Point) = setindex!(vals, val, get_coords(p))
Base.getindex(vals::LatticeVal, p::Point) = getindex(vals, get_coords(p))
Base.get!(vals::LatticeVal, k, v) = get!(vals.dict, k, v)

function Base.show(io::IO, v::LatticeVal{N, T}) where {N, T}
    print(io, length(v), "-element LatticeVal{$N, $T}")
end

abstract type AbstractLatticeVar end

get_num_sites_visited(v::AbstractLatticeVar) = get_num_sites_visited(get_vals(v))

for f in (:length, :eltype)
    @eval Base.$(f)(v::AbstractLatticeVar) = ($f)(get_vals(v))
end

init!(v::AbstractLatticeVar) = (init!(get_vals(v)); nothing)
# This should not be exported to avoid method squatting
init!(d::Distribution) = nothing

###
### LatticeVar
###

"""
    LatticeVar{N, Dist, Vals}

Type representing a quenched random variable on an `N`-dimensional lattice.
Field `dist::Dist` is typically a `Distribution`.

Field `vals::Vals` is typically storage for samples indexed by lattice site,
with `Vals <: LatticeVal`.
"""
struct LatticeVar{N, Dist, Vals} <: AbstractLatticeVar
    vals::Vals
    dist::Dist
end

get_vals(lv::LatticeVar) = lv.vals
get_dist(lv::LatticeVar) = lv.dist

"""
    LatticeVar{N=1}(dist::Distribution)

Construct a sample of a set of i.i.d. random variables with distribution `dist`
indexed by site on an `N`-dimensional hyper-cubic lattice.

### Examples
Exponentially distributed random variables on a square lattice.
The value for a site is generated when it is first accessed.
Subsequent access to the same site return the same value.
```julia-repl
julia> lv =  LatticeVar{2}(Exponential());

julia> lv[0,0]
1.3586837894898518

julia> lv[0,0]
1.3586837894898518

julia> lv[0,1]
0.1412406531803109
```
"""
function LatticeVar{N}(dist::Distribution) where N
    vals = LatticeVal{N, eltype(dist)}()
    return LatticeVar{N, typeof(dist), typeof(vals)}(vals, dist)
end
LatticeVar(dist::Distribution) = LatticeVar{1}(dist)

Base.getindex(var::LatticeVar, inds::Int...) = get!(var.vals, inds, rand(var.dist))
Base.getindex(var::LatticeVar, inds::Tuple) = get!(var.vals, inds, rand(var.dist))
Base.getindex(var::LatticeVar, k::Point) = getindex(var, get_coords(k))
Base.setindex!(var::LatticeVar, val, inds...) = (var.vals[inds...] = val)

function Base.show(io::IO, v::LatticeVar{N}) where {N}
    print(io, length(v), "-element LatticeVar{$N, ")
    show(io, v.dist)
    print(io, "}")
end

###
### LatticeVarParam
###

"""
    LatticeVarParam{N, ParamVar, SampleDist} <: AbstractLatticeVar

Type representing a random variable  with a quenched random parameter on an `N`-dimensional lattice.
More precisely, a quenched parameter is associated with each lattice site.
The value at that site is sampled each time it is queried from a distribution
depending on the quenched parameter.

By default, the sampling distribution is `Exponential`.
"""
struct LatticeVarParam{N, ParamVar, SampleDist} <: AbstractLatticeVar
    param_var::ParamVar
    sample_dist::SampleDist
end

function LatticeVarParam{N}(param_dist::Distribution, sample_dist = Distributions.Exponential) where N
    param_var = LatticeVar{N}(param_dist)
    return LatticeVarParam{N, typeof(param_var), typeof(sample_dist)}(param_var, sample_dist)
end

LatticeVarParam(param_dist::Distribution, sample_dist = Distributions.Exponential) = LatticeVarParam{1}(param_dist, sample_dist)

get_vals(lvp::LatticeVarParam) = get_vals(lvp.param_var)
get_param(lvp::LatticeVarParam) = lvp.param_var
get_param_dist(lvp::LatticeVarParam) = get_dist(get_param(lvp))
get_sample_dist(lvp::LatticeVarParam) = lvp.sample_dist

function Base.show(io::IO, v::LatticeVarParam{N}) where {N}
    print(io, length(v), "-element LatticeVarParam{$N, ")
    show(io, get_param(v).dist)
    print(io, ", ")
    show(io, v.sample_dist)
    print(io, "}")
end

Base.getindex(var::LatticeVarParam, inds::Int...) = rand(var.sample_dist(get_param(var)[inds...]))
Base.getindex(var::LatticeVarParam, inds::Tuple) = rand(var.sample_dist(get_param(var)[inds]))
Base.getindex(var::LatticeVarParam, k::Point) = rand(var.sample_dist(get_param(var)[k]))

end # module LatticeVars


text:
# Find some reasonable values for tick marks.
#
# This is basically Wilkinson's ad-hoc scoring method that tries to balance
# tight fit around the data, optimal number of ticks, and simple numbers.
#
# Args:
#   x_min: minimum value occuring in the data.
#   x_max: maximum value occuring in the data.
#
# Returns:
#   A Float64 vector containing tick marks.
#
function optimize_ticks(x_min::Float64, x_max::Float64)
    if x_min == x_max
        return [x_min]
    end

    # TODO: these should perhaps be part of the theme
    const Q = {(1,1), (5, 0.9), (2, 0.7), (25, 0.5), (3, 0.2)}
    const n = length(Q)
    const k_min   = 2
    const k_max   = 10
    const k_ideal = 5

    xspan = x_max - x_min
    z = ceil(log10(xspan))

    high_score = -Inf
    z_best = 0.0
    k_best = 0.0
    r_best = 0.0
    q_best = 0.0

    while k_max * 10^(z+1) > xspan
        for k in k_min:k_max
            for (q, qscore) in Q
                span = (k - 1) * q * 10^z
                if span < xspan
                    continue
                end

                r = ceil((x_max - span) / (q*10^z))
                while r*q*10^z <= x_min
                    has_zero = r <= 0 && abs(r) < k

                    # simplicity
                    s = has_zero ? 1.0 : 0.0

                    # granularity
                    g = 0 < k < 2k_ideal ? 1 - abs(k - k_ideal) / k_ideal : 0.0

                    # coverage
                    c = 1.5 * xspan/span

                    score = (1/4)g + (1/6)s + (1/3)c + (1/4)qscore

                    # strict limits on coverage
                    if span >= 2*xspan || span < xspan
                        score -= 1000
                    end

                    if score > high_score
                        (q_best, r_best, k_best, z_best) = (q, r, k, z)
                        high_score = score
                    end
                    r += 1
                end
            end
        end
        z -= 1
    end

    S = Array(Float64, int(k_best))
    for i in 0:(k_best - 1)
        S[i+1] = (r_best + i) * q_best * 10^z_best
    end

    S
end


text:
using SpecialMatrices
using LinearAlgebra
using SparseArrays

"""Ankara.jl contain the functions to develop the description of the
discrete and finite harmonic oscillator in terms of the Harper equations."""

# Δ1(j) gives the circulant matrix of size d_{j} × d_{j} where d_{j} = 2j + 1
function Δ1(j)
	Δ1 = zeros(Float64, 2j + 1)
	Δ1[1] = -2
	Δ1[2] = 1
	Δ1[end] = 1
	Δ1 = SpecialMatrices.Circulant(Δ1)
end

# Δ2(j) gives the diagonal matrix of size d_{j} × d_{j} whose diagonal entries
# are oscillating terms
function Δ2(j)
	Δ2 = zeros(Float64, 2j + 1)
	for i in -j:j
		Δ2[i + 1 + j] = -4 * (sin(pi * i / (2j)))^2
	end
	return Δ2 = LinearAlgebra.diagm(Δ2)
end

# If necessary, compose the Hamiltonian (Harper's equation) H from who the
# eigenvectors and eigenvalues are obtained
H(j) = -(1/2.0)*(Δ1(j) + Δ2(j))

# Eigenvectors of H
h(j) = LinearAlgebra.eigvecs(H(j));

# Eigenvalues of H. Altough, this eigvenvalues are not strictly used because
# of the symmetry importation below
k(j) = LinearAlgebra.eigvals(H(j));

# Lineal monotonic increasing values who serves as the (new) eigenvalues
# in the symmetry importation to construct the evolution operator
η(j) = [l for l in 0:(2j + 1)];

# Kernel of the evolution
ker(j, t) = exp.(-1im * pi * t * η(j)/2.0)


text:
module TaskMaster

abstract type Learner end
ask!(learner::Learner) = error("For runner learner needs to implement ask! and tell! methods")
tell!(learner::Learner,message) = error("For runner learner needs to implement ask! and tell! methods")

abstract type Master end

include("abstractions.jl")
include("workmaster.jl")
include("evaluate.jl")

export captureslave!, releaseslave!, releaseall!
export HistoryMaster, WorkMaster, Loop, IgnorantLearner, ask!, tell!, evaluate!

export play!, replay!, PlayRes

end


text:
"""
    get_dose_slopes(X, conditionVar, concentrationVar)

Converts categorical variables in a DataFrame to dosage slopes

# Arguments

- X = DataFrame with at least two categorical variables: conditions, specified
  by `conditionVar`, and concentrations, specified by `concentrationVar`
- conditionVar = Symbol for the categorical condition variable in `X` to use
  for dosage slopes.
- concentrationVar = Symbol for the categorical concentration variable in `X`
  to use for dosage slopes.

# Value

DataFrame

# Some notes

`get_dose_slopes` assumes that within each condition, the concentration levels
are ordered.

"""
function get_dose_slopes(X::DataFrames.DataFrame,
                         conditionVar::Symbol, concentrationVar::Symbol)

    # Ensure that conditionVar is a variable in X
    if !in(conditionVar, propertynames(X))
        error(string(conditionVar, " is not a variable in the DataFrame"))
    end

    # Ensure that concentrationVar is a variable in X
    if !in(concentrationVar, propertynames(X))
        error(string(concentrationVar, " is not a variable in the DataFrame"))
    end

    # Pull out every condition
    allConds = unique(X[:,conditionVar])
    # For each condition, pull out all of the concentrations
    concs = [X[:,concentrationVar][cond .== X[:,conditionVar]]
             for cond in allConds]
    # Get the number of unique concentration levels found for each condition
    numLevels = [length(unique(conc)) for conc in concs]

    # Initialize array of slopes
    XDos = zeros(size(X, 1), length(allConds))

    # Iterate through conditions and assign slopes for concentrations
    for i in 1:length(allConds)
        for j in 1:length(unique(concs[i]))
            XDos[(X[:,conditionVar] .== allConds[i]) .&
                 (X[:,concentrationVar] .== unique(concs[i])[j]), i] .= j
        end
    end

    # Convert results to a DataFrame and rename columns
    XDosDf = convert(DataFrame, XDos)
    rename!(XDosDf, [Symbol(cond) for cond in allConds])

    # Initialize new DataFrame
    newX = DataFrame()

	# Iterate through variables in X
    for var in propertynames(X)
        if !in(var, [conditionVar, concentrationVar])
            # Add the other variables to the new DataFrame
            newX[!,var] = X[:,var]
        else
            # Add the DataFrame of dosage slopes
            if var==conditionVar
                for cond in propertynames(XDosDf)
                    newX[!,cond] = XDosDf[:,cond]
                end
            end
        end
    end

    return newX
end


text:
#=

  Secret Santa problem in Julia ConstraintSolver.jl

  From Ruby Quiz Secret Santa
  http://www.rubyquiz.com/quiz2.html
  """
  Honoring a long standing tradition started by my wife's dad, my friends
  all play a Secret Santa game around Christmas time. We draw names and
  spend a week sneaking that person gifts and clues to our identity. On the
  last night of the game, we get together, have dinner, share stories, and,
  most importantly, try to guess who our Secret Santa was. It's a crazily
  fun way to enjoy each other's company during the holidays.

  To choose Santas, we use to draw names out of a hat. This system was
  tedious, prone to many "Wait, I got myself..." problems. This year, we
  made a change to the rules that further complicated picking and we knew
  the hat draw would not stand up to the challenge. Naturally, to solve
  this problem, I scripted the process. Since that turned out to be more
  interesting than I had expected, I decided to share.

  This weeks Ruby Quiz is to implement a Secret Santa selection script.

  Your script will be fed a list of names on STDIN.
  ...
  Your script should then choose a Secret Santa for every name in the list.
  Obviously, a person cannot be their own Secret Santa. In addition, my friends
  no longer allow people in the same family to be Santas for each other and your
  script should take this into account.
  """

  Comment: Well, this model skips the file input and mail parts. We
           assume that the friends are identified with a number from 1..n,
           and the families is identified with a number 1..num_families.

           And we also do some random family configurations below.

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
using Distributions
const CS = ConstraintSolver
include("constraints_utils.jl")

function secret_santa(family,print_solutions=true,all_solutions=false)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>true,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    n = length(family)
    @variable(model, 1 <= x[1:n] <= n, Int)

    # Everyone gives and receives a Secret Santa
    @constraint(model, x in CS.AllDifferent())

    # Can't be one own's Secret Santa
    for i in 1:n
        @constraint(model, x[i] != i)
    end

    # No Secret Santa to a person in the same family
    @variable(model,1 <= fxi[1:n] <= n, Int)
    for i in 1:n
        my_element(model, x[i],family,fxi[i])
        @constraint(model,family[i] != fxi[i])
    end

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        if all_solutions
            println("num_sols:$num_sols\n")
        end

        if print_solutions
            for sol in 1:num_sols
                # println("solution #$sol")
                xx = convert.(Integer,JuMP.value.(x; result=sol))
                println("x:$xx")

            end
        end
        if all_solutions
            println("num_sols:$num_sols\n")
        end
    else
        println("status:$status")
    end

    return status
end

#
# Testing
#

# The families
secret_santa_problems =  Dict(
:1 =>  [1,1,2,2, 3, 4,4], # Ruby Quiz example. 672 solutions
:2 =>  [1,1,1,1, 2, 3,3,3,3,3, 4,4], # From my Picat model
)
print_solutions=false
all_solutions=true
# Number of solutions for the simple instance
@time secret_santa(secret_santa_problems[:1],print_solutions,all_solutions)

# Generate random families

function random_families(n=10)
    d = round(Int,n/2)
    rand(DiscreteUniform(1,d),n)|>sort
end

print_solutions=true
all_solutions=false

# Test larger instances
for n in [i*10 for i in 1:10]
    println("\nn:$(n)")
    family = random_families(n)
    println("family:$family")
    c = 0
    timeout = false
    while c < 10
        @time status = secret_santa(family,print_solutions,all_solutions)
        c += 1
        if status == MOI.TIME_LIMIT
            timeout = true
        end
        if status == MOI.OPTIMAL || timeout || c >= 10
            break
        end
    end
    if timeout
        break
    end
end


text:
"""
Formatting of incoming Redis Replies
"""
function getline(s::TCPSocket)
    l = chomp(readline(s))
    length(l) > 1 || throw(ProtocolException("Invalid response received: $l"))
    return l
end

convert_reply(reply::Array{UInt8}) = String(reply)
convert_reply(reply::Array) = [convert_reply(r) for r in reply]
convert_reply(x) = x

function read_reply(conn::RedisConnectionBase)
    l = getline(conn.socket)
    reply = parseline(l, conn.socket)
    convert_reply(reply)
end

parse_error(l::AbstractString) = throw(ServerException(l))

function parse_bulk_string(s::TCPSocket, slen::Int)
    b = read(s, slen+2) # add crlf
    if length(b) != slen + 2
        throw(ProtocolException(
            "Bulk string read error: expected $len bytes; received $(length(b))"
        ))
    else
        return resize!(b, slen)
    end
end

function parse_array(s::TCPSocket, slen::Int)
    a = Array{Any, 1}(undef, slen)
    for i = 1:slen
        l = getline(s)
        r = parseline(l, s)
        a[i] = r
    end
    return a
end

function parseline(l::AbstractString, s::TCPSocket)
    reply_type = l[1]
    reply_token = l[2:end]
    if reply_type == '+'
        return reply_token
    elseif reply_type == ':'
        parse(Int, reply_token)
    elseif reply_type == '$'
        slen = parse(Int, reply_token)
        if slen == -1
            nothing
        else
            parse_bulk_string(s, slen)
        end
    elseif reply_type == '*'
        slen = parse(Int, reply_token)
        if slen == -1
            nothing
        else
            parse_array(s, slen)
        end
    elseif reply_type == '-'
        parse_error(reply_token)
    end
end

write_token(io::IO, token::Union{AbstractString, Vector{UInt8}}) =
  write(io, "\$", string(sizeof(token)), "\r\n", token, "\r\n")

function write_token(io::IO, token)
  s = string(token)
  write(io, "\$", string(sizeof(s)), "\r\n", s, "\r\n")
end

"""
Formatting of outgoing commands to the Redis server
"""
function pack_command(io::IO, command::Vector)
    b = write(io, "*$(length(command))\r\n")

    for token in command
        b += write_token(io, token)
    end
    b
end

function execute_command_without_reply(conn::RedisConnectionBase, command)
    is_connected(conn) || throw(ConnectionException("Socket is disconnected"))
    pack_command(conn.socket, command)
end

function execute_command(conn::RedisConnectionBase, command::Vector)
    execute_command_without_reply(conn, command)
    read_reply(conn)
end

baremodule SubscriptionMessageType
    const Message = 0
    const Pmessage = 1
    const Subscribe = 2
    const Unsubscribe = 3
    const Other = 4
end

struct SubscriptionMessage
    message_type
    channel::AbstractString
    message::AbstractString

    function SubscriptionMessage(reply::AbstractArray)
        notification = reply
        message_type = notification[1]
        if message_type == "message"
            new(SubscriptionMessageType.Message, notification[2], notification[3])
        elseif message_type == "pmessage"
            new(SubscriptionMessageType.Pmessage, notification[2], notification[4])
        elseif message_type == "subscribe"
            new(SubscriptionMessageType.Subscribe, notification[2], "")
        elseif message_type == "unsubscribe"
            new(SubscriptionMessageType.Unsubscribe, notification[2], "")
        else
            new(SubscriptionMessageType.Other, "", "")
        end
    end
end


text:
# Define abstract interface
"""
$(TYPEDEF)

Abstract supertype for magnetic moments.
"""
abstract type MagneticMoments{T} end

"""
$(SIGNATURES)

Return the mean magnetic moment [Am²] as `SVector{3,T}` for a given magnetic
field H [T/μ₀] (μ₀ being the vacuum permeability), location r [m] and time t [s].
"""
function meanMagneticMoment(magneticMoment::MagneticMoments{T},H::SVector{3,T},r::SVector{3,T},t::T) where {T}
    error("meanMagneticMoment is not implemented for magneticMoment::$(typeof(magneticMoment)).")
end

# Define specific subtypes
## Langevin model
struct Langevin{T} <: MagneticMoments{T}
    msat::T
    beta::T
end

"""
$(SIGNATURES)

Initialize `Langevin` type `MagneticMoments` with diameter `d` [m], temperature
`T` [K], relative saturation magnetization [T] `Mrel` = Msat*μ₀ (μ₀ being
the vacuum permeability) and numeric `elementType`.
"""
function Langevin(;elementType=Float64,d=25e-9,T=294,Mrel=0.6)
    k = 1.380650424e-23 #Boltzman constant
    μ₀ = 4*π*1e-7 #vacuum permeability
    M = Mrel / μ₀ #saturation magnetization of magnetic material the core is made of
    msat = M * π/6*d^3 #saturation magnetic moment of a single nanoparticle
    beta = msat /(k*T) #H measured in T/μ₀

    Langevin{elementType}(msat,beta)
end

function meanMagneticMoment(magneticMoment::Langevin{T},H::SVector{3,T},r::SVector{3,T},t::T) where {T}
    if norm(H)!=0
        x = magneticMoment.beta*norm(H)
        return magneticMoment.msat*(coth(x) - 1/x)*normalize(H)
    else
        return zero(SVector{3,T})
    end
end