text:
using LearnBase
using Base.Test

# Test if types are exported properly
@test Cost <: Any
@test Penalty <: Cost
@test Loss <: Cost
@test UnsupervisedLoss <: Loss
@test SupervisedLoss <: Loss
@test MarginLoss <: SupervisedLoss
@test DistanceLoss <: SupervisedLoss

@test Minimizable <: Any
@test Transformation <: Any
@test StochasticTransformation <: Transformation

# Test if functions are exported properly
@test typeof(transform) <: Function
@test typeof(transform!) <: Function
@test typeof(getobs) <: Function
@test typeof(getobs!) <: Function
@test typeof(learn) <: Function
@test typeof(learn!) <: Function
@test typeof(update) <: Function
@test typeof(update!) <: Function

@test typeof(value) <: Function
@test typeof(value!) <: Function
@test typeof(meanvalue) <: Function
@test typeof(sumvalue) <: Function
@test typeof(meanderiv) <: Function
@test typeof(sumderiv) <: Function
@test typeof(deriv) <: Function
@test typeof(deriv2) <: Function
@test typeof(deriv!) <: Function
@test typeof(value_deriv) <: Function
@test typeof(grad) <: Function
@test typeof(grad!) <: Function
@test typeof(addgrad!) <: Function
@test typeof(value_grad) <: Function
@test typeof(value_grad!) <: Function
@test typeof(prox) <: Function
@test typeof(prox!) <: Function

@test typeof(isminimizable) <: Function
@test typeof(isdifferentiable) <: Function
@test typeof(istwicedifferentiable) <: Function
@test typeof(isconvex) <: Function
@test typeof(isstrictlyconvex) <: Function
@test typeof(isstronglyconvex) <: Function
@test typeof(isnemitski) <: Function
@test typeof(isunivfishercons) <: Function
@test typeof(isfishercons) <: Function
@test typeof(islipschitzcont) <: Function
@test typeof(islocallylipschitzcont) <: Function
@test typeof(islipschitzcont_deriv) <: Function
@test typeof(isclipable) <: Function
@test typeof(ismarginbased) <: Function
@test typeof(isclasscalibrated) <: Function
@test typeof(isdistancebased) <: Function

@test typeof(issymmetric) <: Function

@test typeof(fit) <: Function
@test typeof(fit!) <: Function
@test typeof(nobs) <: Function

# Test that LearnBase reuses StatsBase functions
using StatsBase
@test StatsBase.fit  == LearnBase.fit
@test StatsBase.fit! == LearnBase.fit!
@test StatsBase.nobs == LearnBase.nobs

# Test superset fallbacks
immutable MyStronglyConvexType end
LearnBase.isstronglyconvex(::MyStronglyConvexType) = true
LearnBase.islipschitzcont(::MyStronglyConvexType) = true
@test isstronglyconvex(MyStronglyConvexType())
@test isstrictlyconvex(MyStronglyConvexType())
@test isconvex(MyStronglyConvexType())
@test islipschitzcont(MyStronglyConvexType())
@test islocallylipschitzcont(MyStronglyConvexType())

# IntervalSet
let s = IntervalSet(-1,1)
    @test typeof(s) == IntervalSet{Int}
    @test typeof(s) <: AbstractSet
    for x in (-1,0,0.5,1,1.0)
        @test x in s
    end
    for x in (-1-1e-10, 1+1e-10, -Inf, Inf, 2, NaN)
        @test !(x in s)
    end
    for i=1:10
        x = rand(s)
        @test typeof(x) == Float64
        @test x in s
    end
    xs = rand(s, 10)
    @test typeof(xs) == Vector{Float64}
    for x in xs
        @test typeof(x) == Float64
        @test x in s
    end
    @test LearnBase.randtype(s) == Float64
    # @show s LearnBase.randtype(s)
end
let s = IntervalSet(-1,1.0)
    @test typeof(s) == IntervalSet{Float64}
    @test typeof(s) <: AbstractSet
    @test 1 in s
    # @show s LearnBase.randtype(s)
end

# DiscreteSet
let s = DiscreteSet([-1,1])
    @test typeof(s) == DiscreteSet{Vector{Int}}
    @test typeof(s) <: AbstractSet
    for x in (-1, 1, -1.0, 1.0)
        @test x in s
    end
    for x in (0, Inf, -Inf, NaN)
        @test !(x in s)
    end
    for i=1:10
        x = rand(s)
        @test typeof(x) == Int
        @test x in s
    end
    xs = rand(s, 10)
    @test typeof(xs) == Vector{Int}
    for x in xs
        @test typeof(x) == Int
        @test x in s
    end
    @test LearnBase.randtype(s) == Int
    @test length(s) == 2
    @test s[1] == -1
    # @show s LearnBase.randtype(s)
end
let s = DiscreteSet([-1,1.0])
    @test typeof(s) == DiscreteSet{Vector{Float64}}
    @test typeof(s) <: AbstractSet
    @test typeof(rand(s)) == Float64
    @test typeof(rand(s, 2)) == Vector{Float64}
    # @show s LearnBase.randtype(s)
end

# TupleSet
let s = TupleSet(IntervalSet(0,1), DiscreteSet([0,1]))
    @test typeof(s) == TupleSet{Tuple{IntervalSet{Int}, DiscreteSet{Vector{Int}}}}
    @test typeof(s) <: AbstractSet
    for x in ([0,0], [0.0,0.0], [0.5,1.0])
        @test x in s
    end
    for x in ([0,0.5], [-1,0])
        @test !(x in s)
    end
    @test typeof(rand(s)) == Vector{Float64}
    @test typeof(rand(s, 2)) == Vector{Vector{Float64}}
    @test typeof(rand(s, Tuple)) == Tuple{Float64,Int}
    @test typeof(rand(s, Tuple, 2)) == Vector{Tuple{Float64,Int}}
    @test LearnBase.randtype(s) == Vector{Float64}
    # @show s LearnBase.randtype(s)
end

# arrays of sets
let s = [IntervalSet(0,1), DiscreteSet([0,1])]
    @test typeof(s) == Vector{AbstractSet}
    for x in ([0,0], [0.0,0.0], [0.5,1.0])
        @test x in s
    end
    for x in ([0,0.5], [-1,0])
        @test !(x in s)
    end
    @test typeof(rand(s)) == Vector{Float64}
    @test typeof(rand(s, 2)) == Vector{Vector{Float64}}
    # @show s LearnBase.randtype(s)
end


text:
test("Utility path intersection, xor, union, difference") do
	pv=[(0.,0.),(2.,0.),(0.,2.)]
	cv=[(0.,0.),(1.,0.),(1.,1.),(0.,1.)]

	success,intersection=pathintersection(pv,cv)
	@test success
	@test intersection[1] == cv[[3,4,1,2]] # Is this always true or is the permutation system dependent?

	success,exclusiveor=pathxor(pv,cv)
	@test success
	@test exclusiveor[1] == [(1.0,0.0),(2.0,0.0),(0.0,2.0),(0.0,1.0),(1.0,1.0)]

	success,union=pathunion(pv,cv)
	@test success
	@test union[1] == pv[[3,1,2]] # Again, is this permutation always true?

	success,difference=pathdifference(pv,cv)
	@test success
	@test difference[1] == exclusiveor[1] # a quirk of this path/clip combination

end

ext: jl


text:
_can_read_file(file_name::String) = return isfile(file_name) && try isreadable(open(file_name, "r")) catch e; false; end
_begins_with_ec2(file_name::String) = return uppercase(String(read(file_name, 3))) == "EC2"
_ends_with_ec2(file_name::String) = return endswith(strip(uppercase(read(file_name, String))), "EC2")

"""
Generate a valid role session name. Currently only ensures that the session name is
64-characters or less.
"""
function _role_session_name(prefix, name, suffix)
    b = IOBuffer()
    write(b, prefix, name)
    truncate(b, min(64 - length(suffix), b.size))  # Assumes ASCII
    write(b, suffix)
    return String(take!(b))
end

"""
Get the value for `key` in the `ini` file for a given `profile`.
"""
function _get_ini_value(
    ini::Inifile, profile::AbstractString, key::AbstractString;
    default_value=nothing
)
    value = get(ini, "profile $profile", key)
    value === :notfound && (value = get(ini, profile, key))
    value === :notfound && (value = default_value)

    return value
end

function _aws_profile_config(ini::Inifile, profile::AbstractString)
    if profile != "default"
        profile = "profile $profile"
    end

    return get(sections(ini), profile, Dict())
end

function _aws_profile_config(ini::Inifile, profile::Nothing)
    _aws_profile_config(ini, _aws_get_profile())
end

function _aws_profile_config(config_file::AbstractString, profile)
    isfile(config_file) || return Dict()
    return _aws_profile_config(read(Inifile(), config_file), profile)
end

function _aws_profile_config(config_file::Nothing, profile)
    return _aws_profile_config(dot_aws_config_file(), profile)
end

"""
Retrieve the `AWSCredentials` for a given role from Security Token Services (STS).
"""
function _aws_get_role(role::AbstractString, ini::Inifile)
    source_profile, role_arn = aws_get_role_details(role, ini)
    source_profile === nothing && return nothing
    credentials = nothing

    for f in (dot_aws_credentials, dot_aws_config)
        credentials = f(source_profile)
        credentials === nothing || break
    end

    credentials === nothing && return nothing
    config = AWSConfig(creds=credentials, region=aws_get_region(config=ini, profile=source_profile))

    # RoleSessionName Documentation
    # https://docs.aws.amazon.com/STS/latest/APIReference/API_AssumeRole.html
    role = AWSServices.sts(
        "AssumeRole",
        LittleDict("RoleArn" => role_arn, "RoleSessionName" => replace(role, r"[^\w+=,.@-]" => s"-"));
        aws_config=config
    )

    role_creds = role["AssumeRoleResult"]["Credentials"]

    return AWSCredentials(
        role_creds["AccessKeyId"],
        role_creds["SecretAccessKey"],
        role_creds["SessionToken"];
        expiry=DateTime(rstrip(role_creds["Expiration"], 'Z'))
    )
end

"""
Get `AWSCredentials` for the specified `profile` from the `inifile`. If targeting the
`~/.aws/config` file, with a non-default `profile`, you must specify `config=true` otherwise
the default credentials will be returned.
"""
function _aws_get_credential_details(profile::AbstractString, ini::Inifile)
    access_key = _get_ini_value(ini, profile, "aws_access_key_id")
    secret_key = _get_ini_value(ini, profile, "aws_secret_access_key")
    token = _get_ini_value(ini, profile, "aws_session_token"; default_value="")

    return (access_key, secret_key, token)
end

"""
Get the default AWS profile
"""
function _aws_get_profile(; default="default")
    @something(
        get(ENV, "AWS_PROFILE", nothing),
        get(ENV, "AWS_DEFAULT_PROFILE", nothing),
        Some(default),
    )
end

"""
Check if credentials will expire within 5 minutes
"""
function _will_expire(aws_creds::AWSCredentials)
    return aws_creds.expiry - now(UTC) <= Minute(5)
end

"""
Retrieve the EC2 meta data from the local AWS endpoint. Return the EC2 metadata request
body, or `nothing` if not running on an EC2 instance.
"""
function _ec2_metadata(metadata_endpoint::String)
    try
        request = @mock HTTP.request("GET", "http://169.254.169.254/latest/meta-data/$metadata_endpoint")

        return request === nothing ? nothing : String(request.body)
    catch e
        e isa HTTP.IOError || e isa HTTP.StatusError && e.status == 404 || rethrow(e)
    end

    return nothing
end


text:
module U2Net
using Flux
using Images

include("utils.jl")
include("model_layers.jl")
include("u2net_model.jl")
include("u2net_s_model.jl")

end # module


text:
s = open("input.txt") do file
    read(file, String);
end

parts = split(s, "\n\n");

rules = split(parts[1],"\n")[1:end];
myticket = split(parts[2],"\n")[2];
otherTickets = split(parts[3],"\n")[2:end];

function GetFields(rules)
    fields = Dict();

    for rule in rules
        m = match(r"([a-z ]+): (\d+)-(\d+) or (\d+)-(\d+)", rule);
        name = m.captures[1];
        nums = map(x->parse(Int,x), m.captures[2:end]);

        fields[name] = nums;
    end

    return fields;
end

function GetValidRanges(rules)
    validRanges = [];

    for rule in rules
        m = match(r": (\d+)-(\d+) or (\d+)-(\d+)", rule);
        nums = map(x->parse(Int,x), m.captures);
        push!(validRanges, nums[1:2])
        push!(validRanges, nums[3:4])
    end

    return validRanges;
end

function CheckValidTickets(ranges, tickets)
    validTickets = falses(length(tickets));

    for i = 1:length(tickets)
        ticket = tickets[i]
        if length(ticket) == 0
            continue
        end

        vals = map(x->parse(Int,x), split(ticket,","));
        validTicket = true;
        for val in vals
            valid = false;
            for range in ranges
                if range[1] <= val <= range[2]
                    valid = true;
                    break;
                end
            end

            if !valid
                validTicket = false;
                break;
            end
        end

        validTickets[i] = validTicket;
    end

    return validTickets;
end

function CouldBeField(x, ranges)
    return ranges[1] <= x <= ranges[2] || ranges[3] <= x <= ranges[4]
end

function GetFieldOrder(tickets, fields)
    ticketNums = vcat(map(x->map(y->parse(Int,y),split(x,",")), tickets)...);
    ticketNums = reshape(ticketNums, (length(fields),length(tickets)));

    fieldOrder = [];

    for i = 1:length(fields)
        col = ticketNums[i,:];

        possibles = [];
        for (field, ranges) in fields
            if all(map(x->CouldBeField(x, ranges), col))
                push!(possibles, field)
            end
        end

        push!(fieldOrder, possibles);
    end

    println(fieldOrder)
    println(map(x->length(x), fieldOrder))
    return fieldOrder;
