text:
module Inflation

export convert_01tuples
export rows_and_cols
export get_M_matrix
export get_b_array
export get_v_array_Variable
export integrate_pfunc

export pfunc1
export pfunc2
export get_pfunct

using Convex
using GLPK
using Printf

export InflKind
export cut
export spiral

@enum InflKind cut spiral

function convert_01tuples(str::String, char::Char)::Vector{String}
	# accepts a string with a certain character 'char'
	# returns all possible combinations where 'char' is
	# replaced with either 0 or 1

	foundall  = false
	newcombos = Vector{String}([str])
	combos    = Vector{String}([])
	while !foundall
		foundall  = true
		combos    = newcombos
		newcombos = []
		for combo in combos
			pos = findfirst(char, combo)
			if pos != nothing
				combo0 = replace(combo, char=>'0', count=1)
				combo1 = replace(combo, char=>'1', count=1)
				push!(newcombos, combo0)
				push!(newcombos, combo1)
				foundall = false
			end
		end
	end

	return combos
end

function rows_and_cols(infl_strings::Vector{String})::Tuple{Int64, Int64}
	# returns the number of rows and columns of the M-vector
	if isempty(infl_strings)
		error("'infl_strings' must have at least one element")
	end

	str1 = infl_strings[1]
	infl_strings_y = convert_01tuples(str1, 'y')

	nrows = length(infl_strings)*length(infl_strings_y)
	ncols = 2^length(str1)

	return (nrows, ncols)
end

function get_M_matrix(infl_strings::Vector{String})::Array{Int64, 2}
	# returns the M-vector needed for the inflation algorithm
	# infl_strings : array of strings that determine the inflation combinations
	# EXAMPLE: infl_strings = ['yyx', 'yxy', 'xyy']

	Nrows, Ncols = rows_and_cols(infl_strings)
	Mvec = zeros(Int64, Nrows, Ncols)

	for (it, str_yx) in enumerate(infl_strings)
		infl_strings_x = convert_01tuples(str_yx, 'y')
		N_x            = length(infl_strings_x)
		for (ix, str_x) in enumerate(infl_strings_x)
			i = (ix-1) + N_x*(it-1)
			fill_row_Mvec(i, str_x, Mvec)
		end
	end

	return Mvec
end

function fill_row_Mvec(i::Int64, str_x::String, Mvec::Array{Int64, 2})::Nothing
	# fills the i^th row of Mvec using the combinations given by str_x
	for str in convert_01tuples(str_x, 'x')
		ic = parse(Int64, str, base=2)
		Mvec[i+1, ic+1] = 1
	end
end

function get_b_array(infl_strings::Vector{String}, pfunc::Function, kind::InflKind)::Array{Float64, 1}
	# returns the b-vector for the Cut Inflation of the triangle causal model
	[@assert only_contains("xy", str) for str in infl_strings]

	Nrows_per_key = length(convert_01tuples(infl_strings[1], 'y'))
	Nkeys         = length(infl_strings)
	barr          = Vector{Float64}(undef, Nrows_per_key*Nkeys)

	for (i, str_yx) in enumerate(infl_strings)
		infl_strings_x = convert_01tuples(str_yx, 'y')
		for (j, str_x) in enumerate(infl_strings_x)
			#index = (i-1) + Nkeys*(j-1) + 1
			index = (j-1) + Nrows_per_key*(i-1) + 1
			barr[index] = inflation_RHS(str_x, pfunc, kind)
		end
	end

	return barr
end

function get_v_array_Variable(infl_strings::Vector{String})
	# gives the Variable instance with the correct size for the v-array
	Nv_size = 2^length(infl_strings[1])
	return Variable(Nv_size)
end

function inflation_RHS_Cut(str_x::String, pfunc::Function)::Float64
	# performs the RHS integration for the Cut inflation

	# check if 'str_x' is valid for the Cut inflation integration
	Nx = sum([char == 'x' for char in str_x])
	@assert Nx            == 1
	@assert length(str_x) == 3

	total = 0.0
	if str_x[3] == 'x'
		# eg) '00x', or '01x', ...
		strA = str_x[1] * "xx"
		strB = "x" * str_x[2] * "x"
		totalA = integrate_pfunc(pfunc, strA, 'x')
		totalB = integrate_pfunc(pfunc, strB, 'x')
		total += totalA*totalB
	else
		# eg) 'x00', '0x1', ...
		total += integrate_pfunc(pfunc, str_x, 'x')
	end

	return total
end

function inflation_RHS_Spiral(str_x::String, pfunc::Function)::Float64
	# performs the RHS integration for the Spiral inflation

	# check if 'str_x' is valid for the Spiral inflation integration
	Nx = sum([char == 'x' for char in str_x])
	@assert Nx            == 3
	@assert length(str_x) == 6

	total = 0.0
	if     match_str(str_x, [2, 4, 6], 'x')
		# 'yxyxyx'
		str = str_x[1] * str_x[3] * str_x[5]
		total += pfunc(str)
	elseif match_str(str_x, [2, 3, 5], 'x')
		# 'yxxyxy'
		str   = str_x[1] * str_x[4] * str_x[6]
		strAB = str[1:2] * "x"
		strC  = "xx" * str[3]

		totalAB = integrate_pfunc(pfunc, strAB, 'x')
		totalC  = integrate_pfunc(pfunc, strC,  'x')
		total  += totalAB*totalC
	elseif match_str(str_x, [1, 4, 5], 'x')
		# 'xyyxxy'
		str   = str_x[2] * str_x[3] * str_x[6]
		strBC = "x" * str[2:3]
		strA  = str[1] * "xx"

		totalBC = integrate_pfunc(pfunc, strBC, 'x')
		totalA  = integrate_pfunc(pfunc, strA,  'x')
		total  += totalA*totalBC
	elseif match_str(str_x, [1, 3, 6], 'x')
		# 'xyxyyx'
		str   = str_x[2] * str_x[4] * str_x[5]
		strAC = str[1] * "x" * str[3]
		strB  = "x" * str[2] * "x"

		totalAC = integrate_pfunc(pfunc, strAC, 'x')
		totalB  = integrate_pfunc(pfunc, strB,  'x')
		total  += totalAC*totalB

	elseif match_str(str_x, [1, 3, 5], 'x')
		# 'xyxyxy'
		str  = str_x[2] * str_x[4] * str_x[6]
		strA = str[1] * "xx"
		strB = "x" * str[2] * "x"
		strC = "xx" * str[3]

		totalA = integrate_pfunc(pfunc, strA, 'x')
		totalB = integrate_pfunc(pfunc, strB, 'x')
		totalC = integrate_pfunc(pfunc, strC, 'x')
		total += totalA*totalB*totalC
	else
		error("Invalid inflation string.")
	end

	return total
end

function match_str(str_x::String, indices::Vector{Int64}, char::Char)::Bool
	# check if 'str_x' has 'char' at certain indices
	matches = true
	for i in indices
		if str_x[i] != char
			matches = false
			break
		end
	end

	return matches
end

function inflation_RHS(str_x::String, pfunc::Function, kind::InflKind)::Float64
	# select which kind of RHS inflation to use
	if     kind == cut
		return inflation_RHS_Cut(str_x, pfunc)
	elseif kind == spiral
		return inflation_RHS_Spiral(str_x, pfunc)
	else
		error("Only the inflations for 'Cut' and 'Spiral' have been implemented.")
	end
end

function integrate_pfunc(pfunc::Function, str_x::String, char::Char)::Float64
	# pfunc : the probability distribution for the triangle causal model
	# str_x : a string of 0's, 1's, and 'char's
	#       : the 'char's will be replaced with 0's and 1's
	#       : and the combination of all these terms will be added together

	# check if str_x only has 0's, 1's, and char's
	@assert only_contains("01"*char, str_x)

	infl_strings_01 = convert_01tuples(str_x, char)
	total = sum([pfunc(str) for str in infl_strings_01])

	return total
end

function get_pfunct(t::Float64)::Function
	# for problem 2, return pfunc(a, b, c; t)
	@assert 0.0 <= t <= 1.0

	function pfunct(str::String)::Float64
		# returns the 1st probability distribution for the triangle causal model
		@assert (a = tryparse(Int64, string(str[1]))) != nothing
		@assert (b = tryparse(Int64, string(str[2]))) != nothing
		@assert (c = tryparse(Int64, string(str[3]))) != nothing
		@assert (a==0 || a==1)
		@assert (b==0 || b==1)
		@assert (c==0 || c==1)

		if (a == b == c)
			return t*0.5
		else
			return (1.0 - t)/6.0
		end
	end

	return pfunct
end

function pfunc1(str::String)::Float64
	@assert length(str) == 3
	@assert (a = tryparse(Int64, string(str[1]))) != nothing
	@assert (b = tryparse(Int64, string(str[2]))) != nothing
	@assert (c = tryparse(Int64, string(str[3]))) != nothing
	return pfunc1(a, b, c)
end

function pfunc1(a::Int64, b::Int64, c::Int64)::Float64
	# returns the 1st probability distribution for the triangle causal model
	@assert (a==0 || a==1)
	@assert (b==0 || b==1)
	@assert (c==0 || c==1)

	if (a == b == c)
		return 0.5
	else
		return 0.0
	end
end

function pfunc2(str::String)::Float64
	@assert length(str) == 3
	@assert (a = tryparse(Int64, string(str[1]))) != nothing
	@assert (b = tryparse(Int64, string(str[2]))) != nothing
	@assert (c = tryparse(Int64, string(str[3]))) != nothing
	return pfunc2(a, b, c)
end

function pfunc2(a::Int64, b::Int64, c::Int64)::Float64
	# returns the 2nd probability distribution for the triangle causal model
	@assert (a==0 || a==1)
	@assert (b==0 || b==1)
	@assert (c==0 || c==1)

	if (a + b + c == 1)
		return 1.0/3.0
	else
		return 0.0
	end
end

function only_contains(ok_chars::String, str::String)::Bool
	# Checks that 'str' only contains elements from 'ok_chars'
	return all([c in ok_chars for c in str])
end

end # module


text:
export Bisection, Newton, Krawczyk

"""
    𝒩(f, f′, X, α)

Single-variable Newton operator.

The symbol for the operator is accessed with `\\scrN<tab>`.
"""
function 𝒩(f, f′, X::Interval{T}, α) where {T}
    m = Interval(mid(X, α))

    return m - (f(m) / f′(X))
end

"""
    𝒩(f, jacobian, X, α)

Multi-variable Newton operator.
"""
function 𝒩(f::Function, jacobian::Function, X::IntervalBox, α)  # multidimensional Newton operator
    m = IntervalBox(Interval.(mid(X, α)))
    J = jacobian(X)

    return IntervalBox(m .- (J \ f(m)))
end

"""
    𝒦(f, f′, X, α)

Single-variable Krawczyk operator.

The symbol for the operator is accessed with `\\scrK<tab>`.
"""
function 𝒦(f, f′, X::Interval{T}, α) where {T}
    m = Interval(mid(X, α))
    Y = 1 / f′(m)

    return m - Y*f(m) + (1 - Y*f′(X)) * (X - m)
end

"""
    𝒦(f, jacobian, X, α)

Multi-variable Krawczyk operator.
"""
function 𝒦(f, jacobian, X::IntervalBox{T}, α) where {T}
    m = mid(X, α)
    mm = IntervalBox(m)
    J = jacobian(X)
    Y = mid.(inv(jacobian(mm)))

    return m - Y*f(mm) + (I - Y*J) * (X.v - m)
end

"""
    Contractor{F}

    Abstract type for contractors.
"""
abstract type Contractor{F} end

"""
    Bisection{F} <: Contractor{F}

Contractor type for the bisection method.
"""
struct Bisection{F} <: Contractor{F}
    f::F
end

function (contractor::Bisection)(X, tol)
    image = (contractor.f)(X)

    if !(contains_zero(image))
        return :empty, X
    end

    return :unknown, X
end

for (Method, 𝒪) in ((:Newton, 𝒩), (:Krawczyk, 𝒦))
    doc = """
        $Method{F, FP} <: Contractor{F}

    Contractor type for the $Method method.

    # Fields
        - `f::F`: function whose roots are searched
        - `f::FP`: derivative or jacobian of `f`

    -----

        (C::$Method)(X, tol, α=where_bisect)

    Contract an interval `X` using $Method operator and return the
    contracted interval together with its status.

    # Inputs
        - `X`: Interval to contract.
        - `tol`: Precision to which unique solutions are refined.
        - `α`: Point of bisection of intervals.
    """

    @eval begin
        struct $Method{F, FP} <: Contractor{F}
            f::F
            f′::FP   # use \prime<TAB> for ′
        end

        function (C::$Method)(X, tol, α=where_bisect)
            op = x -> $𝒪(C.f, C.f′, x, α)
            rt = determine_region_status(op, C.f, X)
            return refine(op, rt, tol)
        end

        @doc $doc $Method
    end
