text:
#############################################################################
# operatornorm.jl
# Handles matrix operator norm (the maximum singular value of a matrix)
# and creates the alias sigmamax
# All expressions and atoms are subtypes of AbstractExpr.
# Please read expressions.jl first.
#############################################################################

### Operator norm

struct OperatorNormAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr}
    size::Tuple{Int,Int}

    function OperatorNormAtom(x::AbstractExpr)
        children = (x,)
        return new(:opnorm, hash(children), children, (1, 1))
    end
end

function sign(x::OperatorNormAtom)
    return Positive()
end

# The monotonicity
function monotonicity(x::OperatorNormAtom)
    return (NoMonotonicity(),)
end

function curvature(x::OperatorNormAtom)
    return ConvexVexity()
end

# in julia, `norm` on matrices is the operator norm
function evaluate(x::OperatorNormAtom)
    return opnorm(evaluate(x.children[1]), 2)
end

sigmamax(x::AbstractExpr) = OperatorNormAtom(x)

function LinearAlgebra.opnorm(x::AbstractExpr, p::Real = 2)
    if length(size(x)) <= 1 || minimum(size(x)) == 1
        throw(ArgumentError("argument to `opnorm` must be a matrix"))
    end
    if p == 1
        return maximum(sum(abs(x), dims = 1))
    elseif p == 2
        return OperatorNormAtom(x)
    elseif p == Inf
        return maximum(sum(abs(x), dims = 2))
    else
        throw(
            ArgumentError("matrix p-norms only defined for p = 1, 2, and Inf"),
        )
    end
end