end

function DetermineFields!(allFieldPossibilities)
    lens = map(x->length(x), allFieldPossibilities);
    smallToBig = sortperm(lens);

    taken = Set()
    for i = 1:length(smallToBig)
        index = smallToBig[i];
        options = allFieldPossibilities[index];

        trueOptions = []
        for opt in options
            if !in(opt, taken)
                push!(trueOptions, opt);
            end
        end

        if length(trueOptions) == 1
            allFieldPossibilities[index] = trueOptions[1];
            push!(taken, trueOptions[1]);
        end

    end

    println(allFieldPossibilities)
    return allFieldPossibilities
end

validRanges = GetValidRanges(rules);
fields = GetFields(rules);

validTickets = CheckValidTickets(validRanges, otherTickets);

trueTickets = otherTickets[validTickets];

possibles = GetFieldOrder(trueTickets, fields);

realOrder = DetermineFields!(possibles);

myTicketVals = map(x->parse(Int,x), split(myticket,","));

departures = map(x->occursin("departure", x), realOrder);
answer = *(myTicketVals[departures]...)


text:
using StudentHousing, StatsBase, Distributions, JuMP, Gurobi

srand(32)

# =============================================================================
# Data
# =============================================================================
nbedrooms_range = collect(1:2)
nbedrooms_frequency = Weights([0.5, 0.5])
nbathrooms_range = collect(1:2)
nbathrooms_frequency = Weights([0.5, 0.5])
prices_range_pp = [800.0, 1000.0]
area_ranges = [0.0]
market_data = StudentHousing.MarketData(nbedrooms_range, nbedrooms_frequency,
    nbathrooms_range, nbathrooms_frequency, prices_range_pp, area_ranges)
budget = 1e6
problem_data = StudentHousingData(market_data, nhouses = 2, budget = budget,
                  demand_distribution = Uniform(0.9, 1.1))

d = problem_data
nhouses = length(d.houses)

# =============================================================================
# Solve problem
# =============================================================================
m, V_generated, λ_generated, house_choice = solve_house_generation(d)

# =============================================================================
# Look at solution
# =============================================================================
getobjectivevalue(m)

# If not integer...
for i = 1:nhouses
    setcategory(m[:λ][i], :Bin)
end
for λnew in λ_generated
    setcategory(λnew, :Bin)
end
solve(m)
getobjectivevalue(m)

# Recover the solution

# In case one of our initial columns was in the optimal basis
results = zeros(nhouses, length(d.patterns_allow))
for k = 1:nhouses
    if getvalue(m[:λ][k]) ≈ 1.0
        results[k, StudentHousing.house_allowedby(k, d)[1]] = 1.0
    end
end
# Check all the other columns we added
for k = 1:length(V_generated)
    if getvalue(λ_generated[k]) ≈ 1.0
        # Find the assignment that we picked for our house
        patterns = find(V_generated[k] .≈ 1.0)
        # Check the house index
        h = house_choice[k]
        # Mark the assignment
        results[h, patterns] .= 1.0
    end
end

println("House 1 assigned to:", find(results[1,:] .≈ 1.0))
# House 1 assigned to:[7, 11]

# sanity check
StudentHousing.house_allowedby(1, d)
# 5-element Array{Int64,1}:
#   6
#   7
#  11
#  12
#  13

println("House 2 assigned to:", find(results[2,:] .≈ 1.0))
# House 2 assigned to:[4]
# sanity check
StudentHousing.house_allowedby(2, d)
# 11-element Array{Int64,1}:
#   3
#   4
#   5
#   6
#   7
#   8
#   9
#  10
#  11
#  12
#  13

# =============================================================================
# Compare with previous MIP results for medium problems
# =============================================================================
# Run a sanity check by solving the GAP with Gurobi
# This is different to just setting a very big budget and demand=0 in the other
# model because of the way we defined shortage to be nonnegative.
m_master = Model()

function bedrooms_scale(i::Int)
    # Number of possible bedrooms = 1 -> i
    nbedrooms_range = collect(1:i)
    # Some synthetic sampling probabilities
    weights = 2.0.^(-collect(1:i))
    weights[1] = 1 - sum(weights[2:end])
    nbedrooms_frequency = Weights(weights)
    # Return the bedroom data
    nbedrooms_range, nbedrooms_frequency
end

srand(1) # houses constructed randomly
nbathrooms_range = collect(1:3)
nbathrooms_frequency = Weights([0.5; 0.4; 0.1])
prices_range_pp = [800.0, 1000.0]
area_ranges = [0.0, 700.0, 850.0]
budget = 1000.0

nhouses = 50

for i = 1:2
    nbedrooms_range, nbedrooms_frequency = bedrooms_scale(i)
    @time market_data = StudentHousing.MarketData(nbedrooms_range,
        nbedrooms_frequency, nbathrooms_range, nbathrooms_frequency,
        prices_range_pp, area_ranges)
    @time problem_data = StudentHousingData(market_data, nhouses = nhouses, budget = budget,
        demand_distribution = Uniform(0.9, 1.1))
    d = problem_data
    P = StudentHousing.get_npatterns(problem_data)
    println("Using $P patterns, solving MIP:")
    tic()
    m_master, V_generated, λ_generated, house_choice = solve_house_generation(problem_data)
    toc()
    println("CG objective = ", getobjectivevalue(m_master))
end

# This represents a solution to GAP with no budget constraints, i.e. can invest
# in as many houses as we like and try to assign as many people as we can

# House 1 assigned to:[7, 11]
# House 2 assigned to:[3]
#   0.000015 seconds (1 allocation: 64 bytes)
#   0.000225 seconds (327 allocations: 319.641 KiB)

# Using 9 patterns, solving MIP:
# 0.000004 seconds (1 allocation: 64 bytes)
# 0.000132 seconds (327 allocations: 319.641 KiB)
# Using 9 patterns, solving MIP:
# elapsed time: 1.867214478 seconds
# CG objective = 4.0
# 0.000003 seconds (1 allocation: 64 bytes)
# 1.012688 seconds (1.05 M allocations: 1.084 GiB, 26.25% gc time)
# Using 329 patterns, solving MIP:
# elapsed time: 102.900935712 seconds
# CG objective = 40.0

results = zeros(nhouses, length(d.patterns_allow))
for k = 1:nhouses
    if getvalue(m_master[:λ][k]) ≈ 1.0
        results[k, StudentHousing.house_allowedby(k, d)[1]] = 1.0
    end
end
# Check all the other columns we added
for k = 1:length(V_generated)
    if getvalue(λ_generated[k]) ≈ 1.0
        # Find the assignment that we picked for our house
        patterns = find(V_generated[k] .≈ 1.0)
        # Check the house index
        h = house_choice[k]
        # Mark the assignment
        results[h, patterns] .= 1.0
    end
end
for h = 1:nhouses
    println("House $h assigned to:", find(results[h,:] .≈ 1.0))
end

# Sanity check, solve original problem with very large budget and demand=1
srand(1)
gap_model = Model()
for i = 1:2
    nbedrooms_range, nbedrooms_frequency = bedrooms_scale(i)
    @time market_data = StudentHousing.MarketData(nbedrooms_range,
        nbedrooms_frequency, nbathrooms_range, nbathrooms_frequency,
        prices_range_pp, area_ranges)
    @time problem_data = StudentHousingData(market_data, nhouses = nhouses, budget = budget,
        demand_distribution = Uniform(0.9, 1.1))
    d = problem_data
    P = StudentHousing.get_npatterns(problem_data)
    println("Using $P patterns, solving MIP:")
    @time gap_model = StudentHousing.gap_model(problem_data)
    tic()
    @assert solve(gap_model) == :Optimal
    toc()
    println("MIP objective = ", getobjectivevalue(gap_model))
end

# Using 9 patterns, solving MIP:
#   0.000535 seconds (3.66 k allocations: 209.766 KiB)
#   0.000484 seconds (6.19 k allocations: 368.250 KiB)
# elapsed time: 0.005137502 seconds
# MIP objective = 4.0
#   0.000004 seconds (1 allocation: 64 bytes)
#   0.938846 seconds (1.05 M allocations: 1.084 GiB, 27.23% gc time)
# Using 329 patterns, solving MIP:
#   0.016782 seconds (72.50 k allocations: 3.942 MiB, 33.63% gc time)
#   0.015228 seconds (234.89 k allocations: 14.217 MiB)
# elapsed time: 0.087777319 seconds
# MIP objective = 40.0

for h = 1:nhouses
    println("House $h assigned to:", find(getvalue.(gap_model[:assignment])[h,:] .≈ 1.0))
end


text:
# Functions like read and readline are blocking, but they can be made nonblocking with Tasks.

while true

    # Read one character. Enter is required after the character.
    # @async begin
    #     a = read(STDIN, Char)
    #     println("got $a")
    # end

    # Read one line. String includes the \n.
    @async begin
        a = readline(STDIN)
        println("got $a (length $(length(a)))")
    end

    # println("doing things") # Not blocked!
    sleep(0.1)

end

max_stars_repo_name: andrewadare/realtime

max_issues_repo_name: andrewadare/realtime

max_forks_repo_name: andrewadare/realtime


text:
function simulate_copy_number(chr, pos_in_chr, cnv_list)
    for cnv in cnv_list
        if chr == cnv["chrom"]
            # check if it is intersect
            if cnv["start"] < pos_in_chr < cnv["end"]
                return cnv["copy"]
            end
        end
    end
    return 1.0
end


text:
# module quantum

export Bosons, Fermions, CavityMode, MultimodeSystem,
        An, Bn, basis, particledensity, Hamiltonian, Jump_operators

using QuantumOptics
using ..potentials
# using ..multiparticlebasis

# Define system

abstract Particles

type Bosons <: Particles
    Nparticles::Int
    Nlevels::Int
    E0::Float64
    potential::Potential
end

type Fermions <: Particles
    Nparticles::Int
    Nlevels::Int
    E0::Float64
    potential::Potential
end

type CavityMode
    index::Int
    Nmodes::Int
    delta::Float64
    eta::Float64
    kappa::Float64
    U0::Float64
end

type MultimodeSystem
    scaling::Float64
    particles::Particles
    modes::Vector{CavityMode}
end

# Define basis of Hilbert space

basis(s::Bosons) = BosonicNParticleBasis(s.Nparticles, s.Nlevels)
basis(s::Fermions) = FermionicNParticleBasis(s.Nparticles, s.Nlevels)
basis(s::CavityMode) = FockBasis(s.Nmodes)
basis(s::MultimodeSystem) = CompositeBasis(basis(s.particles), [basis(mode) for mode=s.modes]...)

# Define interaction between cavity modes and particles

function An(scaling, n, p)
    M = [p.potential.A(scaling, n, i, j) for i=1:p.Nlevels, j=1:p.Nlevels]
    op = SparseOperator(GenericBasis(p.Nlevels), complex(M))
    nparticleoperator_1(basis(p), op)
end

function Bn(scaling, n, p)
    M = [p.potential.B(scaling, n, i, j) for i=1:p.Nlevels, j=1:p.Nlevels]
    op = SparseOperator(GenericBasis(p.Nlevels), complex(M))
    nparticleoperator_1(basis(p), op)
end

An(submode::Int, system::MultimodeSystem) = An(system.scaling, system.modes[submode].index, system.particles)
Bn(submode::Int, system::MultimodeSystem) = Bn(system.scaling, system.modes[submode].index, system.particles)

# Define Hamiltonian
function Hamiltonian(p::Particles)
    M = [p.potential.E(p.E0, j) for j=1:p.Nlevels]
    op = SparseOperator(GenericBasis(p.Nlevels), spdiagm(complex(M)))
    nparticleoperator_1(basis(p), op)
end

Hamiltonian(mode::CavityMode) = (create(basis(mode))*destroy(basis(mode))*(-mode.delta))

function Hamiltonian(system::MultimodeSystem)
    b = basis(system)
    H = embed(b, 1, Hamiltonian(system.particles))

    for i=1:length(system.modes)
        mode = system.modes[i]
        modebasis = basis(mode)
        An_ = An(system.scaling, mode.index, system.particles)
        Bn_ = Bn(system.scaling, mode.index, system.particles)
        H += embed(b, [1,i+1], [An_, number(modebasis)*mode.U0])
        H += embed(b, [1,i+1], [Bn_,(destroy(modebasis)+create(modebasis))*mode.eta])
        H += embed(b, i+1, Hamiltonian(mode))
    end
    return H
end

function Jump_operators(system::MultimodeSystem)
    b = basis(system)
    J = Operator[]
    for i=1:length(system.modes)
        mode = system.modes[i]
        a = destroy(basis(mode))
        j = embed(b, i+1, sqrt(2*mode.kappa)*a)
        push!(J, j)
    end
    return J
end

function particledensity(basis_functions::Vector, rho)
    N_positionpoints = length(basis_functions[1])
    n = zeros(Float64, N_positionpoints)
    for i=1:N_positionpoints
        f(idx_create::Int, idx_destroy::Int) = basis_functions[idx_create][i]*basis_functions[idx_destroy][i]
        n_op = weighted_cdaggeri_cj(basis(rho), f)
        n[i] = real(expect(n_op, rho))
    end
    return n
end

function particledensity(system::Particles, x::Vector{Float64}, rho::Operator)
    @assert basis(system)==rho.basis_l
    @assert basis(system)==rho.basis_r
    basis_functions = Vector{Float64}[]
    for i=1:system.Nlevels
        push!(basis_functions, system.potential.basis_function(i, x))
    end
    return particledensity(basis_functions, rho)
end

# end # module