end

"""
    safe_isempty(X)

Similar to `isempty` function for `IntervalBox`, but also works for `SVector`
of `Interval`.
"""
safe_isempty(X) = isempty(IntervalBox(X))

"""
    determine_region_status(contract, f, X)

Contraction operation for contractors using the first derivative of the
function.

Currently `Newton` and `Krawczyk` contractors use this.
"""
function determine_region_status(op, f, X)
    imX = f(X)

    !(contains_zero(imX)) && return :empty, X

    safe_isempty(imX) && return :empty, X  # X is fully outside of the domain of f

    contracted_X = op(X)

    # Only happens if X is partially out of the domain of f
    safe_isempty(contracted_X) && return :unknown, X  # force bisection

    # given that have the Jacobian, can also do mean value form
    NX = contracted_X ∩ X

    isinf(X) && return :unknown, NX  # force bisection
    safe_isempty(NX) && return :empty, X

    NX ⪽ X  && return :unique, NX  # isinterior; know there's a unique root inside

    return :unknown, NX
end

"""
    refine(op, X::Region, tol)

Generic refine operation for Krawczyk and Newton.
This function assumes that it is already known that `X` contains a unique root.
"""
function refine(op, X::Region, tol)
    while diam(X) > tol  # avoid problem with tiny floating-point numbers if 0 is a root
        NX = op(X) ∩ X
        NX == X && break  # reached limit of precision
        X = NX
    end

    return X
end

"""
    refine(op, X::Tuple{Symbol, Region}, tol)

Wrap the refine method to leave unchanged intervals that are not guaranteed to
contain an unique solution.
"""
function refine(op, rt::Tuple{Symbol, Region}, tol)
    rt[1] != :unique && return rt
    return :unique, refine(op, rt[2], tol)
end

ext: jl


text:
using Documenter, COBREXA
using Literate

# Note: required to deploy the doc from Gitlab CI instead of Travis
ENV["TRAVIS_REPO_SLUG"] = "LCSB-BioCore/COBREXA.jl"
ENV["TRAVIS_BRANCH"] = "master"

# set the merge/pull request ID
if "CI_EXTERNAL_PULL_REQUEST_IID" in keys(ENV)
    ENV["TRAVIS_PULL_REQUEST"] = ENV["CI_EXTERNAL_PULL_REQUEST_IID"]
else
    ENV["TRAVIS_PULL_REQUEST"] = "false"
end

# generate notebooks
notebooks_path = joinpath(@__DIR__, "src", "notebooks")
notebooks =
    joinpath.(notebooks_path, filter(x -> endswith(x, ".jl"), readdir(notebooks_path)))
notebooks_outdir = joinpath(@__DIR__, "src", "notebooks")

# only temporary - will be removed once properly tagged and released
folder = "dev"
branch = "gh-pages"

for notebook in notebooks
    Literate.markdown(
        notebook,
        notebooks_outdir;
        repo_root_url = "https://github.com/$(ENV["TRAVIS_REPO_SLUG"])/blob/master",
        nbviewer_root_url = "https://nbviewer.jupyter.org/github/$(ENV["TRAVIS_REPO_SLUG"])/blob/gh-pages/$(folder)",
        binder_root_url = "https://mybinder.org/v2/gh/$(ENV["TRAVIS_REPO_SLUG"])/$(branch)?filepath=$(folder)",
    )
    Literate.notebook(notebook, notebooks_outdir)
end

# generate index.md from .template and the quickstart in README.md
quickstart = match(
    r"<!--quickstart_begin-->\n([^\0]*)<!--quickstart_end-->",
    open(f -> read(f, String), joinpath(@__DIR__, "..", "README.md")),
).captures[1]
index_md = replace(
    open(f -> read(f, String), joinpath(@__DIR__, "src", "index-template.md")),
    "<!--insert_quickstart-->\n" => quickstart,
)
open(f -> write(f, index_md), joinpath(@__DIR__, "src", "index.md"), "w")

# build the docs
makedocs(
    modules = [COBREXA],
    clean = false,
    sitename = "COBREXA.jl",
    format = Documenter.HTML(
        # Use clean URLs, unless built as a "local" build
        prettyurls = !("local" in ARGS),
        assets = ["assets/favicon.ico"],
        highlights = ["yaml"],
    ),
    authors = "The developers of COBREXA.jl",
    linkcheck = !("skiplinks" in ARGS),
    pages = [
        "Home" => "index.md",
        "Tutorials" => "tutorials.md",
        "Examples and notebooks" => "notebooks.md",
        "Function reference" => "functions.md",
        "How to contribute" => "howToContribute.md",
    ],
)

deploydocs(
    repo = "github.com/$(ENV["TRAVIS_REPO_SLUG"]).git",
    target = "build",
    branch = "gh-pages",
    push_preview = true,
)


text:
"""
GradientGrassmann is an optimisation methdod that keeps the MPS in left-canonical form, and
treats the tensors as points on Grassmann manifolds. It then applies one of the standard
gradient optimisation methods, e.g. conjugate gradient, to the MPS, making use of the
Riemannian manifold structure. A preconditioner is used, so that effectively the metric used
on the manifold is that given by the Hilbert space inner product.

The arguments to the constructor are
method = OptimKit.ConjugateGradient
    The gradient optimisation method to be used. Should either be an instance or a subtype
    of `OptimKit.OptimizationAlgorithm`. If it's an instance, this `method` is simply used
    to do the optimisation. If it's a subtype, then an instance is constructed as
    `method(; maxiter=maxiter, verbosity=verbosity, gradtol=tol)`

finalize! = OptimKit._finalize!
    A function that gets called once each iteration. See OptimKit for details.

tol = Defaults.tol
maxiter = Defaults.maxiter
verbosity = 2
    Arguments passed to the `method` constructor. If `method` is an instance of
    `OptimKit.OptimizationAlgorithm`, these argument are ignored.

In other words, by default conjugate gradient is used. One can easily set `tol`, `maxiter`
and `verbosity` for it, or switch to LBFGS or gradient descent by setting `method`. If more
control is wanted over things like specifics of the linesearch, CG flavor or the `m`
parameter of LBFGS, then the user should create the `OptimKit.OptimizationAlgorithm`
instance manually and pass it as `method`.
"""
struct GradientGrassmann <: Algorithm
    method::OptimKit.OptimizationAlgorithm
    finalize!::Function

    function GradientGrassmann(; method = ConjugateGradient,
                               finalize! = OptimKit._finalize!,
                               tol = Defaults.tol,
                               maxiter = Defaults.maxiter,
                               verbosity = 2)
        if isa(method, OptimKit.OptimizationAlgorithm)
            # We were given an optimisation method, just use it.
            m = method
        elseif method <: OptimKit.OptimizationAlgorithm
            # We were given an optimisation method type, construct an instance of it.
            m = method(; maxiter=maxiter, verbosity=verbosity, gradtol=tol)
        else
            msg = "method should be either an instance or a subtype of OptimKit.OptimizationAlgorithm."
            throw(ArgumentError(msg))
        end
        return new(m, finalize!)
    end
end

function find_groundstate(state::S, H::HT, alg::GradientGrassmann,
                          envs::P=environments(state, H))::Tuple{S,P,Float64} where {S,HT<:Hamiltonian,P}

    !isa(state,FiniteMPS) || dim(state.CR[end]) == 1 || @warn "This is not fully supported - split the mps up in a sum of mps's and optimize seperately"
    normalize!(state)
    res = optimize(GrassmannMPS.fg, (state, envs), alg.method;
                   transport! = GrassmannMPS.transport!,
                   retract = GrassmannMPS.retract,
                   inner = GrassmannMPS.inner,
                   scale! = GrassmannMPS.scale!,
                   add! = GrassmannMPS.add!,
                   finalize! = alg.finalize!,
                   precondition = GrassmannMPS.precondition,
                   isometrictransport = true)
    (x, fx, gx, numfg, normgradhistory) = res
    (state, envs) = x
    return state, envs, normgradhistory[end]
end


text:
#evolution equation SGPE
function phi_step_sgpe(phi, V, t, dt, lambda, m_b, mu, g, stencil_bc, xstep, T)
    noise = sqrt(T * dt) * sampler_complex(size(phi))
    d = length(size(phi))
    out = similar(phi)
    for x in stencil_bc
        out[x[1]] = (
            phi[x[1]] -
            1.0im *
            dt *
            (
                -1.0 / (2.0 * m_b * xstep^2) *
                (sum(phi[i] for i in x[2]) - 2 * d * phi[x[1]]) +
                g * conj(phi[x[1]]) * phi[x[1]] - mu + lambda * V[1][x[1]]
            ) *
            phi[x[1]] + noise[x[1]]
        )
    end
    @. phi = out
end

#if abspath(PROGRAM_FILE) == @__FILE__3
#    include("helper.jl")
#    using BenchmarkTools
#    phi = ones(100,100)+0im*ones(100,100)
#    V = [ones(100,100)+0im*ones(100,100),ones(100,100)+0im*ones(100,100),
#    ones(100,100)+0im*ones(100,100)]
#    xarray = [i for i in range(-5 / 2.0, 5 / 2.0, length = 100)]
#    V = V_fct(V, [0,0], 0, xarray, 1)
#    stencil_bc = tuple_stencil_bc(phi)
#    @btime phi_step_sgpe(phi,V, 0,0.1, 1, 1, 1, 1, stencil_bc, 0.1)
#end

max_stars_repo_name: phyjonas/ImpuritySGPE

max_issues_repo_name: phyjonas/ImpuritySGPE

max_forks_repo_name: phyjonas/ImpuritySGPE


text:
# Test file for Deferred Acceptance algorithm routine

const file_name = "deferred_acceptance.jl"
const function_name = "deferred_acceptance"

include(file_name)
const fn = getfield(Main, Symbol(function_name))

if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

@testset "Testing $file_name" begin

    @testset "$function_name" begin
        # Males' preference orders over females [1, 2, 3] and unmatched 0
        m_prefs = [1 3 2 3;
                   2 1 3 1;
                   3 2 1 2;
                   0 0 0 0]
        # Females' preference orders over males [1, 2, 3, 4] and unmatched 0
        f_prefs = [3 1 3;
                   1 2 0;
                   2 3 2;
                   4 4 1;
                   0 0 4]

        # Unique stable matching
        m_matched_expected = [1, 2, 3, 0]
        f_matched_expected = [1, 2, 3]

        # Male proposal
        m_matched_computed, f_matched_computed = fn(m_prefs, f_prefs)
        @test m_matched_computed == m_matched_expected
        @test f_matched_computed == f_matched_expected

        # Female proposal
        f_matched_computed, m_matched_computed = fn(f_prefs, m_prefs)
        @test m_matched_computed == m_matched_expected
        @test f_matched_computed == f_matched_expected
    end

end


text:
module ModelHousehold

#
# This module provides a slightly abstracted view of a single Scottish household. The fields
# are predominantly based on the FRS, but could in principle come from other sources
# such as SHS, LCF and so on, or perhaps from a web form.
#
using Base: String
using Dates
import Base.show

using ScottishTaxBenefitModel
using .Definitions
using .Utils:
    has_non_z,
    todays_date,
    to_md_table

using .EquivalenceScales:
    EQ_P_Type,
    EQ_Person,
    EQScales,
    get_equivalence_scales,
    eq_dependent_child,
    eq_head,
    eq_other_adult,
    eq_spouse_of_head

import .EquivalenceScales.get_age
import .EquivalenceScales.eq_rel_to_hoh

import Base.isequal
import Base.hash

using .Uprating: uprate, UPRATE_MAPPINGS

export
    Household,
    BenefitUnit,
    BenefitUnits,
    BUAllocation,
    Person,
    People_Dict,

    adult_pids,
    age_then,
    between_ages,
    child_pids,
    default_bu_allocation,
    empl_status_in,
    ge_age,
    get_benefit_units,
    get_head,
    get_spouse,
    has_carer_member,
    has_children,
    has_disabled_member,
    has_income,
    interview_date,
    is_child,
    is_head,
    is_lone_parent,
    is_severe_disability,
    is_single,
    is_spouse,
    isless,
    le_age,
    make_benefit_unit,
    make_eq_scales!,
    num_adults,
    num_adults,
    num_carers,
    num_children,
    num_people,
    oldest_person,
    on_mt_benefits,
    pers_is_carer,
    pers_is_disabled,
    printpids,
    search,
    set_wage!,
    to_string,
    uprate!