# Create the equivalent conic problem:
#   minimize t
#   subject to
#            [tI_m A; A' tI_n] is positive semidefinite
# see eg Recht, Fazel, Parillo 2008 "Guaranteed Minimum-Rank Solutions of Linear Matrix Equations via Nuclear Norm Minimization"
# http://arxiv.org/pdf/0706.4138v1.pdf
function conic_form!(x::OperatorNormAtom, unique_conic_forms)
    if !has_conic_form(unique_conic_forms, x)
        A = x.children[1]
        m, n = size(A)
        t = Variable()
        p = minimize(t, [t*sparse(1.0I, m, m) A; A' t*sparse(1.0I, n, n)] ⪰ 0)
        cache_conic_form!(unique_conic_forms, x, p)
    end
    return get_conic_form(unique_conic_forms, x)
end


text:
const mpi_init_hooks = Any[]

macro mpi_handle(def, mpiname=nothing, extrafields...)
    name = def isa Expr ? def.args[1] : def
    _name = Symbol(:_, name)
    if mpiname == nothing
        mpiname = Symbol(:MPI_, name)
    end
    mpiname_f2c = Symbol(mpiname,:_f2c)
    quote
        Base.@__doc__ mutable struct $(esc(def))
            val::$mpiname
            $(extrafields...)

            $(esc(name))(val::$mpiname, $(extrafields...)) = new(val, $(extrafields...))
        end

        # const initializer
        function $(esc(_name))(ival::Cint, extraargs...)
            if $mpiname == Cint
                return $name(ival, extraargs...)
            else
                x = $name(C_NULL, extraargs...)
                push!(mpi_init_hooks, () -> x.val = $mpiname_f2c(ival))
                return x
            end
        end

        function Base.:(==)(a::$name, b::$name)
            a.val == b.val
        end

        # allows us to pass XX objects directly into MPI_XX ccall signatures
        function Base.cconvert(::Type{$mpiname}, obj::$name)
            obj.val
        end
        # allows us to pass XX objects directly into Ptr{MPI_XX} ccall signatures
        function Base.unsafe_convert(::Type{Ptr{$mpiname}}, obj::$name)
            convert(Ptr{$mpiname}, pointer_from_objref(obj))
        end
    end
end

"""
    MPI.free(obj)

Free the MPI object handle `obj`. This is typically used as the finalizer, and so need not be called directly unless otherwise noted.
"""
function free
end

ext: jl


text:
using ROS
ROS.@include actionlib_tutorials: FibonacciAction, FibonacciFeedback, FibonacciResult, FibonacciGoal

action_server = nothing

function init()
    ROS.init("test_actions")
    nh = ROS.NodeHandle()

    global action_server = ROS.actionlib.SimpleActionServer("test_action_server", ROS.actionlib_tutorials_FibonacciAction)
    ROS.actionlib.start(action_server)
    ROS.actionlib.registerGoalCallback(action_server, execute_cb)

    action_client = ROS.actionlib.SimpleActionClient("test_action_server", ROS.actionlib_tutorials_FibonacciAction)
    ROS.actionlib.sendGoal(action_client, ROS.actionlib_tutorials_FibonacciGoal())
    ROS.actionlib.waitForResult(action_client, ROS.Duration(3))

    action_result_called = ROS.actionlib.getResult(action_client)

    while ROS.ok()
        state = ROS.actionlib.getState(action_client)
        println(state)
        ROS.spinOnce()
    end

function execute_cb()
    goal = ROS.actionlib.acceptNewGoal(action_server)
    action_feedback = ROS.actionlib_tutorials_FibonacciFeedback()
    f::Vector{Number} = []
    for i in 1:100
        push!(f,i)
        action_feedback.sequence = f
        ROS.actionlib.publishFeedback(action_server, action_feedback)
    end
    action_result = ROS.actionlib_tutorials_FibonacciResult()
    action_result.sequence = f
    ROS.actionlib.setSucceeded(action_server, action_result, ":)")
end

init()


text:
#jl #! format: off
# # Operations problems with [PowerSimulations.jl](https://github.com/NREL-SIIP/PowerSimulations.jl)

# **Originally Contributed by**: Clayton Barrows

# ## Introduction

# PowerSimulations.jl supports the construction and solution of optimal power system
# scheduling problems (Operations Problems). Operations problems form the fundamental
# building blocks for [sequential simulations](https://nbviewer.jupyter.org/github/NREL-SIIP/SIIPExamples.jl/blob/master/notebook/3_PowerSimulations_examples/02_sequential_simulations.ipynb).
# This example shows how to specify and customize a the mathematics that will be applied to the data with
# an `ProblemTemplate`, build and execute an `DecisionModel`, and access the results.

# ## Dependencies
using SIIPExamples

# ### Modeling Packages
using PowerSystems
using PowerSimulations
const PSI = PowerSimulations
using PowerSystemCaseBuilder

# ### Optimization packages
using HiGHS # solver

# ### Data
# This data depends upon the [RTS-GMLC](https://github.com/gridmod/rts-gmlc) dataset. Let's
# use [PowerSystemCaseBuilder.jl](https://nbviewer.jupyter.org/github/NREL-SIIP/SIIPExamples.jl/blob/master/notebook/3_PowerSimulations_examples/10_PowerSystemCaseBuilder.ipynb) to download and build a `System`.

sys = build_system(PSITestSystems, "modified_RTS_GMLC_DA_sys")

# ## Define a problem specification with an `ProblemTemplate`
# You can create an empty template with:
template_uc = ProblemTemplate()

# Now, you can add a `DeviceModel` for each device type to create an assignment between PowerSystems device types
# and the subtypes of `AbstractDeviceFormulation`. PowerSimulations has a variety of different
# `AbstractDeviceFormulation` subtypes that can be applied to different PowerSystems device types,
# each dispatching to different methods for populating optimization problem objectives, variables,
# and constraints.

print_tree(PSI.AbstractDeviceFormulation)

# ### Branch Formulations
# Here is an example of relatively standard branch formulations. Other formulations allow
# for selective enforcement of transmission limits and greater control on transformer settings.
set_device_model!(template_uc, Line, StaticBranch)
set_device_model!(template_uc, Transformer2W, StaticBranch)
set_device_model!(template_uc, TapTransformer, StaticBranch)

# ### Injection Device Formulations
# Here we define template entries for all devices that inject or withdraw power on the
# network. For each device type, we can define a distinct `AbstractDeviceFormulation`. In
# this case, we're defining a basic unit commitment model for thermal generators,
# curtailable renewable generators, and fixed dispatch (net-load reduction) formulations
# for `HydroDispatch` and `RenewableFix` devices.

set_device_model!(template_uc, ThermalStandard, ThermalStandardUnitCommitment)
set_device_model!(template_uc, RenewableDispatch, RenewableFullDispatch)
set_device_model!(template_uc, PowerLoad, StaticPowerLoad)
set_device_model!(template_uc, HydroDispatch, FixedOutput)
set_device_model!(template_uc, HydroEnergyReservoir, HydroDispatchRunOfRiver)
set_device_model!(template_uc, RenewableFix, FixedOutput)

# ### Service Formulations
# We have two `VariableReserve` types, parameterized by their direction. So, similar to
# creating `DeviceModel`s, we can create `ServiceModel`s. The primary difference being
# that `DeviceModel` objects define how constraints get created, while `ServiceModel` objects
# define how constraints get modified.
set_service_model!(template_uc, VariableReserve{ReserveUp}, RangeReserve)
set_service_model!(template_uc, VariableReserve{ReserveDown}, RangeReserve)

# ### Network Formulations
# Finally, we can define the transmission network specification that we'd like to model. For simplicity, we'll
# choose a copper plate formulation. But there are dozens of specifications available through
# an integration with [PowerModels.jl](https://github.com/lanl-ansi/powermodels.jl). *Note that
# many formulations will require appropriate data and may be computationally intractable*
set_network_model!(template_uc, NetworkModel(CopperPlatePowerModel))

# ## `DecisionModel`
# Now that we have a `System` and an `ProblemTemplate`, we can put the two together
# to create an `DecisionModel` that we solve.

# ### Optimizer
# It's most convenient to define an optimizer instance upfront and pass it into the
# `DecisionModel` constructor. For this example, we can use the free HiGHS solver with a
# relatively relaxed MIP gap (`ratioGap`) setting to improve speed.

solver = optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.5)

# ### Build an `DecisionModel`
# The construction of an `DecisionModel` essentially applies an `ProblemTemplate`
# to `System` data to create a JuMP model.

problem = DecisionModel(template_uc, sys; optimizer = solver, horizon = 24)

build!(problem, output_dir = mktempdir())

#nb # The principal component of the `DecisionModel` is the JuMP model. For small problems,
#nb # you can inspect it by simply printing it to the screen:
#nb # ```julia
#nb # PSI.get_jump_model(problem)
#nb # ```
#nb #
#nb # For anything of reasonable size, that will be unmanageable. But you can print to a file:
#nb # ```julia
#nb # f = open("testmodel.txt","w"); print(f,PSI.get_jump_model(problem)); close(f)
#nb # ```
#nb #
# In addition to the JuMP model, an `DecisionModel` keeps track of a bunch of metadata
# about the problem and some references to pretty names for constraints and variables.
# All of these details are contained within the `optimization_container` field.

print_struct(typeof(PSI.get_optimization_container(problem)))

# ### Solve an `DecisionModel`

solve!(problem)

# ## Results Inspection
# PowerSimulations collects the `DecisionModel` results into a struct:

print_struct(PSI.ProblemResults)

res = ProblemResults(problem);

# ### Optimizer Stats
# The optimizer summary is included

get_optimizer_stats(res)

# ### Objective Function Value

get_objective_value(res)

# ### Variable, Parameter, Auxillary Variable, Dual, and Expression Values
# The solution value data frames for variables, parameters, auxillary variables, duals and
# expressions can be accessed using the `read_` methods:

read_variables(res)

# Or, you can read a single parameter values for parameters that exist in the results.
list_parameter_names(res)
read_parameter(res, "ActivePowerTimeSeriesParameter__RenewableDispatch")

# ## Plotting
# Take a look at the examples in [the plotting folder.](../../notebook/3_PowerSimulations_examples/Plotting)

max_stars_repo_name: NREL-SIIP/Examples

max_issues_repo_name: NREL-SIIP/Examples

max_forks_repo_name: NREL-SIIP/Examples


text:
export Cubic, Spline, npatches, funk

# see http://mathworld.wolfram.com/CubicSpline.html

import Base: getindex, show

"""
`Cubic(a,b,c,d)` is a 3rd degree polynomial `a+bx+cx^2+dx^3`.
Use `f(x)` to evalue `f` at the value `x`.
"""
struct Cubic
    a::Number
    b::Number
    c::Number
    d::Number
end

(f::Cubic)(x::Number) = f.a + x*(f.b + x*(f.c+x*f.d))

"""
`f'` where `f` is a `Cubic` or `Spline` is the derivative
of `f`.
"""
Base.adjoint(f::Cubic) = Cubic(f.b,2f.c,3f.d,0)

struct Spline
    patches::Array{Cubic,1}
    closed::Bool
end

is_closed(S::Spline) = S.closed

function Base.adjoint(S::Spline)
    plist = adjoint.(S.patches)
    return Spline(plist,S.closed)
end

"""
`Spline(vals,kind)` returns a cubic spline based on the values in `vals`.
The resulting spline `S` will have the property that `S(1)==y[1]`,
`S(2)==y[2]`, and so on up to `S(n)==y[n]` where `n` is the length of `y`.

+ If `kind` is `:open` then the second derivatives at the end points will be zero. (This is the default.)
+ If `kind` is `:closed` then we assume that we are interpolating a periodic function where `S(n+1)==S(1)`.
"""
function Spline(y::Array{T,1},kind::Symbol=:open)::Spline where T<:Number
    n = length(y)
    if kind == :open
        @assert n>2 "Open splines must have at least three points"
        return open_spline(y)
    end
    if kind == :closed
        @assert n>3 "Closed splines must have at least four points"
        return closed_spline(y)
    end

    error("Spline type must be :open or :closed, not $kind")
end

function npatches(S::Spline)
    return length(S.patches)
end

function show(io::IO, S::Spline)
    adjective = is_closed(S) ? "Closed" : "Open"
    print(io, "$adjective spline with $(npatches(S)) patches")
end

getindex(S::Spline, idx::Int) = S.patches[idx]

"""
`funk(S)` converts the Spline `S` into a callable function
(e.g., that can be passed to `plot`).
"""
function funk(S)::Function
    return x -> S(x)
end

function (S::Spline)(x::Real)
    np = npatches(S)

    if is_closed(S)

        p = Int(floor(x))
        x = mod(x,np)
        p = Int(floor(x))

        if p==0
            p = np
            x += np
        end
        f = S[p]
        return f(x-p)
    end

    # open spline
    p = Int(floor(x))
    if x < 1
        p = 1
    end

    if p > np
        p = np
    end
    f = S[p]
    return f(x-p)
end

function open_spline(y::Array{T,1})::Spline where T<:Number
    n = length(y)
    M = zeros(n,n)
    for i=1:n
        M[i,i] = 4
    end
    M[1,1] = 2
    M[n,n] = 2

    for i=1:n-1
        M[i,i+1] = 1
        M[i+1,i] = 1
    end

    rhs = zeros(T,n)
    rhs[1] = 3*(y[2]-y[1])
    for k=2:n-1
        rhs[k] = 3*(y[k+1]-y[k-1])
    end
    rhs[n] = 3*(y[n]-y[n-1])

    D = M\rhs

    a = zeros(Number,n-1)
    b = zeros(Number,n-1)
    c = zeros(Number,n-1)
    d = zeros(Number,n-1)

    for j=1:n-1
        a[j] = y[j]
        b[j] = D[j]
        c[j] = 3*(y[j+1]-y[j])-2D[j]-D[j+1]
        d[j] = 2*(y[j]-y[j+1])+D[j]+D[j+1]
    end
    return Spline( [ Cubic(a[i],b[i],c[i],d[i]) for i=1:n-1 ] , false)
end

function closed_spline(y::Array{T,1})::Spline where T <: Number
    yy = copy(y)
    n = length(y)
    prepend!(yy,y[end-3:end])
    append!(yy,y[1:4])

    S = open_spline(yy)

    plist = S.patches[5:end-3]
    return Spline(plist,true)

end


text:
# ---
# title: 462. Minimum Moves to Equal Array Elements II
# id: problem462
# author: Indigo
# date: 2021-02-04
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/minimum-moves-to-equal-array-elements-ii/description/>
# hidden: true
# ---
#
# Given a **non-empty** integer array, find the minimum number of moves required
# to make all array elements equal, where a move is incrementing a selected
# element by 1 or decrementing a selected element by 1.
#
# You may assume the array's length is at most 10,000.
#
# **Example:**
#
#
#
#     Input:
#     [1,2,3]
#
#     Output:
#     2
#
#     Explanation:
#     Only two moves are needed (remember each move increments or decrements one element):
#
#     [1,2,3]  =>  [2,2,3]  =>  [2,2,2]
#
#
#
## @lc code=start
using LeetCode

function min_moves2(nums::Vector{Int})
    median = partialsort!(nums, (length(nums) + 1) ÷ 2)
    return sum(abs, num - median for num in nums)
end
## @lc code=end


text:
df = DataFrame(CSV.File("../datasets/Ailerons/ailerons.data", header=0))
data = convert(Matrix, df[:,1:40])
target = convert(Array, df[:, 41])

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using ForwardDiff, DiffBase
using BenchmarkTools

include(joinpath(dirname(dirname(@__FILE__)), "test", "utils.jl"))

name(f) = last(split(string(f), '.'))

const SUITE = BenchmarkGroup()

const vecs = map(n -> rand(MersenneTwister(1), n), (10, 100, 1000))
const mats = map(n -> rand(MersenneTwister(1), n, n), (5, 16, 32))

const value_group = addgroup!(SUITE, "value")
const derivative_group = addgroup!(SUITE, "derivative")
const gradient_group = addgroup!(SUITE, "gradient")
const jacobian_group = addgroup!(SUITE, "jacobian")
const hessian_group = addgroup!(SUITE, "hessian")

for f in (DiffBase.NUMBER_TO_NUMBER_FUNCS..., DiffBase.NUMBER_TO_ARRAY_FUNCS...)
    x = 1.0
    y = f(x)

    value_group[name(f)] = @benchmarkable $(f)($x)

    out = isa(y, Number) ? DiffBase.DiffResult(y, y) : DiffBase.DiffResult(similar(y), similar(y))
    derivative_group[name(f)] = @benchmarkable ForwardDiff.derivative!($out, $f, $x)
end

for f in (DiffBase.VECTOR_TO_NUMBER_FUNCS..., DiffBase.MATRIX_TO_NUMBER_FUNCS...)
    fval = addgroup!(value_group, name(f))
    fgrad = addgroup!(gradient_group, name(f))
    fhess = addgroup!(hessian_group, name(f))
    arrs = in(f, DiffBase.VECTOR_TO_NUMBER_FUNCS) ? vecs : mats
    for x in arrs
        y = f(x)

        fval[length(x)] = @benchmarkable $(f)($x)

        gout = DiffBase.DiffResult(y, similar(x, typeof(y)))
        gcfg = ForwardDiff.GradientConfig(nothing, x)
        fgrad[length(x)] = @benchmarkable ForwardDiff.gradient!($gout, $f, $x, $gcfg)

        hout = DiffBase.DiffResult(y, similar(x, typeof(y)), similar(x, typeof(y), length(x), length(x)))
        hcfg = ForwardDiff.HessianConfig(nothing, hout, x)
        fhess[length(x)] = @benchmarkable ForwardDiff.hessian!($hout, $f, $x, $hcfg)
    end
end

for f in DiffBase.ARRAY_TO_ARRAY_FUNCS
    fval = addgroup!(value_group, name(f))
    fjac = addgroup!(jacobian_group, name(f))
    for x in mats
        y = f(x)
        fval[length(x)] = @benchmarkable $(f)($x)

        out = DiffBase.JacobianResult(y, x)
        cfg = ForwardDiff.JacobianConfig(nothing, y, x)
        fjac[length(x)] = @benchmarkable ForwardDiff.jacobian!($out, $f, $y, $x, $cfg)
    end
end


text:
export add_type_E!, remove_type_E!

"""If left kink and right_kink are type-E-entangled it returns a tuple of the two kinks whos orbs
are sorted in a way that k an j of both kinks are the common orbitals.
Otherwise it returns false.
This does not check wether the two kinks are neighbouring"""
function is_type_E(left_kink::T4, right_kink::T4)
    c_orb1 = intersect(creators(left_kink), annihilators(right_kink))
    c_orb2 = intersect(annihilators(left_kink), creators(right_kink))
    if (length(c_orb1) == 1) & (length(c_orb2) == 1)
        c_orb1 = first(c_orb1)
        c_orb2 = first(c_orb2)
        noncommon_orb_leftk_left = first(setdiff(annihilators(left_kink), creators(right_kink)))
        noncommon_orb_rightk_right = first(setdiff(creators(right_kink), annihilators(left_kink)))
        noncommon_orb_leftk_right = first(setdiff(creators(left_kink), annihilators(right_kink)))
        noncommon_orb_rightk_left = first(setdiff(annihilators(right_kink), creators(left_kink)))
        @assert length(unique( (noncommon_orb_leftk_left, noncommon_orb_rightk_right, noncommon_orb_leftk_right, noncommon_orb_rightk_left ) ) ) == 4 " The given orbitals are not pairwise distinct. "
        @assert noncommon_orb_leftk_right.vec - noncommon_orb_leftk_left.vec == noncommon_orb_rightk_left.vec - noncommon_orb_rightk_right.vec
        return T4(noncommon_orb_leftk_right, c_orb1, c_orb2, noncommon_orb_leftk_left), T4(noncommon_orb_rightk_right, c_orb2, c_orb1, noncommon_orb_rightk_left)
    else
        return false
    end
end

"""Return a Set of tuples of pairs of two imaginary times and kinks of 'neighbouring' kinks that are Type-E entangled AND removable.
'neighbouring' refers to that only Tuples of Kinks that are the closest Kink to act on an orbital of the
other kink in the corresponding direktion are looked at.
The Tuples are always arranged in a way that the Kink who gets neighboured by
the opther stands first.(vice versa does not have to be the case)
The Set consists of the pairs where the Type-E-entanglement is oriented
#to the left of the first τ."""
function left_type_E_removable_pairs(ck)
    pairs_left = Set{Tuple{Pair{ImgTime,T4},Pair{ImgTime,T4}}}()
    for (i,(τ,kink)) in enumerate(ck)
        kink_orbs = orbs(kink)
        i_left = prev_affecting(ck, kink_orbs, τ)
        sorted_kinks = is_type_E(last(ck[i_left]), kink)
        if sorted_kinks == false
            continue
        end
        i_right = next_affecting(ck, kink_orbs, τ)# TODO: this is only called to do the assertions
        @assert i_left != i_right
        @assert i_right != 0
        if norm(last(sorted_kinks).i.vec - last(sorted_kinks).k.vec) <= ex_radius
            push!(pairs_left, (τ => last(sorted_kinks), first(ck[i_left]) => first(sorted_kinks)))# TODO: how do we now that the first kink can be removed?
        end
    end
    return pairs_left
end

"""Return a Tuple of 2 imaginary times of 'neighbouring' Kinks that are Type-E entangled AND removable.
'neighbouring' refers to that only Tuples of Kinks that are the closest Kink to act on an orbital of the
other kink in the corresponding direktion are looked at.
The Tuples are always arranged in a way that the Kink who gets neighboured by
the opther stands first.(vice versa does not have to be the case)
The Set consists of the pairs where the Type-E-entanglement is oriented
#to the right of the first τ."""
function right_type_E_removable_pairs(ck)
    pairs_right = Set{Tuple{Pair{ImgTime,T4},Pair{ImgTime,T4}}}()
    for (i,(τ,kink)) in enumerate(ck)
        kink_orbs = orbs(kink)
        i_right = next_affecting(ck, kink_orbs, τ)

        sorted_kinks = is_type_E(kink, last(ck[i_right]))
        if sorted_kinks == false
            continue
        end
        i_left = prev_affecting(ck, kink_orbs, τ)# TODO: this is only called to do the assertions
        @assert i_left != i_right
        @assert i_right != 0
        if norm(first(sorted_kinks).i.vec - first(sorted_kinks).k.vec) <= ex_radius
            push!(pairs_right, (τ => first(sorted_kinks), first(ck[i_right]) => last(sorted_kinks)))
        end
    end
    return pairs_right
end

"""
    possible_new_kink_new_occ_orb(occs, new_kink_old_creator, changed_kink_old_creator, new_kink_old_annihilator, changed_kink_old_annihilator)
This function will find possibilites for the choice of the new Orbital that is occupied in occs and should be part of both Kinks after an add_type_E!-update,
when which orbital of the old-Kink is new_kink_old_creator, changed_kink_old_creator, new_kink_old_annihilator and changed_kink_old_annihilator
is already selected.
To check possibility this will in particular check conservation laws and the occupation of the resulting fourth orb.
"""
function possible_new_kink_new_occ_orb(occs, new_kink_old_creator, changed_kink_old_creator, new_kink_old_annihilator, changed_kink_old_annihilator)
    possibilities = filter(new_kink_new_annihilator -> (!in(find_fourth_orb_for_kink(new_kink_old_creator, new_kink_old_annihilator, new_kink_new_annihilator),occs)
                                                        && (!in(find_fourth_orb_for_kink(new_kink_old_creator, new_kink_old_annihilator, new_kink_new_annihilator), (new_kink_old_annihilator, changed_kink_old_annihilator)))),
                    intersect!(union!(sphere_with_same_spin(new_kink_old_creator, dk = ex_radius),
                            sphere_with_same_spin(PlaneWave(new_kink_old_creator.vec, flip(new_kink_old_annihilator.spin)), dk = ex_radius)), occs))
    return setdiff!(possibilities, (new_kink_old_creator, changed_kink_old_creator, new_kink_old_annihilator, changed_kink_old_annihilator) )
end

"""
    possible_new_kink_new_occ_orb(occs, new_kink_old_creator, changed_kink_old_creator, new_kink_old_annihilator, changed_kink_old_annihilator)
This function will find possibilites for the choice of the new Orbital that is unoccupied in occs and should be part of both Kinks after an add_type_E!-update,
when which orbital of the old-Kink is new_kink_old_creator, changed_kink_old_creator, new_kink_old_annihilator and changed_kink_old_annihilator
is already selected.
To check possibility this will in particular check conservation laws and the occupation of the resulting fourth orb.
"""
function possible_new_kink_new_unocc_orb(occs, new_kink_old_creator, changed_kink_old_creator, new_kink_old_annihilator, changed_kink_old_annihilator)
    possibilities = filter(new_kink_new_annihilator -> (in(find_fourth_orb_for_kink(new_kink_old_creator, new_kink_old_annihilator, new_kink_new_annihilator),occs)
                                                        && (!in(find_fourth_orb_for_kink(new_kink_old_creator, new_kink_old_annihilator, new_kink_new_annihilator), (new_kink_old_creator, changed_kink_old_creator)))),
                    setdiff!(union!(sphere_with_same_spin(new_kink_old_creator, dk = ex_radius),
                            sphere_with_same_spin(PlaneWave(new_kink_old_creator.vec, flip(new_kink_old_annihilator.spin)), dk = ex_radius)), occs))
    return setdiff!(possibilities, (new_kink_old_creator, changed_kink_old_creator, new_kink_old_annihilator, changed_kink_old_annihilator) )
end

function add_type_E!(m::Model, e::Ensemble, c::Configuration) :: Tuple{Float64,Step}
    #After the updates the i and l komponents of both kinks will contain the old kink, while the j and k components contain the old orbitals
    prop_prob = 1.0
    if isempty(c.kinks)
        return 1.0, Step()
    end
    old_kink_index = rand(1:length(c.kinks))# TODO: use this in order to pass the index of the removed kink to Step()
    old_kink = c.kinks[ old_kink_index ]
    prop_prob *= 1.0/length(c.kinks)
    occs = occupations_at(c, first(old_kink))
    prop_prob *= 0.5 # left or right
    direction = rand([:left, :right])
    #now choose orbitals of the old kinks that shell also be part of the new kink
    prop_prob *= 0.5
    #choose exitation to of old kink to happen in first kink after update and exitation to happen second
    if rand() > 0.5
        new_kink_old_creator = last(old_kink).i
        changed_kink_old_creator = last(old_kink).j
    else
        new_kink_old_creator = last(old_kink).j
        changed_kink_old_creator = last(old_kink).i
    end
    prop_prob *= 0.5
    if rand() > 0.5
        new_kink_old_annihilator = last(old_kink).k
        changed_kink_old_annihilator = last(old_kink).l
    else
        new_kink_old_annihilator = last(old_kink).l
        changed_kink_old_annihilator = last(old_kink).k
    end

    #find occupied orb for creation of Type_E
    if direction == :left
        opportunities_new_kink_new_annihilator = possible_new_kink_new_occ_orb(occs, new_kink_old_creator, changed_kink_old_creator, new_kink_old_annihilator, changed_kink_old_annihilator)
    else
        opportunities_new_kink_new_annihilator = possible_new_kink_new_unocc_orb(occs, new_kink_old_creator, changed_kink_old_creator, new_kink_old_annihilator, changed_kink_old_annihilator)
    end
    if isempty(opportunities_new_kink_new_annihilator)
        return 1.0, Step()
    end
    new_kink_new_annihilator = rand(opportunities_new_kink_new_annihilator)
    prop_prob *= 1.0/length(opportunities_new_kink_new_annihilator)
    #calculate new creator

    new_kink_new_creator = find_fourth_orb_for_kink(new_kink_old_creator, new_kink_old_annihilator, new_kink_new_annihilator)

    @assert(!in(new_kink_new_creator, orbs(last(old_kink))))

    if direction == :left
        @assert(!in(new_kink_new_creator, occs))
        τ_Intervall = first(old_kink) - τ_prev_affecting( c.kinks, (new_kink_old_creator, new_kink_new_creator, new_kink_new_annihilator, new_kink_old_annihilator), first(old_kink) )
    else
        @assert(in(new_kink_new_creator, occs))
        τ_Intervall = τ_next_affecting(c.kinks, (new_kink_old_creator, new_kink_new_creator, new_kink_new_annihilator, new_kink_old_annihilator), first(old_kink)) - first(old_kink)
    end

    if τ_Intervall < 0
        τ_Intervall +=1
    end
    @assert τ_Intervall > 0
    delta_τ = ImgTime(rand()*τ_Intervall)
    if direction == :left
        τ_new_kink = first(old_kink) - delta_τ
    else
        τ_new_kink = first(old_kink) + delta_τ
    end
    if τ_new_kink > 1
        τ_new_kink -= 1
        change_orbs = true
    elseif τ_new_kink < 0
        τ_new_kink += 1
        change_orbs = true
    else
        change_orbs = false
    end
    #no 2 kinks at same τ
    while haskey(c.kinks, τ_new_kink)
        delta_τ = ImgTime(rand()*τ_Intervall)
        if direction == :left
            τ_new_kink = first(old_kink) - delta_τ
        else
            τ_new_kink = first(old_kink) + delta_τ
        end
        if τ_new_kink > 1
            τ_new_kink -= 1
            change_orbs = true
        elseif τ_new_kink < 0
            τ_new_kink += 1
            change_orbs = true
        else
            change_orbs = false
        end
    end

    prop_prob *= 1.0/float(τ_Intervall)

    # see if c.occupations change
    if change_orbs
        if direction == :left
            orbs_drop = (new_kink_new_annihilator, new_kink_old_annihilator)
            orbs_add = (new_kink_old_creator,new_kink_new_creator)
        else
            orbs_drop = (new_kink_old_creator, new_kink_new_creator)
            orbs_add = (new_kink_new_annihilator, new_kink_old_annihilator)
        end
    else
        orbs_drop = nothing
        orbs_add = nothing
    end

    # shuffle indices
    prop_prob *= 1.0/16.0
    # shuffle changed kink: shuffle creators and shuffle annihilators
    add_kink1 = first(old_kink) => T4( random_shuffle( new_kink_new_annihilator, changed_kink_old_creator )..., random_shuffle( changed_kink_old_annihilator, new_kink_new_creator )... )
    # shuffle new kink: shuffle creators and shuffle annihilators
    add_kink2 = τ_new_kink => T4( random_shuffle( new_kink_new_creator, new_kink_old_creator )..., random_shuffle( new_kink_old_annihilator, new_kink_new_annihilator )... )

    # MC Step generated by this update
    Δ = Step(orbs_drop, (old_kink,), orbs_add, (add_kink1, add_kink2))

    if direction == :left
        delta_di = ΔWdiag_element(m, e, c, new_kink_old_creator, new_kink_new_creator, new_kink_new_annihilator, new_kink_old_annihilator, τ_new_kink, first(old_kink))
    else
        delta_di = ΔWdiag_element(m, e, c, new_kink_new_annihilator, new_kink_old_annihilator, new_kink_old_creator, new_kink_new_creator, first(old_kink), τ_new_kink)
    end

    apply_step!(c,Δ)
    @assert (is_type_E(last(add_kink2), last(add_kink1)) != false)
    # calculate weight difference
    dw_off_diag = abs( Woffdiag_element(m, e, last(add_kink2)) * Woffdiag_element(m, e, last(add_kink1)) / Woffdiag_element(m, e, last(old_kink)) )
    if direction == :left
        dw = e.β * dw_off_diag* exp(-(e.β * delta_τ*(energy(m, last(add_kink2).i) + energy(m, last(add_kink2).j)- energy(m, last(add_kink2).k) - energy(m, last(add_kink2).l)) + e.β * delta_di))
        inverse_prop_prob = 0.125 / length( right_type_E_removable_pairs( c.kinks ) )# TODO: is there a faster way than to explicitly calculate all these pairs in order to get their number?
    else
        dw = e.β * dw_off_diag* exp(-(e.β * delta_τ*(energy(m, last(add_kink2).k) + energy(m, last(add_kink2).l) - energy(m, last(add_kink2).i) - energy(m, last(add_kink2).j)) + e.β * delta_di))
        inverse_prop_prob = 0.125 / length( left_type_E_removable_pairs( c.kinks ) )# TODO: is there a straighter way than to explicitly calculate all these pairs in order to get their number ?
    end
    @assert !isinf(inverse_prop_prob)
    @assert !iszero(prop_prob)
    @assert !iszero(inverse_prop_prob)
    @assert (  0 <= (inverse_prop_prob/prop_prob)*dw < Inf)
    return (inverse_prop_prob/prop_prob)*dw, Δ
end

function remove_type_E!(m::Model, e::Ensemble, c::Configuration) :: Tuple{Float64,Step}
    prop_prob = 0.5
    direction = rand([:left, :right])
    #removed kink left of changed kink
    if direction == :left
        opportunities = right_type_E_removable_pairs(c.kinks)
    else
        opportunities = left_type_E_removable_pairs(c.kinks)
    end
    if isempty(opportunities)
        return 1.0, Step()
    end
    removed_kink_pair, changed_kink_pair = rand(opportunities)
    removed_kink_τ =  first(removed_kink_pair)
    changed_kink_τ =  first(changed_kink_pair)
    prop_prob *= 1.0/length(opportunities)

    # safe those for later
    removed_kink = last(removed_kink_pair)
    changed_kink_old = last(changed_kink_pair)
    # shuffle Indices
    prop_prob *= 0.25
    # shuffle creators and shuffle annihilators
    add_kink = changed_kink_τ => T4( random_shuffle( changed_kink_old.i, removed_kink.i )..., random_shuffle( changed_kink_old.l, removed_kink.l )...)
    if direction == :left
        τ_left = removed_kink_τ
        τ_right = changed_kink_τ
    else
        τ_left = changed_kink_τ
        τ_right = removed_kink_τ
    end
    # see if c.occupations change
    if τ_left > τ_right
        if direction == :left
            orbs_drop = (removed_kink.i,removed_kink.j)
            orbs_add = (removed_kink.k,removed_kink.l)
        else
            orbs_drop = (removed_kink.k,removed_kink.l)
            orbs_add = (removed_kink.i, removed_kink.j)
        end
    else
        orbs_drop = nothing
        orbs_add = nothing
    end

    kinks_drop = (changed_kink_τ => c.kinks[changed_kink_τ], removed_kink_τ => c.kinks[removed_kink_τ])

    # MC Step generated by this update
    Δ = Step(orbs_drop, kinks_drop, orbs_add, (add_kink,))

    apply_step!(c,Δ) # get proposed (new) configuration for inverse proposal probability

    # calculate reverse_prop_prob
    occs = occupations_at(c, changed_kink_τ)
    if direction == :left
        opportunities_new_kink_new_annihilator = possible_new_kink_new_occ_orb(occs, removed_kink.i, changed_kink_old.i, removed_kink.l, changed_kink_old.l)
        @assert(in(removed_kink.k,opportunities_new_kink_new_annihilator))
        τ_Intervall = changed_kink_τ - τ_prev_affecting( c.kinks, (removed_kink.k, removed_kink.l,removed_kink.i, removed_kink.j), changed_kink_τ )
    else
        opportunities_new_kink_new_annihilator = possible_new_kink_new_unocc_orb(occs, removed_kink.i, changed_kink_old.i, removed_kink.l, changed_kink_old.l)
        @assert(in(removed_kink.k,opportunities_new_kink_new_annihilator))
        τ_Intervall =  τ_next_affecting( c.kinks, (removed_kink.i, removed_kink.j,removed_kink.k, removed_kink.l), changed_kink_τ ) - changed_kink_τ
    end

    if τ_Intervall < 0
        τ_Intervall +=1
    end

    # the number of kinks is reduced by one with this update
    # TODO: explain all factors
    inverse_prop_prob = 0.5 / (length(c.kinks)) / length(opportunities_new_kink_new_annihilator) / float(τ_Intervall) / 4 / 16

    # calculate weight change
    dw_off_diag = abs( Woffdiag_element(m, e, removed_kink) * Woffdiag_element(m, e, changed_kink_old) / Woffdiag_element(m, e, last(add_kink)) )

    delta_τ = float(τ_right - τ_left)
    if delta_τ < 0
        delta_τ +=1
    end
    if direction == :left
        delta_di = ΔWdiag_element(m, e, c, removed_kink.i, removed_kink.j, removed_kink.k, removed_kink.l, removed_kink_τ, changed_kink_τ)
        dw = (1.0/e.β)* (1.0/dw_off_diag) * exp(e.β * delta_τ * (energy(m,removed_kink.i) + energy(m,removed_kink.j) - energy(m,removed_kink.k) - energy(m,removed_kink.l)) + e.β * delta_di)
    else
        delta_di = ΔWdiag_element(m, e, c, removed_kink.k, removed_kink.l, removed_kink.i,removed_kink.j, changed_kink_τ, removed_kink_τ)
        dw = (1.0/e.β)*(1.0/dw_off_diag) * exp(e.β * delta_τ*(energy(m, removed_kink.k) + energy(m, removed_kink.l) - energy(m, removed_kink.i) - energy(m, removed_kink.j)) + e.β * delta_di)
    end
    @assert (  0 <= (inverse_prop_prob/prop_prob)*dw < Inf)
    return (inverse_prop_prob/prop_prob) * dw, Δ
end

function isuseful(c::Configuration, up::typeof(add_type_E!))
    if isempty(c.kinks)
        return false
    else
        return true
    end
end

function isuseful(c::Configuration, up::typeof(remove_type_E!))
    if length(c.kinks) < 3
        return false
    else
        return true
    end
end


text:
module NRES

# TODO TemporalType kan også includeres ved "using SAT" inne i NRES-module!
# -> Vi kan mao. gå ut fra activate!(TemporalType) , deactivate!, is_active, osv..

"""
    abstract type Representation end

Representation interface:
    - function map_to_SAT   ->   returns SAT that corresponds to coordinate.
"""
abstract type Representation end
# TODO: PLAN:
#   - treng mdl-variabel the_active_trait::T
#   - treng funktion activate_in_muex!(map::Representation, trait::T)
#       -> activates trait and deactivates current trait (if any). Save as prev. active (for learning)?
#   - function active_trait_in(map::Representation)

#function the_active_trait_in_muex(map::Representation)

"""
    abstract type TemporalType
Type for representing SAT and other elements that can partake in a temporal sequence.
    - function is_active
    - function activate!
    - function deactivate!
"""
abstract type TemporalType end

include("tile_map.jl")
include("composite_map.jl")

# export Node, SAT, MuExS
end # module

max_stars_repo_name: leikanger/NRES

max_issues_repo_name: leikanger/NRES

max_forks_repo_name: leikanger/NRES


text:
"""
    RobertsBL1D(delta::FT, h::FT) where FT

The β for a given boundary layer
as described in section 5.6 (page 333) of
Computational Fluid Mechanics and Heat Transfer,
2nd edition, J. Tannehill et al.

# Input
```
hmin     = wall boundary (minimum value)
hmax     = wall boundary (maximum value)
delta    = thickness of boundary layer
```

The `if` statement protects against the case when
β = infinity.
"""
function RobertsBL1D(delta::FT, h::FT) where FT
  if delta < h*FT(99/100)
    β = (FT(1) - delta/h)^(-FT(1/2))
  else
    β = FT(1000)
  end
  return β
end

RobertsBL(delta::FT, hmin::FT, hmax::FT) where FT = RobertsBL1D(delta, hmax-hmin)

function Hartmann_BL_1D(Ha::FT, hmin::FT, hmax::FT) where FT
  @assert Ha>eps(FT)
  return RobertsBL((hmax-hmin)/Ha, hmin, hmax)
end

function Reynolds_BL_1D(Re::FT, hmin::FT, hmax::FT) where FT
  @assert Re>eps(FT)
  β = RobertsBL((hmax-hmin)/sqrt(Re), hmin, hmax)
end

function HartmannBL(Ha::FT, hmin::Vector{FT}, hmax::Vector{FT}) where FT
  @assert Ha>eps(FT)
  @assert length(hmin)==length(hmax)
  return [RobertsBL((hmax[i]-hmin[i])/Ha, hmin[i],hmax[i]) for i in 1:3]
end

function Re_Ha_BL_1D(Re::FT, Ha::FT, hmin::FT, hmax::FT) where FT
 temp1 = Reynolds_BL_1D(Re, hmin, hmax)
 temp2 = Hartmann_BL_1D(Ha, hmin, hmax)
 return min(temp1, temp2)
end

function ReynoldsBL(Re::FT, hmin::Vector{FT}, hmax::Vector{FT}) where FT
  @assert Re>eps(FT)
  β = [RobertsBL((hmax[i] - hmin[i])/sqrt(Re), hmin[i], hmax[i]) for i in 1:3]
end

function Re_Ha_BL(Re::FT, Ha::FT, hmin::Vector{FT}, hmax::Vector{FT}) where FT
  @assert length(hmin)==length(hmax)
  temp1 = ReynoldsBL(Re, hmin, hmax)
  temp2 = HartmannBL(Ha, hmin, hmax)
  β = [min(temp1[i],temp2[i]) for i in 1:length(hmin)]
end


text:
export RandomVariable, randomvar, simplerandomvar

## Random variable implementation

mutable struct RandomVariable <: AbstractVariable
    name                :: Symbol
    inputmsgs           :: Vector{LazyObservable{AbstractMessage}}
    marginal            :: Union{Nothing, MarginalObservable}
    pipeline            :: AbstractPipelineStage
    local_constraint
    prod_constraint
    prod_strategy
    form_constraint
    form_check_strategy
end

function randomvar(name::Symbol; pipeline = EmptyPipelineStage(), local_constraint = Marginalisation(), prod_constraint = ProdAnalytical(), prod_strategy = FoldLeftProdStrategy(), form_constraint = UnspecifiedFormConstraint(), form_check_strategy = FormConstraintCheckPickDefault())
    return RandomVariable(name, Vector{LazyObservable{AbstractMessage}}(), nothing, pipeline, local_constraint, prod_constraint, prod_strategy, form_constraint, form_check_strategy)
end

function randomvar(name::Symbol, dims::Tuple; pipeline = EmptyPipelineStage(), local_constraint = Marginalisation(), prod_constraint = ProdAnalytical(), prod_strategy = FoldLeftProdStrategy(), form_constraint = UnspecifiedFormConstraint(), form_check_strategy = FormConstraintCheckPickDefault())
    return randomvar(name, dims...; pipeline = pipeline, local_constraint = local_constraint, prod_constraint = prod_constraint, prod_strategy = prod_strategy, form_constraint = form_constraint, form_check_strategy = form_check_strategy)
end

function randomvar(name::Symbol, dims::Vararg{Int}; pipeline = EmptyPipelineStage(), local_constraint = Marginalisation(), prod_constraint = ProdAnalytical(), prod_strategy = FoldLeftProdStrategy(), form_constraint = UnspecifiedFormConstraint(), form_check_strategy = FormConstraintCheckPickDefault())
    vars = Array{RandomVariable}(undef, dims)
    for i in CartesianIndices(axes(vars))
        @inbounds vars[i] = randomvar(Symbol(name, :_, Symbol(join(i.I, :_))); pipeline = pipeline, local_constraint = local_constraint, prod_constraint = prod_constraint, prod_strategy = prod_strategy, form_constraint = form_constraint, form_check_strategy = form_check_strategy)
    end
    return vars
end

degree(randomvar::RandomVariable)              = length(randomvar.inputmsgs)
name(randomvar::RandomVariable)                = randomvar.name
local_constraint(randomvar::RandomVariable)    = randomvar.local_constraint
prod_constraint(randomvar::RandomVariable)     = randomvar.prod_constraint
prod_strategy(randomvar::RandomVariable)       = randomvar.prod_strategy
form_constraint(randomvar::RandomVariable)     = randomvar.form_constraint
form_check_strategy(randomvar::RandomVariable) = _form_check_strategy(randomvar.form_check_strategy, randomvar)

_form_check_strategy(::FormConstraintCheckPickDefault, randomvar::RandomVariable) = default_form_check_strategy(form_constraint(randomvar))
_form_check_strategy(form_check_strategy, randomvar::RandomVariable)              = form_check_strategy

messages_prod_fn(randomvar::RandomVariable) = messages_prod_fn(prod_strategy(randomvar), prod_constraint(randomvar), UnspecifiedFormConstraint(), default_form_check_strategy(UnspecifiedFormConstraint()))
marginal_prod_fn(randomvar::RandomVariable) = marginal_prod_fn(prod_strategy(randomvar), prod_constraint(randomvar), form_constraint(randomvar), form_check_strategy(randomvar))

getlastindex(randomvar::RandomVariable) = degree(randomvar) + 1

messagein(randomvar::RandomVariable, index::Int)  = @inbounds randomvar.inputmsgs[index]
messageout(randomvar::RandomVariable, index::Int) = collectLatest(AbstractMessage, Message, skipindex(randomvar.inputmsgs, index), messages_prod_fn(randomvar))

get_pipeline_stages(randomvar::RandomVariable)        = randomvar.pipeline
add_pipeline_stage!(randomvar::RandomVariable, stage) = randomvar.pipeline = (randomvar.pipeline + stage)

_getmarginal(randomvar::RandomVariable)                                = randomvar.marginal
_setmarginal!(randomvar::RandomVariable, marginal::MarginalObservable) = randomvar.marginal = marginal
_makemarginal(randomvar::RandomVariable)                               = collectLatest(AbstractMessage, Marginal, randomvar.inputmsgs, marginal_prod_fn(randomvar))

function setmessagein!(randomvar::RandomVariable, index::Int, messagein)
    if index === degree(randomvar) + 1
        push!(randomvar.inputmsgs, messagein)
    else
        throw_setmessagein!(randomvar, index)
    end
end

# https://github.com/JuliaLang/julia/issues/29688
throw_setmessagein!(randomvar, index) = throw("Inconsistent state in setmessagein! function for random variable $(randomvar). `index` should be equal to `degree(randomvar) + 1 = $(degree(randomvar) + 1)`, $(index) is given instead")


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modeling language for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# deprecated.jl
# List of deprecated macros and functions, mostly from the "big renaming"
macro deprecate_macro(old,new)
    oldmac = Symbol(string("@",old))
    newmac = Symbol(string("@",new))
    s = string(oldmac," is deprecated, use ", newmac, " instead.")
    if VERSION > v"0.5-"
        # backtraces are ok on 0.5
        depwarn = :(Base.depwarn($s,$(quot(oldmac))))
    else
        # backtraces are junk on 0.4
        depwarn = :(Base.warn_once($s))
    end
    @eval macro $old(args...)
        return Expr(:block, $depwarn, Expr(:macrocall, $(quot(newmac)), [esc(x) for x in args]...))
    end
    eval(Expr(:export,oldmac))
    return
end

@deprecate_macro defVar variable
@deprecate_macro defVars variables
@deprecate_macro setObjective objective
@deprecate_macro setNLObjective NLobjective
@deprecate_macro addConstraint constraint
@deprecate_macro addConstraints constraints
@deprecate_macro addNLConstraint NLconstraint
@deprecate_macro addNLConstraints NLconstraints
@deprecate_macro addSDPConstraint SDconstraint
@deprecate_macro defExpr expression
@deprecate_macro defNLExpr NLexpression
@deprecate_macro defNLParam NLparameter
@deprecate_macro defConstrRef constraintref

# callback macros
@deprecate_macro addLazyConstraint lazyconstraint
@deprecate_macro addusercut usercut

@deprecate JuMPNLPEvaluator JuMP.NLPEvaluator
@deprecate addConstraint JuMP.addconstraint
@deprecate addCutCallback addcutcallback
@deprecate addHeuristicCallback addheuristiccallback
@deprecate addInfoCallback addinfocallback
@deprecate addLazyCallback addlazycallback
@deprecate addLazyConstraint JuMP.addlazyconstraint
@deprecate addSolution addsolution
@deprecate addToExpression JuMP.addtoexpr
@deprecate addUserCut JuMP.addusercut
@deprecate affToStr string
@deprecate buildInternalModel JuMP.build
@deprecate chgConstrRHS JuMP.setRHS
@deprecate conToStr string
@deprecate exprToStr string
@deprecate getCategory getcategory
@deprecate getConstraintBounds JuMP.constraintbounds
@deprecate getDual getdual
@deprecate getInternalModel internalmodel
@deprecate getLinearIndex linearindex
@deprecate getLower getlowerbound
@deprecate getName getname
@deprecate getObjective getobjective
@deprecate getObjectiveSense getobjectivesense
@deprecate getObjectiveValue getobjectivevalue
@deprecate getUpper getupperbound
@deprecate getValue getvalue
@deprecate getVar getvariable
@deprecate quadToStr string
@deprecate setCategory setcategory
@deprecate setLower setlowerbound
@deprecate setName setname
@deprecate setObjective JuMP.setobjective
@deprecate setObjectiveSense setobjectivesense
@deprecate setPrintHook JuMP.setprinthook
@deprecate setSolutionValue! setsolutionvalue
@deprecate setSolveHook JuMP.setsolvehook
@deprecate setSolver setsolver
@deprecate setUpper setupperbound
@deprecate setValue setvalue

function registerNLFunction(args...;autodiff=false)
    Base.depwarn("registerNLFunction is deprecated, use JuMP.register instead.",:registerNLFunction)
    JuMP.register(args...,autodiff=autodiff)
end
export registerNLFunction


text:
# We use nHidden as a vector, containing the number of hidden units in each layer

# Function that returns total number of parameters
function NeuralNet_nParams(d,nHidden)

	# Connections from inputs to first hidden layer
	nParams = d*nHidden[1]

	# Connections between hidden layers
	for h in 2:length(nHidden)
		nParams += nHidden[h-1]*nHidden[h]
	end

	# Connections from last hidden layer to output
	nParams += nHidden[end]

end

# Compute squared error and gradient
# for a single training example (x,y)
# (x is assumed to be a column-vector)
function NeuralNet_backprop(bigW,x,y,nHidden)

	if typeof(x)<:Array{Float64,1}
		d=length(x)
		n=1
	else
		n,d = size(x)
	end
	nLayers = length(nHidden)

	#### Reshape 'bigW' into vectors/matrices
	W1 = reshape(bigW[1:nHidden[1]*d],nHidden[1],d)
	ind = nHidden[1]*d
	Wm = Matrix{Float64}[]
	for layer in 2:nLayers
		push!(Wm, reshape(bigW[ind+1:ind+nHidden[layer]*nHidden[layer-1]],nHidden[layer],nHidden[layer-1]))
		ind += nHidden[layer]*nHidden[layer-1]
	end
	w = bigW[ind+1:end]

	#### Define activation function and its derivative
	h(z) = (tanh.(z))
	dh(z) = (sech.(z)).^2

	#### Forward propagation
	z = Matrix{Float64}[]

	push!(z, x*W1')

	for layer in 2:nLayers
		push!(z, h(z[layer-1]*Wm[layer-1]'))
	end

	yhat = (h(z[end])*w)

	r = yhat-y
	f = (1/2)r.^2

	#### Backpropagation
	dr = r
	err = dr

	# Output weights
	Gout = err'*h(z[end])

	Gm = Matrix{Float64}[]
	if nLayers > 1
		# Last Layer of Hidden Weights
		#print(size(err), size(dh(z[end])),size(w))
		backprop = (err.*(dh(z[end]).*w'))
		#print(size(backprop),size(h(z[end-1])))
		push!(Gm,backprop'*h(z[end-1]))

		# Other Hidden Layers
		for layer in nLayers-2:-1:1
			#print(size(Wm[layer+1]), size(backprop),size(dh(z[layer+1])))
			backprop = (backprop*Wm[layer+1]).*dh(z[layer+1])
			#print(size(backprop),size(h(z[layer])))
			push!(Gm,backprop'*h(z[layer]))
		end
		GM=reverse!(Gm)
		# Input Weights
		#print(size(Wm[1]), size(backprop),size(dh(z[1])))
		backprop = (backprop*Wm[1]).*dh(z[1])
		G1 = backprop'*x
		#print(size(G1))
	else
		# Input weights
		G1 = x'*(err.*(dh(z[1]).*w'))
	end

	#### Put gradients into vector
	g = zeros(size(bigW))
	g[1:nHidden[1]*d] = G1
	ind = nHidden[1]*d
	for layer in 2:nLayers
		g[ind+1:ind+nHidden[layer]*nHidden[layer-1]] = Gm[layer-1]
		ind += nHidden[layer]*nHidden[layer-1]
	end
	g[ind+1:end] = Gout

	return (f,g/n)
end

# Computes predictions for a set of examples X
function NeuralNet_predict(bigW,Xhat,nHidden)
	(t,d) = size(Xhat)
	nLayers = length(nHidden)

	#### Reshape 'bigW' into vectors/matrices
	W1 = reshape(bigW[1:nHidden[1]*d],nHidden[1],d)
	ind = nHidden[1]*d
	Wm = Matrix{Float64}[]
	for layer in 2:nLayers
		push!(Wm, reshape(bigW[ind+1:ind+nHidden[layer]*nHidden[layer-1]],nHidden[layer],nHidden[layer-1]))
		ind += nHidden[layer]*nHidden[layer-1]
	end
	w = bigW[ind+1:end]

	#### Define activation function and its derivative
	h(z) = tanh.(z)
	dh(z) = (sech.(z)).^2

	#### Forward propagation on each example to make predictions
	yhat = zeros(t,1)
	for i in 1:t
		# Forward propagation
		z = fill([],(nLayers))
		z[1] = W1*Xhat[i,:]
		for layer in 2:nLayers
			z[layer] = Wm[layer-1]*h(z[layer-1])
		end
		yhat[i] = w'*h(z[end])
	end
	return yhat
end


text:
using RNLA, Test, LinearAlgebra

@testset "rnmf_test1" begin
    m = 50
    k = 10
    n = 30
    W = rand(1:100, (m, k))
    H = rand(1:100, (k, n))
    X = W * H

    rW, rH = rnmf(X, k)
    @test norm(rW - W) <= sqrt(eps())
    @test norm(rH - H) <= sqrt(eps())
end

@testset "rnmf_test2" begin
    m = 100
    k = 5
    n = 15
    W = rand(1:100, (m, k))
    H = rand(1:100, (k, n))
    X = W*H

    rW, rH = rnmf(X, k)
    @test norm(rW - W) <= sqrt(eps())
    @test norm(rH - H) <= sqrt(eps())
end


text:
import Revise
import GLRenderer as GL
import ThreeDP3 as T
import Images as I
import MiniGSG as S
import PoseComposition: Pose
using Gen
import MeshCatViz

import MeshCatViz

MeshCatViz.setup_visualizer()

YCB_DIR = "/home/nishadg/mcs/ThreeDVision.jl/data/ycbv2"
world_scaling_factor = 100.0
id_to_cloud, id_to_shift, id_to_box  = T.load_ycbv_models_adjusted(YCB_DIR, world_scaling_factor);
all_ids = sort(collect(keys(id_to_cloud)));

IDX =100
@show T.get_ycb_scene_frame_id_from_idx(YCB_DIR,IDX)
gt_poses, ids, rgb_image, gt_depth_image, cam_pose, original_camera = T.load_ycbv_scene_adjusted(
    YCB_DIR, IDX, world_scaling_factor, id_to_shift
);
camera = T.scale_down_camera(original_camera, 4)
@show camera
img = I.colorview(I.Gray, gt_depth_image ./ maximum(gt_depth_image))
I.colorview(I.RGB, permutedims(Float64.(rgb_image)./255.0, (3,1,2)))

# +
resolution = 0.5

renderer = GL.setup_renderer(camera, GL.DepthMode())
for id in all_ids
    cloud = id_to_cloud[id]
    v,n,f = GL.mesh_from_voxelized_cloud(GL.voxelize(cloud, resolution), resolution)
    GL.load_object!(renderer, v, f)
end
depth_image = GL.gl_render(
    renderer, ids, gt_poses, T.IDENTITY_POSE)
depth_image[depth_image .== 50000.0] .= 200.0

function get_cloud_from_ids_and_poses(ids, poses, camera_pose)
    depth_image = GL.gl_render(renderer, ids, poses, camera_pose)
    cloud = GL.flatten_point_cloud(GL.depth_image_to_point_cloud(depth_image, camera))
    cloud = cloud[:,cloud[3,:] .< (camera.far - 1.0)]
    if size(cloud)[2] == 0
        cloud = zeros(3,1)
    end
    cloud
end

img = I.colorview(I.Gray, depth_image ./ maximum(depth_image))

# +

hypers = T.Hyperparams(slack_dir_conc=300.0, slack_offset_var=0.5, p_outlier=0.01,
    noise=0.2, resolution=0.5, parent_face_mixture_prob=0.99, floating_position_bounds=(-1000.0, 1000.0, -1000.0,1000.0,-1000.0,1000.0))

params = T.SceneModelParameters(
    boxes=vcat([id_to_box[id] for id in ids],[S.Box(40.0,40.0,0.1)]),
    get_cloud_from_poses_and_idx=
    (poses, idx, p) -> get_cloud_from_ids_and_poses(ids, poses[1:end-1], cam_pose),
    camera_pose=cam_pose,
    hyperparams=hypers, N=1)

table_pose = Pose(0.0, 0.0, -0.1)
num_obj = length(params.boxes)
g = T.graph_with_edges(num_obj, [])
constraints = Gen.choicemap((T.structure_addr()) => g, T.floating_pose_addr(num_obj) => table_pose)
for i=1:num_obj-1
   constraints[T.floating_pose_addr(i)] = T.get_c_relative_to_a(cam_pose, gt_poses[i])
end

obs_cloud = GL.flatten_point_cloud(GL.depth_image_to_point_cloud(gt_depth_image, original_camera))
obs_cloud = GL.move_points_to_frame_b(obs_cloud, cam_pose)
obs_cloud =obs_cloud[:, 1.2 .< obs_cloud[3, :] .< 40.0];
obs_cloud =obs_cloud[:, -25.0 .< obs_cloud[1, :] .< 25.0];
obs_cloud =obs_cloud[:, -25.0 .< obs_cloud[2, :] .< 25.0];
constraints[T.obs_addr()] = T.voxelize(obs_cloud, params.hyperparams.resolution)
MeshCatViz.reset_visualizer()
MeshCatViz.viz(constraints[T.obs_addr()])

trace, _ = generate(T.scene, (params,), constraints);
# -

MeshCatViz.reset_visualizer()
MeshCatViz.viz(T.get_obs_cloud(trace); color=I.colorant"black", channel_name=:obs_cloud)
MeshCatViz.viz(T.get_gen_cloud(trace); color=I.colorant"red", channel_name=:gen_cloud)

traces = [trace];

for i in 1:length(ids)
    if !T.isFloating(T.get_structure(trace),i)
        continue
    end
    pose_addr = T.floating_pose_addr(i)
    for _ in 1:50
        trace, acc = T.drift_move(trace, pose_addr, 0.5, 10.0)
        if acc traces = push!(traces, trace) end

        trace, acc = T.drift_move(trace, pose_addr, 1.0, 100.0)
        if acc traces = push!(traces, trace) end
        trace, acc = T.drift_move(trace, pose_addr, 1.5, 1000.0)
        if acc traces = push!(traces, trace) end
        trace, acc = T.drift_move(trace, pose_addr, 0.1, 1000.0)
        if acc traces = push!(traces, trace) end
        trace, acc = T.drift_move(trace, pose_addr, 0.1, 100.0)
        if acc traces = push!(traces, trace) end
        trace, acc = T.drift_move(trace, pose_addr, 0.5, 5000.0)
        if acc traces = push!(traces, trace) end
        trace, acc = T.pose_flip_move(trace, pose_addr, 1, 1000.0)
        if acc traces = push!(traces, trace) end
        trace, acc = T.pose_flip_move(trace, pose_addr, 2, 1000.0)
        if acc traces = push!(traces, trace) end
        trace, acc = T.pose_flip_move(trace, pose_addr, 3, 1000.0)
        if acc traces = push!(traces, trace) end
    end
end

CRAZY_POSE = Pose(-100.0 * ones(3), T.IDENTITY_ORIENTATION)
for i in 1:length(ids)
    if !T.isFloating(T.get_structure(trace),i)
        continue
    end
    addr = T.floating_pose_addr(i)

    t, = Gen.update(trace, Gen.choicemap(T.floating_pose_addr(i) => CRAZY_POSE))
    valid = T.get_unexplained_obs_cloud(t)

    p = T.icp_object_pose(trace[addr], trace[:obs][:,valid],
                          p-> get_cloud_from_ids_and_poses(
                            [ids[i]], [p], cam_pose), false)
    for _ in 1:3
        trace, acc = T.pose_mixture_move(trace, addr, [trace[addr], p], [0.5, 0.5], 0.001, 5000.0)
        if acc traces = push!(traces, trace) end
    end
end

MeshCatViz.reset_visualizer()
MeshCatViz.viz(T.get_obs_cloud(trace); color=I.colorant"black", channel_name=:obs_cloud)
MeshCatViz.viz(T.get_gen_cloud(trace); color=I.colorant"red", channel_name=:gen_cloud)

depth_image = GL.gl_render(
    renderer, ids, T.get_poses(trace)[1:end-1], cam_pose)
depth_image[depth_image .== 50000.0] .= 200.0
img = I.colorview(I.Gray, depth_image ./ maximum(depth_image))

# +
renderer_texture = GL.setup_renderer(original_camera, GL.TextureMode())
obj_paths = T.load_ycb_model_obj_file_paths(YCB_DIR)
texture_paths = T.load_ycb_model_texture_file_paths(YCB_DIR)
for id in all_ids
    v,n,f,t = renderer_texture.gl_instance.load_obj_parameters(
        obj_paths[id]
    )
    v = v * world_scaling_factor
    v .-= id_to_shift[id]'

    GL.load_object!(renderer_texture, v, n, f, t,
        texture_paths[id]
    )
end

# -

rgb_image, depth_image = GL.gl_render(
    renderer_texture, ids, T.get_poses(trace)[1:end-1], cam_pose)
I.colorview(I.RGBA, permutedims(rgb_image,(3,1,2)))


text:
function load(path)
  groups = []
  open(path) do group_data
    group = []
    while !eof(group_data)
      line = readline(group_data)
      if line == ""
        push!(groups, join(group, ','))
        empty!(group)
      else
        push!(group, strip(line))
      end
    end

    push!(groups, join(group, ','))
  end

  groups
end

function count_answers(data)
  answer_counts :: Array{Int, 1} = []
  total_answers = 0
  total_everyone = 0
  for group in data
    num_members = 1
    answer_counts = zeros(26)
    for c in group
      if c != ','
        answer_counts[convert(Int, c) - 96] += 1
      else
        num_members += 1
      end
    end

    new_answers = 0
    new_everyone = 0
    for c in answer_counts
      if c > 0
        new_answers += 1
        if c == num_members
          new_everyone += 1
        end
      end
    end

    total_answers += new_answers
    total_everyone += new_everyone
  end

  total_answers, total_everyone
end

function solve(data)
  println(@time count_answers(data))
end

data = load("input.txt")
solve(data)

max_stars_repo_name: wbthomason/advent

max_issues_repo_name: wbthomason/advent

max_forks_repo_name: wbthomason/advent


text:
@testset "ZafarData" begin
	X, y = genZafarData()
	@test scitype(y) == AbstractArray{Multiclass{2}, 1}
	@test string.(names(X)) == ["x1", "x2", "z"]
end

@testset "ZafarData2" begin
  X, y = genZafarData2()
  @test scitype(y) == AbstractArray{Multiclass{2}, 1}
  @test string.(names(X)) == ["x1", "x2", "z"]
end

@testset "SubgroupData" begin
	X, y = genSubgroupData(;setting = "B1")
	@test scitype(y) <: AbstractVector{<:Finite}
	@test string.(names(X)) == ["x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8", "x9", "x10", "z"]
end

@testset "BiasedSampleData" begin
  X, y = genBiasedSampleData()
  @test scitype(y) == AbstractArray{Multiclass{2}, 1}
  @test string.(names(X)) == ["x1", "x2", "x3", "x4", "z"]
end

@testset "Law School" begin
  A = law_school(10)
  @test length(A)==6
end


text:
"""
Strs package

This pulls together a number of different packages from JuliaString.org

It has two main goals, first, to be a drop-in replacement for the built-in `String` and `Char`
types, adding types that are both faster and easier to use, that are also using for interfacing
with other languages, and are safer to use and second,
to provide a better option than the built-in string literal syntax and the `@printf`/`@sprintf`
macros for formatted output.
See StrBase, StrRegex, StrLiterals, StrEntities, StrFormat and Format for more documentation.

Copyright 2017-2020 Gandalf Software, Inc., Scott P. Jones
Licensed under MIT License, see LICENSE.md
"""
module Strs

using ModuleInterfaceTools
export @api

using PCRE2

import StrTables: lookupname, matchchar, matches, longestmatches, completions
export lookupname, matchchar, matches, longestmatches, completions

@api extend! StrRegex, StrLiterals, Format
@api use! StrFormat, StrEntities

using InternedStrings: intern
export @i_str, @I_str, intern

macro i_str(str) ; Expr(:call, :intern, interpolated_parse(str, StrBase._str)) ; end
macro I_str(str) ; Expr(:call, :intern, interpolated_parse(str, StrBase._str, true)) ; end

@api reexport StrRegex, StrLiterals

@api freeze

end # module Strs


text:
# const AbstractScalar{T} = Union{T,AbstractArray{T,0},Ref{T},Tuple{T}}
# using Preferences
# function get_experiment()
#     default_set = "off"
#     set = @load_preference("experiment", default_set)
#     if set ∉ ("on", "off")
#         @error("Invalid setting \"$(provider)\"; valid settings include [\"on\", \"off\"], defaulting to \"off\"")
#         set = default_set
#     end
#     return set
# end

# # Read in preferences, see if any users have requested a particular backend
# const experiment = get_experiment()

# function set_experiment!(set; export_prefs::Bool = false)
#     if set !== nothing && set !== missing && set ∉ ("on", "off")
#         throw(ArgumentError("Invalid set '$(set)'"))
#     end
#     set_preferences!(@__MODULE__, "experiment" => set; export_prefs, force = true)
#     if set != experiment
#         # Re-fetch to get default values in the event that `nothing` or `missing` was passed in.
#         set = get_experiment()
#         @info("Experiment function changed; restart Julia for this change to take effect", set)
#     end
# end

# @static if experiment == "on"
#     @warn "Experiment on, some base behaviour changed"
#     include("nonlazyscalar.jl")
# end
# include("lazygenerator.jl")


text:
# Plot results

export
    plotpcs,
    plotfluctuations

using Gadfly
using Colors
import Cairo

"Dictionary of default colours for plotting."
const default_colours = Dict{String, String}(
    "exprose"=> "#3498DB", # Blue
    "input_1"=> "#E74C3C", # Red
    "input_2"=> "#2ECC71", # Green
    "mod"    => "#FF9900", # Orange
    "extra"  => "#000000", # Black
)

"`Theme` for plots."
plottheme() = Theme(
    background_color=colorant"white",
    panel_stroke=colorant"black",
    major_label_color=colorant"black",
    minor_label_color=colorant"black",
    major_label_font_size=20pt,
    minor_label_font_size=16pt,
    grid_line_width=0mm,
)

"`Theme` for plot points."
pointtheme(colour::AbstractString) = Theme(
    default_color=parse(Colorant, colour),
    point_size=1.5mm,
    highlight_width=0mm,
)

"Plot the projections of an ensemble onto the principal components of the motion."
function plotpcs(out_prefix::AbstractString,
                    pcs::Array{Float64};
                    pcs_ref_one::Array{Float64}=Float64[],
                    pcs_ref_two::Array{Float64}=Float64[],
                    pcs_mod::Array{Float64}=Float64[],
                    pcs_extra::Array{Float64}=Float64[],
                    comps_to_plot::Array{Tuple{Int,Int},1}=[(1,2), (1,3), (2,3)])
    checkfilepath(out_prefix)
    for (i, j) in comps_to_plot
        out_filepath = "$(out_prefix)_$(i)_$j.png"
        @assert length(pcs) > 0 "No projections in array"
        @assert max(i, j) <= size(pcs, 1) "Chose a PC number higher than those available"
        layers = layer(
            x=pcs[i,:],
            y=pcs[j,:],
            Geom.point,
            pointtheme(default_colours["exprose"]),
        )
        pcsaddlayer!(layers, pcs_ref_one, (i, j), default_colours["input_1"])
        pcsaddlayer!(layers, pcs_ref_two, (i, j), default_colours["input_2"])
        pcsaddlayer!(layers, pcs_mod, (i, j), default_colours["mod"])
        pcsaddlayer!(layers, pcs_extra, (i, j), default_colours["extra"])
        p = plot(
            layers...,
            Guide.xlabel("Projections onto PC $i / Å"),
            Guide.ylabel("Projections onto PC $j / Å"),
            plottheme,
        )
        draw(PNG(out_filepath, 20cm, 20cm), p)
    end
    println("Plotted principal components to file(s) \"$(out_prefix)_x_y.png\"")
end

"Adds a layer of points to an existing layer for PC plotting."
function pcsaddlayer!(layers::Array{Layer,1}, pcs::Array{Float64}, comps::Tuple{Int,Int}, colour::AbstractString)
    if length(pcs) > 0
        i, j = comps
        @assert max(i, j) <= size(pcs, 1) "Chose a PC number higher than those available"
        push!(layers, layer(
            x=pcs[i,:],
            y=pcs[j,:],
            Geom.point,
            pointtheme(colour),
        )[1])
    end
end

"""
Plot root mean square fluctuations of an aligned ensemble.
Arguments are the output filepath and the array of fluctuations.
"""
function plotfluctuations(out_filepath::AbstractString,
                    flucs::Array{Float64,1};
                    flucs_mod::Array{Float64,1}=Float64[])
    checkfilepath(out_filepath)
    n_res = length(flucs)
    @assert n_res > 0 "Fluctuation array is empty"
    layers = layer(
        x=collect(1:n_res),
        y=flucs,
        Geom.line,
        pointtheme(default_colours["exprose"]),
    )
    if length(flucs_mod) > 0
        @assert length(flucs_mod) == n_res "Number of fluctuations differ"
        push!(layers, layer(
            x=collect(1:n_res),
            y=flucs_mod,
            Geom.line,
            pointtheme(default_colours["mod"]),
        )[1])
    end
    p = plot(
        layers...,
        Guide.xlabel("Residue index"),
        Guide.ylabel("RMSF / Å"),
        plottheme,
    )
    draw(PNG(out_filepath, 20cm, 15cm), p)
    println("Plotted RMSFs to file \"$out_filepath\"")
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule extxyz_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("extxyz")
JLLWrappers.@generate_main_file("extxyz", UUID("6ecdc6fc-93a8-5528-aee3-ac7ae1c60be7"))
end  # module extxyz_jll


text:
#!/usr/bin/env julia
# Copyright (c) 2017 Kevin Murray <kdmfoss@gmail.com>
# This Source Code Form is subject to the terms of the Mozilla Public License,
# v. 2.0. If a copy of the MPL was not distributed with this file, You can
# obtain one at http://mozilla.org/MPL/2.0/.

module BarcodeFactory

import Combinatorics
import Iterators
using Bio.Seq
import DataStructures
import YAML

immutable Barcode
	name::String
	r1barcode::DNASequence
	r2barcode::Nullable{DNASequence}
	Barcode(n, s1::DNASequence, s2::DNASequence) = new(n, s1, Nullable{DNASequence}(s2))
	Barcode(n, s) = new(n, s, Nullable{DNASequence}())
end

abstract Demuxer
immutable Axe <: Demuxer end
immutable AdapterRemoval <: Demuxer end
immutable Fastx <: Demuxer end
immutable Flexbar <: Demuxer end

function Base.write(io::IO, ::Type{Axe}, bcds::Vector{Barcode})
	for bcd in bcds
		r2bcd = isnull(bcd.r2barcode) ? "" : "\t$(get(bcd.r2barcode))"
		println(io, "$(bcd.r1barcode)$r2bcd\t$(bcd.name)")
	end
end

function Base.write{T<:Union{AdapterRemoval, Fastx}}(io::IO, ::Type{T},
                                                     bcds::Vector{Barcode})
	for bcd in bcds
		r2bcd = isnull(bcd.r2barcode) ? "" : "\t$(get(bcd.r2barcode))"
		println(io, "$(bcd.name)\t$(bcd.r1barcode)$r2bcd")
	end
end

function Base.write(io::IO, ::Type{Flexbar}, bcds::Vector{Barcode})
	for bcd in bcds
		println(io, ">$(bcd.name)\n$(bcd.r1barcode)")
	end
end

Base.write{T<:Demuxer}(::Type{T}, bcds) = write(STDOUT, T, bcds)

"""
Computes the cloud of words, w \in A^k, such that d(w_i, w_j) >= dist for all i
and j
"""
function hamming_space(k, dist)
    accepted = Vector{DNASequence}()
    for candidate_nt in Iterators.product(repeated(ACGT, k)...)
        candidate = DNASequence(DNANucleotide[candidate_nt...])
        acceptable = true
        for j in accepted
            if mismatches(candidate, j) < dist
                acceptable = false
                break
            end
        end
        if acceptable
            push!(accepted, candidate)
        end
    end
    return accepted
end

"Generate a list of alphabetic codes of length >= n"
function alphacode(n)
    A = collect('A':'Z')
    len = Int(ceil(log(n)/log(26)))
    codes = String[]
    for code in Iterators.product(repeated(A, len)...)
        push!(codes, join(code))
    end
    return codes
end

function generate_index_set(num_indices::Vector{Int}, idx_length::Int,
                            distance::Int=3; maxlength::Int=0)
    barcodes = Barcode[]
    all_index_seqs = hamming_space(idx_length, distance)
    shuffle!(all_index_seqs)

    totaln = prod(num_indices)
    maxn = maximum(num_indices)
    if maxn > length(all_index_seqs)
        error("More indicies requested than available with "*
              "length $idx_length and distance $distance")
    end

    idx_names = alphacode(totaln)

    if maxlength > idx_length
        for i in 1:length(all_index_seqs)
            additional = i % (maxlength - idx_length)
            all_index_seqs[i] *= dna"A"^additional
        end
    end

    # For each pair in all R1 & R2 index seqs
    for (i, index) in enumerate(Iterators.product(map(x -> (1:x), num_indices)...))
        idx_seqs = collect(map(x -> all_index_seqs[x], index))
        push!(barcodes, Barcode(idx_names[i], idx_seqs...))
    end
    return barcodes
end

const DEMUXERS = Dict{String, Type}(
        "axe" => Axe,
        "ar" => AdapterRemoval,
        "flexbar" => Flexbar,
        "fastx" => Fastx,
)

function generate_from_yaml(yamlfile::String, root::String; seed::Int=3301)
    sets = YAML.load_file(yamlfile)
    for set in sets
        srand(seed)
        setname = set["name"]
        indices = generate_index_set(set["num_indices"], set["length"],
                                     set["dist"], maxlength=get(set, "maxlength", 0))
        for demuxer in set["demuxers"]
            ext = ".$demuxer"
            if demuxer == "flexbar"
                ext = "_flexbar.fasta"
            end
            open("$root/$setname$ext", "w") do outf
                write(outf, DEMUXERS[demuxer], indices)
            end
        end
    end
end

end # module BarcodeFactory

using ArgParse

function main()
    ap = ArgParseSettings()
    @add_arg_table ap begin
        "settings"
            help="YAML-encoding settings file"
            required=true
            arg_type=String
        "rootdir"
            help="Output directory"
            default="."
            required=false
            arg_type=String
        "--seed"
            help="Random seed"
            default=3301
            arg_type=Int
    end
    args = parse_args(ap)

    BarcodeFactory.generate_from_yaml(args["settings"], args["rootdir"],
                                      seed=args["seed"])
end

main()


text: include("Misfits.jl")


text:
function part1(input)
    expenses = falses(2021)
    for expense_string in readlines(input)
        expense = parse(Int, expense_string)
        if expenses[1 + 2020 - expense]
            return expense * (2020 - expense)
        end
        expenses[1 + expense] = true
    end
    error("No solution found.")
end

function part2(input)
    expenses = Int[]
    two_expenses = Dict{Int, Int}()
    for expense_string in readlines(input)
        expense = parse(Int, expense_string)
        if haskey(two_expenses, 2020 - expense)
            return two_expenses[2020 - expense] * expense
        end
        for expense2 in expenses
            two_expenses[expense + expense2] = expense * expense2
        end
        push!(expenses, expense)
    end
    error("No solution found.")
end


text:
using DifferentialEquations
using SimpleDiffEq
using DiffEqCallbacks
using Random
using Distributions
using DiffEqParamEstim
using DataFrames
using StatsPlots
using BenchmarkTools

function sir_ode!(du,u,p,t)
    (S,I,R,Y) = u
    (β,c,γ) = p
    N = S+I+R
    infection = β*c*I/N*S
    recovery = γ*I
    @inbounds begin
        du[1] = -infection
        du[2] = infection - recovery
        du[3] = recovery
        du[4] = infection
    end
    nothing
end;

δt = 0.1
tmax = 40.0
tspan = (0.0,tmax)
t = 0.0:δt:tmax
obstimes = 0:1.0:tmax;

u0 = [990.0,10.0,0.0,0.0]; # S,I.R,Y

p = [0.05,10.0,0.25]; # β,c,γ

affect!(integrator) = integrator.u[4] = 0.0
cb_zero = PresetTimeCallback(obstimes,affect!)

prob_ode = ODEProblem(sir_ode!,u0,tspan,p)

sol_ode = solve(prob_ode,callback=cb_zero);

df_ode = DataFrame(sol_ode(obstimes)')
df_ode[!,:t] = obstimes;

@df df_ode plot(:t,
    [:x1 :x2 :x3 :x4],
    label=["S" "I" "R" "Y"],
    xlabel="Time",
    ylabel="Number")

Random.seed!(1234);

data = rand.(Poisson.(df_ode[!,:x4]))

plot(obstimes,data)
plot!(obstimes,df_ode[!,:x4])

using Optim

function ss1(β)
    prob = remake(prob_ode,u0=[990.0,10.0,0.0,0.0],p=[β,10.0,0.25])
    sol = solve(prob,Tsit5(),callback=cb_zero,saveat=obstimes)
    sol_data = sol(obstimes)[4,:]
    return(sum((sol_data - data) .^2))
end

function nll1(β)
    prob = remake(prob_ode,u0=[990.0,10.0,0.0,0.0],p=[β,10.0,0.25])
    sol = solve(prob,Tsit5(),callback=cb_zero,saveat=obstimes)
    sol_data = sol(obstimes)[4,:]
    -sum(logpdf.(Poisson.(sol_data),data))
end

lower1 = 0.0
upper1 = 1.0
initial_x1 = 0.1

opt1_ss = Optim.optimize(ss1,lower1,upper1)

opt1_nll = Optim.optimize(nll1,lower1,upper1)

function ss2(x)
    (i0,β) = x
    I = i0*1000.0
    prob = remake(prob_ode,u0=[1000-I,I,0.0,0.0],p=[β,10.0,0.25])
    sol = solve(prob,Tsit5(),callback=cb_zero,saveat=obstimes)
    sol_data = sol(obstimes)[4,:]
    return(sum((sol_data - data) .^2))
end

function nll2(x)
    (i0,β) = x
    I = i0*1000.0
    prob = remake(prob_ode,u0=[1000-I,I,0.0,0.0],p=[β,10.0,0.25])
    sol = solve(prob,Tsit5(),callback=cb_zero,saveat=obstimes)
    sol_data = sol(obstimes)[4,:]
    -sum(logpdf.(Poisson.(sol_data),data))
end

lower2 = [0.0,0.0]
upper2 = [1.0,1.0]
initial_x2 = [0.01,0.1]

opt2_ss = Optim.optimize(ss2,lower2,upper2,initial_x2)

opt2_nll = Optim.optimize(nll2,lower2,upper2,initial_x2)

function loss_function(sol)
    sol_data = DataFrame(sol(obstimes)')[!,:x4]
    -sum(logpdf.(Poisson.(sol_data),data))
end

prob_generator = (prob,q) -> remake(prob,
    u0=[1000-(q[1]*1000),q[1]*1000,0.0,0.0],
    p=[q[2],10.0,0.25])

cost_function = build_loss_objective(prob_ode,
    Tsit5(),
    loss_function,
    prob_generator = prob_generator,
    maxiters=10000,
    verbose=false,
    callback=cb_zero)

opt_pe1 = Optim.optimize(cost_function,lower2,upper2,initial_x2)

using NLopt
opt = Opt(:LD_MMA, 2)
opt.lower_bounds = lower2
opt.upper_bounds = upper2
opt.min_objective = cost_function
(minf,minx,ret) = NLopt.optimize(opt,initial_x2)

using BlackBoxOptim
bound1 = Tuple{Float64, Float64}[(0.0,1.0),(0.0, 1.0)]
result = bboptimize(cost_function;SearchRange = bound1, MaxSteps = 110e3)

include(joinpath(@__DIR__,"tutorials","appendix.jl"))
appendix()

max_stars_repo_name: FSUcilab/sir-julia

max_issues_repo_name: FSUcilab/sir-julia

max_forks_repo_name: FSUcilab/sir-julia


text:
"""
    FilePaths.@compat fn

Generates a compatibility method for handling paths as strings based on the function
definition passed in.

Wrapper method properties:

- Any arguments that accepted `P <: AbstractPath` will accept `Union{String, P}`
- All string path inputs will be converted to path types
- If a path return type was specified then the path result will be converted back to a string.

# Examples

```
julia> using FilePathsBase; using FilePathsBase: /; using FilePaths

julia> FilePaths.@compat function myrelative(x::AbstractPath, y::AbstractPath)
           return relative(x, y)
       end
myrelative (generic function with 2 methods)

julia> FilePaths.@compat function myjoin(x::P, y::String)::P where P <: AbstractPath
           return x / y
       end
myjoin (generic function with 2 methods)

julia> myrelative(cwd(), home())
p"repos/FilePaths.jl"

julia> myrelative(pwd(), homedir())
p"repos/FilePaths.jl"

julia> myjoin(parent(cwd()), "FilePaths.jl")
p"/Users/rory/repos/FilePaths.jl"

julia> myjoin("/Users/rory/repos", "FilePaths.jl")
"/Users/rory/repos/FilePaths.jl"
```
"""
macro compat(ex)
    mod::Module = QuoteNode(__module__).value
    new_ex = compat_exp(mod, deepcopy(ex))

    return quote
        $ex
        $new_ex
    end |> esc
end

function compat_exp(mod::Module, ex::Expr)
    fdef = splitdef(ex)
    args = Symbol[]
    kwargs = Expr[]
    convert_vars = Symbol[]
    body = Expr[]
    params = Symbol[]

    # A function that identifies args that need to be converted
    function parse_arg(a)
        if a.args[2] in params
            push!(convert_vars, a.args[1])
        elseif _ispath(mod, a.args[2])
            # Modify the arg type declaration from P<:AbstractPath to Union{String, P}
            # NOTE: If the variable is parameterized then we should have already updated the
            # parameterized type in the where clause
            a.args[2] = :(Union{AbstractString, $(a.args[2])})
            push!(convert_vars, a.args[1])
        end
    end

    # Identify any where params that are a subtype of AbstractPath
    if haskey(fdef, :whereparams)
        for (i, p) in enumerate(fdef[:whereparams])
            # If the param is a subtype of AbstracPath
            # then we store the lookup symbol in the params array
            if _ispath(mod, p.args[2])
                # Modify parameterized where clause from P<:AbstractPath to Union{String, P}
                p.args[2] = :(Union{AbstractString, $(p.args[2])})
                push!(params, p.args[1])
            end
        end
    end

    # Identify args that need to be converted
    if haskey(fdef, :args)
        for (i, a) in enumerate(fdef[:args])
            # An arg can be an expression or a symbol (no type information)
            if isa(a, Expr)
                if (a.head) === :kw
                    # Optional arguments show up as `kw` and need to be parsed as such
                    T = a.args[1]
                    parse_arg(T)
                    push!(args, T.args[1])
                else
                    parse_arg(a)
                    push!(args, a.args[1])
                end
            elseif isa(a, Symbol)
                push!(args, a)
            else
                throw(ArgumentError("Uknown argument type $a"))
            end
        end
    end

    # Identify kwargs that need to be converted
    if haskey(fdef, :kwargs)
        for (i, k) in enumerate(fdef[:kwargs])
            T = k.args[1]
            parse_arg(T)
            # Rewrite the kwarg expression to pass them along
            push!(kwargs, :($(T.args[1])=$(T.args[1])))
        end
    end

    # Insert our convert statements for the appropriate variables
    for v in convert_vars
        push!(body, :($v = Path($v)))
    end

    # Push our paths method call into the body
    push!(body, :(result = $(fdef[:name])($(args...); $(kwargs...))))

    # If we have a return type and it's a path then convert the result back to a string
    if haskey(fdef, :rtype) && (fdef[:rtype] in params || _ispath(mod, fdef[:rtype]))
        push!(body, :(result = string(result)))
        # Set the return type to String to avoid an incorrect conversion back to a path
        fdef[:rtype] = :String
    end

    # Finally, insert a return statement
    push!(body, :(return result))

    # Update our definition with the new body
    fdef[:body].args = body

    # Combine the modified definition back into an expression
    return MacroTools.combinedef(fdef)
end

function _ispath(mod::Module, t::Symbol)
    T = getfield(mod, t)
    return T <: AbstractPath
end