text:
import Base: copy, copyto!, ==
import Random: seed!
import RandomNumbers: gen_seed, seed_type, unsafe_copyto!, unsafe_compare

@inline threefry_rotl(x::UInt64, N) = (x << (N & 63)) | (x >> ((64-N) & 63))
@inline threefry_rotl(x::UInt32, N) = (x << (N & 31)) | (x >> ((32-N) & 31))

@inline SKEIN_KS_PARITY(::Type{UInt32}) = 0x1BD11BDA
@inline SKEIN_KS_PARITY(::Type{UInt64}) = 0x1BD11BDAA9FC1A22

@inline Rx4_0_0(::Type{UInt64}) = 14
@inline Rx4_0_1(::Type{UInt64}) = 16
@inline Rx4_1_0(::Type{UInt64}) = 52
@inline Rx4_1_1(::Type{UInt64}) = 57
@inline Rx4_2_0(::Type{UInt64}) = 23
@inline Rx4_2_1(::Type{UInt64}) = 40
@inline Rx4_3_0(::Type{UInt64}) =  5
@inline Rx4_3_1(::Type{UInt64}) = 37
@inline Rx4_4_0(::Type{UInt64}) = 25
@inline Rx4_4_1(::Type{UInt64}) = 33
@inline Rx4_5_0(::Type{UInt64}) = 46
@inline Rx4_5_1(::Type{UInt64}) = 12
@inline Rx4_6_0(::Type{UInt64}) = 58
@inline Rx4_6_1(::Type{UInt64}) = 22
@inline Rx4_7_0(::Type{UInt64}) = 32
@inline Rx4_7_1(::Type{UInt64}) = 32

@inline Rx2_0_0(::Type{UInt64}) = 16
@inline Rx2_1_0(::Type{UInt64}) = 42
@inline Rx2_2_0(::Type{UInt64}) = 12
@inline Rx2_3_0(::Type{UInt64}) = 31
@inline Rx2_4_0(::Type{UInt64}) = 16
@inline Rx2_5_0(::Type{UInt64}) = 32
@inline Rx2_6_0(::Type{UInt64}) = 24
@inline Rx2_7_0(::Type{UInt64}) = 21

@inline Rx4_0_0(::Type{UInt32}) = 10
@inline Rx4_0_1(::Type{UInt32}) = 26
@inline Rx4_1_0(::Type{UInt32}) = 11
@inline Rx4_1_1(::Type{UInt32}) = 21
@inline Rx4_2_0(::Type{UInt32}) = 13
@inline Rx4_2_1(::Type{UInt32}) = 27
@inline Rx4_3_0(::Type{UInt32}) = 23
@inline Rx4_3_1(::Type{UInt32}) =  5
@inline Rx4_4_0(::Type{UInt32}) =  6
@inline Rx4_4_1(::Type{UInt32}) = 20
@inline Rx4_5_0(::Type{UInt32}) = 17
@inline Rx4_5_1(::Type{UInt32}) = 11
@inline Rx4_6_0(::Type{UInt32}) = 25
@inline Rx4_6_1(::Type{UInt32}) = 10
@inline Rx4_7_0(::Type{UInt32}) = 18
@inline Rx4_7_1(::Type{UInt32}) = 20
@inline Rx2_0_0(::Type{UInt32}) = 13
@inline Rx2_1_0(::Type{UInt32}) = 15
@inline Rx2_2_0(::Type{UInt32}) = 26
@inline Rx2_3_0(::Type{UInt32}) =  6
@inline Rx2_4_0(::Type{UInt32}) = 17
@inline Rx2_5_0(::Type{UInt32}) = 29
@inline Rx2_6_0(::Type{UInt32}) = 16
@inline Rx2_7_0(::Type{UInt32}) = 24

"""
```julia
Threefry2x{T, R} <: R123Generator2x{T}
Threefry2x([seed, R])
Threefry2x(T[, seed, R])
```

Threefry2x is one kind of Threefry Counter-Based RNGs. It generates two numbers at a time.

`T` is `UInt32` or `UInt64`(default).

`seed` is a `Tuple` of two `Integer`s which will both be automatically converted to `T`.

`R` denotes to the Rounds which must be at least 1 and no more than 32. With 20 rounds (by default), it has a
considerable safety margin over the minimum number of rounds with no known statistical flaws, but still has
excellent performance.
"""
mutable struct Threefry2x{T<:Union{UInt32, UInt64}, R} <: R123Generator2x{T}
    x1::T
    x2::T
    key1::T
    key2::T
    ctr1::T
    ctr2::T
    p::Int
end

function Threefry2x(::Type{T}=UInt64, seed::NTuple{2, Integer}=gen_seed(T, 2), R::Integer=20) where
        T <: Union{UInt32, UInt64}
    @assert 1 <= R <= 32
    r = Threefry2x{T, Int(R)}(0, 0, 0, 0, 0, 0, 0)
    seed!(r, seed)
end
Threefry2x(seed::NTuple{2, Integer}, R::Integer=20) = Threefry2x(UInt64, seed, R)

function seed!(r::Threefry2x{T}, seed::NTuple{2, Integer}=gen_seed(T, 2)) where T <: Union{UInt32, UInt64}
    r.x1 = r.x2 = 0
    r.key1 = seed[1] % T
    r.key2 = seed[2] % T
    r.ctr1 = r.ctr2 = 0
    random123_r(r)
    r
end

@inline seed_type(::Type{Threefry2x{T, R}}) where {T, R} = NTuple{2, T}

function copyto!(dest::Threefry2x{T, R}, src::Threefry2x{T, R}) where {T, R}
    unsafe_copyto!(dest, src, T, 6)
    dest.p = src.p
    dest
end

copy(src::Threefry2x{T, R}) where {T, R} = Threefry2x{T, R}(src.x1, src.x2, src.key1, src.key2,
                                                            src.ctr1, src.ctr2, src.p)

==(r1::Threefry2x{T, R}, r2::Threefry2x{T, R}) where {T, R} = unsafe_compare(r1, r2, T, 6) && r1.p == r2.p

@inline function random123_r(r::Threefry2x{T, R}) where {T <: Union{UInt32, UInt64}, R}
    ks2 = SKEIN_KS_PARITY(T)
    ks0 = r.key1
    x0 = r.ctr1
    ks2 ⊻= ks0
    ks1 = r.key2
    x1 = r.ctr2
    ks2 ⊻= ks1
    x0 += ks0
    x1 += ks1

    if R > 0 x0 += x1; x1 = threefry_rotl(x1, Rx2_0_0(T)); x1 ⊻= x0; end
    if R > 1 x0 += x1; x1 = threefry_rotl(x1, Rx2_1_0(T)); x1 ⊻= x0; end
    if R > 2 x0 += x1; x1 = threefry_rotl(x1, Rx2_2_0(T)); x1 ⊻= x0; end
    if R > 3 x0 += x1; x1 = threefry_rotl(x1, Rx2_3_0(T)); x1 ⊻= x0; end
    if R > 3
        x0 += ks1; x1 += ks2;
        x1 += 1 % T;
    end
    if R > 4 x0 += x1; x1 = threefry_rotl(x1, Rx2_4_0(T)); x1 ⊻= x0; end
    if R > 5 x0 += x1; x1 = threefry_rotl(x1, Rx2_5_0(T)); x1 ⊻= x0; end
    if R > 6 x0 += x1; x1 = threefry_rotl(x1, Rx2_6_0(T)); x1 ⊻= x0; end
    if R > 7 x0 += x1; x1 = threefry_rotl(x1, Rx2_7_0(T)); x1 ⊻= x0; end
    if R > 7
        x0 += ks2; x1 += ks0;
        x1 += 2 % T;
    end
    if R > 8 x0 += x1; x1 = threefry_rotl(x1, Rx2_0_0(T)); x1 ⊻= x0; end
    if R > 9 x0 += x1; x1 = threefry_rotl(x1, Rx2_1_0(T)); x1 ⊻= x0; end
    if R > 10 x0 += x1; x1 = threefry_rotl(x1, Rx2_2_0(T)); x1 ⊻= x0; end
    if R > 11 x0 += x1; x1 = threefry_rotl(x1, Rx2_3_0(T)); x1 ⊻= x0; end
    if R > 11
        x0 += ks0; x1 += ks1;
        x1 += 3 % T;
    end
    if R > 12 x0 += x1; x1 = threefry_rotl(x1, Rx2_4_0(T)); x1 ⊻= x0; end
    if R > 13 x0 += x1; x1 = threefry_rotl(x1, Rx2_5_0(T)); x1 ⊻= x0; end
    if R > 14 x0 += x1; x1 = threefry_rotl(x1, Rx2_6_0(T)); x1 ⊻= x0; end
    if R > 15 x0 += x1; x1 = threefry_rotl(x1, Rx2_7_0(T)); x1 ⊻= x0; end
    if R > 15
        x0 += ks1; x1 += ks2;
        x1 += 4 % T;
    end
    if R > 16 x0 += x1; x1 = threefry_rotl(x1, Rx2_0_0(T)); x1 ⊻= x0; end
    if R > 17 x0 += x1; x1 = threefry_rotl(x1, Rx2_1_0(T)); x1 ⊻= x0; end
    if R > 18 x0 += x1; x1 = threefry_rotl(x1, Rx2_2_0(T)); x1 ⊻= x0; end
    if R > 19 x0 += x1; x1 = threefry_rotl(x1, Rx2_3_0(T)); x1 ⊻= x0; end
    if R > 19
        x0 += ks2; x1 += ks0;
        x1 += 5 % T;
    end
    if R > 20 x0 += x1; x1 = threefry_rotl(x1, Rx2_4_0(T)); x1 ⊻= x0; end
    if R > 21 x0 += x1; x1 = threefry_rotl(x1, Rx2_5_0(T)); x1 ⊻= x0; end
    if R > 22 x0 += x1; x1 = threefry_rotl(x1, Rx2_6_0(T)); x1 ⊻= x0; end
    if R > 23 x0 += x1; x1 = threefry_rotl(x1, Rx2_7_0(T)); x1 ⊻= x0; end
    if R > 23
        x0 += ks0; x1 += ks1;
        x1 += 6 % T;
    end
    if R > 24 x0 += x1; x1 = threefry_rotl(x1, Rx2_0_0(T)); x1 ⊻= x0; end
    if R > 25 x0 += x1; x1 = threefry_rotl(x1, Rx2_1_0(T)); x1 ⊻= x0; end
    if R > 26 x0 += x1; x1 = threefry_rotl(x1, Rx2_2_0(T)); x1 ⊻= x0; end
    if R > 27 x0 += x1; x1 = threefry_rotl(x1, Rx2_3_0(T)); x1 ⊻= x0; end
    if R > 27
        x0 += ks1; x1 += ks2;
        x1 += 7 % T;
    end
    if R > 28 x0 += x1; x1 = threefry_rotl(x1, Rx2_4_0(T)); x1 ⊻= x0; end
    if R > 29 x0 += x1; x1 = threefry_rotl(x1, Rx2_5_0(T)); x1 ⊻= x0; end
    if R > 30 x0 += x1; x1 = threefry_rotl(x1, Rx2_6_0(T)); x1 ⊻= x0; end
    if R > 31 x0 += x1; x1 = threefry_rotl(x1, Rx2_7_0(T)); x1 ⊻= x0; end
    if R > 31
        x0 += ks2; x1 += ks0;
        x1 += 8 % T;
    end
    r.x1, r.x2 = x0, x1
end

"""
```julia
Threefry4x{T, R} <: R123Generator4x{T}
Threefry4x([seed, R])
Threefry4x(T[, seed, R])
```

Threefry2x is one kind of Threefry Counter-Based RNGs. It generates four numbers at a time.

`T` is `UInt32` or `UInt64`(default).

`seed` is a `Tuple` of four `Integer`s which will all be automatically converted to `T`.

`R` denotes to the Rounds which must be at least 1 and no more than 32. With 20 rounds (by default), it has a
considerable safety margin over the minimum number of rounds with no known statistical flaws, but still has
excellent performance.
"""
mutable struct Threefry4x{T<:Union{UInt32, UInt64}, R} <: R123Generator4x{T}
    x1::T
    x2::T
    x3::T
    x4::T
    key1::T
    key2::T
    key3::T
    key4::T
    ctr1::T
    ctr2::T
    ctr3::T
    ctr4::T
    p::Int
end