mutable struct Person{RT<:Real}
    hid::BigInt # == sernum
    pid::BigInt # == unique id (year * 100000)+
    pno:: Int # person number in household
    is_hrp :: Bool
    default_benefit_unit:: Int
    is_standard_child :: Bool
    age:: Int

    sex::Sex
    ethnic_group::Ethnic_Group
    marital_status::Marital_Status
    highest_qualification::Qualification_Type

    sic::SIC_2007
    occupational_classification::Standard_Occupational_Classification
    public_or_private :: Employment_Sector
    principal_employment_type :: Employment_Type

    socio_economic_grouping::Socio_Economic_Group
    age_completed_full_time_education:: Int
    years_in_full_time_work:: Int
    employment_status::ILO_Employment
    actual_hours_worked::RT
    usual_hours_worked::RT

    age_started_first_job :: Int

    income::Incomes_Dict{RT}
    benefit_ratios :: Incomes_Dict{RT}

    jsa_type :: JSAType # FIXME change this name
    esa_type :: JSAType
    dla_self_care_type :: LowMiddleHigh
    dla_mobility_type :: LowMiddleHigh
    attendance_allowance_type :: LowMiddleHigh
    pip_daily_living_type :: PIPType
    pip_mobility_type ::  PIPType

    bereavement_type :: BereavementType
    had_children_when_bereaved :: Bool

    assets::Asset_Dict{RT}
    over_20_k_saving :: Bool
    pay_includes ::Included_In_Pay_Dict{Bool}

    # contracted_out_of_serps::Bool

    registered_blind::Bool
    registered_partially_sighted::Bool
    registered_deaf::Bool

    disabilities::Disability_Dict{Bool}

    health_status::Health_Status

    has_long_standing_illness :: Bool
    adls_are_reduced :: ADLS_Inhibited
    how_long_adls_reduced :: Illness_Length

    relationships::Relationship_Dict
    relationship_to_hoh :: Relationship
    is_informal_carer::Bool
    receives_informal_care_from_non_householder::Bool
    hours_of_care_received::RT
    hours_of_care_given::RT
    #
    # Childcare fields; assigned to children
    #
    hours_of_childcare :: RT
    cost_of_childcare :: RT
    childcare_type :: Child_Care_Type
    employer_provides_child_care :: Bool

    company_car_fuel_type :: Fuel_Type
    company_car_value :: RT
    company_car_contribution :: RT
    fuel_supplied :: RT
    onerand :: String
end

People_Dict = Dict{BigInt,Person{T}} where T<:Real
Pid_Array = Vector{BigInt}

mutable struct Household{RT<:Real}
    sequence:: Int # position in current generated dataset
    hid::BigInt
    data_year :: Int
    interview_year:: Int
    interview_month:: Int
    quarter:: Int
    tenure::Tenure_Type
    region::Standard_Region
    ct_band::CT_Band
    council_tax::RT
    water_and_sewerage ::RT
    mortgage_payment::RT
    mortgage_interest::RT
    years_outstanding_on_mortgage:: Int
    mortgage_outstanding::RT
    year_house_bought:: Int
    gross_rent::RT # rentg Gross rent including Housing Benefit  or rent Net amount of last rent payment
    rent_includes_water_and_sewerage::Bool
    other_housing_charges::RT
    gross_housing_costs::RT
    total_income::RT
    total_wealth::RT
    house_value::RT
    weight::RT
    council :: Symbol
    nhs_board :: Symbol
    bedrooms :: Int
    head_of_household :: BigInt
    people::People_Dict{RT}
    onerand :: String
    equivalence_scales :: EQScales{RT}
end

function to_string( hh :: Household ) :: String
    s = to_md_table( hh, exclude=[:people,:onerand])
    for pid in sort(collect(keys(hh.people)))
        s *= to_string( hh.people[pid] )
    end
    return s
end

function to_string( pers :: Person ) :: String
    s = """
    ### Person $(pers.pno)
    """
    s *= to_md_table( pers, depth=1,exclude=[:onerand])
    return s
end

"""
for equivalence scale implicit interface
"""
function get_age( p::Person ) :: Int
    return p.age;
end

"""
for the eq scale implicit interface
"""
function eq_rel_to_hoh( p :: Person ) :: EQ_P_Type
    if p.relationship_to_hoh == This_Person
        return eq_head
    elseif p.is_standard_child
        return eq_dependent_child
    elseif p.relationship_to_hoh in [Spouse,Cohabitee]
        return eq_spouse_of_head
    else
        return eq_other_adult
    end
end

function make_eq_scales!( hh :: Household{T} ) where T
    pl = collect(values(hh.people))
    # sz = size( pl )[1]
    # print( "num people $sz" )
    hh.equivalence_scales = get_equivalence_scales( T, pl )
end

function interview_date( hh :: Household ) :: Date
    return Date( hh.interview_year, hh.interview_month, 15 )
end

function uprate!( pid :: BigInt, year::Integer, quarter::Integer, person :: Person )
    for (t,inc) in person.income
            person.income[t] = uprate( inc, year, quarter, UPRATE_MAPPINGS[t])
    end
    for (a,ass) in person.assets
            person.assets[a] = uprate( ass, year, quarter, UPRATE_MAPPINGS[a])
    end
    person.cost_of_childcare = uprate( person.cost_of_childcare, year, quarter, upr_earnings )
end

function uprate!( hh :: Household )

    hh.water_and_sewerage  = uprate( hh.water_and_sewerage , hh.interview_year, hh.quarter, upr_housing_rents )
    hh.mortgage_payment = uprate( hh.mortgage_payment, hh.interview_year, hh.quarter, upr_housing_oo )
    hh.mortgage_interest = uprate( hh.mortgage_interest, hh.interview_year, hh.quarter, upr_housing_oo )
    hh.mortgage_outstanding = uprate( hh.mortgage_outstanding, hh.interview_year, hh.quarter, upr_housing_oo )
    hh.gross_rent = uprate( hh.gross_rent, hh.interview_year, hh.quarter, upr_housing_rents )
    hh.other_housing_charges = uprate( hh.other_housing_charges, hh.interview_year, hh.quarter, upr_nominal_gdp )
    hh.gross_housing_costs = uprate( hh.gross_housing_costs, hh.interview_year, hh.quarter, upr_nominal_gdp )
    hh.total_income = uprate( hh.total_income, hh.interview_year, hh.quarter, upr_nominal_gdp )
    hh.total_wealth = uprate( hh.total_wealth, hh.interview_year, hh.quarter, upr_nominal_gdp )
    hh.house_value = uprate( hh.house_value, hh.interview_year, hh.quarter, upr_housing_oo )
    for (pid,person) in hh.people
        uprate!( pid, hh.interview_year, hh.quarter, person )
    end

end

function oldest_person( people :: People_Dict ) :: NamedTuple
    oldest = ( age=-999, pid=BigInt(0))
    for person in people
        if person.age > oldest.age
            oldest.age = person.age
            oldest.pid = person.pid
        end
    end
    oldest
end

PeopleArray = Vector{Person}

struct BenefitUnit{T<:Real}
    people :: People_Dict{T}
    head :: BigInt
    spouse :: BigInt
    adults :: Pid_Array
    children :: Pid_Array
end

BenefitUnits = Vector{BenefitUnit}
BUAllocation = Vector{PeopleArray}

function num_children( bu :: BenefitUnit ) :: Integer
    size( bu.children )[1]
end

function num_children( hh :: Household ) :: Integer
    count( hh, is_child )
end

"""
For Budget constraints: set the wage and make all the needed
corrections to employment status, hours worked, pensions(? pensions not done yet)
"""
function set_wage!( pers :: Person, gross :: Real, wage :: Real )
    pers.income[wages] = gross
    h = gross/wage
    pers.usual_hours_worked = h
    pers.employment_status = if h < 5
        Unemployed
    elseif h < 30
        Part_time_Employee
    else
        Full_time_Employee
    end
    # println( "made usual_hours_worked = $h gross=$gross wage=$wage pers.employment_status=$(pers.employment_status)")
end

#
# This creates a array of references to each person in the houshold, broken into
# benefit units using the default FRS/EFS benefit unit number.
#
function default_bu_allocation( hh :: Household ) :: BUAllocation
    bua = BUAllocation()
    nbus = 0
    # how many bus
    for (pid,person) in hh.people
        nbus = max( nbus, person.default_benefit_unit )
    end
    ## create benefit units
    for buno in 1:nbus
        push!( bua, PeopleArray())
    end
    sz = size( bua )[1]
    ## allocate peobuallocationple to them
    for (pid,person) in hh.people
        pp = bua[person.default_benefit_unit]
        push!( pp, person )
    end
    ## sort people in each by pid
    for buno in 1:nbus
        sort!( bua[buno], lt=(left,right)->isless(right.age,left.age))
    end
    bua
end

"""

Create a benefit unit from an array of people
really only needed for some unit tests.

FIXME: don't use outside of tests! likely mess up children and spouses.

"""
function make_benefit_unit(
    people :: PeopleArray,
    head :: BigInt,
    spouse:: BigInt = -1 ) :: BenefitUnit
    npeople = size( people )[1]
    T = typeof( people[1].hours_of_care_received )
    pd = People_Dict{T}()
    children = Pid_Array()
    adults = Pid_Array()
    push!( adults, head )
    if spouse > 0
        push!( adults, spouse )
    end
    palloc = spouse <= 0 ? 2 : 3
    for n in palloc:npeople
        push!(children, people[n].pid)
        @assert people[n].age <= 21 # vague idiot check
    end
    for pers in people
        pd[pers.pid] = pers
    end
    return BenefitUnit( pd, head, spouse, adults, children )
end

function child_pids( hh :: Household ) :: Pid_Array
    chlds = Pid_Array()
    for (pid,pers) in hh.people
        if is_child( pers )
            push!( chlds, pid )
        end
    end
    return chlds
end

#
# This creates a array of references to each person in the houshold, broken into
# benefit units using the default FRS/EFS benefit unit number.
#
function allocate_to_bus( T::Type, bua :: BUAllocation ) :: BenefitUnits
    nbus = size(bua)[1]
    bus = BenefitUnits(undef, nbus)
    for i in 1:nbus
        people = People_Dict{T}()
        head_pid :: BigInt = -1
        spouse_pid :: BigInt = -1
        children = Pid_Array()
        adults = Pid_Array()
        npeople = size( bua[i])[1]
        for p in 1:npeople
            person = bua[i][p]
            people[person.pid] = person
            if p == 1
                head_pid = person.pid
                push!( adults, head_pid )
            else
                # println( "on bu $i person $p relationships $(person.relationships)")
                reltohead = person.relationships[head_pid]
                if reltohead in [Spouse,Cohabitee,Civil_Partner]
                    spouse_pid = person.pid
                    push!( adults, spouse_pid )
                    # FIXME we need to remove these checks if
                    # we're using a non-default allocation to bus
                    # @assert person.age >= 16
                else
                    # @assert person.age <= 19
                    push!( children, person.pid )
                end
            end
        end
        new_bu = BenefitUnit( people, head_pid, spouse_pid, adults, children )
        bus[i] = new_bu
    end
    bus
end

function get_benefit_units(
    hh :: Household{T},
    allocator :: Function=default_bu_allocation ) :: BenefitUnits where T
    allocate_to_bus( T, allocator(hh))
end

function num_people( bu :: BenefitUnit )::Integer
    length( bu.people )
end

function num_adults( bu :: BenefitUnit )::Integer
    size( bu.adults )[1]
end

function adult_pids( hh :: Household )::Pid_Array
    pids = Pid_Array()
    for (pid,pers) in hh.people
        if ! is_child( pers )
            push!( pids, pid )
        end
    end
    return pids
end

function num_adults( hh )
    return size( adult_pids(hh))[1]
end

function num_people( hh :: Household ) :: Integer
    length( hh.people )
end

function get_head( bu :: BenefitUnit )::Person
    bu.people[bu.head]
end

function get_head( hh :: Household ) :: Person
    return hh.people[hh.head_of_household]
end

function is_head( bu :: BenefitUnit, pers :: Person  ) :: Bool
    return bu.head == pers.pid
end

function is_spouse( unit, pers :: Person  ) :: Bool
    head = get_head(unit)
    p = partner_of( head )
    if p === nothing
        return false
    end
    return p == pers.pid
end

function partner_of( pers :: Person ) :: Union{Nothing,BigInt}
    for (pid,rel) in pers.relationships
        if rel in [Spouse,Cohabitee]
            return pid
        end
    end
    return nothing
end

function get_spouse( hh :: Household ) ::Union{Nothing,Person}
    head = get_head( hh )
    p = partner_of( head )
    if p !== nothing
        return hh.people[p]
    end
    return nothing
end

function get_spouse( bu :: BenefitUnit )::Union{Nothing,Person}
    if bu.spouse <= 0
        return nothing
    end
    bu.people[bu.spouse]
end

function has_children( bu :: BenefitUnit ) :: Bool
    num_children( bu ) > 0
end

"""
relies on `is_child` being sensible
"""
function has_children( hh :: Household )::Bool
    for (pid, pers ) in hh.people
        if is_child( pers )
            return true
        end
    end
    return false
end

function is_lone_parent( bu :: BenefitUnit ) :: Bool
    return bu.spouse < 0 && has_children( bu )
end

function is_child( pers :: Person )
    # FIXME crude
    pers.is_standard_child
    # below won't work for hh  since
    # wouldn't be child if in own BU & we can't check that with 1 person
    #if pers.age <= 15
    #    return true
    #end
    #(pers.age <= 19) && (pers.employment_status in [Student])
end

#
# fixme just count people???
#
function is_lone_parent( hh :: Household,
    allocator :: Function=default_bu_allocation ) :: Bool

    bus = get_benefit_units( hh, allocator )
    if size(bus)[1] > 1
        return false
    end
    return bus[1].spouse < 0 && size( bus[1].children )[1] > 0
end

function is_single( bu :: BenefitUnit ) :: Bool
    num_people( bu ) == 1
end

function is_single( hh :: Household ) :: Bool
    num_people( hh ) == 1
end

"""
FIXME we're going to do this solely on benefit receipt/is_informal_carer flag
for now until we get the regressions done
we use historic benefits here since this uses actual data
"""
function pers_is_carer( pers :: Person, params ... ) :: Bool
    if has_non_z( pers.income, carers_allowance )
        return true
    end
    # FIXME check/parameterise this
    if (pers.is_informal_carer) && (pers.hours_of_care_given > 0)
        return true
    end
    return false
end

const DISABLE_BENEFITS = [
    severe_disability_allowance,
    attendance_allowance,
    incapacity_benefit,
    dlaself_care,
    dlamobility,
    personal_independence_payment_daily_living,
    personal_independence_payment_mobility]

"""

"""
function is_severe_disability( pers :: Person )
    return pers_is_disabled( pers ) ||
           pers.dla_self_care_type in [mid, high] ||
           pers.dla_mobility_type == high ||
           pers.attendance_allowance_type == high ||
           pers.pip_daily_living_type == enhanced_pip ||
           pers.pip_mobility_type == enhanced_pip ||
           pers.adls_are_reduced == reduced_a_lot ||
           pers.health_status == Very_Bad

end

"""
FIXME we're going to do this solely on benefit receipt
for now until we get the regressions done
we use historic benefits here since this uses actual data
"""
function pers_is_disabled( pers :: Person, params ... ) :: Bool
    if pers.registered_blind
        return true
    end
    if pers.employment_status in [Permanently_sick_or_disabled]
        return true
    end
    if pers.adls_are_reduced == reduced_a_lot
        return true
    end
    for dis in [mobility,
        dexterity,
        learning,
        memory,
        mental_health]
        if haskey(pers.disabilities,dis) && pers.disabilities[dis]
            return true
        end
    end

    for k in DISABLE_BENEFITS
        if has_non_z( pers.income, k )
            return true
        end
    end # loop
    return false
end

function empl_status_in( pers :: Person, statuses ...)
    return pers.employment_status in statuses
end

le_age( pers :: Person, age ... ) = pers.age <= age[1]

ge_age( pers :: Person, age ... ) = pers.age >= age[1]

between_ages( pers :: Person, age ... ) = age[1] >= pers.age <= age[2]

has_income( pers::Person, which :: Incomes_Type ) = haskey( pers.income, which )

function has_income( pers::Person, which ... )::Bool
    for k in which
        if has_income( pers, k )
            return true
        end
    end
    return false
end

function search( people :: People_Dict, func :: Function, params... ) :: Bool
    for ( pid,pers ) in people
        if func( pers, params... )
            return true
        end
    end
    return false
end

function count( people :: People_Dict, func :: Function, params... ) :: Int
    n = 0
    for ( pid,pers ) in people
        if func( pers, params... )
            n += 1
        end
    end
    return n
end

function search( bu :: BenefitUnit, func :: Function, params... ) :: Bool
    return search( bu.people, func, params... )
end

function search( hh :: Household, func :: Function, params ... ) :: Bool
    return search( hh.people, func, params... )
end

function count( bu :: BenefitUnit, func :: Function, params... ) :: Integer
     return count( bu.people, func, params... )
end

function count( hh :: Household, func :: Function, params ... ) :: Integer
    return count( hh.people, func, params... )
end

function on_mt_benefits( pers :: Person ) :: Bool
    if has_income( pers, [income_support, working_tax_credit, child_tax_credit, housing_benefit, universal_credit]...)
        return true
    end
    if pers.jsa_type == income_related_jsa
        return true
    end
    if pers.esa_type == income_related_jsa
        return true
    end
    return false
end

function on_mt_benefits( bu :: BenefitUnit ) :: Bool
    return search( bu, on_mt_benefits )
end

function on_mt_benefits( hh :: Household ) :: Bool
    return search( hh, on_mt_benefits )
end

function has_disabled_member( bu :: BenefitUnit ) :: Bool
    search( bu.people, pers_is_disabled )
end

function has_disabled_member( hh :: Household ) :: Bool
    search( hh.people, pers_is_disabled )
end

function has_carer_member( bu :: BenefitUnit ) :: Bool
    search( bu.people, pers_is_carer )
end

function has_carer_member( hh :: Household ) :: Bool
    search( hh.people, pers_is_carer )
end

function num_carers( bu :: BenefitUnit ) :: Int
    count( bu.people, pers_is_carer )
end

function carers( hh :: Household ) :: Int
    count( hh.people, pers_is_carer )
end

# simple diagnostic prints for testing allocation

function printpids( pers :: Person )
    println( "pid: $(pers.pid) age $(pers.age) sex $(pers.sex) default_benefit_unit $(pers.default_benefit_unit) ")
end

function printpids( bu::BenefitUnit)
      head = get_head( bu )
      print( "HEAD: ")
      printpids( head )
      spouse = get_spouse( bu )
      if spouse !== nothing
            print( "SPOUSE: ")
            printpids( spouse )

      end
      for chno in bu.children
            print( "CHILD_$(chno):")
            child = bu.people[chno]
            printpids( child )
      end
end

function printpids( people :: People_Dict )
    ks = sort(collect(keys(people)))
    for k in ks
        print( "PERSON_$k")
        printpids( people[k] )
    end
end

function printpids( buas::BUAllocation )
    nbus = size( buas )[1]
    for i in 1:nbus
        npeople = size( buas[i] )[1]
        for p in 1:npeople
            print( "PERS_$(p): ")
            printpids( buas[i][p])
        end
    end
end

function age_then( pers :: Person, when :: Int ) :: Int
    Utils.age_then( pers.age, when );
end

function age_then( pers :: Person, when :: Date = todays_date() ) :: Int
    Utils.age_then( pers.age, Dates.year( when ))
end

end # module

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
function test_mapped_dicts()
    d = UnitCircle()
    m = parameterization(d)
    b = Fourier(20)
    c = mapped_dict(b, m)
    f1(x,y) = exp(x+im*y)
    u = approximate(c, f1)
    @test in_support(c, 1, SVector(1.0, 0.0))
    @test abs(u(SVector(cos(10),sin(10)))-f1(cos(10),sin(10))) < 1e-6
end


text:
using Revise
using ADSeismic
using ADCME
using JLD2
using Test

param = ElasticPropagatorParams(NX=10, NY=10, NSTEP=10, DELTAT=1e-4,
    DELTAX=1.0, DELTAY=1.0)
nv = 10

idx1 = (1:5)
idx2 = (6:10)
V = zeros(param.NSTEP+1, param.NX, param.NY)
@load "data/dat0-4.jld2" Vx param
param00 = param
V[:, idx1, idx1] = Vx
@load "data/dat1-4.jld2" Vx param
param01 = param
V[:, idx1, idx2] = Vx
@load "data/dat2-4.jld2" Vx param
param10 = param
V[:, idx2, idx1] = Vx
@load "data/dat3-4.jld2" Vx param
param11 = param
V[:, idx2, idx2] = Vx

V = V[nv,:,:]
@load "data/dat0-1.jld2" Vx param
Vx = Vx[nv,:,:]

@test V≈Vx

V = zeros(param.NSTEP+1, param.NX, param.NY)
@load "data/dat0-4.jld2" Vy
V[:, idx1, idx1] = Vy
@load "data/dat1-4.jld2" Vy
V[:, idx1, idx2] = Vy
@load "data/dat2-4.jld2" Vy
V[:, idx2, idx1] = Vy
@load "data/dat3-4.jld2" Vy
V[:, idx2, idx2] = Vy
V = V[nv,:,:]
@load "data/dat0-1.jld2" Vy
Vy = Vy[nv,:,:]
@test V≈Vy

V = zeros(param.NSTEP+1, param.NX, param.NY)
@load "data/dat0-4.jld2" Sxx
V[:, idx1, idx1] = Sxx
@load "data/dat1-4.jld2" Sxx
V[:, idx1, idx2] = Sxx
@load "data/dat2-4.jld2" Sxx
V[:, idx2, idx1] = Sxx
@load "data/dat3-4.jld2" Sxx
V[:, idx2, idx2] = Sxx
V = V[nv,:,:]
@load "data/dat0-1.jld2" Sxx
Sxx = Sxx[nv,:,:]
@test V≈Sxx

V = zeros(param.NSTEP+1, param.NX, param.NY)
@load "data/dat0-4.jld2" Syy
V[:, idx1, idx1] = Syy
@load "data/dat1-4.jld2" Syy
V[:, idx1, idx2] = Syy
@load "data/dat2-4.jld2" Syy
V[:, idx2, idx1] = Syy
@load "data/dat3-4.jld2" Syy
V[:, idx2, idx2] = Syy
V = V[nv,:,:]
@load "data/dat0-1.jld2" Syy
Syy = Syy[nv,:,:]
@test V≈Syy

V = zeros(param.NSTEP+1, param.NX, param.NY)
@load "data/dat0-4.jld2" Sxy
V[:, idx1, idx1] = Sxy
@load "data/dat1-4.jld2" Sxy
V[:, idx1, idx2] = Sxy
@load "data/dat2-4.jld2" Sxy
V[:, idx2, idx1] = Sxy
@load "data/dat3-4.jld2" Sxy
V[:, idx2, idx2] = Sxy
V = V[nv,:,:]
@load "data/dat0-1.jld2" Sxy
Sxy = Sxy[nv,:,:]
@test V≈Sxy

# close("all")
# visualize_wavefield(V, param)


text:
include("board.jl")
module dParse

using BM
global DB
using SQLite

function Parse(pathToDatabase)

  #captures = readdml("captures.txt")
  db = SQLite.DB(pathToDatabase) #Opens the database gamefile
  #calculates all black piece positions
  res = SQLite.query(db,"SELECT key,value FROM meta where key = 'type';")
  global gameType = get(res[1,2])
  res = SQLite.query(db,"SELECT key,value FROM meta where key = 'seed';")
  global seed = parse(get(res[1,2]))