function Threefry4x(::Type{T}=UInt64, seed::NTuple{4, Integer}=gen_seed(T, 4), R::Integer=20) where
        T <: Union{UInt32, UInt64}
    @assert 1 <= R <= 72
    r = Threefry4x{T, Int(R)}(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
    seed!(r, seed)
end
Threefry4x(seed::NTuple{4, Integer}, R::Integer=20) = Threefry4x(UInt64, seed, R)

function seed!(r::Threefry4x{T}, seed::NTuple{4, Integer}=gen_seed(T, 4)) where T <: Union{UInt32, UInt64}
    r.x1 = r.x2 = r.x3 = r.x4 = 0
    r.key1, r.key2, r.key3, r.key4 = seed[1] % T, seed[2] % T, seed[3] % T, seed[4] % T
    r.ctr1 = r.ctr2 = r.ctr3 = r.ctr4 = 0
    random123_r(r)
    r
end

@inline seed_type(::Type{Threefry4x{T, R}}) where {T, R} = NTuple{4, T}

function copyto!(dest::Threefry4x{T, R}, src::Threefry4x{T, R}) where {T, R}
    unsafe_copyto!(dest, src, T, 12)
    dest.p = src.p
    dest
end

copy(src::Threefry4x{T, R}) where {T, R} = Threefry4x{T, R}(src.x1, src.x2, src.x3, src.x4,
                                                     src.key1, src.key2, src.key3, src.key4,
                                                     src.ctr1, src.ctr2, src.ctr3, src.ctr4, src.p)

==(r1::Threefry4x{T, R}, r2::Threefry4x{T, R}) where {T, R} = unsafe_compare(r1, r2, T, 12) && r1.p == r2.p

@inline function random123_r(r::Threefry4x{T, R}) where {T <: Union{UInt32, UInt64}, R}
    ks4 = SKEIN_KS_PARITY(T)
    ks0 = r.key1
    x0 = r.ctr1
    ks4 ⊻= ks0
    ks1 = r.key2
    x1 = r.ctr2
    ks4 ⊻= ks1
    ks2 = r.key3
    x2 = r.ctr3
    ks4 ⊻= ks2
    ks3 = r.key4
    x3 = r.ctr4
    ks4 ⊻= ks3
    x0 += ks0; x1 += ks1; x2 += ks2; x3 += ks3;

    if R > 0
        x0 += x1; x1 = threefry_rotl(x1, Rx4_0_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_0_1(T)); x3 ⊻= x2;
    end
    if R > 1
        x0 += x3; x3 = threefry_rotl(x3, Rx4_1_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_1_1(T)); x1 ⊻= x2;
    end
    if R > 2
        x0 += x1; x1 = threefry_rotl(x1, Rx4_2_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_2_1(T)); x3 ⊻= x2;
    end
    if R > 3
        x0 += x3; x3 = threefry_rotl(x3, Rx4_3_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_3_1(T)); x1 ⊻= x2;
    end
    if R > 3
        x0 += ks1; x1 += ks2; x2 += ks3; x3 += ks4;
        x3 += 1 % T;
    end
    if R > 4
        x0 += x1; x1 = threefry_rotl(x1, Rx4_4_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_4_1(T)); x3 ⊻= x2;
    end
    if R > 5
        x0 += x3; x3 = threefry_rotl(x3, Rx4_5_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_5_1(T)); x1 ⊻= x2;
    end
    if R > 6
        x0 += x1; x1 = threefry_rotl(x1, Rx4_6_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_6_1(T)); x3 ⊻= x2;
    end
    if R > 7
        x0 += x3; x3 = threefry_rotl(x3, Rx4_7_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_7_1(T)); x1 ⊻= x2;
    end
    if R > 7
        x0 += ks2; x1 += ks3; x2 += ks4; x3 += ks0;
        x3 += 2 % T;
    end
    if R > 8
        x0 += x1; x1 = threefry_rotl(x1, Rx4_0_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_0_1(T)); x3 ⊻= x2;
    end
    if R > 9
        x0 += x3; x3 = threefry_rotl(x3, Rx4_1_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_1_1(T)); x1 ⊻= x2;
    end
    if R > 10
        x0 += x1; x1 = threefry_rotl(x1, Rx4_2_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_2_1(T)); x3 ⊻= x2;
    end
    if R > 11
        x0 += x3; x3 = threefry_rotl(x3, Rx4_3_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_3_1(T)); x1 ⊻= x2;
    end
    if R > 11
        x0 += ks3; x1 += ks4; x2 += ks0; x3 += ks1;
        x3 += 3 % T;
    end
    if R > 12
        x0 += x1; x1 = threefry_rotl(x1, Rx4_4_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_4_1(T)); x3 ⊻= x2;
    end
    if R > 13
        x0 += x3; x3 = threefry_rotl(x3, Rx4_5_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_5_1(T)); x1 ⊻= x2;
    end
    if R > 14
        x0 += x1; x1 = threefry_rotl(x1, Rx4_6_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_6_1(T)); x3 ⊻= x2;
    end
    if R > 15
        x0 += x3; x3 = threefry_rotl(x3, Rx4_7_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_7_1(T)); x1 ⊻= x2;
    end
    if R > 15
        x0 += ks4; x1 += ks0; x2 += ks1; x3 += ks2;
        x3 += 4 % T;
    end
    if R > 16
        x0 += x1; x1 = threefry_rotl(x1, Rx4_0_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_0_1(T)); x3 ⊻= x2;
    end
    if R > 17
        x0 += x3; x3 = threefry_rotl(x3, Rx4_1_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_1_1(T)); x1 ⊻= x2;
    end
    if R > 18
        x0 += x1; x1 = threefry_rotl(x1, Rx4_2_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_2_1(T)); x3 ⊻= x2;
    end
    if R > 19
        x0 += x3; x3 = threefry_rotl(x3, Rx4_3_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_3_1(T)); x1 ⊻= x2;
    end
    if R > 19
        x0 += ks0; x1 += ks1; x2 += ks2; x3 += ks3;
        x3 += 5 % T;
    end
    if R > 20
        x0 += x1; x1 = threefry_rotl(x1, Rx4_4_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_4_1(T)); x3 ⊻= x2;
    end
    if R > 21
        x0 += x3; x3 = threefry_rotl(x3, Rx4_5_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_5_1(T)); x1 ⊻= x2;
    end
    if R > 22
        x0 += x1; x1 = threefry_rotl(x1, Rx4_6_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_6_1(T)); x3 ⊻= x2;
    end
    if R > 23
        x0 += x3; x3 = threefry_rotl(x3, Rx4_7_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_7_1(T)); x1 ⊻= x2;
    end
    if R > 23
        x0 += ks1; x1 += ks2; x2 += ks3; x3 += ks4;
        x3 += 6 % T;
    end
    if R > 24
        x0 += x1; x1 = threefry_rotl(x1, Rx4_0_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_0_1(T)); x3 ⊻= x2;
    end
    if R > 25
        x0 += x3; x3 = threefry_rotl(x3, Rx4_1_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_1_1(T)); x1 ⊻= x2;
    end
    if R > 26
        x0 += x1; x1 = threefry_rotl(x1, Rx4_2_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_2_1(T)); x3 ⊻= x2;
    end
    if R > 27
        x0 += x3; x3 = threefry_rotl(x3, Rx4_3_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_3_1(T)); x1 ⊻= x2;
    end
    if R > 27
        x0 += ks2; x1 += ks3; x2 += ks4; x3 += ks0;
        x3 += 7 % T;
    end
    if R > 28
        x0 += x1; x1 = threefry_rotl(x1, Rx4_4_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_4_1(T)); x3 ⊻= x2;
    end
    if R > 29
        x0 += x3; x3 = threefry_rotl(x3, Rx4_5_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_5_1(T)); x1 ⊻= x2;
    end
    if R > 30
        x0 += x1; x1 = threefry_rotl(x1, Rx4_6_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_6_1(T)); x3 ⊻= x2;
    end
    if R > 31
        x0 += x3; x3 = threefry_rotl(x3, Rx4_7_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_7_1(T)); x1 ⊻= x2;
    end
    if R > 31
        x0 += ks3; x1 += ks4; x2 += ks0; x3 += ks1;
        x3 += 8 % T;
    end
    if R > 32
        x0 += x1; x1 = threefry_rotl(x1, Rx4_0_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_0_1(T)); x3 ⊻= x2;
    end
    if R > 33
        x0 += x3; x3 = threefry_rotl(x3, Rx4_1_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_1_1(T)); x1 ⊻= x2;
    end
    if R > 34
        x0 += x1; x1 = threefry_rotl(x1, Rx4_2_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_2_1(T)); x3 ⊻= x2;
    end
    if R > 35
        x0 += x3; x3 = threefry_rotl(x3, Rx4_3_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_3_1(T)); x1 ⊻= x2;
    end
    if R > 35
        x0 += ks4; x1 += ks0; x2 += ks1; x3 += ks2;
        x3 += 9 % T;
    end
    if R > 36
        x0 += x1; x1 = threefry_rotl(x1, Rx4_4_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_4_1(T)); x3 ⊻= x2;
    end
    if R > 37
        x0 += x3; x3 = threefry_rotl(x3, Rx4_5_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_5_1(T)); x1 ⊻= x2;
    end
    if R > 38
        x0 += x1; x1 = threefry_rotl(x1, Rx4_6_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_6_1(T)); x3 ⊻= x2;
    end
    if R > 39
        x0 += x3; x3 = threefry_rotl(x3, Rx4_7_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_7_1(T)); x1 ⊻= x2;
    end
    if R > 39
        x0 += ks0; x1 += ks1; x2 += ks2; x3 += ks3;
        x3 += 10 % T;
    end
    if R > 40
        x0 += x1; x1 = threefry_rotl(x1, Rx4_0_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_0_1(T)); x3 ⊻= x2;
    end
    if R > 41
        x0 += x3; x3 = threefry_rotl(x3, Rx4_1_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_1_1(T)); x1 ⊻= x2;
    end
    if R > 42
        x0 += x1; x1 = threefry_rotl(x1, Rx4_2_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_2_1(T)); x3 ⊻= x2;
    end
    if R > 43
        x0 += x3; x3 = threefry_rotl(x3, Rx4_3_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_3_1(T)); x1 ⊻= x2;
    end
    if R > 43
        x0 += ks1; x1 += ks2; x2 += ks3; x3 += ks4;
        x3 += 11 % T;
    end
    if R > 44
        x0 += x1; x1 = threefry_rotl(x1, Rx4_4_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_4_1(T)); x3 ⊻= x2;
    end
    if R > 45
        x0 += x3; x3 = threefry_rotl(x3, Rx4_5_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_5_1(T)); x1 ⊻= x2;
    end
    if R > 46
        x0 += x1; x1 = threefry_rotl(x1, Rx4_6_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_6_1(T)); x3 ⊻= x2;
    end
    if R > 47
        x0 += x3; x3 = threefry_rotl(x3, Rx4_7_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_7_1(T)); x1 ⊻= x2;
    end
    if R > 47
        x0 += ks2; x1 += ks3; x2 += ks4; x3 += ks0;
        x3 += 12 % T;
    end
    if R > 48
        x0 += x1; x1 = threefry_rotl(x1, Rx4_0_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_0_1(T)); x3 ⊻= x2;
    end
    if R > 49
        x0 += x3; x3 = threefry_rotl(x3, Rx4_1_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_1_1(T)); x1 ⊻= x2;
    end
    if R > 50
        x0 += x1; x1 = threefry_rotl(x1, Rx4_2_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_2_1(T)); x3 ⊻= x2;
    end
    if R > 51
        x0 += x3; x3 = threefry_rotl(x3, Rx4_3_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_3_1(T)); x1 ⊻= x2;
    end
    if R > 51
        x0 += ks3; x1 += ks4; x2 += ks0; x3 += ks1;
        x3 += 13 % T;
    end
    if R > 52
        x0 += x1; x1 = threefry_rotl(x1, Rx4_4_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_4_1(T)); x3 ⊻= x2;
    end
    if R > 53
        x0 += x3; x3 = threefry_rotl(x3, Rx4_5_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_5_1(T)); x1 ⊻= x2;
    end
    if R > 54
        x0 += x1; x1 = threefry_rotl(x1, Rx4_6_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_6_1(T)); x3 ⊻= x2;
    end
    if R > 55
        x0 += x3; x3 = threefry_rotl(x3, Rx4_7_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_7_1(T)); x1 ⊻= x2;
    end
    if R > 55
        x0 += ks4; x1 += ks0; x2 += ks1; x3 += ks2;
        x3 += 14 % T;
    end
    if R > 56
        x0 += x1; x1 = threefry_rotl(x1, Rx4_0_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_0_1(T)); x3 ⊻= x2;
    end
    if R > 57
        x0 += x3; x3 = threefry_rotl(x3, Rx4_1_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_1_1(T)); x1 ⊻= x2;
    end
    if R > 58
        x0 += x1; x1 = threefry_rotl(x1, Rx4_2_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_2_1(T)); x3 ⊻= x2;
    end
    if R > 59
        x0 += x3; x3 = threefry_rotl(x3, Rx4_3_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_3_1(T)); x1 ⊻= x2;
    end
    if R > 59
        x0 += ks0; x1 += ks1; x2 += ks2; x3 += ks3;
        x3 += 15 % T;
    end
    if R > 60
        x0 += x1; x1 = threefry_rotl(x1, Rx4_4_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_4_1(T)); x3 ⊻= x2;
    end
    if R > 61
        x0 += x3; x3 = threefry_rotl(x3, Rx4_5_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_5_1(T)); x1 ⊻= x2;
    end
    if R > 62
        x0 += x1; x1 = threefry_rotl(x1, Rx4_6_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_6_1(T)); x3 ⊻= x2;
    end
    if R > 63
        x0 += x3; x3 = threefry_rotl(x3, Rx4_7_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_7_1(T)); x1 ⊻= x2;
    end
    if R > 63
        x0 += ks1; x1 += ks2; x2 += ks3; x3 += ks4;
        x3 += 16 % T;
    end
    if R > 64
        x0 += x1; x1 = threefry_rotl(x1, Rx4_0_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_0_1(T)); x3 ⊻= x2;
    end
    if R > 65
        x0 += x3; x3 = threefry_rotl(x3, Rx4_1_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_1_1(T)); x1 ⊻= x2;
    end
    if R > 66
        x0 += x1; x1 = threefry_rotl(x1, Rx4_2_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_2_1(T)); x3 ⊻= x2;
    end
    if R > 67
        x0 += x3; x3 = threefry_rotl(x3, Rx4_3_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_3_1(T)); x1 ⊻= x2;
    end
    if R > 67
        x0 += ks2; x1 += ks3; x2 += ks4; x3 += ks0;
        x3 += 17 % T;
    end
    if R > 68
        x0 += x1; x1 = threefry_rotl(x1, Rx4_4_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_4_1(T)); x3 ⊻= x2;
    end
    if R > 69
        x0 += x3; x3 = threefry_rotl(x3, Rx4_5_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_5_1(T)); x1 ⊻= x2;
    end
    if R > 70
        x0 += x1; x1 = threefry_rotl(x1, Rx4_6_0(T)); x1 ⊻= x0;
        x2 += x3; x3 = threefry_rotl(x3, Rx4_6_1(T)); x3 ⊻= x2;
    end
    if R > 71
        x0 += x3; x3 = threefry_rotl(x3, Rx4_7_0(T)); x3 ⊻= x0;
        x2 += x1; x1 = threefry_rotl(x1, Rx4_7_1(T)); x1 ⊻= x2;
    end
    if R > 71
        x0 += ks3; x1 += ks4; x2 += ks0; x3 += ks1;
        x3 += 18 % T;
    end
    r.x1, r.x2, r.x3, r.x4 = x0, x1, x2, x3
end


text:
@testset "$M{$K,$T}" for M in [RefStdMap, GCStdMap, STL.SharedStdMap],
                         K in sort!(collect(STL.StdMap_keys); by=string),
                         T in sort!(collect(STL.StdMap_types); by=string)
    # These segfault; converting SharedStdString to Ptr{StdString} doesn't work
    T ≡ RefStdVector{STL.SharedStdString} && continue
    T ≡ STL.SharedStdVector{STL.SharedStdString} && continue
    K ≡ STL.SharedStdString && continue
    # @show M K T

    # TODO: Need to free the generatedp RefStdString objects
    mkK(i) = K ≡ Bool ? i % Bool : K <: Ptr ? K(i) : K ≡ RefStdString ? K(string(i)) : K ≡ STL.SharedStdString ? K(string(i)) : i
    function mkT(i)
        return T ≡ Bool ? i % Bool :
               T <: Ptr ? T(i) :
               T ≡ RefStdString ? T(string(i)) :
               T ≡ STL.SharedStdString ? T(string(i)) :
               T ≡ RefStdVector{RefStdString} ? (x = T(1); x[0] = RefStdString(string(i)); x) :
               T ≡ RefStdVector{STL.SharedStdString} ? (x = T(1); x[0] = STL.SharedStdString(string(i)); x) :
               T ≡ STL.SharedStdVector{RefStdString} ? (x = T(1); x[0] = RefStdString(string(i)); x) :
               T ≡ STL.SharedStdVector{STL.SharedStdString} ? (x = T(1); x[0] = STL.SharedStdString(string(i)); x) : i
    end

    # @show 1
    map = M{K,T}()
    @test keytype(map) ≡ K
    @test eltype(map) ≡ T

    # @show 2
    @test length(map) == 0
    @test isempty(map)

    # @show 3
    N = K ≡ Bool ? 2 : 10
    # @show 31
    for i in 0:(N - 1)
        map[mkK(i)] = mkT(i)
    end
    # @show 32
    @test length(map) == N
    # @show 33
    for i in 0:(N - 1)
        # @show 331
        @test haskey(map, mkK(i))
        # @show 332
        k = mkK(i)
        # @show 3321
        t = mkT(i)
        # @show 3322
        q = map[mkK(i)]
        # @show 3323
        q == mkT(i)
        @test map[mkK(i)] == mkT(i)
        # @show 333
        if K ≢ Bool
            # @show 334
            @test !haskey(map, mkK(i + 20))
        end
        # @show 335
    end

    # @show 4
    elts = Dict([elt for elt in map])
    @test typeof(elts) == Dict{K,T}
    @test elts == Dict([mkK(i) => mkT(i) for i in 0:(N - 1)])

    # @show 5
    for i in 0:(N - 1)
        @test haskey(map, mkK(i))
        delete!(map, mkK(i))
        @test !haskey(map, mkK(i))
    end

    # @show 6
    @test isempty(map)

    # @show 7
    if M == RefStdMap
        free(map)
    end
    GC.gc(true)
    # @show 8
end


text:
#using Random
# using Interpolations

U_func_time_trivial(t) = 0.0;
U_func_time_trivial(t, u) = zeros(size(u));
U_func_space_trivial(u) = zeros(size(u));

#####
##### Spike-like controllers
#####
# FROM: https://discourse.julialang.org/t/diffeqflux-with-time-as-additional-input-to-neural-ode/26456/3
function U_func_time(t, u, U_width, U_starts;
                        F_dim=2, amplitude=1.0)
    F = zeros(size(u))
    if any(i -> i < U_width, abs.(U_starts .- t))
        F[F_dim] = amplitude
    end
    return F
end

function U_func_time_multivariate(t, u, U_widths, U_starts;
                                F_dim=[2], amplitudes=nothing)
    if amplitudes == nothing
        amplitudes = ones(size(F_dim))
    end
    F = zeros(size(u))
    # Each forced dimension can have different timing
    for i_dim in 1:length(F_dim)
        if any(i -> i < U_widths[i_dim], abs.(U_starts[i_dim,:] .- t))
            F[i_dim] = amplitudes[i_dim]
        end
    end
    return F
end

function spring_forcing_example(tspan)
    Random.seed!(11)
    num_ctr = 2;
    U_starts = rand(num_ctr,1) .* tspan[2]
    U_width = 0.05;
    U_func(t, u) = U_func_time(t, u, U_width, U_starts)

    return U_func
end

#####
##### Piecewise constant controllers
#####

# Returns the function with signature:
#   f = generate_U_func_piecewise(t, u)
#   f(t, u) # control at that time
function generate_U_func_piecewise(tspan, u_size;
                                    num_jumps = 10,
                                    U_min=-1, U_max=1,
                                    F_dim=2)
    jump_locations = sort(rand(num_jumps) .* tspan[end])
    jump_locations = [jump_locations; tspan[end]]
    U_series = zeros(u_size)
    my_rand() = (U_max-U_min)*rand() + U_min
    let jump_ind=1, current_val=my_rand()
        for i in 1:length(tspan)
            if tspan[i] > jump_locations[jump_ind]
                jump_ind += 1
                current_val = my_rand()
            end
            U_series[F_dim, i] = current_val
        end
    end
    # Generates a function that can be called like f(t)
    # return U_series
    U_func = generate_map_rows_function(
                (x)->CubicSplineInterpolation(tspan, x), U_series)
    return (t, tmp)->U_func(t) # Needs the right signature
end

#####
##### Coordinate (spatial) controllers
#####
function U_func_hard_wall(X, x_dim=1, v_dim=2, wall_strength=100.0)
    F = zeros(size(X))
    X[x_dim]<0.0 ? F[v_dim]=wall_strength : return F
    # F = [wall_strength*σ(-5.0(x+0.0)), 0.0, 0.0]
    return F
end

function U_func_soft_wall(X, x_dim=1, v_dim=2, wall_strength=100.0)
    F = zeros(size(X))
    F[y_dim] = wall_strength*σ(-10.0*X[x_dim])
    return F
end

function U_func_spring(X; k=1000.0, r=0.1, damping=1.0,
                       r_dim=1, F_dim=2, damp_dim=2)
    F = zeros(size(X))
    # Damped spring at the 0 coordinate
    if r - X[r_dim] > 0.0
        Δ = r - X[r_dim]
        F[F_dim] = k*Δ + -damping*abs(X[damp_dim])
    end
    return F
end

function U_func_spring_sphere(X; k=1000.0,
                              r=0.1, location=nothing, damping=1.0,
                              F_dim=2, damp_dim=2)
    F = zeros(size(X))
    if location == nothing
        # Default is the origin
        location = zeros(size(X))
    end
    Δ = sum(abs2,location .- X)
    if Δ < r
        F[F_dim] = k*Δ + -damping*abs(X[damp_dim])
    end
    return F
end

#####
##### Export
#####
export U_func_time, U_func_time_multivariate,
        spring_forcing_example, U_func_spring_sphere,
        U_func_hard_wall, U_func_soft_wall, U_func_spring,
        U_func_time_trivial, U_func_space_trivial,
        generate_U_func_piecewise


text:
import Oceananigans.Fields: Field, location

function Field((LX, LY, LZ)::Tuple, grid::DistributedGrid, data, old_bcs, indices::Tuple, op, status)
    arch = architecture(grid)
    new_bcs = inject_halo_communication_boundary_conditions(old_bcs, arch.local_rank, arch.connectivity)
    return Field{LX, LY, LZ}(grid, data, new_bcs, indices, op, status)
end

const DistributedField      = Field{<:Any, <:Any, <:Any, <:Any, <:DistributedGrid}
const DistributedFieldTuple = NamedTuple{S, <:NTuple{N, DistributedField}} where {S, N}


text:
module JuliaSets

export escape_time, JuliaSet, call

type JuliaSet
   Q::Function
end

function call{T<:Real}(J::JuliaSet, xs::Array{T}, ys::Array{T}, n_iter::Int64=100, escape_tol::T=100.0)
  @assert size(xs) == size(ys)
  ans = xs + 1im * ys
  for i = 1:n_iter
    ans = J.Q(ans)
  end
  return abs(ans) .< T(escape_tol)
end

function escape_time{T <: Real}(Q::Function, xs::Array{T,1}, ys::Array{T,1}, n_iter=100, escape_tol=100)
  A = zeros(length(xs), length(ys))
  for i in 1:length(xs)
    for j in 1:length(ys)
      z = xs[i] + 1im * ys[j]
      for k = 1:n_iter
        z = Q(z)
        if abs(z) > escape_tol # diverged
          A[i,j] = k
          break
        end
      end
      if A[i,j] == 0 # didn't diverge
        A[i,j] = escape_tol + 1
      end
    end
  end
  return A
end

end # module


text: return (42 + true)

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using Optim
using PrettyTables
using LinearAlgebra
using Tables

struct MultinomialLogitModel <: LogitModel
    coefnames::Vector{Symbol}
    coefs::Vector{Float64}
    vcov::Matrix{Float64}
    init_ll::Float64
    final_ll::Float64
    # TODO log likelihood at constants
end

extract_val(::Val{T}) where T = T
# why does this get wrapped in an extra type?
extract_val(::Type{Val{T}}) where T = T

function extract_namedtuple_bool(nt, ::Val{key})::Bool where key
    nt[key]::Bool
end

#=
Much work has gone into optimizing this to have zero allocations. Key optimizations:
- chosen_col is passed as a Val, so the column name is dispatched on. This means that the compiler knows which
  column will be used for chosen at compile time, and since it knows column types (from the type of NamedTuple row)
  can infer that chosen will always be an Int64 and avoid allocations
- Similarly, a barrier function with a val argument is used to prevent type instability when reading availability columns.
- params is received as a Vector{T} not an AbstractVector{T} - for some reason this saves two allocations
- FunctionWrappers are used to indicate to the compiler that ufunc will always return the same type (tested, the wrapper is necessary)

- Note that this was tested from a script calling multinomial_logit not inside a function, some of these optimizations may not be
  necessary inside a function.
=#
function mnl_ll_row(row, params::Vector{T}, utility_functions, ::Val{chosen_col}, avail_cols)::T where {T <: Number, chosen_col}
    util_sum = zero(T)

    chosen = row[chosen_col]
    local chosen_exputil::T
    for (choiceidx, ufunc) in enumerate(utility_functions)
        exp_util = if isnothing(avail_cols) || extract_namedtuple_bool(row, @inbounds Val(avail_cols[choiceidx]))
            # choice is available, either implicitly or explicitly
            exp(ufunc(params, row))
        else
            zero(T)
        end

        if choiceidx == chosen
            chosen_exputil = exp_util
        end
        util_sum += exp_util
    end

    log(chosen_exputil / util_sum)
end

function multinomial_logit_log_likelihood(utility_functions, chosen_col, avail_cols, data, parameters::Vector{T}) where T
    @debug "objective called with params" parameters
    R = rowtype(data)
    U = typeof(utility_functions)
    C = typeof(chosen_col)
    A = typeof(avail_cols)
    rowwise_loglik(
        FunctionWrapper{T, Tuple{R, Vector{T}, U, C, A}}(mnl_ll_row),
        data, parameters, utility_functions, chosen_col, avail_cols)
end

function multinomial_logit(
    utility,
    chosen,
    data;
    availability::Union{Nothing, AbstractVector{<:Pair{<:Any, <:Any}}}=nothing,
    method=BFGS(),
    se=true,
    verbose=:no,
    iterations=1_000
    )

    # if data isa JuliaDB.AbstractIndexedTable
    #     check_perfect_prediction(data, chosen, [utility.columnnames...])
    # end

    data, choice_col, avail_cols = prepare_data(data, chosen, utility.alt_numbers, availability)
    row_type = rowtype(data)
    obj(p::AbstractVector{T}) where T = -multinomial_logit_log_likelihood(
        FunctionWrapper{T, Tuple{Vector{T}, row_type}}.(utility.utility_functions),
        Val(choice_col), avail_cols, data, p)::T
    init_ll = -obj(utility.starting_values)

    @info "Log-likelihood at starting values $(init_ll)"

    results = optimize(
        TwiceDifferentiable(obj, utility.starting_values, autodiff=:forward),
        copy(utility.starting_values),
        method,
        Optim.Options(show_trace=verbose == :medium || verbose == :high, extended_trace=verbose==:high, iterations=iterations)
    )

    if !Optim.converged(results)
        @error "Failed to converge!"
        #throw(ConvergenceException(Optim.iterations(results)))
    else
        @info "Optimization converged successfully after $(Optim.iterations(results)) iterations"
    end

    @info """
    Using method $(Optim.summary(results)),
    $(Optim.f_calls(results)) function evaluations, $(Optim.g_calls(results)) gradient evaluations
    """

    final_ll = -Optim.minimum(results)
    params = Optim.minimizer(results)

    # put any fixed parameters back into the data
    final_coefnames = [utility.coefnames..., keys(utility.fixed_coefs)...]
    final_coefs = [params..., values(utility.fixed_coefs)...]

    if se
        @info "Calculating and inverting Hessian"

        # compute standard errors
        hess = ForwardDiff.hessian(obj, params)
        local inv_hess
        try
            inv_hess = inv(hess)
        catch e
            !(e isa LinearAlgebra.SingularException) && rethrow()
            @error "Hessian is singular. Not reporting standard errors, and you should probably be suspicious of point estimates."
            se = false
        end

        if se
            vcov = similar(inv_hess, length(final_coefs), length(final_coefs))
            vcov[:, :] .= convert(eltype(vcov), NaN)
            vcov[1:length(params), 1:length(params)] = inv_hess
        end
    end

    if !se
        vcov = fill(NaN, length(final_coefs), length(final_coefs))
    end

    return MultinomialLogitModel(final_coefnames, final_coefs, vcov, init_ll, final_ll)
end

function Base.summary(res::MultinomialLogitModel)
    mcfadden = 1 - res.final_ll / res.init_ll
    header = """
Multinomial logit model
Initial log-likelhood (at starting values): $(res.init_ll)
Final log-likelihood: $(res.final_ll)
McFadden's pseudo-R2 (relative to starting values): $mcfadden
"""

    vc = vcov(res)
    # nan variance in fixed params
    if !all(filter(x->!isnan(x), diag(vc)) .≥ 0)
        @error "Some estimated variances are negative, not showing std. errors!"
        ses = diag(vc)
        pval = fill(NaN, length(ses))
        selab = "Var"
    else
        ses = sqrt.(diag(vc))
        pval = coef(res) ./ ses
        selab = "Std. Err."
    end

    data = hcat(
        coefnames(res),
        coef(res),
        ses,
        pval
    )

    table = pretty_table(String, data, header=["", "Coef", selab, "Z-stat"],
        header_crayon=crayon"yellow bold", formatters=ft_printf("%.5f", 2:4))

    return header * table
end

multinomial_logit(NamedTuple) = error("Not enough arguments. Make sure arguments to @utility are enclosed in parens")


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

# Some image manipulation functions

module Image

export downSample, interpolateFlow, warpMotion, readImage, writeFlo

# Down sample an image to a smaller size
#   We can think of the image as dividing up into a grid of equal sized pixels
#   The pixels of the desired smaller image overlap the pixels of the larger original image
#   We want the value of the output pixel to be the average of the input pixels overlapped weighted by
#   area of overlap
# Pre: size(oi)>=(nw,nh)
function downSample(oi::Matrix{Float32}, nw::Int, nh::Int)
# Returns Matrix{Float32} size=(nh,nw)
  (oh, ow) = size(oi)
  if ow==nw && oh==nh
    return oi
  end
  fx = Float32(ow)/Float32(nw)   # width of old pixels per new pixel
  fy = Float32(oh)/Float32(nh)   # height of old pixels per new pixel
  ni = zeros(Float32, nh, nw)
  for y in 1:nh
    for x in 1:nw
      r = 0.0f0
      fi = Float32(y-1)*fy     # start vertical in old pixels of new pixel
      fli = Float32(y)*fy      # end vertical in old pixels of new pixel
      li = round(Int, floor(fli))+1  # end vertical index in old pixels of new pixel
      fj = Float32(x-1)*fx     # start horozontal in old pixels of new pixel
      flj = Float32(x)*fx      # end horozontal in old pixels of new pixel
      lj = round(Int, floor(flj))+1  # end horozontal index in old pixels of new pixel
      for i in round(Int, floor(fi))+1:min(li, oh)
        for j in round(Int, floor(fj))+1:min(lj, ow)
          c = 1.0f0  # weight each full overlapped old pixel by one, adjusted below
          if i-fi<1  # adjust for start vertical overlap
            c = i-fi
          end
          if j-fj<1  # adjust for start horozontal overlap
            c *= j-fj
          end
          if i==li   # adjust for end vertical overlap
            c *= fli-li+1
          end
          if j==lj   # adjust for end horozontal overlap
            c *= flj-lj+1
          end
          r += c*oi[i,j]
        end
      end
      ni[y,x] = r/(fx*fy)  # adjust for total area of old pixels overlapped by new pixel
    end
  end
  return ni
end

# Linearly interpolate between array values
# Pre: 1<=x<=size(a,2) && 1<=y<=size(a,1) && 1<=j<=size(a,3)
function interpolate(a::Array{Float32, 3}, x::Float32, y::Float32, j)
# Returns Float32
  xx = round(Int, floor(x))
  yy = round(Int, floor(y))
  alpha = x-xx
  beta = y-yy
  if alpha>0.0f0
    if beta>0.0f0
      return        beta *(alpha*a[yy+1,xx+1,j]+(1.0f0-alpha)*a[yy+1,xx,j])+
             (1.0f0-beta)*(alpha*a[yy,  xx+1,j]+(1.0f0-alpha)*a[yy,  xx,j])
    else
      return alpha*a[yy,xx+1,j]+(1.0f0-alpha)*a[yy,xx,j]
    end
  else
    if beta>0.0f0
      return beta*a[yy+1,xx,j]+(1.0f0-beta)*a[yy,xx,j]
    else
      return a[yy,xx,j]
    end
  end
end

# Linearly interpolate between array values
# Pre: 1<=x<=size(a,2) && 1<=y<=size(a,1)
function interpolate(a::Matrix{Float32}, x::Float32, y::Float32)
# Returns Float32
  xx = round(Int, floor(x))
  yy = round(Int, floor(y))
  alpha = x-xx
  beta = y-yy
  if alpha>0.0f0
    if beta>0.0f0
      return        beta *(alpha*a[yy+1,xx+1]+(1.0f0-alpha)*a[yy+1,xx])+
             (1.0f0-beta)*(alpha*a[yy,  xx+1]+(1.0f0-alpha)*a[yy,  xx])
    else
      return alpha*a[yy,xx+1]+(1.0f0-alpha)*a[yy,xx]
    end
  else
    if beta>0.0f0
      return beta*a[yy+1,xx]+(1.0f0-beta)*a[yy,xx]
    else
      return a[yy,xx]
    end
  end
end

# Given a flow at a coarser grain, interpolate it to a finer grain
# Pre: size(ou,1)<=nh && size(ou,2)<=nw && size(ou,3)=2
function interpolateFlow(ou::Array{Float32, 3}, nw::Int, nh::Int)
# Returns Array{Float32, 3} size=(nh,nw,2)
  oh = size(ou, 1)
  ow = size(ou, 2)
  if ow==nw && oh==nh
    return ou
  end
  nu = zeros(Float32, nh, nw, 2)
  s = Float32(Float64(nw)/Float64(ow))
  for y in 1:nh
    for x in 1:nw
      nx = Float32(x-1)/Float32(nw)*Float32(ow)+1
      ny = Float32(y-1)/Float32(nh)*Float32(oh)+1
      if nx>=ow
        nx = Float32(ow)
      end
      if ny>=oh
        ny = Float32(oh)
      end
      nu[y,x,1] = interpolate(ou, nx, ny, 1)/s
      nu[y,x,2] = interpolate(ou, nx, ny, 2)/s
    end
  end
  return nu
end

# Given a flow and an image, move the image according to the flow
# Pre: size(i,1)==size(u,1)==size(ii,1) && size(i,2)=size(u,2)==size(ii,2) && size(u,3)=2
function warpMotion(i::Matrix{Float32}, u::Array{Float32, 3}, ii::Matrix{Float32})
# Returns Matrix{Float32} size=size(i)
  (h, w) = size(i)
  wi = zeros(Float32, h, w)
  for y in 1:h
    for x in 1:w
      nx = x+u[y,x,1]
      ny = y+u[y,x,2]
      wi[y,x] = (nx<1 || nx>w || ny<1 || ny>h ? ii[y,x] : interpolate(i, nx, ny))
    end
  end
  return wi
end

function readImage(fn::AbstractString)
  f = open(fn, "r")
  w = read(f, UInt32)
  h = read(f, UInt32)
  img = zeros(Float32, h, w)
  for i in 1:h
    for j in 1:w
      img[i,j] = Float32(read(f, UInt8))/255.0
    end
  end
  close(f)
  return img
end

function writeFlo(u::Array{Float32, 3}, fn::AbstractString)
  h = size(u, 1)
  w = size(u, 2)
  f = open(fn, "w")
  write(f, w)
  write(f, h)
  for i in 1:h
    for j in 1:w
      write(f, u[i, j, 1])
      write(f, u[i, j, 2])
    end
  end
  close(f)
end

end


text:
# Example from Nocedal & Wright, p. 281
# Used to test all the different algorithms
@testset "2by2" begin

function f_2by2!(F, x)
    F[1] = (x[1]+3)*(x[2]^3-7)+18
    F[2] = sin(x[2]*exp(x[1])-1)
end

function g_2by2!(J, x)
    J[1, 1] = x[2]^3-7
    J[1, 2] = 3*x[2]^2*(x[1]+3)
    u = exp(x[1])*cos(x[2]*exp(x[1])-1)
    J[2, 1] = x[2]*u
    J[2, 2] = u
end

df = OnceDifferentiable(f_2by2!, g_2by2!, [ -0.5; 1.4], [ -0.5; 1.4])

# Test trust region
r = nlsolve(df, [ -0.5; 1.4], method = :trust_region, autoscale = true)
@test converged(r)
@test norm(r.zero - [ 0; 1]) < 1e-7
r = nlsolve(df, [ -0.5; 1.4], method = :trust_region, autoscale = false)
@test converged(r)
@test norm(r.zero - [ 0; 1]) < 1e-7

df32 = OnceDifferentiable(f_2by2!, g_2by2!, [ -0.5f0; 1.4f0], [ -0.5f0; 1.4f0])
r = nlsolve(df32, [ -0.5f0; 1.4f0], method = :trust_region, autoscale = true)
@test eltype(r.zero) == Float32
@test converged(r)
@test norm(r.zero - [ 0; 1]) < 1e-7
r = nlsolve(df32, [ -0.5f0; 1.4f0], method = :trust_region, autoscale = false)
@test eltype(r.zero) == Float32
@test converged(r)
@test norm(r.zero - [ 0; 1]) < 1e-7

# Test Newton
r = nlsolve(df, [ -0.5; 1.4], method = :newton, linesearch = LineSearches.BackTracking(), ftol = 1e-6)
@test converged(r)
@test norm(r.zero - [ 0; 1]) < 1e-6
r = nlsolve(df32, [ -0.5f0; 1.4f0], method = :newton, linesearch = LineSearches.BackTracking(), ftol = 1e-3)
@test eltype(r.zero) == Float32
@test converged(r)
@test norm(r.zero - [ 0; 1]) < 1e-6
r = nlsolve(df, [ -0.5; 1.4], method = :newton, linesearch = LineSearches.HagerZhang(), ftol = 1e-6)
@test converged(r)
@test norm(r.zero - [ 0; 1]) < 1e-6
r = nlsolve(df, [ -0.5; 1.4], method = :newton, linesearch = LineSearches.StrongWolfe(), ftol = 1e-6)
@test converged(r)
@test norm(r.zero - [ 0; 1]) < 1e-6

# test local convergence of Anderson: close to a fixed-point and with
# a small beta, f should be almost affine, in which case Anderson is
# equivalent to GMRES and should converge
r = nlsolve(df, [ 0.01; .99], method = :anderson, m = 10, beta=.01)
@test converged(r)
@test norm(r.zero - [ 0; 1]) < 1e-8
end


text:
# scoring.jl
# 2019-12-28
# David Merrell
#
# Methods for scoring performance in network inference tasks.
# In particular, for the simulation study (where we score
# edge prediction as a probabilistic classification task)

module Scoring

using JSON, CSV, ArgParse

export aucroc, aucpr, average_precision, julia_main

"""
Compute AUCROC from the given predictions,
with respect to the given ground-truth.
Integral is computed via trapezoidal rule.

Has n*log(n) complexity
(superlinear expense incurred by sorting;
 the AUCROC loop has linear expense).
"""
function aucroc(preds::Vector, truths::Vector{Bool}; return_curve::Bool=false)

    if return_curve
        curve = Vector{Vector{Float64}}()
        push!(curve, [0.0;0.0])
    end

    # Sort the predictions (and corresponding ground truth)
    # by confidence level
    srt_inds = sortperm(preds, order=Base.Order.Reverse)
    srt_preds = preds[srt_inds]
    srt_truths = truths[srt_inds]

    N = length(srt_preds)
    npos = sum(truths)
    nneg = N - npos

    area = 0.0
    tp = 0
    fp = 0
    prev_fp = 0
    prev_tp = 0

    # Loop through the predictions
    # (in order of decreasing confidence level)
    for i=1:N
        if (i > 1) && (srt_preds[i] != srt_preds[i-1])
            area += 0.5*(prev_tp + tp)*(fp - prev_fp)/(npos*nneg)
            if return_curve
                push!(curve, [fp/nneg; tp/npos])
            end
            prev_tp = tp
            prev_fp = fp
        end

        if srt_truths[i]
            tp += 1
        else
            fp += 1
        end
    end
    area += 0.5*(prev_tp + tp)*(fp - prev_fp)/(npos*nneg)

    if return_curve
        push!(curve, [1.0; 1.0])
        return area, transpose(hcat(curve...))
    else
        return area
    end

end

"""
Compute AUCPR from the given predictions,
with respect to the given ground-truth.
Integral is computed via trapezoidal rule.

Has n*log(n) complexity
(superlinear expense incurred by sorting;
 the AUCPR loop has linear expense).
"""
function aucpr(preds::Vector, truths::Vector{Bool}; return_curve::Bool=false)

    if return_curve
        curve = Vector{Vector{Float64}}()
        push!(curve, [0.0;1.0])
    end

    # Sort the predictions (and corresponding ground truth)
    # by confidence level
    srt_inds = sortperm(preds, order=Base.Order.Reverse)
    srt_preds = preds[srt_inds]
    srt_truths = truths[srt_inds]

    area = 0.0
    tp = 0
    fp = 0
    prev_p = 1.0
    prev_r = 0.0
    N = length(srt_truths)
    npos = sum(srt_truths)
    nneg = N - npos

    for i=1:N
        if (i>1) && (srt_preds[i] < srt_preds[i-1])

            recall = tp / npos
            precision = tp / (tp + fp)
            if return_curve
                push!(curve, [recall; precision])
            end
            area += 0.5*(precision + prev_p)*(recall - prev_r)

            prev_r = recall
            prev_p = precision
        end

        if srt_truths[i]
            tp += 1
        else
            fp += 1
        end
    end

    recall = tp / npos
    precision = tp / (tp + fp)
    area += 0.5*(precision + prev_p)*(recall - prev_r)

    if return_curve
        push!(curve, [recall; precision])
        return area, transpose(hcat(curve...))
    end

    return area
end

"""
Compute average precision from the given predictions,
with respect to the given ground-truth.

Has n*log(n) complexity
(superlinear expense incurred by sorting)

return_curve: return the PR curve for these predictions
"""
function average_precision(preds::Vector, truths::Vector{Bool}; return_curve::Bool=false)

    if return_curve
        curve = Vector{Vector{Float64}}()
        push!(curve, [0.0;1.0])
    end

    # Sort the predictions (and corresponding ground truth)
    # by confidence level
    srt_inds = sortperm(preds, order=Base.Order.Reverse)
    srt_preds = preds[srt_inds]
    srt_truths = truths[srt_inds]

    s = 0.0
    tp = 0
    fp = 0
    prev_r = 0.0
    N = length(srt_truths)
    npos = sum(srt_truths)
    nneg = N - npos

    for i=1:N
        if (i>1) && (srt_preds[i] < srt_preds[i-1])

            recall = tp / npos
            precision = tp / (tp + fp)
            if return_curve
                push!(curve, [recall; precision])
            end
            s += precision*(recall - prev_r)

            prev_r = recall
        end

        if srt_truths[i]
            tp += 1
        else
            fp += 1
        end
    end

    recall = tp / npos
    precision = tp / (tp + fp)
    s += precision*(recall - prev_r)

    if return_curve
        push!(curve, [recall; precision])
        return s, transpose(hcat(curve...))
    end

    return s
end

"""
Take a "vector of vectors" representation of the
network edges and return a single vector in canonical order.
"""
function ps_to_vec(parent_sets::Vector)
    return [e for ps in parent_sets for e in ps]
end

"""
Read a JSON file containing edge predictions.
We assume there's a "edge_conf_key" field specifying
where to access the predicted parent sets.

Return a single vector of confidence levels for all edges.
"""
function load_ps_predictions(pred_fname)

    f = open(pred_fname, "r")
    js_d = JSON.parse(read(f, String))
    close(f)

    key = js_d["edge_conf_key"]

    parent_sets = js_d[key]

    return ps_to_vec(parent_sets)
end

"""
Read a CSV file containing a true adjacency matrix.
We assume the file is comma-delimited with no header
or index column; the (i,j)-th entry denotes the jth
parent of the ith vertex (that is: the columns are
parent sets). Each entry must be 0 or 1.
"""
function load_ps_truths(truths_fname)

    df = CSV.read(truths_fname; header=false)
    adj_mat = convert(Matrix{Bool}, df)

    parent_sets = [adj_mat[:,i] for i=1:size(adj_mat,2)]

    return ps_to_vec(parent_sets)
end

"""
Given
* a vector of edge prediction filenames
* a vector of edge truth filenames
* a scoring function;
return a vector of the scores
"""
function edge_score(pred_fname::String,
                    truth_fname::String,
                    score_fn::Function)

    preds = load_ps_predictions(pred_fname)
    truths = load_ps_truths(truth_fname)
    return score_fn(preds, truths)

end

function get_args(args::Vector{String})

    s = ArgParseSettings()
    @add_arg_table! s begin
        "--pred-file"
            help="The path to a JSON file containing edge predictions"
            required=true
            arg_type=String
        "--truth-file"
            help="The path to a CSV file containing true adjacency matrice(s)."
            required=true
            arg_type=String
        "--output-file"
            help="The path of an output JSON file"
            required=true
            arg_type=String
        "--skip-roc"
            help="Do NOT compute ROC or AUCROC for the predictions."
            action=:store_true
    end

    args = parse_args(args, s)
    return args
end

# Main function -- for purposes of static compilation
function julia_main()

    args = get_args(ARGS)

    pred_fname = args["pred-file"]
    truth_fname = args["truth-file"]
    output_fname = args["output-file"]
    skip_roc = args["skip-roc"]

    result = Dict{String,Any}()

    # Compute mean precision and PR curve
    pr_score_fn = (pred, truth) -> average_precision(pred, truth; return_curve=true)
    pr_stuff = edge_score(pred_fname, truth_fname, pr_score_fn)
    result["aucpr"] = pr_stuff[1]
    result["pr_curve"] = pr_stuff[2]

    if !skip_roc
        # Compute AUCROC and ROC curve for each set of predictions
        roc_score_fn = (pred, truth) -> aucroc(pred, truth; return_curve=true)
        roc_stuff = edge_score(pred_fname, truth_fname, roc_score_fn)
        result["aucroc"] = roc_stuff[1]
        result["roc_curve"] = roc_stuff[2]
    end

    f = open(output_fname, "w")
    JSON.print(f, result)
    close(f)

    return 0
end

# END MODULE
end

using .Scoring

julia_main()

max_stars_repo_name: gitter-lab/ssps

max_issues_repo_name: gitter-lab/ssps

max_forks_repo_name: gitter-lab/ssps


text:
config = MOI.Test.Config()
optimize_bridged!(mock) = MOIU.mock_optimize!(mock, [4.0, 1.0])
optimize_cached!(mock) = MOIU.mock_optimize!(mock, [1.0, 4.0])
for mock in [bridged_mock(optimize_bridged!), cached_mock(optimize_cached!)]
    Tests.dsos_quartic_constant_test(mock, config)
    Tests.sdsos_quartic_constant_test(mock, config)
    Tests.sos_quartic_constant_test(mock, config)
end

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
using Dispersal
using Test
using DynamicGrids

@testset "Mortality without resistance grid" begin
    init = (
        exposure   = [100.0 0.0;
                      100.0 0.0],
        population = [100.0 100.0;
                      0.0   0.0],
    )

    @testset "LoglogisticMortality" begin
        loglogmortalityrule = LoglogisticMortality{Tuple{:population,:exposure},:population}(
            median=20, hillcoefficient=0.5, timestep=1
        )
        output = ArrayOutput(init; tspan=1:3)
        sim!(output, loglogmortalityrule)
        @test output[1].exposure == output[2].exposure == output[3].exposure
        @test output[2].population ≈ [30.9017 100.0; 0.0 0.0] atol=1e-4
        @test output[3].population ≈ [9.54915 100.0; 0.0 0.0] atol=1e-4
    end

    @testset "ExponentialMortality low threshold" begin
        expmortalityrule = ExponentialMortality{Tuple{:population,:exposure},:population}(
            rate=0.01, threshold = 20.0, timestep=1
        )
        output = ArrayOutput(init; tspan=1:3)
        sim!(output, expmortalityrule)
        @test output[1].exposure == output[2].exposure == output[3].exposure
        @test output[2].population ≈ [44.9329 100.0; 0.0 0.0] atol=1e-4
        @test output[3].population ≈ [20.1897 100.0; 0.0 0.0] atol=1e-4
    end
    @testset "ExponentialMortality high threshold" begin
        expmortalityrule = ExponentialMortality{Tuple{:population,:exposure},:population}(
            rate=0.01, threshold=200.0, timestep=1
        )
        output = ArrayOutput(init; tspan=1:3)
        sim!(output, expmortalityrule)
        @test output[1].exposure == output[2].exposure == output[3].exposure
        @test output[1].population ==  output[2].population ==  output[3].population
    end
end

@testset "Mortality with Resistance grid" begin

    phenotyperule = Cell{:phenotype, :phenotype}() do data, state, I
        return state * 10
    end

    @testset "LoglogisticMortality with Grid median" begin

        init = (
            exposure   = [100.0 0.0;
                          100.0 0.0],
            population = [100.0 100.0;
                          0.0   0.0],
            phenotype  = [20 20;
                          20 20],
        )
        loglogmortalityrule = LoglogisticMortality{Tuple{:population,:exposure},:population}(
            median=Grid(:phenotype), hillcoefficient=0.5, timestep=1
        )
        output = ArrayOutput(init; tspan=1:3)
        sim!(output, Ruleset(loglogmortalityrule, phenotyperule))
        @test output[1].exposure == output[2].exposure == output[3].exposure
        @test output[2].population ≈ [30.9016 100.0; 0.0 0.0] atol=1e-4
        @test output[3].population ≈ [18.1017 100.0; 0.0 0.0] atol=1e-4
    end

    @testset "ExponentialMortality with Grid rate" begin
        init = (
            exposure   = [100.0 0.0;
                          100.0 0.0],
            population = [100.0 100.0;
                          0.0   0.0],
            phenotype  = [0.01 0.01;
                          0.01 0.01],
        )

        expmortalityrule = ExponentialMortality{Tuple{:population,:exposure},:population}(
            rate=Grid(:phenotype), threshold=20.0, timestep=1
        )
        output = ArrayOutput(init; tspan=1:3)
        sim!(output, Ruleset(expmortalityrule, phenotyperule))
        @test output[1].exposure == output[2].exposure == output[3].exposure
        @test output[2].population ≈ [44.9328 100.0; 0.0 0.0] atol=1e-4
        @test output[3].population ≈ [0.01507 100.0; 0.0 0.0] atol=1e-4
    end
end


text:
abstract type Component end
# `f` is function: (state, val) --> new_state
# `state` is String, Vector{String}, or AllKeys (which parts of the state the component can change).
# `children`
# `id` is required to know which server function `my_component.f` to call.

struct AllKeys end

#-----------------------------------------------------------------------------# Button
struct Button <: Component
    f
    state
    children
    id::String
    Button(f; state = AllKeys(), children) = new(f, state, children, randstring(20))
end

function Base.show(io::IO, M::MIME"text/html", o::Button)
    onClick = "onClick=\${e => this.action(\"$(o.id)\", e.target.value)}"
    node = h.button(o.children, type="button", __VERBATIM__ = onClick)
    show(io, M, node)
end

# #-----------------------------------------------------------------------------# Handler
# # How we set things like: onClick={e => handle(e.target.value)}
# struct Handler
#     key::String
#     function_string::String
# end
# function Base.show(io::IO, ::MIME"text/html", h::Handler)

# end

# #-----------------------------------------------------------------------------# Component
# struct Component
#     f::Function  # (state, value) → new_state
#     node::Hyperscript.Node
#     id::String  # use same `id` in html as well as in Julia
# end

# function button(f, args...; kw...)
#     id = randstring(20)
#     Component(f, m("button", type="button", args...; id, kw...), id)
# end

# #-----------------------------------------------------------------------------# Button
# struct Button <: Component
#     server_function  # (state,val) -> mutated state
#     children
#     active::Bool
#     id::String
# end
# Button(f, children; active::Bool=false, id=randstring(10)) = Button(f, children, active ,id)

# function js(o::Button)
#     """
#     function Button({children, active}) {
#         const c = active ? "text-white bg-indigo-600 hover:bg-indigo-700" : "text-indigo-700 bg-indigo-100 hover:bg-indigo-200"
#         return html`
#             <button type="button"

#                 class={`inline-flex items-center m-2 px-4 py-2 border border-transparent text-sm font-medium rounded-md shadow-sm \${c} focus:outline-none focus:ring-2 focus:ring-offset-2 focus:ring-indigo-500`}

#                 onClick={e => rountTrip($(o.id), e.target.value)}
#             >
#                 {children}
#             </button>
#         `
#       }
#     """
# end


text:
function parsestring(s)
    s = replace(s, "long" => "") # In the 2017 do file they declared some variables as long
    tokens = String[]
    n = 0
    i = 1
    while i <= lastindex(s)
        if s[i] == ' '
            i = i + 1
        else
            nx = findnext("  ", s, i)
            if nx === nothing
                push!(tokens, s[i:end])
                break
            else
                push!(tokens, s[i:(first(nx) - 1) ])
                i = first(nx)
            end
        end
    end
    tokens
end

function process_tok(tok)
    alltokens = String[]
    ready = false
    for i = 1:length(tok)
        w = tok[i]
        if "infix" ⊆ w
            ready = true
            continue
        end
        "using" ⊆ w && break
        ready && append!(alltokens, parsestring(w))
    end
    alltokens
end

function str2range(s)
    myr = r"\d+"
    d = eachmatch(myr, s)
    c = collect(d)
    parse(Int, first(c).match):parse(Int, last(c).match)
end

function read_fixedwidth(data, toks)
    names = String[tok[1] for tok in toks]
    dat = zeros(length(data), length(toks))
    for i in eachindex(data)
        line = data[i]
        for j in eachindex(toks)
            tok = toks[j]
            r = tok[2]
            dat[i, j] = Parsers.parse(Float64, line[r])
        end
    end
    DataFrame(dat, Symbol.(names))
end

function readPSID(zipname)
    zipp = "$zipname.zip"
    ZIPNAME = uppercase(zipname)
    #t = mktempdir()
    t = "datafiles/$zipname"
    #run(`unzip $zipname -d $t`)
    isdir(t) || run(DataDeps.unpack_cmd("$zipp", "$t", ".zip", ""))
    tok = readlines("$t/$ZIPNAME.do")
    alltokens = process_tok(tok)
    toks = [(alltokens[i], str2range(alltokens[i+1])) for i in 1:2:length(alltokens)]
    data = readlines("$t/$ZIPNAME.txt")
    out = read_fixedwidth(data, toks)
end

function unzip_data()
    years = [collect(1968:1997); collect(1999:2:2017)]
    filenames = [year <= 1993 ? "fam$year" : "fam$(year)er" for year in years]
    datas =  SortedDict(year => readPSID(filename) for (year, filename) in zip(years, filenames))
    inddata = readPSID("ind2017er")
    (famdatas = datas, inddata = inddata)
end


text:
export HomotopyWithCache

struct HomotopyWithCache{H<:AbstractHomotopy,C<:AbstractHomotopyCache} <: AbstractHomotopy
    homotopy::H
    cache::C
end

HomotopyWithCache(H::AbstractHomotopy, x, t) = HomotopyWithCache(H, cache(H, x, t))

cache(::HomotopyWithCache, x, t) = HomotopyNullCache()

(H::HomotopyWithCache)(x, t) = evaluate(H, x, t)

Base.size(H::HomotopyWithCache) = size(H.homotopy)
evaluate(H::HomotopyWithCache, x, t, ::HNC = HNC()) = evaluate(H.homotopy, x, t, H.cache)
@propagate_inbounds evaluate!(u, H::HomotopyWithCache, x, t, ::HNC = HNC()) =
    evaluate!(u, H.homotopy, x, t, H.cache)
jacobian(H::HomotopyWithCache, x, t, ::HNC = HNC()) = jacobian(H.homotopy, x, t, H.cache)
@propagate_inbounds jacobian!(U, H::HomotopyWithCache, x, t, ::HNC = HNC()) =
    jacobian!(U, H.homotopy, x, t, H.cache)
dt(H::HomotopyWithCache, x, t, ::HNC = HNC()) = dt(H.homotopy, x, t, H.cache)
@propagate_inbounds dt!(u, H::HomotopyWithCache, x, t, ::HNC = HNC()) =
    dt!(u, H.homotopy, x, t, H.cache)
jacobian_and_dt(H::HomotopyWithCache, x, t, ::HNC = HNC()) =
    jacobian_and_dt(H.homotopy, x, t, H.cache)
@propagate_inbounds jacobian_and_dt!(U, u, H::HomotopyWithCache, x, t, ::HNC = HNC()) =
    jacobian_and_dt!(U, u, H.homotopy, x, t, H.cache)
evaluate_and_jacobian(H::HomotopyWithCache, x, t, ::HNC = HNC()) =
    evaluate_and_jacobian(H.homotopy, x, t, H.cache)
@propagate_inbounds evaluate_and_jacobian!(
    u,
    U,
    H::HomotopyWithCache,
    x,
    t,
    ::HNC = HNC(),
) = evaluate_and_jacobian!(u, U, H.homotopy, x, t, H.cache)
basehomotopy(H::HomotopyWithCache) = basehomotopy(H.homotopy)


text:
function init_constraint!(
    com::CS.CoM,
    constraint::AndConstraint,
    fct,
    set::AndSet;
)
    set_impl_functions!(com, constraint.lhs)
    set_impl_functions!(com, constraint.rhs)
    if constraint.lhs.impl.init
        !init_constraint!(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set) && return false
    end
    if constraint.rhs.impl.init
        !init_constraint!(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set) && return false
    end
    return true
end

"""
    is_constraint_solved(
        constraint::AndConstraint,
        fct,
        set::AndSet,
        values::Vector{Int},
    )

Checks if both inner constraints are solved
"""
function is_constraint_solved(
    constraint::AndConstraint,
    fct,
    set::AndSet,
    values::Vector{Int},
)
    lhs_num_vars = get_num_vars(constraint.lhs.fct)
    lhs_solved = is_constraint_solved(constraint.lhs, constraint.lhs.fct, constraint.lhs.set, values[1:lhs_num_vars])
    !lhs_solved && return false
    rhs_num_vars = get_num_vars(constraint.rhs.fct)
    rhs_solved = is_constraint_solved(constraint.rhs, constraint.rhs.fct, constraint.rhs.set, values[end-rhs_num_vars+1:end])
    return rhs_solved
end

"""
    function is_constraint_violated(
        com::CoM,
        constraint::AndConstraint,
        fct,
        set::AndSet,
    )

Check if either of the inner constraints are violated already
"""
function is_constraint_violated(
    com::CoM,
    constraint::AndConstraint,
    fct,
    set::AndSet,
)
    lhs_violated = is_constraint_violated(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set)
    lhs_violated && return true
    rhs_violated = is_constraint_violated(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set)
    return lhs_violated || rhs_violated
end

"""
    still_feasible(com::CoM, constraint::AndConstraint, fct, set::AndSet, vidx::Int, value::Int)

Return whether the constraint can be still fulfilled when setting a variable with index `vidx` to `value`.
**Attention:** This assumes that it isn't violated before.
"""
function still_feasible(
    com::CoM,
    constraint::AndConstraint,
    fct,
    set::AndSet,
    vidx::Int,
    value::Int,
)
    lhs_indices = constraint.lhs.indices
    for i in 1:length(lhs_indices)
        if lhs_indices[i] == vidx
            !still_feasible(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set, vidx, value) && return false
        end
    end
    rhs_indices = constraint.rhs.indices
    for i in 1:length(rhs_indices)
        if rhs_indices[i] == vidx
            !still_feasible(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set, vidx, value) && return false
        end
    end
    return true
end

"""
    prune_constraint!(com::CS.CoM, constraint::AndConstraint, fct, set::AndSet; logs = true)

Reduce the number of possibilities given the `AndConstraint` by pruning both parts
Return whether still feasible
"""
function prune_constraint!(
    com::CS.CoM,
    constraint::AndConstraint,
    fct,
    set::AndSet;
    logs = true,
)
    !activate_lhs!(com, constraint) && return false
    !prune_constraint!(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set; logs=logs) && return false
    !activate_rhs!(com, constraint) && return false
    feasible = prune_constraint!(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set; logs=logs)
    return feasible
end


text:
type ColorScheme{Object}
    x::Dict
end

"""
Tree structure, to allow for hierarchical coloring.
Each key holds either a key to a Color, or a new key into the Dict.
Keys are symbols, so values must not be symbols, as they will then be used for lookup.
"""
immutable ColorScheme
    container::Dict{Symbol, Any}
end
function ColorScheme(values::Vector)
    # we're mega hip and make every value a signal,
    # so that one can interactively change the colorscheme
    dict = Dict{Symbo, Any}(map(values) do pair
        key, value = pair
        if !isa(value, Symbol) && !isa(value, Signal)
            return key => Signal(value)
        end
        key=>value
    end)
    ColorScheme(dict)
end
function Base.setindex!(x::ColorScheme, key::Symbol, value)
    x.container[key] = value
end
function Base.get(td::ColorScheme, key::Symbol, default)
    circle_check = Set()
    # values are either new keys, or values.
    # Keys are always symbols, values never. So we iterate as long as the value is a Symbol
    while !isa(key, Symbol)
        haskey(td, key) || return default # if doesn't contain key, return default
        if value in circle_check # check if we visited key already
            error("Circle detected for value: $value, in tree dict: $td")
        end
        push!(circle_check, key)
        key = td.container[key]
    end
    # found non symbol, must be a value!
    key
end

function color_scheme(code::Code)
    code.language != "Julia" && error("only Julia code supported. Found: $code.language")
    ColorScheme([
        :Any => RGB(0,0,0),
        :Keyword => RGB(0,0,0),
        :Scope => RGB(0,0,0),
        :Comment => RGB(0.9, 0.9, 0.9),
        :Number => RGB(0.9, 0.9, 0.9),
        :FloatingPoint => RGB(0.9, 0.9, 0.9),
        :Integer => RGB(0.9, 0.9, 0.9),
        :Operator => RGB(1,0,0),
        :Call => RGB(0,0,1),
        :String => RGB(0.9, 0.9, 0.9),
        :symbol => RGB(0.8, 0.8, 0.1),

        :macro_call => :Call,
        :call => :Call,

        :ControlFlow => :Keyword,
        :function => :Keyword,
        :macro => :Keyword,
        :do => :Keyword,
        :OPKeyword => :Keyword,
        :(::) => :OPKeyword,
        symbol(",") => :OPKeyword,
        symbol(";") => :OPKeyword,
        symbol(":") => :OPKeyword,
        symbol("...") => :OPKeyword,
        :end => :Scope,
        :while => :ControlFlow,
        :continue => :ControlFlow,
        :try => :ControlFlow,
        :catch => :ControlFlow,
        :finally => :ControlFlow,
        :break => :ControlFlow,
        :for => :ControlFlow,
        :if => :ControlFlow,
        :else => :ControlFlow,
        :elseif => :ControlFlow,
        :(&&) => :ControlFlow,
        :(||) => :ControlFlow,
        :ArrayScope => :Scope,
        symbol("(") => :ArrayScope,
        symbol(")") => :ArrayScope,
        symbol("[") => :ArrayScope,
        symbol("]") => :ArrayScope,
        :Curly       => :Scope,
        symbol("{") => :Curly,
        symbol("}") => :Curly
    ])
end

token2symbol(token::Lexer.AbstractToken) = token2symbol(token.val)
token2symbol(token::Char) = symbol(token)
token2symbol(token::Number) = :Number
token2symbol(token::FloatingPoint) = :FloatingPoint
token2symbol(token::Integer) = :Integer
token2symbol(token) = symbol(string(token))

function token2color(colorscheme, token)
    sym = token2symbol(token)
    get(colorscheme, sym, get(colorscheme, :Any, RGB(0.,0.,0.)))
end
function tkstream2color(colorscheme, tokenstream)
    string_started = false
    symbol_started = false
    macro_call_started = false
    last_token = peek(tokenstream)
    while !Lexer.eof(token)
        token2color(token)
    end
end

# Block elements

printmd(io::IO, md::MD) = printmd(io, md.content)

function printmd(io::RichText, content::Vector)
    for md in content
        printmd(io, md)
        println(io)
    end
end

function printmd{l}(io::RichText, header::Header{l})
    withtag(io, "h$l") do
        htmlinline(io, header.text)
    end
end

function printmd(io::RichText, code::Code)
    insert!(io, code)
end

function printmd(io::RichText, md::Paragraph)
    withtag(io, :p) do
        htmlinline(io, md.content)
    end
end

function printmd(io::RichText, md::BlockQuote)
    withtag(io, :blockquote) do
        println(io)
        printmd(io, md.content)
    end
end

function printmd(io::RichText, md::List)
    withtag(io, md.ordered ? :ol : :ul) do
        for item in md.items
            println(io)
            withtag(io, :li) do
                htmlinline(io, item)
            end
        end
        println(io)
    end
end

function printmd(io::RichText, md::HorizontalRule)
    insert!(io, Line(Point2f0(0), Point2f0(area.w, 0)))
end

# Inline elements

function htmlinline(io::IO, content::Vector)
    for x in content
        htmlinline(io, x)
    end
end

function htmlinline(io::IO, code::Code)
    withtag(io, :code) do
        htmlesc(io, code.code)
    end
end

function htmlinline(io::IO, md::Union{Symbol, AbstractString})
    htmlesc(io, md)
end

function htmlinline(io::IO, md::Bold)
    withtag(io, :strong) do
        htmlinline(io, md.text)
    end
end

function htmlinline(io::IO, md::Italic)
    withtag(io, :em) do
        htmlinline(io, md.text)
    end
end

function htmlinline(io::IO, md::Image)
    tag(io, :img, :src=>md.url, :alt=>md.alt)
end

function htmlinline(io::IO, link::Link)
    withtag(io, :a, :href=>link.url) do
        htmlinline(io, link.text)
    end
end

function htmlinline(io::IO, br::LineBreak)
    tag(io, :br)
end


text:
module MathOptInterfaceSCS

export SCSOptimizer

using MathOptInterface
const MOI = MathOptInterface
const CI = MOI.ConstraintIndex
const VI = MOI.VariableIndex

const MOIU = MOI.Utilities

const SF = Union{MOI.SingleVariable, MOI.ScalarAffineFunction{Float64}, MOI.VectorOfVariables, MOI.VectorAffineFunction{Float64}}
const SS = Union{MOI.EqualTo{Float64}, MOI.GreaterThan{Float64}, MOI.LessThan{Float64}, MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, MOI.SecondOrderCone, MOI.ExponentialCone, MOI.PositiveSemidefiniteConeTriangle}

using SCS

include("types.jl")

mutable struct SCSOptimizer <: MOI.AbstractOptimizer
    cone::Cone
    maxsense::Bool
    data::Union{Void, Data} # only non-Void between MOI.copy! and MOI.optimize!
    sol::Solution
    function SCSOptimizer()
        new(Cone(), false, nothing, Solution())
    end
end

function MOI.isempty(optimizer::SCSOptimizer)
    !optimizer.maxsense && optimizer.data === nothing
end
function MOI.empty!(optimizer::SCSOptimizer)
    optimizer.maxsense = false
    optimizer.data = nothing # It should already be nothing except if an error is thrown inside copy!
end

MOI.canaddvariable(optimizer::SCSOptimizer) = false

MOI.supports(::SCSOptimizer, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}) = true
MOI.supportsconstraint(::SCSOptimizer, ::Type{<:SF}, ::Type{<:SS}) = true
MOI.copy!(dest::SCSOptimizer, src::MOI.ModelLike; copynames=true) = MOIU.allocateload!(dest, src, copynames)

# Implements optimize! : translate data to SCSData and call SCS_solve
include("solve.jl")

# Implements getter for result value and statuses
include("attributes.jl")

end # module