board = BM.startGame(gameType)
  #=----Gets last move number----=#
  res = SQLite.query(db,"SELECT MAX(move_number) FROM moves;") #Finds the last played move (maximum move_number)
  lastMoveIDNullable = res[1,1] #SQL query with max move_number (POSSIBLY "NULL" if no moves have been made)

  if (!isnull(lastMoveIDNullable)) #Checks that lastMoveID was not NULL
    lastMoveID = parse(get(res[1,1])) #Parses the last move move_number as an Int

  else #lastMoveID is NULL
    lastMoveID = 0 #move_number "0" is unsused. This implies no moves have been made

  end
  currentMoveID=1
    #=----Replays the game until move_id = lastMoveID----=#
    while currentMoveID<=lastMoveID
      res = SQLite.query(db,"SELECT sourcex,sourcey,targetx,targety,move_type,option FROM moves WHERE move_number = $(currentMoveID);")
      sourcexNullable = res[1][1]
      sourceyNullable = res[1][2]
      targetxNullable = res[1][3]
      targetyNullable = res[1][4]
      move_type = get(res[1][5])
      optionNullable = res[1][6]

      if move_type == "move" #Regular Move
        targetx = get(targetxNullable)
        targety = get(targetyNullable)
        sourcex = get(sourcexNullable)
        sourcey = get(sourceyNullable)
        if ( (targetx,targety) in keys(board))# Non-empty piece at target: capture
          #push(captures,board[targetx][targety])
        end

        board = BM.nextState(board,(sourcex,sourcey,targetx,targety))

      elseif move_type == "drop"
        option = get(optionNullable)
        try
          #deleteat!(captures,findfirst(x->x.piece==option))
        end
        targetx = get(targetxNullable)
        targety = get(targetyNullable)
        board[targetx][targety].piece = option
        if iseven(currentMoveID)
         board[(targetx,targety)] = (option,'w')
        else
          board[(targetx,targety)] = (option,'b')
        end
        newBoard.currentPlayer = (newBoard.currentPlayer == 'b' ? 'w' : 'b')
      elseif move_type == "resign"
        #Do nothing
      end

    end

  return (board,gameType,seed,lastMoveID)
end #End function

#End module
end

max_stars_repo_name: kinoroy/Shogi

max_issues_repo_name: kinoroy/Shulia

max_forks_repo_name: kinoroy/Shulia


text:
# Test file for Deferred Acceptance algorithm routine

const file_name = "my_Gale_Shap.jl"
const function_name = "my_Gale_Shap"

include(file_name)
const fn = getfield(Main, Symbol(function_name))

if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

@testset "Testing $file_name" begin

    @testset "$function_name: one-to-one" begin
        # Males' preference orders over females [1, 2, 3] and unmatched 0
        m_prefs = [1 3 2 3;
                   2 1 3 1;
                   3 2 1 2;
                   0 0 0 0]
        # Females' preference orders over males [1, 2, 3, 4] and unmatched 0
        f_prefs = [3 1 3;
                   1 2 0;
                   2 3 2;
                   4 4 1;
                   0 0 4]

        # Unique stable matching
        m_matched_expected = [1, 2, 3, 0]
        f_matched_expected = [1, 2, 3]

        # Male proposal
        m_matched_computed, f_matched_computed = fn(m_prefs, f_prefs)
        @test m_matched_computed == m_matched_expected
        @test f_matched_computed == f_matched_expected

        # Female proposal
        f_matched_computed, m_matched_computed = fn(f_prefs, m_prefs)
        @test m_matched_computed == m_matched_expected
        @test f_matched_computed == f_matched_expected
    end

    @testset "$function_name: many-to-one with caps 1" begin
        # Males' preference orders over females [1, 2, 3] and unmatched 0
        m_prefs = [1 3 2 3;
                   2 1 3 1;
                   3 2 1 2;
                   0 0 0 0]
        # Females' preference orders over males [1, 2, 3, 4] and unmatched 0
        f_prefs = [3 1 3;
                   1 2 0;
                   2 3 2;
                   4 4 1;
                   0 0 4]

        # Capacities for females
        caps = [1, 1, 1]
        indptr_expected = [1, 2, 3, 4]

        # Unique stable matching
        m_matched_expected = [1, 2, 3, 0]
        f_matched_expected = [1, 2, 3]

        # Male proposal
        m_matched_computed, f_matched_computed, indptr_computed =
            fn(m_prefs, f_prefs, caps)
        @test m_matched_computed == m_matched_expected
        @test f_matched_computed == f_matched_expected
        @test indptr_computed == indptr_expected
    end

    @testset "$function_name: many-to-one" begin
        # From http://www.columbia.edu/~js1353/pubs/qst-many-to-one.pdf
        # Originally from Gusfield and Irving (1989, Section 1.6.5)

        m = 11  # Number of students
        n = 5   # Number of colleges

        # Students' preference orders over colleges 1, ..., 5 and unmatched 0
        s_prefs = [3, 1, 5, 4, 0, 2,
                   1, 3, 4, 2, 5, 0,
                   4, 5, 3, 1, 2, 0,
                   3, 4, 1, 5, 0, 2,
                   1, 4, 2, 0, 3, 5,
                   4, 3, 2, 1, 5, 0,
                   2, 5, 1, 3, 0, 4,
                   1, 3, 2, 5, 4, 0,
                   4, 1, 5, 0, 2, 3,
                   3, 1, 5, 2, 4, 0,
                   5, 4, 1, 3, 2, 0]
        s_prefs = reshape(s_prefs, n+1, m)

        # Colleges' preference orders over students 1, ..., 11 and unmatched 0
        c_prefs = [3, 7, 9, 11, 5, 4, 10, 8, 6, 1,
                   2, 0, 5, 7, 10, 6, 8, 2, 3, 11,
                   0, 1, 4, 9, 11, 6, 8, 3, 2, 4,
                   7, 1, 10, 0, 5, 9, 10, 1, 2, 11,
                   4, 9, 5, 3, 6, 8, 0, 7, 2, 4,
                   10, 7, 6, 1, 8, 3, 11, 9, 0, 5]
        c_prefs = reshape(c_prefs, m+1, n)

        # Capacities for colleges
        caps = [4, 1, 3, 2, 1]
        indptr_expected = [1, 5, 6, 9, 11, 12]

        # Unique stable matching
        s_matched_expected = [3, 1, 4, 3, 1, 3, 2, 1, 4, 1, 5]
        c_matched_expected = [2, 5, 8, 10, 7, 1, 4, 6, 3, 9, 11]

        # Male proposal
        s_matched_computed, c_matched_computed, indptr_computed =
            fn(s_prefs, c_prefs, caps)
        @test s_matched_computed == s_matched_expected
        @test indptr_computed == indptr_expected

        # Sort matched students for each college
        for j in 1:n
            sort!(sub(c_matched_computed,
                      indptr_expected[j]:indptr_expected[j+1]-1)
            )
        end
        @test c_matched_computed == c_matched_expected
    end

end

max_stars_repo_name: SUZUKITAISHI/matching

max_issues_repo_name: SUZUKITAISHI/matching

max_forks_repo_name: SUZUKITAISHI/matching


text:
function AlignmentDetector(arg0::Orbit, arg1::PVCoordinatesProvider, arg2::jdouble)
    return AlignmentDetector((Orbit, PVCoordinatesProvider, jdouble), arg0, arg1, arg2)
end

function AlignmentDetector(arg0::jdouble, arg1::Orbit, arg2::PVCoordinatesProvider, arg3::jdouble)
    return AlignmentDetector((jdouble, Orbit, PVCoordinatesProvider, jdouble), arg0, arg1, arg2, arg3)
end

function AlignmentDetector(arg0::jdouble, arg1::jdouble, arg2::PVCoordinatesProvider, arg3::jdouble)
    return AlignmentDetector((jdouble, jdouble, PVCoordinatesProvider, jdouble), arg0, arg1, arg2, arg3)
end

function g(obj::AlignmentDetector, arg0::SpacecraftState)
    return jcall(obj, "g", jdouble, (SpacecraftState,), arg0)
end

function get_align_angle(obj::AlignmentDetector)
    return jcall(obj, "getAlignAngle", jdouble, ())
end

function get_pv_coordinates_provider(obj::AlignmentDetector)
    return jcall(obj, "getPVCoordinatesProvider", PVCoordinatesProvider, ())
end


text:
using ConfParser
using Base.Test

conf = ConfParse("confs/config.ini")
parse_conf!(conf)

@test retrieve(conf, "database", "user") == "dbuser"
@test retrieve(conf, "database", "password") == "abc123"
@test commit!(conf, "default", "database", "newuser") == true
@test erase!(conf, "foobarness") == true
save!(conf, "confs/out.conf")

# validate new conf write
conf = ConfParse("confs/out.conf")
parse_conf!(conf)
@test retrieve(conf, "database", "user") == "dbuser"
@test retrieve(conf, "default", "database") == "newuser"


text:
using DataStructures, OPFInstances
using JuMP, KNITRO
using BenchmarkTools#, ProfileView

!isdefined(:MPC) && (const MPC = MathProgComplex)

function main()
    tic()
    # problem_c, pt = MPC.import_from_dat(getinstancepath("Matpower", "QCQP", "case1354pegase"))
    problem_c, pt = MPC.import_from_dat(getinstancepath("Matpower", "QCQP", "case13659pegase"))
    #problem_c, pt = import_from_dat(getinstancepath("Matpower", "QCQP", "case9"))
    println("import_from_dat : ", toq())
    tic()
    problem = MPC.pb_cplx2real(problem_c)
    # @profile MPC.pb_cplx2real(problem_c)
    # open("prof.txt", "w") do s
    #      Profile.print(IOContext(s, :displaysize => (24, 500)), format=:flat)
    #      # Profile.print(IOContext(s))
    #  end
    println("pb_cplx2real : ", toq())
    return

    # problem_c, pt = MPC.import_from_dat(getinstancepath("Matpower", "QCQP", "case1354pegase"))
    # #problem_c, pt = import_from_dat(getinstancepath("Matpower", "QCQP", "case9"))
    # println("import_from_dat : ")
    # @btim²e MPC.import_from_dat($(getinstancepath("Matpower", "QCQP", "case1354pegase")))

    # # @profile MPC.pb_cplx2real(problem_c)
    # # ProfileView.view()
    # # open("prof.txt", "w") do s
    # #      Profile.print(IOContext(s, :displaysize => (24, 500)), format=:flat)
    # #      # Profile.print(IOContext(s))
    # #  end

    # problem = MPC.pb_cplx2real(problem_c)
    # println("pb_cplx2real : ")
    # @btime MPC.pb_cplx2real($problem_c)
    # println()
    # tic()
    #
    # solver = KnitroSolver(KTR_PARAM_OUTLEV=3,
    #                     KTR_PARAM_MAXIT=600,
    #                     KTR_PARAM_SCALE=0,
    #                     KTR_PARAM_FEASTOL=1.0,
    #                     KTR_PARAM_OPTTOL=1.0,
    #                     KTR_PARAM_FEASTOLABS=1e-6,
    #                     KTR_PARAM_OPTTOLABS=1e-3,
    #                     KTR_PARAM_BAR_INITPT=2,
    #                     KTR_PARAM_PRESOLVE=0,
    #                     KTR_PARAM_HONORBNDS=0)
    #
    # println("KnitroSolver : ", toq())
    # tic()
    # m, JuMPvars = get_JuMP_cartesian_model(problem, solver)
    # println("get_JuMP_cartesian_model : ", toq())
    # tic()
    # solve(m)
    # println("solve : ", toq())
    #
    # println("Objective value : $(getobjectivevalue(m))\n")

   # for (varname, var) in JuMPvars
   #   value = getvalue(var)
   #   println("$varname   $value")
   # end
   #
   #
   #  amplexportpath = joinpath(Pkg.dir("MathProgComplex"), "Knitro_runs")
   #  point = Point()
   #
   #  export_to_dat(problem, amplexportpath)
   #  run_knitro(amplexportpath)

end

main()

# Timings, case9
# Without Manuel's update to order.jl, Exponent.jl

# import_from_dat time:
#   9.342 ms (30396 allocations: 1.27 MiB)
# pb_cplx2real time:
#   38.879 ms (370727 allocations: 16.47 MiB)
# get_JuMP_cartesian_model time:
#   892.305 μs (8573 allocations: 429.27 KiB)

# import_from_dat time:
#   7.991 ms (18790 allocations: 756.09 KiB)
# pb_cplx2real time:
#   16.707 ms (136490 allocations: 5.73 MiB)
# get_JuMP_cartesian_model time:
#   874.254 μs (8574 allocations: 429.52 KiB)

#   case300

#   import_from_dat time:
#   331.475 ms (764052 allocations: 29.46 MiB)
# pb_cplx2real time:
#   930.399 ms (6244658 allocations: 255.18 MiB)
# get_JuMP_cartesian_model time:
#   52.450 ms (395145 allocations: 19.43 MiB)

#   import_from_dat time:
#   457.236 ms (1841742 allocations: 79.81 MiB)
# pb_cplx2real time:
#   2.460 s (20014923 allocations: 889.81 MiB)
# get_JuMP_cartesian_model time:
#   54.524 ms (395145 allocations: 19.43 MiB)


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    # Make sure to exercise the interpreter
    interp = ASTInterpreter.enter_call_expr(nothing,Expr(:call, Base.gcd, 21,  7))
    ASTInterpreter.finish!(interp, recursive=true)
    precompile(ASTInterpreter.RunDebugREPL,(ASTInterpreter.Interpreter,))
    nothing
end


text:
#=
# The BUGS Book
# Example 2.1.2, page 17:
# For
#   Y ~ Binomial(0.5, 8)
# we want to know Pr(Y <= 2)
#

  See ~/jags/bugs_book_2_1_2.jags
       Mean     SD Naive SE Time-series SE
  P2 0.1442 0.3513 0.001434       0.001434
  Y  4.0025 1.4177 0.005788       0.005839

  ~/webppl/bugs_book_2_1_2.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function bugs_book_2_1_2()

    y ~ Binomial(8,0.5)
    p2 ~ y <= 2 ? flip(1.0) : flip(0.0);

    return y, p2

end

model = bugs_book_2_1_2()

num_chains = 4
# chains = sample(model, Prior(), MCMCThreads(), 10_000, num_chains)

# chains = sample(model, MH(), MCMCThreads(), 40_000, num_chains)

# chains = sample(model, PG(15), MCMCThreads(), 10_000, num_chains)

chains = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chains)

# chains = sample(model, IS(), MCMCThreads(), 10_000, num_chains)

display(chains)
# display(plot(chains))

gen = generated_quantities(model, chains)
show_var_dist_pct(gen, 120)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
export POSTagger

using ..Units
using TextModels
using BSON

MISSING_TAG = "-MISSING-"

struct POSTagger <: AbstractPOSTagger
  pos_model::TextModels.PerceptronTagger

  function POSTagger(model_file)
    tagger = TextModels.PerceptronTagger();
    pretrained = BSON.load(model_file)

    tagger.model.weights = pretrained[:weights]
    tagger.tagdict = pretrained[:tagdict]
    tagger.classes = tagger.model.classes = Set(pretrained[:classes])
    println("loaded successfully")

    new(tagger)
  end
end

function (tagger::POSTagger)(sentence_tokens::Vector{String})
  map(tagger.pos_model(sentence_tokens)) do (token, tag)
    tag = tag === missing ? MISSING_TAG : tag

    (token, tag)
  end |> vec -> convert(Vector{Tuple{String, String}}, vec)
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    LeaveBallOut(ball; loss=Dict())

Leave-`ball`-out (a.k.a. spatial leave-one-out) validation.
Optionally, specify `loss` function from `LossFunctions.jl`
for some of the variables.

    LeaveBallOut(radius; loss=Dict())

By default, use Euclidean ball of given `radius` in space.

## References

* Le Rest et al. 2014. [Spatial leave-one-out cross-validation
  for variable selection in the presence of spatial autocorrelation]
  (https://onlinelibrary.wiley.com/doi/full/10.1111/geb.12161)
"""
struct LeaveBallOut{B<:BallNeighborhood} <: AbstractErrorEstimator
  ball::B
  loss::Dict{Symbol,SupervisedLoss}
end

LeaveBallOut(ball::BallNeighborhood; loss=Dict()) =
  LeaveBallOut{typeof(ball)}(ball, loss)

LeaveBallOut(radius::Number; loss=Dict()) =
  LeaveBallOut(BallNeighborhood(radius), loss=loss)

function error(solver::AbstractLearningSolver,
               problem::LearningProblem,
               eestimator::LeaveBallOut)
  sdata = sourcedata(problem)
  ovars = outputvars(task(problem))
  ball  = eestimator.ball
  loss  = eestimator.loss
  for var in ovars
    if var ∉ keys(loss)
      loss[var] = defaultloss(sdata[var][1])
    end
  end

  # efficient neighborhood search
  searcher = NeighborhoodSearcher(sdata, ball)

  # number of folds
  nfolds = nelms(sdata)

  # error for a ball b
  function ε(b)
    # points inside and outside ball
    coords  = coordinates(sdata, b)
    inside  = search(coords, searcher)
    outside = setdiff(1:nelms(sdata), inside)

    # source and target indices
    sinds = outside
    tinds = [b]

    # source and target data
    train = view(sdata, sinds)
    hold  = view(sdata, tinds)

    # setup and solve sub-problem
    subproblem = LearningProblem(train, hold, task(problem))
    solution   = solve(subproblem, solver)

    # loss for each variable
    losses = map(ovars) do var
      y = sdata[var][b]
      ŷ = solution[var][1]
      var => value(loss[var], y, ŷ)
    end

    Dict(losses)
  end

  # compute error for each fold in parallel
  εs = foldxt(vcat, Map(ε), 1:nfolds)

  # combine error from different folds
  Dict(var => mean(get.(εs, var, 0)) for var in ovars)
end


text:
import Lazy
using Base.Meta
using Calculus

#Operator overloading function for AD types
const B_OP_START = 1
const OP = Array{Symbol,1}([:+,:-,:*,:/,:^])
const B_OP_END = length(OP)
const U_OP_START = B_OP_END + 1

const S_TO_OC = Dict{Symbol,Int}()
const S_TO_DIFF = Dict{Symbol,Array{Any,1}}()
const S_TO_DIFF_FLAG = Dict{Symbol,Array{Bool,1}}()

##########################################################################################
#
# Differetitate rule building using Calculus
#   https://github.com/johnmyleswhite/Calculus.jl/blob/master/src/differentiate.jl
#
##########################################################################################
import Calculus.differentiate

function differentiate(x::SymbolParameter{:abs},args,wrt)
    x = args[1]
    xp = differentiate(x,wrt)
    return :($x>0?$xp:-$xp)
end

for sym in OP
    (dx,dy) = differentiate("x $(sym) y",[:x,:y])
    # @show dx, dy
    dxx = differentiate(dx,:x)
    dxxi = isa(dxx,Number)
    # @show dxx
    dxy = differentiate(dx,:y)
    dxyi = isa(dxy,Number)
    # @show dxy
    dyy = differentiate(dy,:y)
    dyyi = isa(dyy,Number)
    # @show dyy
    # @show sym,dxx,dxy,dyy
    # @show sym,dxxi,dxyi,dyyi
    S_TO_DIFF[sym] = [dx,dy,dxx,dxy,dyy]
    S_TO_DIFF_FLAG[sym] = [dxxi,dxyi,dyyi]
end

for (sym, dx) = symbolic_derivatives_1arg()
    push!(OP,sym)
    try
        dxx = differentiate(dx,:x)
        dxxi = isa(dxx,Number)
        S_TO_DIFF[sym] = [dx,dxx]
        S_TO_DIFF_FLAG[sym] = [dxxi]
    catch e
        pop!(OP)
        warn(e)
    end
end

const U_OP_END = length(OP)

##########################################################################################
#
# Operator overload for types to build tape
#   tape statistic is not initialized
#
##########################################################################################
for oc = B_OP_START:B_OP_END
    o = OP[oc]
    S_TO_OC[o] = oc
    @eval   ($o){I}(l::AD{I},r::AD{I}) =
            begin
                return AD_O($(quot(o)),l,r)
            end
end

for oc = U_OP_START:U_OP_END
    o = OP[oc]
    S_TO_OC[o] = oc
    @eval   ($o){I}(l::AD{I}) =
            begin
                return AD_O($(quot(o)),l)
            end
end

(+){I}(l::AD{I}, args...) = AD_O(:+,tuple(l,args...))
(*){I}(l::AD{I}, args...) = AD_O(:*,tuple(l,args...))

##########################################################################################

function replace_sym(s::Symbol,ex::AbstractString,r::AbstractString)
    return replace_sym(s,parse(ex),parse(r))
end
function replace_sym(s::Symbol,ex::Expr,r::AbstractString)
    return replace_sym(s,ex,parse(r))
end

function replace_sym(s::Symbol,ex::Number,r)
    return ex
end
function replace_sym(s::Symbol,ex::Symbol,r)
    if(ex==s) return r
    else return ex
    end
end
function replace_sym(s::Symbol,ex::Expr, r)
    for i = 1:length(ex.args)
        ex.args[i] = replace_sym(s,ex.args[i],r)
    end
    return ex
end

##########################################################################################
#
# 0ord eval
#
##########################################################################################
## one argument functions
switchblock = Expr(:block)
for i = U_OP_START:U_OP_END
    o = OP[i]
    ex = :(return $(o)(v))
    push!(switchblock.args,quot(o),ex)
end
switchexpr = Expr(:macrocall, Expr(:.,:Lazy,quot(Symbol("@switch"))), :s,switchblock)
@eval @inline function eval_0ord{V}(s::Symbol, v::V)
    if s==:-
        return -v
    else
        $switchexpr
    end
end

# 2 argument functions
switchblock = Expr(:block)

for i = B_OP_START:B_OP_END
    o = OP[i]
    ex = Expr(:block)
    push!(ex.args,parse("return $(o)(l,r)"))
    # ex = :(@inbounds return $(o)(l,r))
    push!(switchblock.args,quot(o),ex)
end

switchexpr = Expr(:macrocall, Expr(:.,:Lazy,quot(Symbol("@switch"))), :s,switchblock)
@eval @inline function eval_0ord{V}(s::Symbol,l::V,r::V)
    $switchexpr
end

# 2+ argument functions
@eval @inline function eval_0ord{I,V}(s::Symbol, v::Vector{V}, i::I, e::I)
    # @show s,v
    ret = v[i]
    if s==:+
        @simd for j = i+1:e
            @inbounds ret +=  v[j]
        end
    else
        @simd for j = i+1:e
            @inbounds ret *=  v[j]
        end
    end
    return ret
end

##########################################################################################
#
# 1ord eval
#
##########################################################################################
# one argument functions
switchblock = Expr(:block)
for i = U_OP_START:U_OP_END
    o = OP[i]
    dx = S_TO_DIFF[o][1]
    dx = replace_sym(:x,dx,"v")
    # @show dx
    ex = Expr(:block)
    push!(ex.args,parse("@inbounds imm[1]=$(dx)"))
    push!(switchblock.args,quot(o),ex)
end
switchexpr = Expr(:macrocall,Expr(:.,:Lazy,quot(Symbol("@switch"))),:s,switchblock)
@eval @inline function eval_1ord{V}(s::Symbol,v::V,imm::Vector{V})
    if s==:-
        @inbounds imm[1] = -one(V)
    else
        $switchexpr
    end
    nothing
end

# 2 argument functions
switchblock = Expr(:block)

for i = B_OP_START:B_OP_END
    o = OP[i]
    ex = Expr(:block)

    (dx,dy) = S_TO_DIFF[o]
    dx = replace_sym(:y,replace_sym(:x,dx,"l"),"r")
    dy = replace_sym(:y,replace_sym(:x,dy,"l"),"r")
    # @show dx, dy
    ex = Expr(:block)
    push!(ex.args,parse("@inbounds imm[1]=$(dx)"))
    push!(ex.args,parse("@inbounds imm[2]=$(dy)"))
    push!(switchblock.args,quot(o),ex)
end

switchexpr = Expr(:macrocall, Expr(:.,:Lazy,quot(Symbol("@switch"))), :s,switchblock)
# @show switchexpr
@eval @inline function eval_1ord{V}(s::Symbol,l::V, r::V, imm::Vector{V})
    if s==:^
        exponent::V = r
        base::V = l
        t0 = l<0.0?0.0:log(l)
        t = l^r
        if exponent == 2.0
            @inbounds imm[1] = 2.0*l
            @inbounds imm[2] = t*t0
        else
            @inbounds imm[1] = r*l^(r-1.0)
            @inbounds imm[2] = t*t0
        end
    else
        $switchexpr
    end
    nothing
end

@eval @inline function eval_1ord{I,V}(s::Symbol, v::Vector{V},i::I,e::I,imm::Vector{V})
    if s==:+
        # @inbounds imm[1:(e-i)+1] = one(V)
    else
        imm_i = zero(I)
        @simd for j0=i:e
            c=one(V)
            @simd for j1=i:e
                if j0!=j1
                    @inbounds c *= v[j1]
                end
            end
            imm_i += 1
            @inbounds imm[imm_i] = c
        end
    end
    nothing
end

##########################################################################################
#
# 2ord eval
#
##########################################################################################
#one argument funtion
switchblock = Expr(:block)
for i = U_OP_START:U_OP_END
    o = OP[i]
    (dx,dxx) = S_TO_DIFF[o]
    dx = replace_sym(:x,dx,"v")
    dxx = replace_sym(:x,dxx,"v")
    # @show dx, dxx
    ex = Expr(:block)
    push!(ex.args,parse("@inbounds imm[1] = $(dx)"))
    push!(ex.args,parse("@inbounds imm[2] = $(dxx)"))
    push!(switchblock.args,quot(o),ex)
end
switchexpr = Expr(:macrocall,Expr(:.,:Lazy,quot(Symbol("@switch"))),:s,switchblock)
@eval @inline function eval_2ord{V}(s::Symbol,v::V,imm::Vector{V})
    if s==:-
        @inbounds imm[1] = -one(V)
        @inbounds imm[2] = zero(V)
    else
        $switchexpr
    end
    nothing
end

# 2 argument function
switchblock = Expr(:block)
for i = B_OP_START:B_OP_END
    o = OP[i]
    ex = Expr(:block)
    (dx,dy,dxx,dxy,dyy) = S_TO_DIFF[o]
    dx = replace_sym(:y,replace_sym(:x,dx,"l"),"r")
    dy = replace_sym(:y,replace_sym(:x,dy,"l"),"r")
    dxx = replace_sym(:y,replace_sym(:x,dxx,"l"),"r")
    dxy = replace_sym(:y,replace_sym(:x,dxy,"l"),"r")
    dyy = replace_sym(:y,replace_sym(:x,dyy,"l"),"r")
    # @show dx
    # @show dy
    # @show dxx
    # @show dyy
    # @show dxy
    ex = Expr(:block)
    push!(ex.args,parse("@inbounds imm[1] = $(dx)"))
    push!(ex.args,parse("@inbounds imm[2] = $(dy)"))
    push!(ex.args,parse("@inbounds imm[3] = $(dxx)"))
    push!(ex.args,parse("@inbounds imm[4] = $(dxy)"))
    push!(ex.args,parse("@inbounds imm[5] = $(dyy)"))
    push!(switchblock.args,quot(o),ex)
end
switchexpr = Expr(:macrocall, Expr(:.,:Lazy,quot(Symbol("@switch"))), :s,switchblock)
@eval @inline function eval_2ord{V}(s::Symbol,l::V, r::V, imm::Vector{V})
    if s == :^
        exponent = r
        base = l
        if exponent == 2.0
            t = base*base
            t2 = base<0.0?0.0:log(base)
            @inbounds imm[1] = 2.0*base
            @inbounds imm[2] = t*t2
            @inbounds imm[3] = 2.0
            @inbounds imm[4] = base + 2.0*base*t2
            @inbounds imm[5] = t*t2*t2
        else
            t = base^exponent
            t2 = base<0.0?0.0:log(base)
            t3 = base^(exponent-1.0)
            @inbounds imm[1] = exponent*base^t3
            @inbounds imm[2] = t*t2
            @inbounds imm[3] = exponent*(exponent-1.0)*base^(exponent-2.0)
            @inbounds imm[4] = t3 + exponent*t3*t2
            @inbounds imm[5] = t*t2*t2
        end
    else
        $switchexpr
    end
    nothing
end

@eval @inline function eval_2ord{I,V}(s::Symbol,v::Vector{V},i::I,e::I,imm::Vector{V})
    if s==:+
        # @inbounds imm[1:(e-i)+1] = one(V)
    else
        #first order
        imm_i = zero(I)
        for j0=i:e
            f_ord = one(V)
            for j1=i:j0-1
                @inbounds f_ord *= v[j1]
            end
            for j1=j0+1:e
                @inbounds f_ord *= v[j1]
            end
            imm_i += 1
            @inbounds imm[imm_i] = f_ord
        end
        #second order
        for j0=i:e
            for j1=j0+1:e
                s_ord = one(V)
                for j3=i:j0-1
                    @inbounds s_ord *= v[j3]
                end
                for j3=j0+1:j1-1
                    @inbounds s_ord *= v[j3]
                end
                for j3=j1+1:e
                    @inbounds s_ord *= v[j3]
                end
                imm_i+=1
                @inbounds imm[imm_i] = s_ord
            end
        end
    end
    nothing
end
##########################################################################################


text:
using PMDPs
using PMDPs.LP
using BSON
using DataFrames
using Random, RandomNumbers.Xorshifts

@testset "hindsight_policy.jl" begin
    pp_params = Dict{Symbol,Any}(:T=>1000,:c=>40,:res_budget_μ=>5.0,:nᵣ=>10,:expected_res=>1200.0)
    pp = PMDPs.linear_pp(;pp_params...)

    mg = PMDPs.PMDPg(pp)
    rnd = Xorshift128Plus(1)
    traces = [PMDPs.simulate_trace(mg, rnd) for i in 1:3]

    # get results by applying MILP based hindsight policy
    results = PMDPs.hindsight(pp, traces, MersenneTwister(1))

    # compare the result of applying the MILP policy with value of the optimization criterion from MILP
    for (i,tr) in enumerate(traces)
        (r, u, alloc, action_seq, requests) = PMDPs.LP.MILP_hindsight_pricing(mg, tr)
        @test r==results.r[i]
    end

end


text:
using SplineFEEC
using Test

@testset "SplineFEEC.jl" begin
    # Write your tests here.
end


text:
using MathProgBase.SolverInterface
#export FOSAlgorithm

abstract type AbstractSolution end

mutable struct Solution <: AbstractSolution
    x::Array{Float64, 1}
    y::Array{Float64, 1}
    s::Array{Float64, 1}
    status::Symbol
end

abstract type FOSAlgorithm <: AbstractMathProgSolver end
abstract type FOSSolverData end
struct FOSSolverDataPlaceholder <: FOSSolverData end

abstract type FOSModel end

abstract type AbstractStatus end
mutable struct NoStatus <: AbstractStatus
    status::Symbol
end
# Define Solver for interface
# struct FOSSolver <: AbstractMathProgSolver
#     options
# end
# FOSSolver(;kwargs...) = FOSSolver(kwargs)
#

mutable struct FOSMathProgModel <: AbstractConicModel
    input_numconstr::Int64            # Only needed for interface?
    input_numvar::Int64               # Only needed for interface?
    K1::ConeProduct
    K2::ConeProduct
    A::SparseMatrixCSC{Float64,Int}   # The A matrix (equalities)
    b::Vector{Float64}                # RHS
    c::Vector{Float64}                # The objective coeffs (always min)
    alg::FOSAlgorithm                 #This descides which algorithm to call
    data::FOSSolverData               # Solver specific data is stored here
#    orig_sense::Symbol               # Original objective sense
    # Post-solve
    solve_stat::Symbol
    obj_val::Float64
    primal_sol::Vector{Float64}
    dual_sol::Vector{Float64}
    slack::Vector{Float64}
    options::Dict{Symbol, Any}
    enditr::Int64
    status_generator::Function
    init_duration::UInt64           # In nano-seconds
    history::ValueHistories.MVHistory
end

function FOSMathProgModel(s::FOSAlgorithm; kwargs...)
    FOSMathProgModel(0, 0, ConeProduct(), ConeProduct(), spzeros(0, 0),
                     Float64[], Float64[], s, FOSSolverDataPlaceholder(), :NotSolved,
                     0.0, Float64[], Float64[], Float64[], Dict{Symbol,Any}(kwargs), -1,
                     x -> (@error "No status generator defined"),
                     UInt64(1), ValueHistories.MVHistory())
end

# To handle both custom types and
const AbstractFOSModel = Union{FOSMathProgModel, FOSModel}


text:
module OpenTelemetrySDK

using OpenTelemetryAPI

const API = OpenTelemetryAPI

include("instrumentation.jl")
include("resource.jl")
include("trace/sampling.jl")
include("trace/id_generator.jl")
include("trace/span_processor.jl")
include("trace/trace.jl")

end # module


text:
abstract type AbstractRefinery end

needs_refinement(refinery::AbstractRefinery, cell) = error("Your refinery should implement this function")
refine_data(refinery::AbstractRefinery, cell, indices) = error("Your refinery should implement this function")

function adaptivesampling!(root::Cell, refinery::AbstractRefinery)
    refinement_queue = [root]
    refine_function =
        (cell, indices) -> refine_data(refinery, cell, indices)
    while !isempty(refinement_queue)
        cell = pop!(refinement_queue)
        if needs_refinement(refinery, cell)
            split!(cell, refine_function)
            append!(refinement_queue, children(cell))
        end
    end
    root
end


text:
abstract type AbstractScenarioProblems{S <: AbstractScenario} end
struct ScenarioProblems{S <: AbstractScenario} <: AbstractScenarioProblems{S}
    scenarios::Vector{S}
    problems::Vector{JuMP.Model}
    parent::JuMP.Model

    function ScenarioProblems(::Type{S}) where S <: AbstractScenario
        return new{S}(Vector{S}(), Vector{JuMP.Model}(), JuMP.Model())
    end

    function ScenarioProblems(scenarios::Vector{S}) where S <: AbstractScenario
        return new{S}(scenarios, Vector{JuMP.Model}(), JuMP.Model())
    end
end
ScenarioProblemChannel{S} = RemoteChannel{Channel{ScenarioProblems{S}}}
struct DScenarioProblems{S <: AbstractScenario} <: AbstractScenarioProblems{S}
    parent::JuMP.Model
    scenario_distribution::Vector{Int}
    scenarioproblems::Vector{ScenarioProblemChannel{S}}

    function DScenarioProblems(scenario_distribution::Vector{Int}, scenarioproblems::Vector{ScenarioProblemChannel{S}}) where S <: AbstractScenario
        return new{S}(JuMP.Model(), scenario_distribution, scenarioproblems)
    end
end
function ScenarioProblems(::Type{S}, procs::Vector{Int}) where S <: AbstractScenario
    if (length(procs) == 1 || nworkers() == 1) && procs[1] == 1
        return ScenarioProblems(S)
    else
        isempty(procs) && error("No requested procs.")
        length(procs) <= nworkers() || error("Not enough workers to satisfy requested number of procs. There are ", nworkers(), " workers, but ", length(procs), " were requested.")
        scenarioproblems = Vector{ScenarioProblemChannel{S}}(undef, length(procs))
        @sync begin
            for (i,p) in enumerate(procs)
                scenarioproblems[i] = RemoteChannel(() -> Channel{ScenarioProblems{S}}(1), p)
                @async remotecall_fetch((sp,S)->put!(sp, ScenarioProblems(S)),
                                        p,
                                        scenarioproblems[i],
                                        S)
            end
        end
        return DScenarioProblems(zeros(Int, length(procs)), scenarioproblems)
    end
end

function ScenarioProblems(scenarios::Vector{S}, procs::Vector{Int}) where S <: AbstractScenario
    if (length(procs) == 1 || nworkers() == 1) && procs[1] == 1
        return ScenarioProblems(scenarios)
    else
        isempty(procs) && error("No requested procs.")
        length(procs) <= nworkers() || error("Not enough workers to satisfy requested number of procs. There are ", nworkers(), " workers, but ", length(procs), " were requested.")
        scenarioproblems = Vector{ScenarioProblemChannel{S}}(undef, length(procs))
        (nscen, extra) = divrem(length(scenarios), length(procs))
        start = 1
        stop = nscen + (extra > 0)
        scenario_distribution = zeros(Int, length(procs))
        @sync begin
            for (i,p) in enumerate(procs)
                n = nscen + (extra > 0)
                scenarioproblems[i] = RemoteChannel(() -> Channel{ScenarioProblems{S}}(1), p)
                scenario_range = start:stop
                @async remotecall_fetch((sp,scenarios)->put!(sp, ScenarioProblems(scenarios)),
                                        p,
                                        scenarioproblems[i],
                                        scenarios[scenario_range])
                scenario_distribution[i] = n
                start = stop + 1
                stop += n
                stop = min(stop, length(scenarios))
                extra -= 1
            end
        end
        return DScenarioProblems(scenario_distribution, scenarioproblems)
    end
end

# Base overloads #
# ========================== #
Base.getindex(sp::DScenarioProblems, i::Integer) = sp.scenarioproblems[i]
# ========================== #

# Getters #
# ========================== #
function scenario(scenarioproblems::ScenarioProblems{S}, i::Integer) where S <: AbstractScenario
    return scenarioproblems.scenarios[i]
end
function scenario(scenarioproblems::DScenarioProblems{S}, i::Integer) where S <: AbstractScenario
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if i <= n+j
            return remotecall_fetch((sp,i)->fetch(sp).scenarios[i], w, scenarioproblems[w-1], i-j)
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, i))
end
function scenarios(scenarioproblems::ScenarioProblems{S}) where S <: AbstractScenario
    return scenarioproblems.scenarios
end
function scenarios(scenarioproblems::DScenarioProblems{S}) where S <: AbstractScenario
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    partial_scenarios = Vector{Vector{S}}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_scenarios[i] = remotecall_fetch((sp)->fetch(sp).scenarios, w, scenarioproblems[w-1])
        end
    end
    return reduce(vcat, partial_scenarios)
end
function expected(scenarioproblems::ScenarioProblems{S}) where S <: AbstractScenario
    return expected(scenarioproblems.scenarios)
end
function expected(scenarioproblems::DScenarioProblems{S}) where S <: AbstractScenario
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    partial_expecations = Vector{S}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_expecations[i] = remotecall_fetch((sp) -> expected(fetch(sp)).scenario, w, scenarioproblems[w-1])
        end
    end
    return expected(partial_expecations)
end
function scenariotype(scenarioproblems::ScenarioProblems{S}) where S <: AbstractScenario
    return S
end
function scenariotype(scenarioproblems::DScenarioProblems{S}) where S <: AbstractScenario
    return S
end
function subproblem(scenarioproblems::ScenarioProblems{S}, i::Integer) where S <: AbstractScenario
    return scenarioproblems.problems[i]
end
function subproblem(scenarioproblems::DScenarioProblems{S}, i::Integer) where S <: AbstractScenario
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if i <= n+j
            return remotecall_fetch((sp,i)->fetch(sp).problems[i],w,scenarioproblems[w-1],i-j)
        end
        j += n
    end
    throw(BoundsError(scenarioproblems,i))
end
function subproblems(scenarioproblems::ScenarioProblems{S}) where S <: AbstractScenario
    return scenarioproblems.problems
end
function subproblems(scenarioproblems::DScenarioProblems{S}) where S <: AbstractScenario
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    partial_subproblems = Vector{JuMP.Model}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_subproblems[i] = remotecall_fetch((sp)->fetch(sp).problems,w,scenarioproblems[w-1])
        end
    end
    return reduce(vcat, partial_subproblems)
end
function nsubproblems(scenarioproblems::ScenarioProblems{S}) where S <: AbstractScenario
    return length(scenarioproblems.problems)
end
function nsubproblems(scenarioproblems::DScenarioProblems{S}) where S <: AbstractScenario
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    partial_lengths = Vector{Int}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_lengths[i] = remotecall_fetch((sp) -> nsubproblems(fetch(sp)),w,scenarioproblems[w-1])
        end
    end
    return sum(partial_lengths)
end
function parentmodel(scenarioproblems::AbstractScenarioProblems)
    return scenarioproblems.parent
end
function recourse_length(scenarioproblems::ScenarioProblems)
    return scenarioproblems.problems[1].numCols
end
function recourse_length(scenarioproblems::DScenarioProblems)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    return remotecall_fetch((sp)->recourse_length(fetch(sp)), 2, scenarioproblems[1])
end
function probability(scenarioproblems::ScenarioProblems{S}) where S <: AbstractScenario
    return probability(scenarioproblems.scenarios)
end
function probability(scenarioproblems::DScenarioProblems{S}) where S <: AbstractScenario
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    partial_probabilities = Vector{Float64}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_probabilities[i] = remotecall_fetch((sp) -> probability(fetch(sp)), w, scenarioproblems[w-1])
        end
    end
    return sum(partial_probabilities)
end
function nscenarios(scenarioproblems::ScenarioProblems{S}) where S <: AbstractScenario
    return length(scenarioproblems.scenarios)
end
function nscenarios(scenarioproblems::DScenarioProblems{S}) where S <: AbstractScenario
    return sum(scenarioproblems.scenario_distribution)
end
distributed(scenarioproblems::ScenarioProblems) = false
distributed(scenarioproblems::DScenarioProblems) = true
# ========================== #

# Setters
# ========================== #
function set_stage_data!(scenarioproblems::ScenarioProblems{P}, data::P) where P
    scenarioproblems.stage.data = data
    return nothing
end
function set_stage_data!(scenarioproblems::DScenarioProblems{P}, data::P) where P
    scenarioproblems.stage.data = data
    @sync begin
        for w in workers()
            @async remotecall_fetch((sp, data)->set_stage_data(fetch(sp), data), w, scenarioproblems[w-1], data)
        end
    end
    return nothing
end
function add_scenario!(scenarioproblems::ScenarioProblems{S}, scenario::S) where S <: AbstractScenario
    push!(scenarioproblems.scenarios, scenario)
    return nothing
end
function add_scenario!(scenarioproblems::DScenarioProblems{S}, scenario::S) where S <: AbstractScenario
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    _, w = findmin(scenarioproblems.scenario_distribution)
    add_scenario!(scenarioproblems, scenario, w+1)
    return nothing
end
function add_scenario!(scenarioproblems::DScenarioProblems{S}, scenario::S, w::Integer) where S <: AbstractScenario
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    remotecall_fetch((sp, scenario) -> add_scenario!(fetch(sp), scenario),
                     w,
                     scenarioproblems[w-1],
                     scenario)
    scenarioproblems.scenario_distribution[w-1] += 1
    return nothing
end
function add_scenario!(scenariogenerator::Function, scenarioproblems::ScenarioProblems)
    add_scenario!(scenarioproblems, scenariogenerator())
    return nothing
end
function add_scenario!(scenariogenerator::Function, scenarioproblems::DScenarioProblems)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    _, w = findmin(scenarioproblems.scenario_distribution)
    add_scenario!(scenariogenerator, scenarioproblems, w+1)
    return nothing
end
function add_scenario!(scenariogenerator::Function, scenarioproblems::DScenarioProblems, w::Integer)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    remotecall_fetch((sp, generator) -> add_scenario!(fetch(sp), generator()),
                     w,
                     scenarioproblems[w-1],
                     scenariogenerator)
    scenarioproblems.scenario_distribution[w] += 1
    return nothing
end
function add_scenarios!(scenarioproblems::ScenarioProblems{S}, scenarios::Vector{S}) where S <: AbstractScenario
    append!(scenarioproblems.scenarios, scenarios)
    return nothing
end
function add_scenarios!(scenariogenerator::Function, scenarioproblems::ScenarioProblems{S}, n::Integer) where S <: AbstractScenario
    for i = 1:n
        add_scenario!(scenarioproblems) do
            return scenariogenerator()
        end
    end
    return nothing
end
function add_scenarios!(scenarioproblems::DScenarioProblems{S}, scenarios::Vector{S}) where S <: AbstractScenario
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    (nscen, extra) = divrem(length(scenarios), nworkers())
    start = 1
    stop = nscen + (extra > 0)
    @sync begin
        for w in workers()
            n = nscen + (extra > 0)
            scenario_range = start:stop
            @async remotecall_fetch((sp, scenarios) -> add_scenarios!(fetch(sp), scenarios),
                                    w,
                                    scenarioproblems[w-1],
                                    scenarios[scenario_range])
            scenarioproblems.scenario_distribution[w-1] += n
            start = stop + 1
            stop += n
            stop = min(stop, length(scenarios))
            extra -= 1
        end
    end
    return nothing
end
function add_scenarios!(scenarioproblems::DScenarioProblems{S}, scenarios::Vector{S}, w::Integer) where S <: AbstractScenario
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    remotecall_fetch((sp, scenarios) -> add_scenarios!(fetch(sp), scenarios),
                     w,
                     scenarioproblems[w-1],
                     scenarios)
    scenarioproblems.scenario_distribution[w-1] += length(scenarios)
    return nothing
end
function add_scenarios!(scenariogenerator::Function, scenarioproblems::DScenarioProblems, n::Integer)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    (nscen, extra) = divrem(n, nworkers())
    @sync begin
        for w in workers()
            m = nscen + (extra > 0)
            @async remotecall_fetch((sp, gen, n) -> add_scenarios!(gen, fetch(sp), n),
                                    w,
                                    scenarioproblems[w-1],
                                    scenariogenerator,
                                    m)
            scenarioproblems.scenario_distribution[w-1] += m
            extra -= 1
        end
    end
    return nothing
end
function add_scenarios!(scenariogenerator::Function, scenarioproblems::DScenarioProblems{S}, n::Integer, w::Integer) where S <: AbstractScenario
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    remotecall_fetch((sp, gen) -> add_scenarios!(gen, fetch(sp), n),
                     w,
                     scenarioproblems[w-1],
                     scenariogenerator,
                     n)
    scenarioproblems.scenario_distribution[w-1] += n
    return nothing
end
function remove_scenarios!(scenarioproblems::ScenarioProblems)
    empty!(scenarioproblems.scenarios)
    return nothing
end
function remove_scenarios!(scenarioproblems::DScenarioProblems)
    @sync begin
        for w in workers()
            @async remotecall_fetch((sp)->remove_scenarios!(fetch(sp)), w, scenarioproblems[w-1])
            scenarioproblems.scenario_distribution[w-1] = 0
        end
    end
    return nothing
end
function remove_subproblems!(scenarioproblems::ScenarioProblems)
    empty!(scenarioproblems.problems)
    return nothing
end
function remove_subproblems!(scenarioproblems::DScenarioProblems)
    @sync begin
        for w in workers()
            @async remotecall_fetch((sp)->remove_subproblems!(fetch(sp)), w, scenarioproblems[w-1])
        end
    end
    return nothing
end
# ========================== #

# Sampling #
# ========================== #
function sample!(scenarioproblems::ScenarioProblems{S}, sampler::AbstractSampler{S}, n::Integer) where S <: AbstractScenario
    _sample!(scenarioproblems, sampler, n, nscenarios(scenarioproblems), 1/n)
end
function sample!(scenarioproblems::ScenarioProblems{S}, sampler::AbstractSampler{Scenario}, n::Integer) where S <: Scenario
    _sample!(scenarioproblems, sampler, n, nscenarios(scenarioproblems), 1/n)
end
function sample!(scenarioproblems::DScenarioProblems{S}, sampler::AbstractSampler{S}, n::Integer) where S <: AbstractScenario
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    m = nscenarios(scenarioproblems)
    (nscen, extra) = divrem(n, nworkers())
    @sync begin
        for w in workers()
            d = nscen + (extra > 0)
            @async remotecall_fetch((sp,sampler,n,m,π)->_sample!(fetch(sp),sampler,n,m,π), w, scenarioproblems[w-1], sampler, d, m, 1/n)
            scenarioproblems.scenario_distribution[w-1] += d
            extra -= 1
        end
    end
    return nothing
end
function sample!(scenarioproblems::DScenarioProblems{S}, sampler::AbstractSampler{Scenario}, n::Integer) where S <: Scenario
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    m = nscenarios(scenarioproblems)
    (nscen, extra) = divrem(n, nworkers())
    @sync begin
        for w in workers()
            d = nscen + (extra > 0)
            @async remotecall_fetch((sp,sampler,n,m,π)->_sample!(fetch(sp),sampler,n,m,π), w, scenarioproblems[w-1], sampler, d, m, 1/n)
            scenarioproblems.scenario_distribution[w-1] += d
            extra -= 1
        end
    end
    return nothing
end
function _sample!(scenarioproblems::ScenarioProblems, sampler::AbstractSampler, n::Integer, m::Integer, π::AbstractFloat)
    if m > 0
        # Rescale probabilities of existing scenarios
        for scenario in scenarioproblems.scenarios
            p = probability(scenario) * m / (m+n)
            set_probability!(scenario, p)
        end
        π *= n/(m+n)
    end
    for i = 1:n
        add_scenario!(scenarioproblems) do
            return sample(sampler, π)
        end
    end
    return nothing
end
# ========================== #


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

relaxed_tols = (default_tol_relax = 1000,)
insts = OrderedDict()
insts["minimal"] = [
    ((2, 4, true, true),),
    ((2, 4, true, false),),
    ((2, 2, false, true),),
    ((2, 2, false, false),),
    ]
insts["fast"] = [
    ((3, 6, true, true),),
    ((3, 6, true, false),),
    ((5, 5, false, true),),
    ((5, 5, false, false),),
    ((10, 20, true, true),),
    ((10, 20, true, false),),
    ((15, 15, false, true),),
    ((15, 15, false, false),),
    ((30, 30, false, false),),
    ((30, 30, false, true),),
    ]
insts["various"] = [
    ((6, 6, true, true), nothing, relaxed_tols),
    ((6, 6, true, false), nothing, relaxed_tols),
    ((6, 6, false, true),),
    ((6, 6, false, false),),
    ((12, 12, true, true), nothing, relaxed_tols),
    ((12, 12, false, true),),
    ((12, 12, false, false),),
    ((24, 24, true, true), nothing, relaxed_tols),
    ((24, 24, false, true), nothing, relaxed_tols),
    ((24, 24, false, false), nothing, relaxed_tols),
    ]
return (LyapunovStabilityJuMP, insts)