text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using JuMP, MathOptInterface, DataStructures
import JuMP: value, fix, set_name

# Extend some JuMP functions so that decision variables can be safely replaced by
# (constant) floating point numbers.
function value(x::Float64)
    x
end

function fix(x::Float64, v::Float64; force)
    abs(x - v) < 1e-6 || error("Value mismatch: $x != $v")
end

function set_name(x::Float64, n::String)
    # nop
end

mutable struct UnitCommitmentModel
    mip::JuMP.Model
    vars::DotDict
    eqs::DotDict
    exprs::DotDict
    instance::UnitCommitmentInstance
    isf::Array{Float64, 2}
    lodf::Array{Float64, 2}
    obj::AffExpr
end

function build_model(;
                     filename::Union{String, Nothing}=nothing,
                     instance::Union{UnitCommitmentInstance, Nothing}=nothing,
                     isf::Union{Array{Float64,2}, Nothing}=nothing,
                     lodf::Union{Array{Float64,2}, Nothing}=nothing,
                     isf_cutoff::Float64=0.005,
                     lodf_cutoff::Float64=0.001,
                     optimizer=nothing,
                     model=nothing,
                     variable_names::Bool=false,
                    ) :: UnitCommitmentModel

    if (filename === nothing) && (instance === nothing)
        error("Either filename or instance must be specified")
    end

    if filename !== nothing
        @info "Reading: $(filename)"
        time_read = @elapsed begin
            instance = UnitCommitment.read(filename)
        end
        @info @sprintf("Read problem in %.2f seconds", time_read)
    end

    if length(instance.buses) == 1
        isf = zeros(0, 0)
        lodf = zeros(0, 0)
    else
        if isf === nothing
            @info "Computing injection shift factors..."
            time_isf = @elapsed begin
                isf = UnitCommitment.injection_shift_factors(lines=instance.lines,
                                                             buses=instance.buses)
            end
            @info @sprintf("Computed ISF in %.2f seconds", time_isf)

            @info "Computing line outage factors..."
            time_lodf = @elapsed begin
                lodf = UnitCommitment.line_outage_factors(lines=instance.lines,
                                                          buses=instance.buses,
                                                          isf=isf)
            end
            @info @sprintf("Computed LODF in %.2f seconds", time_lodf)

            @info @sprintf("Applying PTDF and LODF cutoffs (%.5f, %.5f)", isf_cutoff, lodf_cutoff)
            isf[abs.(isf) .< isf_cutoff] .= 0
            lodf[abs.(lodf) .< lodf_cutoff] .= 0
        end
    end

    @info "Building model..."
    time_model = @elapsed begin
        if model === nothing
            if optimizer === nothing
                mip = Model()
            else
                mip = Model(optimizer)
            end
        else
            mip = model
        end
        model = UnitCommitmentModel(mip,
                                    DotDict(),  # vars
                                    DotDict(),  # eqs
                                    DotDict(),  # exprs
                                    instance,
                                    isf,
                                    lodf,
                                    AffExpr(),  # obj
                                   )
        for field in [:prod_above, :segprod, :reserve, :is_on, :switch_on, :switch_off,
                      :net_injection, :curtail, :overflow, :loads, :startup]
            setproperty!(model.vars, field, OrderedDict())
        end
        for field in [:startup_choose, :startup_restrict, :segprod_limit, :prod_above_def,
                      :prod_limit, :binary_link, :switch_on_off, :ramp_up, :ramp_down,
                      :startup_limit, :shutdown_limit, :min_uptime, :min_downtime, :power_balance,
                      :net_injection_def, :min_reserve]
            setproperty!(model.eqs, field, OrderedDict())
        end
        for field in [:inj, :reserve, :net_injection]
            setproperty!(model.exprs, field, OrderedDict())
        end
        for lm in instance.lines
            add_transmission_line!(model, lm)
        end
        for b in instance.buses
            add_bus!(model, b)
        end
        for g in instance.units
            add_unit!(model, g)
        end
        for ps in instance.price_sensitive_loads
            add_price_sensitive_load!(model, ps)
        end
        build_net_injection_eqs!(model)
        build_reserve_eqs!(model)
        build_obj_function!(model)
    end
    @info @sprintf("Built model in %.2f seconds", time_model)

    if variable_names
        set_variable_names!(model)
    end

    return model
end

function add_transmission_line!(model, lm)
    vars, obj, T = model.vars, model.obj, model.instance.time
    for t in 1:T
        overflow = vars.overflow[lm.name, t] = @variable(model.mip, lower_bound=0)
        add_to_expression!(obj, overflow, lm.flow_limit_penalty[t])
    end
end

function add_bus!(model::UnitCommitmentModel, b::Bus)
    mip, vars, exprs = model.mip, model.vars, model.exprs
    for t in 1:model.instance.time
        # Fixed load
        exprs.net_injection[b.name, t] = AffExpr(-b.load[t])

        # Reserves
        exprs.reserve[b.name, t] = AffExpr()

         # Load curtailment
        vars.curtail[b.name, t] = @variable(mip, lower_bound=0, upper_bound=b.load[t])
        add_to_expression!(exprs.net_injection[b.name, t], vars.curtail[b.name, t], 1.0)
        add_to_expression!(model.obj,
                           vars.curtail[b.name, t],
                           model.instance.power_balance_penalty[t])
    end
end

function add_price_sensitive_load!(model::UnitCommitmentModel, ps::PriceSensitiveLoad)
    mip, vars = model.mip, model.vars
    for t in 1:model.instance.time
        # Decision variable
        vars.loads[ps.name, t] = @variable(mip, lower_bound=0, upper_bound=ps.demand[t])

        # Objective function terms
        add_to_expression!(model.obj, vars.loads[ps.name, t], -ps.revenue[t])

        # Net injection
        add_to_expression!(model.exprs.net_injection[ps.bus.name, t], vars.loads[ps.name, t], -1.0)
    end
end

function add_unit!(model::UnitCommitmentModel, g::Unit)
    mip, vars, eqs, exprs, T = model.mip, model.vars, model.eqs, model.exprs, model.instance.time
    gi, K, S = g.name, length(g.cost_segments), length(g.startup_categories)

    if !all(g.must_run) && any(g.must_run)
        error("Partially must-run units are not currently supported")
    end

    if g.initial_power === nothing || g.initial_status === nothing
        error("Initial conditions for $(g.name) must be provided")
    end

    is_initially_on = (g.initial_status > 0 ? 1.0 : 0.0)

    # Decision variables
    for t in 1:T
        for k in 1:K
            model.vars.segprod[gi, t, k] = @variable(model.mip, lower_bound=0)
        end
        model.vars.prod_above[gi, t] = @variable(model.mip, lower_bound=0)
        if g.provides_spinning_reserves[t]
            model.vars.reserve[gi, t] = @variable(model.mip, lower_bound=0)
        else
            model.vars.reserve[gi, t] = 0.0
        end
        for s in 1:S
            model.vars.startup[gi, t, s] = @variable(model.mip, binary=true)
        end
        if g.must_run[t]
            model.vars.is_on[gi, t] = 1.0
            model.vars.switch_on[gi, t] = (t == 1 ? 1.0 - is_initially_on : 0.0)
            model.vars.switch_off[gi, t] = 0.0
        else
            model.vars.is_on[gi, t] = @variable(model.mip, binary=true)
            model.vars.switch_on[gi, t] = @variable(model.mip, binary=true)
            model.vars.switch_off[gi, t] = @variable(model.mip, binary=true)
        end
    end

    for t in 1:T
        # Time-dependent start-up costs
        for s in 1:S
            # If unit is switching on, we must choose a startup category
            eqs.startup_choose[gi, t, s] =
                @constraint(mip, vars.switch_on[gi, t] == sum(vars.startup[gi, t, s] for s in 1:S))

            # If unit has not switched off in the last `delay` time periods, startup category is forbidden.
            # The last startup category is always allowed.
            if s < S
                range = (t - g.startup_categories[s + 1].delay + 1):(t - g.startup_categories[s].delay)
                initial_sum = (g.initial_status < 0 && (g.initial_status + 1 in range) ? 1.0 : 0.0)
                eqs.startup_restrict[gi, t, s] =
                    @constraint(mip, vars.startup[gi, t, s]
                                        <= initial_sum + sum(vars.switch_off[gi, i] for i in range if i >= 1))
            end

            # Objective function terms for start-up costs
            add_to_expression!(model.obj,
                               vars.startup[gi, t, s],
                               g.startup_categories[s].cost)
        end

        # Objective function terms for production costs
        add_to_expression!(model.obj, vars.is_on[gi, t], g.min_power_cost[t])
        for k in 1:K
            add_to_expression!(model.obj, vars.segprod[gi, t, k], g.cost_segments[k].cost[t])
        end

        # Production limits (piecewise-linear segments)
        for k in 1:K
            eqs.segprod_limit[gi, t, k] =
                @constraint(mip, vars.segprod[gi, t, k] <= g.cost_segments[k].mw[t] * vars.is_on[gi, t])
        end

        # Definition of production
        eqs.prod_above_def[gi, t] =
            @constraint(mip, vars.prod_above[gi, t] == sum(vars.segprod[gi, t, k] for k in 1:K))

        # Production limit
        eqs.prod_limit[gi, t] =
            @constraint(mip,
                        vars.prod_above[gi, t] + vars.reserve[gi, t]
                                <= (g.max_power[t] - g.min_power[t]) * vars.is_on[gi, t])

        # Binary variable equations for economic units
        if !g.must_run[t]

            # Link binary variables
            if t == 1
                eqs.binary_link[gi, t] =
                    @constraint(mip,
                                vars.is_on[gi, t] - is_initially_on ==
                                    vars.switch_on[gi, t] - vars.switch_off[gi, t])
            else
                eqs.binary_link[gi, t] =
                    @constraint(mip,
                                vars.is_on[gi, t] - vars.is_on[gi, t-1] ==
                                    vars.switch_on[gi, t] - vars.switch_off[gi, t])
            end

            # Cannot switch on and off at the same time
            eqs.switch_on_off[gi, t] =
                @constraint(mip, vars.switch_on[gi, t] + vars.switch_off[gi, t] <= 1)
        end

        # Ramp up limit
        if t == 1
            if is_initially_on == 1
                eqs.ramp_up[gi, t] =
                    @constraint(mip,
                                vars.prod_above[gi, t] + vars.reserve[gi, t] <=
                                        (g.initial_power - g.min_power[t]) + g.ramp_up_limit)
            end
        else
            eqs.ramp_up[gi, t] =
                @constraint(mip,
                            vars.prod_above[gi, t] + vars.reserve[gi, t] <=
                                    vars.prod_above[gi, t-1] + g.ramp_up_limit)
        end

        # Ramp down limit
        if t == 1
            if is_initially_on == 1
                eqs.ramp_down[gi, t] =
                    @constraint(mip,
                                vars.prod_above[gi, t] >=
                                        (g.initial_power - g.min_power[t]) - g.ramp_down_limit)
            end
        else
            eqs.ramp_down[gi, t] =
                @constraint(mip,
                            vars.prod_above[gi, t] >=
                                vars.prod_above[gi, t-1] - g.ramp_down_limit)
        end

        # Startup limit
        eqs.startup_limit[gi, t] =
            @constraint(mip,
                        vars.prod_above[gi, t] + vars.reserve[gi, t] <=
                            (g.max_power[t] - g.min_power[t]) * vars.is_on[gi, t]
                            - max(0, g.max_power[t] - g.startup_limit) * vars.switch_on[gi, t])

        # Shutdown limit
        if g.initial_power > g.shutdown_limit
            eqs.shutdown_limit[gi, 0] =
                @constraint(mip, vars.switch_off[gi, 1] <= 0)
        end
        if t < T
            eqs.shutdown_limit[gi, t] =
                @constraint(mip,
                            vars.prod_above[gi, t] <=
                                (g.max_power[t] - g.min_power[t]) * vars.is_on[gi, t]
                                - max(0, g.max_power[t] - g.shutdown_limit) * vars.switch_off[gi, t+1])
        end

        # Minimum up-time
        eqs.min_uptime[gi, t] =
            @constraint(mip,
                        sum(vars.switch_on[gi, i]
                            for i in (t - g.min_uptime + 1):t if i >= 1
                           ) <= vars.is_on[gi, t])

        # # Minimum down-time
        eqs.min_downtime[gi, t] =
            @constraint(mip,
                        sum(vars.switch_off[gi, i]
                            for i in (t - g.min_downtime + 1):t if i >= 1
                           ) <= 1 - vars.is_on[gi, t])

        # Minimum up/down-time for initial periods
        if t == 1
            if g.initial_status > 0
                eqs.min_uptime[gi, 0] =
                    @constraint(mip, sum(vars.switch_off[gi, i]
                                         for i in 1:(g.min_uptime - g.initial_status) if i <= T) == 0)
            else
                eqs.min_downtime[gi, 0] =
                    @constraint(mip, sum(vars.switch_on[gi, i]
                                         for i in 1:(g.min_downtime + g.initial_status) if i <= T) == 0)
            end
        end

        # Add to net injection expression
        add_to_expression!(exprs.net_injection[g.bus.name, t], vars.prod_above[g.name, t], 1.0)
        add_to_expression!(exprs.net_injection[g.bus.name, t], vars.is_on[g.name, t], g.min_power[t])

        # Add to reserves expression
        add_to_expression!(exprs.reserve[g.bus.name, t], vars.reserve[gi, t], 1.0)
    end
end

function build_obj_function!(model::UnitCommitmentModel)
    @objective(model.mip, Min, model.obj)
end

function build_net_injection_eqs!(model::UnitCommitmentModel)
    T = model.instance.time
    for t in 1:T, b in model.instance.buses
        net = model.vars.net_injection[b.name, t] = @variable(model.mip)
        model.eqs.net_injection_def[t, b.name] =
            @constraint(model.mip, net == model.exprs.net_injection[b.name, t])
    end
    for t in 1:T
        model.eqs.power_balance[t] =
            @constraint(model.mip, sum(model.vars.net_injection[b.name, t]
                                       for b in model.instance.buses) == 0)
    end
end

function build_reserve_eqs!(model::UnitCommitmentModel)
    reserves = model.instance.reserves
    for t in 1:model.instance.time
        model.eqs.min_reserve[t] =
            @constraint(model.mip, sum(model.exprs.reserve[b.name, t]
                                       for b in model.instance.buses) >= reserves.spinning[t])
    end
end

function enforce_transmission(;
                              model::UnitCommitmentModel,
                              violation::Violation,
                              isf::Array{Float64,2},
                              lodf::Array{Float64,2})::Nothing

    instance, mip, vars = model.instance, model.mip, model.vars
    limit::Float64 = 0.0

    if violation.outage_line === nothing
        limit = violation.monitored_line.normal_flow_limit[violation.time]
        @info @sprintf("    %8.3f MW overflow in %-5s time %3d (pre-contingency)",
                       violation.amount,
                       violation.monitored_line.name,
                       violation.time)
    else
        limit = violation.monitored_line.emergency_flow_limit[violation.time]
        @info @sprintf("    %8.3f MW overflow in %-5s time %3d (outage: line %s)",
                       violation.amount,
                       violation.monitored_line.name,
                       violation.time,
                       violation.outage_line.name)
    end

    fm = violation.monitored_line.name
    t = violation.time
    flow = @variable(mip, base_name="flow[$fm,$t]")

    overflow = vars.overflow[violation.monitored_line.name, violation.time]
    @constraint(mip, flow  <= limit + overflow)
    @constraint(mip, -flow <= limit + overflow)

    if violation.outage_line === nothing
        @constraint(mip, flow == sum(vars.net_injection[b.name, violation.time] *
                                         isf[violation.monitored_line.offset, b.offset]
                                     for b in instance.buses
                                     if b.offset > 0))
    else
        @constraint(mip, flow == sum(vars.net_injection[b.name, violation.time] * (
                                         isf[violation.monitored_line.offset, b.offset] + (
                                             lodf[violation.monitored_line.offset, violation.outage_line.offset] *
                                             isf[violation.outage_line.offset, b.offset]
                                         )
                                     )
                                     for b in instance.buses
                                     if b.offset > 0))
    end
    nothing
end

function set_variable_names!(model::UnitCommitmentModel)
    @info "Setting variable and constraint names..."
    time_varnames = @elapsed begin
        set_jump_names!(model.vars)
        set_jump_names!(model.eqs)
    end
    @info @sprintf("Set names in %.2f seconds", time_varnames)
end

function set_jump_names!(dict)
    for name in keys(dict)
        for idx in keys(dict[name])
            idx_str = join(map(string, idx), ",")
            set_name(dict[name][idx], "$name[$idx_str]")
        end
    end
end

function get_solution(model::UnitCommitmentModel)
    instance, T = model.instance, model.instance.time
    function timeseries(vars, collection)
        return OrderedDict(b.name => [round(value(vars[b.name, t]), digits=5) for t in 1:T]
                           for b in collection)
    end
    function production_cost(g)
        return [value(model.vars.is_on[g.name, t]) * g.min_power_cost[t] +
                sum(Float64[value(model.vars.segprod[g.name, t, k]) * g.cost_segments[k].cost[t]
                            for k in 1:length(g.cost_segments)])
                for t in 1:T]
    end
    function production(g)
        return [value(model.vars.is_on[g.name, t]) * g.min_power[t] +
                sum(Float64[value(model.vars.segprod[g.name, t, k])
                            for k in 1:length(g.cost_segments)])
                for t in 1:T]
    end
    function startup_cost(g)
        S = length(g.startup_categories)
        return [sum(g.startup_categories[s].cost * value(model.vars.startup[g.name, t, s])
                    for s in 1:S)
                for t in 1:T]
    end
    sol = OrderedDict()
    sol["Production (MW)"] = OrderedDict(g.name => production(g) for g in instance.units)
    sol["Production cost (\$)"] = OrderedDict(g.name => production_cost(g) for g in instance.units)
    sol["Startup cost (\$)"] = OrderedDict(g.name => startup_cost(g) for g in instance.units)
    sol["Is on"] = timeseries(model.vars.is_on, instance.units)
    sol["Switch on"] = timeseries(model.vars.switch_on, instance.units)
    sol["Switch off"] = timeseries(model.vars.switch_off, instance.units)
    sol["Reserve (MW)"] = timeseries(model.vars.reserve, instance.units)
    sol["Net injection (MW)"] = timeseries(model.vars.net_injection, instance.buses)
    sol["Load curtail (MW)"] = timeseries(model.vars.curtail, instance.buses)
    if !isempty(instance.lines)
        sol["Line overflow (MW)"] = timeseries(model.vars.overflow, instance.lines)
    end
    if !isempty(instance.price_sensitive_loads)
        sol["Price-sensitive loads (MW)"] = timeseries(model.vars.loads, instance.price_sensitive_loads)
    end
    return sol
end

function fix!(model::UnitCommitmentModel, solution)::Nothing
    vars, instance, T = model.vars, model.instance, model.instance.time
    for g in instance.units
        for t in 1:T
            is_on = round(solution["Is on"][g.name][t])
            production = round(solution["Production (MW)"][g.name][t], digits=5)
            reserve = round(solution["Reserve (MW)"][g.name][t], digits=5)
            JuMP.fix(vars.is_on[g.name, t], is_on, force=true)
            JuMP.fix(vars.prod_above[g.name, t], production - is_on * g.min_power[t], force=true)
            JuMP.fix(vars.reserve[g.name, t], reserve, force=true)
        end
    end
end

function set_warm_start!(model::UnitCommitmentModel, solution)::Nothing
    vars, instance, T = model.vars, model.instance, model.instance.time
    for g in instance.units
        for t in 1:T
            JuMP.set_start_value(vars.is_on[g.name, t], solution["Is on"][g.name][t])
            JuMP.set_start_value(vars.switch_on[g.name, t], solution["Switch on"][g.name][t])
            JuMP.set_start_value(vars.switch_off[g.name, t], solution["Switch off"][g.name][t])
        end
    end
end

function optimize!(model::UnitCommitmentModel;
                   time_limit=3600,
                   gap_limit=1e-4,
                   two_phase_gap=true,
                   )::Nothing

    function set_gap(gap)
        try
            JuMP.set_optimizer_attribute(model.mip, "MIPGap", gap)
            @info @sprintf("MIP gap tolerance set to %f", gap)
        catch
            @warn "Could not change MIP gap tolerance"
        end
    end

    instance = model.instance
    initial_time = time()

    large_gap = false
    has_transmission = (length(model.isf) > 0)

    if has_transmission && two_phase_gap
        set_gap(1e-2)
        large_gap = true
    else
        set_gap(gap_limit)
    end

    while true
        time_elapsed = time() - initial_time
        time_remaining = time_limit - time_elapsed
        if time_remaining < 0
            @info "Time limit exceeded"
            break
        end

        @info @sprintf("Setting MILP time limit to %.2f seconds", time_remaining)
        JuMP.set_time_limit_sec(model.mip, time_remaining)

        @info "Solving MILP..."
        JuMP.optimize!(model.mip)

        has_transmission || break

        violations = find_violations(model)
        if isempty(violations)
            @info "No violations found"
            if large_gap
                large_gap = false
                set_gap(gap_limit)
            else
                break
            end
        else
            enforce_transmission(model, violations)
        end
    end

    nothing
end

function find_violations(model::UnitCommitmentModel)
    instance, vars = model.instance, model.vars
    length(instance.buses) > 1 || return []
    violations = []
    @info "Verifying transmission limits..."
    time_screening = @elapsed begin
        non_slack_buses = [b for b in instance.buses if b.offset > 0]
        net_injections = [value(vars.net_injection[b.name, t])
                          for b in non_slack_buses, t in 1:instance.time]
        overflow = [value(vars.overflow[lm.name, t])
                    for lm in instance.lines, t in 1:instance.time]
        violations = UnitCommitment.find_violations(instance=instance,
                                                    net_injections=net_injections,
                                                    overflow=overflow,
                                                    isf=model.isf,
                                                    lodf=model.lodf)
    end
    @info @sprintf("Verified transmission limits in %.2f seconds", time_screening)
    return violations
end

function enforce_transmission(model::UnitCommitmentModel, violations::Array{Violation, 1})
    for v in violations
        enforce_transmission(model=model,
                             violation=v,
                             isf=model.isf,
                             lodf=model.lodf)
    end
end

export UnitCommitmentModel, build_model, get_solution, optimize!


text:
# implements a multi-core, single machine parallel value iteration solver

type ParallelSolver <: Solver

    numProcessors::Int64 # number of processors to use

    stateOrder::Vector # order for processing chunks at each iter, e.g. [[1,1000],[1000,2000]]

    maxIteratons::Int64 # number of iterations

    tolerance::Float64 # Bellman residual

    gaussSeidel::Bool # flag for gauss-seidel value iteration (regular uses x2 more memory)

    includeV::Bool # flag for including the utility array

    includeQ::Bool # flag for including the Q-matrix

    includeA::Bool # flag for including the policy

end

# over-loaded constructor
function ParallelSolver(numProcessors::Int64; stateOrder::Vector=[], maxIterations::Int64=1000,
                        tolerance::Float64=1e-3, gaussSiedel::Bool=true,
                        includeV::Bool=true, includeQ::Bool=true, includeA::Bool=true)
    return ParallelSolver(numProcessors, stateOrder, maxIterations, tolerance, gaussSiedel, includeV, includeQ, includeA)
end

# returns the utility function and the Q-matrix
function solve(solver::ParallelSolver, mdp::DiscreteMDP; verbose::Bool=false)
    nProcs   = solver.numProcessors
    maxProcs = int(CPU_CORES / 2)

    # processor restriction checks
    nProcs < 2 ? error("Less than 2 processors not allowed") : nothing
    nProcs > maxProcs ? error("Number of requested processors is too large, try $maxProcs") : nothing

    # check if a default state ordering is needed
    order = solver.stateOrder
    if isempty(order)
        nStates = numStates(mdp)
        solver.stateOrder = {1:nStates}
    end

    gs = solver.gaussSeidel
    u = []
    q = []

    # check gauss-seidel flag
    if gs
        u, q = solveGS(solver, mdp, verbose=verbose)
    else
        u, q = solveRegular(solver, mdp, verbose=verbose)
    end

    p = []

    # check policy flag
    if solver.includeA
        p = computePolicy(mdp, u)
    end

    policy = DiscretePolicy(V=u, Q=q, P=p)

    return policy
end

function solveGS(solver::ParallelSolver, mdp::DiscreteMDP; verbose::Bool=false)

    nStates  = numStates(mdp)
    nActions = numActions(mdp)

    maxIter = solver.maxIteratons
    tol     = solver.tolerance

    nProcs  = solver.numProcessors

    # number of chunks in each DP iteration
    nChunks = length(solver.stateOrder)
    order   = solver.stateOrder

    # shared utility function and Q-matrix
    util = SharedArray(Float64, (nStates), init = S -> S[localindexes(S)] = 0.0, pids = [1:nProcs])
    valQ  = SharedArray(Float64, (nActions, nStates), init = S -> S[localindexes(S)] = 0.0, pids = [1:nProcs])

    iterTime  = 0.0
    totalTime = 0.0

    for i = 1:maxIter

        tic()
        residual = 0.0

        for c = 1:nChunks

            lst = segment(nProcs, order[c])

            results = pmap(x -> (idxs = x; solveChunk(mdp, util, valQ, idxs)), lst)

            newResidual = maximum(results)
            newResidual > residual ? (residual = newResidual) : (nothing)

        end # chunk loop

        iterTime = toq();
        totalTime += iterTime
        verbose ? (@printf("Iter %i: resid = %.3e, runtime = %.3e, net runtime = %.3e\n", i, residual, iterTime, totalTime)) : nothing

        # terminate if tolerance value is reached
        if residual < tol; break; end

    end # main iteration loop

    return util, valQ

end

function solveRegular(solver::ParallelSolver, mdp::DiscreteMDP; verbose::Bool=false)

    nStates  = numStates(mdp)
    nActions = numActions(mdp)

    maxIter = solver.maxIteratons
    tol     = solver.tolerance

    nProcs  = solver.numProcessors

    # number of chunks in each DP iteration
    nChunks = length(solver.stateOrder)
    order   = solver.stateOrder

    # shared utility function and Q-matrix
    util1 = SharedArray(Float64, (nStates), init = S -> S[localindexes(S)] = 0.0, pids = [1:nProcs])
    util2 = SharedArray(Float64, (nStates), init = S -> S[localindexes(S)] = 0.0, pids = [1:nProcs])
    valQ  = SharedArray(Float64, (nActions, nStates), init = S -> S[localindexes(S)] = 0.0, pids = [1:nProcs])

    iterTime  = 0.0
    totalTime = 0.0

    uCount = 0
    lastIdx = 1
    for i = 1:maxIter
        # residual tolerance
        residual = 0.0
        uIdx = 1
        tic()
        for c = 1:nChunks
            # util array to update: 1 or 2
            uIdx = uCount % 2 + 1
            lst = segment(nProcs, order[c])

            if uIdx == 1
                # returns the residual
                results = pmap(x -> (idxs = x; solveChunk(mdp, util1, util2, valQ, idxs)), lst)
                newResidual = maximum(results)
                newResidual > residual ? (residual = newResidual) : (nothing)
                # update the old utility array in parallel
                results = pmap(x -> (idxs = x; updateChunk(util1, util2, idxs)), lst)
            else
                # returns the residual
                results = pmap(x -> (idxs = x; solveChunk(mdp, util2, util1, valQ, idxs)), lst)
                newResidual = maximum(results)
                newResidual > residual ? (residual = newResidual) : (nothing)
                # update the old utility array, this is computationally costly
                results = pmap(x -> (idxs = x; updateChunk(util2, util1, idxs)), lst)
            end

            uCount += 1
        end # chunk loop

        iterTime = toq();
        totalTime += iterTime
        verbose ? (@printf("Iter %i: resid = %.3e, runtime = %.3e, net runtime = %.3e\n", i, residual, iterTime, totalTime)) : nothing

        # terminate if tolerance value is reached
        if residual < tol; lastIdx = uIdx; break; end

    end # main iteration loop
    lastIdx == 1 ? (return util2, valQ) : (return util1, valQ)
end

# updates the shared array utility and returns the residual
# valOld is used to update, and valNew is the updated value function
function solveChunk(mdp::DiscreteMDP, valOld::SharedArray, valNew::SharedArray, valQ::SharedArray, iter::UnitRange)

    nActions = numActions(mdp)
    residual = 0.0

    for si = iter

        qHi = -Inf

        for ai = 1:nActions

            states, probs = nextStates(mdp, si, ai)
            qNow = reward(mdp, si, ai)

            for sp = 1:length(states)
                spi = states[sp]
                qNow += 0.95 * probs[sp] * valOld[spi]  # gamma = 0.95
            end # sp loop

            valQ[ai, si] = qNow

            if ai == 1 || qNow > qHi
                qHi = qNow
                valNew[si] = qHi
            end

        end # action loop

        newResidual = (valOld[si] - valNew[si])^2
        newResidual > residual ? (residual = newResidual) : (nothing)

    end # state loop
    return residual
end

# updates shared utility and Q-Matrix for gauss-seidel value iteration
function solveChunk(mdp::DiscreteMDP, util::SharedArray, valQ::SharedArray, iter::UnitRange)

    nActions = numActions(mdp)
    residual = 0.0

    for si = iter

        qHi = -Inf
        utilOld = util[si]
        utilNew = utilOld

        for ai = 1:nActions

            states, probs = nextStates(mdp, si, ai)
            qNow = reward(mdp, si, ai)

            for sp = 1:length(states)
                spi = states[sp]
                qNow += 0.95 * probs[sp] * util[spi]  # gamma = 0.95
            end # sp loop

            valQ[ai, si] = qNow

            if ai == 1 || qNow > qHi

                qHi = qNow
                util[si] = qHi

                utilNew = qHi

            end

        end # action loop

        newResidual = (utilOld - utilNew)^2
        newResidual > residual ? (residual = newResidual) : (nothing)

    end # state loop

    return residual

end

# for updating the utility array in parallel
function updateChunk(utilOld::SharedArray, utilNew::SharedArray, iter::UnitRange)
    for i = iter
        utilOld[i] = utilNew[i]
    end
    return iter
end

ext: jl


text:
export ace_high

_ace_high = true

function ace_high(bit::Bool)::Bool
    global _ace_high
    _ace_high = bit
end

"""
    ace_high(bit::Bool)::Bool

Determine if aces are highest value cards in a deck. After `ace_high(true)` aces are greater than kings.
After `ace_high(false)`, aces are lower than twos. This affects the result of `<` comparisons and sorting.

    ace_high()::Bool

Return the current status of aces (`true` means aces are high, `false` means aces are low).
"""
function ace_high()::Bool
    global _ace_high
    return _ace_high
end

import Base: (<), isless

function (<)(C::Card, D::Card)::Bool

    rC = rank(C)
    rD = rank(D)

    if ace_high()
        if rC == 1
            rC += 13
        end
        if rD == 1
            rD += 13
        end
    end

    if rC < rD
        return true
    end

    if rC > rD
        return false
    end

    # compare by suits

    sC = suit(C).s
    sD = suit(D).s

    return sC < sD
end

isless(C::Card, D::Card)::Bool = C < D


text:
using ModernGL, GeometryTypes, GLAbstraction, GLWindow, Images

# Load our texture. See "downloads.jl" to get the images.
img = load(GLAbstraction.dir("tutorials", "images", "kitten.png"))

# Create the window. This sets all the hints and makes the context current.
window = create_glcontext("Textures exercise 2", resolution=(800,600))

vao = glGenVertexArrays()
glBindVertexArray(vao)

# The positions of the vertices in our rectangles
# This is a two-rectangle solution; alternatively, one could modify
# the fragment shader to compute the reflection.
vertex_positions = Point{2,Float32}[(-0.5,  0.5),     # top-left
                                    ( 0.5,  0.5),     # top-right
                                    ( 0.5,  0.0),     # middle-right
                                    (-0.5,  0.0),     # middle-left
                                    ( 0.5, -0.5),     # bottom-right
                                    (-0.5, -0.5)]     # bottom-left
# Duplicate the middle ones
vertex_positions = vertex_positions[[1:4;4;3;5:6]]

# The texture coordinates of each vertex
vertex_texcoords = Vec2f0[(0, 0),
                          (1, 0),
                          (1, 1),
                          (0, 1)]
vertex_texcoords = vertex_texcoords[[1:4;4;3;2;1]]

# Specify how vertices are arranged into faces
elements = Face{3,UInt32,-1}[(0,1,2),
                             (2,3,0),
                             (4,5,6),
                             (6,7,4)]

# The vertex shader---note the `vert` in front of """
vertex_shader = vert"""
#version 150

in vec2 position;
in vec2 texcoord;

out vec2 Texcoord;

void main()
{
    Texcoord = texcoord;
    gl_Position = vec4(position, 0.0, 1.0);
}
"""

# The fragment shader
fragment_shader = frag"""
# version 150

in vec2 Texcoord;

out vec4 outColor;

uniform sampler2D tex;

void main()
{
    outColor = texture(tex, Texcoord);
}
"""

# Link everything together, using the corresponding shader variable as
# the Dict key
bufferdict = Dict(:position=>Buffer(vertex_positions),
                  :texcoord=>Buffer(vertex_texcoords),
                  :tex=>Texture(data(img)),
                  :indexes=>indexbuffer(elements)) # special for element buffers

ro = std_renderobject(bufferdict,
                      LazyShader(vertex_shader, fragment_shader))

# Draw until we receive a close event
while !GLFW.WindowShouldClose(window)
    glClear(GL_COLOR_BUFFER_BIT)
    GLAbstraction.render(ro)
    GLFW.SwapBuffers(window)
    GLFW.PollEvents()
    if GLFW.GetKey(window, GLFW.KEY_ESCAPE) == GLFW.PRESS
        GLFW.SetWindowShouldClose(window, true)
    end
end

GLFW.DestroyWindow(window)  # needed if you're running this from the REPL


text:
export
    geocentric_planet,
    planet_dimension,
    vsop_to_fk5,
    vsop87d_dimension,
    planet_data

#=
    The VSOP87d planetary position model

    L is ecliptic longitude
    B is ecliptic latitude
    R is radius vector from sun

=#

# This global dict will contain any vsop planet data loaded from file...

global planet_data = Dict{AbstractString,Any}()

function load_planet_data(planet)
    # if data for planet hasn't been loaded before, load some from file
    global planet_data, dname

    datapath = join([dname, "/data/"])

    if !haskey(planet_data, planet)
        planet_data[planet] = include("$(datapath)/$(lowercase(planet))_vsop_data.jl")
    end
    return planet_data[planet]
end

"""
    Return one of heliocentric ecliptic longitude, latitude and radius.
        [Meeus-1998: pg 218]

        Parameters:
            jd : Julian Day in dynamical time
            planet_data_dim: a slice of the data for L, B, or R

        Returns:
            longitude in radians, or
            latitude in radians, or
            radius in au

"""
function planet_dimension(jd, planet_data_dim)
    X = 0.0
    tauN = 1.0
    tau = jd_to_jcent(jd) / 10.0
    c = planet_data_dim
    for series in c
        seriesSum = 0
        for s in series
           (A, B, C) = s
           seriesSum = seriesSum + (A * cos(B + C * tau))
        end
        X = X + seriesSum * tauN
        tauN = tauN*tau # last one is wasted
    end

    return X
end

"""
    Return heliocentric ecliptic longitude, latitude and radius.

        Parameters:
            jd : Julian Day in dynamical time
            planet : must be one of ("Mercury", "Venus", "Earth", "Mars",
                "Jupiter", "Saturn", "Uranus", "Neptune")

        Returns:
            longitude in radians
            latitude in radians
            radius in au

"""
function vsop87d_dimension(jd, planet)
    pdata = load_planet_data(planet)
    L = planet_dimension(jd, pdata["L"])
    L = mod2pi(L)
    B = planet_dimension(jd, pdata["B"])
    R = planet_dimension(jd, pdata["R"])
    return (L, B, R)
end

# Constant terms

_k0 = deg2rad(-1.397)
_k1 = deg2rad(-0.00031)
_k2 = deg2rad(dms_to_d(0, 0, -0.09033))
_k3 = deg2rad(dms_to_d(0, 0,  0.03916))

"""
    Convert VSOP to FK5 coordinates.

    This is required only when using the full precision of the
    VSOP model.

    [Meeus-1998: pg 219]

    Parameters:
        jd : Julian Day in dynamical time
        L : longitude in radians
        B : latitude in radians

    Returns:
        corrected longitude in radians
        corrected latitude in radians

"""
function vsop_to_fk5(jd, L, B)
    T = jd_to_jcent(jd)
    L1 = polynomial([L, _k0, _k1], T)
    cosL1 = cos(L1)
    sinL1 = sin(L1)
    deltaL = _k2 + _k3 * (cosL1 + sinL1) * tan(B)
    deltaB = _k3 * (cosL1 - sinL1)
    return (mod2pi(L + deltaL), B + deltaB)
end

"""
    Calculate the equatorial coordinates of a planet

    The results will be geocentric, corrected for light-time and
    aberration.

    Parameters:
        jd : Julian Day in dynamical time
        planet : must be one of ("Mercury", "Venus", "Earth", "Mars",
            "Jupiter", "Saturn", "Uranus", "Neptune")
        deltaPsi : nutation in longitude, in radians
        epsilon : True obliquity (corrected for nutation), in radians
        delta : desired accuracy, in days

    Returns:
        right accension, in radians
        declination, in radians

"""
function geocentric_planet(jd, planet, deltaPsi, epsilon, delta)
    t = jd
    l0 = -100.0 # impossible value
    # We need to iterate to correct for light-time and aberration.
    # At most three passes through the loop always nails it.
    # Note that we move both the Earth and the other planet during
    # the iteration.
    trial = 1
    for bailout = 1:30
    	# heliocentric geometric ecliptic coordinates of the Earth
        (L0, B0, R0) = vsop87d_dimension(t, "Earth")

        # heliocentric geometric ecliptic coordinates of the planet
        (L, B, R) = vsop87d_dimension(t, planet)

        # rectangular offset
        cosB0 = cos(B0)
        cosB = cos(B)
        x = R * cosB * cos(L) - R0 * cosB0 * cos(L0)
        y = R * cosB * sin(L) - R0 * cosB0 * sin(L0)
        z = R * sin(B)        - R0 * sin(B0)

        # geocentric geometric ecliptic coordinates of the planet
        x2 = x*x
        y2 = y*y
        l = atan(y, x)
        b = atan(z, sqrt(x2 + y2))

        # distance to planet in AU
        dist = sqrt(x2 + y2 + z*z)

        # light time in days
        tau = 0.0057755183 * dist

        if abs(diff_angle(l, l0)) < 2pi * delta
            break
        end

        # adjust for light travel time and try again
        l0 = l
        t = jd - tau
    end
    if trial > 19
        error("Failed - too many tries")
    end

    # transform to FK5 ecliptic and equinox
    (l, b) = vsop_to_fk5(jd, l, b)

    # nutation in longitude
    l = l + deltaPsi

    # equatorial coordinates
    (ra, declin) = ecl_to_equ(l, b, epsilon)

    return (ra, declin)
end


text:
import Images
import TestImages
import ImageMetadata
import ImageAxes
import ImageDraw
import ImageView
import ImageTransformations
import ImageFiltering

image2Filter = TestImages.testimage("lighthouse")
#add noise to image
noiseAddition = Images.channelview(image2Filter)
noiseAddition = eltype(noiseAddition).(clamp.(noiseAddition.+ randn.() .* 0.1, 0, 1))
noisyImage = Images.colorview(Images.RGB, noiseAddition)

image2Use = Images.Gray.(noisyImage)

ImageView.imshow(image2Use)
#sigma = 0.8
#gauss = ImageFiltering.KernelFactors.gaussian(sigma)
#gaussKernel = kernelfactors((gauss, gauss))

function bilateralFunction( imageWindow )
    std_spatial = 8.0
    std_pixelVal = 0.4

    normFactor_space = 1 / (2* pi * std_spatial)
    normFactor_pixelVal = 1 /(2 * pi * std_pixelVal)

    filteredValue = 0.0
    filteredNormFactor = 0.0
    centerRowIndex = ceil(Int, length(axes(imageWindow, 1)) / 2)
    centerColIndex = ceil(Int, length(axes(imageWindow, 2)) / 2)

    for rowIndex = axes(imageWindow, 1)
        for colIndex = axes(imageWindow, 2)
            spatialDistance = abs2(rowIndex - centerRowIndex ) + abs2(colIndex - centerColIndex)
            pixelDistance   = abs(imageWindow[rowIndex, colIndex] - imageWindow[centerRowIndex,centerColIndex] )

            gauss_spatial   = normFactor_space * exp(spatialDistance / (2 * abs2(std_spatial)))
            gauss_pixel     = normFactor_pixelVal * exp(pixelDistance / (2 * abs2(std_pixelVal)))

            filteredValue +=  gauss_spatial * gauss_pixel *  imageWindow[rowIndex,colIndex]
            filteredNormFactor += gauss_spatial * gauss_pixel
        end
    end

    filteredValue = filteredValue / (filteredNormFactor)
    return filteredValue
end

prova = [5,5]
filteredImage = ImageFiltering.MapWindow.mapwindow(bilateralFunction, image2Use,prova)

ImageView.imshow(filteredImage)


text:
struct SurfaceGroup{T, S, R} <: AbstractFPGroup
    genus::Int
    boundaries::Int
    gens::Vector{T}
    relations::Vector{<:Pair{S,S}}
    rws::R
end

include("symplectic_twists.jl")

genus(S::SurfaceGroup) = S.genus

function Base.show(io::IO, S::SurfaceGroup)
    print(io, "π₁ of the orientable surface of genus $(genus(S))")
    if S.boundaries > 0
        print(io, " with $(S.boundaries) boundary components")
    end
end

function SurfaceGroup(genus::Integer, boundaries::Integer)
    @assert genus > 1

    # The (confluent) rewriting systems comes from
    # S. Hermiller, Rewriting systems for Coxeter groups
    # Journal of Pure and Applied Algebra
    # Volume 92, Issue 2, 7 March 1994, Pages 137-148
    # https://doi.org/10.1016/0022-4049(94)90019-1
    # Note: the notation is "inverted":
    # a_g of the article becomes A_g here.

    ltrs = String[]
    for i in 1:genus
        subscript = join('₀'+d for d in reverse(digits(i)))
        append!(ltrs, ["A" * subscript, "a" * subscript, "B" * subscript, "b" * subscript])
    end
    Al = Alphabet(reverse!(ltrs))

    for i in 1:genus
        subscript = join('₀'+d for d in reverse(digits(i)))
        KnuthBendix.set_inversion!(Al, "a" * subscript, "A" * subscript)
        KnuthBendix.set_inversion!(Al, "b" * subscript, "B" * subscript)
    end

    if boundaries == 0
        word = Int[]

        for i in reverse(1:genus)
            x = 4 * i
            append!(word, [x, x-2, x-1, x-3])
        end
        comms = Word(word)
        word_rels = [ comms => one(comms) ]

        rws = KnuthBendix.RewritingSystem(word_rels, KnuthBendix.RecursivePathOrder(Al))
        KnuthBendix.knuthbendix!(rws)
    elseif boundaries == 1
        S = typeof(one(Word(Int[])))
        word_rels = Pair{S, S}[]
        rws = RewritingSystem(word_rels, KnuthBendix.LenLex(Al))
    else
        throw("Not Implemented")
    end

    F = FreeGroup(alphabet(rws))
    rels = [F(lhs)=>F(rhs) for (lhs,rhs) in word_rels]

    return SurfaceGroup(genus, boundaries, KnuthBendix.letters(Al)[2:2:end], rels, rws)
end

rewriting(S::SurfaceGroup) = S.rws
KnuthBendix.alphabet(S::SurfaceGroup) = alphabet(rewriting(S))
relations(S::SurfaceGroup) = S.relations

function symplectic_twists(π₁Σ::SurfaceGroup)
    g = genus(π₁Σ)

    saut = SpecialAutomorphismGroup(FreeGroup(2g), maxrules=100)

    Aij  = [SymplecticMappingClass(saut, :A, i, j) for i in 1:g for j in 1:g if i≠j]

    Bij  = [SymplecticMappingClass(saut, :B, i, j) for i in 1:g for j in 1:g if i≠j]

    mBij = [SymplecticMappingClass(saut, :B, i, j, minus=true) for i in 1:g for j in 1:g if i≠j]

    Bii  = [SymplecticMappingClass(saut, :B, i, i) for i in 1:g]

    mBii = [SymplecticMappingClass(saut, :B, i, i, minus=true) for i in 1:g]

    return [Aij; Bij; mBij; Bii; mBii]
end

KnuthBendix.alphabet(G::AutomorphismGroup{<:SurfaceGroup}) = rewriting(G)

function AutomorphismGroup(π₁Σ::SurfaceGroup; kwargs...)
    S = vcat(symplectic_twists(π₁Σ)...)
    A = Alphabet(S)

    # this is to fix the definitions of symplectic twists:
    # with i->gens(π₁Σ, i) the corresponding automorphisms return
    # reversed words
    domain = ntuple(i->inv(gens(π₁Σ, i)), 2genus(π₁Σ))
    return AutomorphismGroup(π₁Σ, S, A, domain)
end


text:
module Kelpie

import EzXML: link!, EzXML

export html
export html_div
export html_element

const HTML_ELEMENTS = [
    :head,
    :link,
    :meta,
    :style,
    :title,
    :body,
    :article,
    :aside,
    :footer,
    :header,
    :h1,
    :h2,
    :h3,
    :h4,
    :h5,
    :h6,
    :main,
    :nav,
    :section,
    :blockquote,
    :dd,
    # :div, Special case: we don't want to override the Base.div function
    :dl,
    :dt,
    :hr,
    :li,
    :ol,
    :p,
    :ul,
    :a,
    :b,
    :br,
    :code,
    :em,
    :i,
    :span,
    :strong,
    :sub,
    :sup,
    :u,
    :img,
    :iframe,
    :script,
    :table,
    :tbody,
    :td,
    :tfoot,
    :th,
    :thead,
    :tr,
    :button,
]

"""
    link_or_text!(node, content...)

Converts each `content` to an `EzXML.TextNode` if it isn't already an `EzXML.Node` and links
it to `node`.

```jldoctest
julia> import EzXML: ElementNode, prettyprint; import Kelpie: link_or_text!

julia> prettyprint(link_or_text!(ElementNode("div"), ElementNode("br")))
<div>
  <br/>
</div>

julia> prettyprint(link_or_text!(ElementNode("h1"), "Kelpie.jl"))
<h1>Kelpie.jl</h1>

julia> prettyprint(link_or_text!(ElementNode("div"), "The end", ElementNode("hr")))
<div>The end<hr/></div>
```
"""
function link_or_text!(node, content...)
    for con in content
        if !isnothing(con)
            if typeof(con) <: EzXML.Node
                link!(node, con)
            else
                link!(node, EzXML.TextNode(string(con)))
            end #if
        end #if
    end #for

    return node
end #function

"""
    html_element(name, content...=nothing; kwargs...)

Creates a new `EzXML.Node` with name `name`, containing `content`, and with attributes
specified by `kwargs`.

# Example

```jldoctest
julia> import EzXML: prettyprint

julia> prettyprint(html_element("img"; src="https://millironx.com/images/charolette.jpg"))
<img src="https://millironx.com/images/charolette.jpg"/>

julia> prettyprint(html_element("span", "MillironX"; class="label-primary"))
<span class="label-primary">MillironX</span>
```
"""
function html_element(name::AbstractString, content...=nothing; kwargs...)
    el = EzXML.ElementNode(name)

    for (key, value) in kwargs
        link!(el, EzXML.AttributeNode(replace(string(key), "_" => "-"), string(value)))
    end #for

    if !isnothing(content)
        link_or_text!(el, content...)
    end #if

    return el
end #function

"""
    html(content...=nothing; kwargs...)

Creates a new HTML document filled with `content`.

# Example

```jldoctest
julia> import EzXML: prettyprint

julia> prettyprint(html())
<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
<!DOCTYPE html SYSTEM "about:legacy-compat">
<html/>
```
"""
function html(content...=nothing; kwargs...)
    doc = EzXML.HTMLDocumentNode("about:legacy-compat", nothing)
    node = html_element("html", content...; kwargs...)
    link!(doc, node)
    return doc
end #function

for symbol in HTML_ELEMENTS
    name = string(symbol)

    @eval function $symbol(content...=nothing; kwargs...)
        return html_element($name, content...; kwargs...)
    end #function

    @eval export $symbol
end

function html_div(content...=nothing; kwargs...)
    return html_element("div", content...; kwargs...)
end #function

end #module


text:
# Load Julia packages (libraries)

using CSV, DataFrames, Statistics
using AxisKeys
using StanSample

df = CSV.read(joinpath(@__DIR__, "..", "..", "..", "data", "chimpanzees.csv"), DataFrame);

# Define the Stan language model

stan10_4 = "
data{
    int N;
    int N_actors;
    int pulled_left[N];
    int prosoc_left[N];
    int condition[N];
    int actor[N];
}
parameters{
    vector[N_actors] a;
    real bp;
    real bpC;
}
model{
    vector[N] p;
    bpC ~ normal( 0 , 10 );
    bp ~ normal( 0 , 10 );
    a ~ normal( 0 , 10 );
    for ( i in 1:504 ) {
        p[i] = a[actor[i]] + (bp + bpC * condition[i]) * prosoc_left[i];
        p[i] = inv_logit(p[i]);
    }
    pulled_left ~ binomial( 1 , p );
}
";

data = (N = size(df, 1), N_actors = length(unique(df.actor)),
    actor = df.actor, pulled_left = df.pulled_left,
    prosoc_left = df.prosoc_left, condition = df.condition);

# Sample using cmdstan

m10_4s = SampleModel("m10.4s", stan10_4)
rc10_4s = stan_sample(m10_4s; data);

if success(rc10_4s)
    chns = read_samples(m10_4s, :keyedarray)

    # Display the chns

    chns |> display
    println()

    # Display the keys

    axiskeys(chns) |> display
    println()

    axiskeys(chns, :param) |> display
    println()

    chns(chain=1) |> display
    println()

    chns[:, 1, 8] |> display
    println()

    chns(chain=1, param=:bp) |> display
    println()

    chns(chain=[1, 3], param=[:bp, :bpC]) |> display
    println()

    # Select all elements starting with 'a'

    chns_a = matrix(chns, :a)
    chns_a |> display
    println()

    mean(chns_a, dims=1) |> display
    println()

    typeof(chns_a.data) |> display
    println()

    ndraws_a, nchains_a, nparams_a = size(chns_a)
    chn_a = reshape(chns_a, ndraws_a*nchains_a, nparams_a)
    println()

    for row in eachrow(chn_a)
        # ...
    end

    # Or use read_samples to only use chains 2 and 4 using the chains kwarg.

    chns2 = read_samples(m10_4s; chains=[2, 4])
    chns2_a = matrix(chns2, :a)
    ndraws2_a, nchains2_a, nparams2_a = size(chns2_a)
    chn2_a = reshape(chns2_a, ndraws2_a*nchains2_a, nparams2_a)
    mean(chns2_a, dims=1) |> display

end


text:
module MCIntegrate

export SMC, IMC

using Distributions, Statistics

function SMC(; f::Function, X₀::Real, X₁::Real, Ŋ::Integer = 100000)
    ΔX = X₁ - X₀
    fₓ = f.(ΔX .* rand(Ŋ) .+ X₀)
    Avg = mean(fₓ)
    ∫fₓdx = Avg * ΔX
    Δ = std(fₓ) / √Ŋ
    return ∫fₓdx, Δ
end

function IMC(; f::Function, g::Function, ∫gₓdx::Real, Đₓ::Function, Ŋ::Integer = 100000)
    fₓ╱gₓ = x -> f(x) / g(x)
    Samples = fₓ╱gₓ.(Đₓ.(rand(Ŋ)))
    Avg = mean(Samples)
    σ = std(Samples)
    ∫fₓdx = Avg * ∫gₓdx
    Δ = ∫gₓdx * σ / √Ŋ
    return ∫fₓdx, Δ
end

end


text:
using PyPlot

# Herein we compare the approximate kernel with the exact kernel.

ntime = 1000
time = collect(linspace(0,10,ntime))

kex = zeros(ntime,ntime)
kap = zeros(ntime,ntime)
beta = 3.0
w = 1.e-9
#nc = 4
#amp = [2.953823582534888+1.5e-3,-0.04961614737314587,-1.6900719772140012,-0.21413545794774125]
#b = [1.589510614489265,3.4729345651560872,1.8698492189074047,6.363451205766626]
#omega = [0.3329119766576696,4.5999414581901386,1.7115890219170917,2.205939216106191]
nc = 3
amp = [2.7518074963750734,-0.3546256365636662,-1.3968344706386124]
b = [1.5312499407621147,4.753504292338284,1.7967431129812566]
omega = [0.34148888178746056,2.042571639590954,1.7641684073510928]
# Matern 3/2 kernel:
#nc = 2
#amp = [-99.0,100.0]
#b = [1.0,0.99]
#omega = [0.0,0.0]
for i=1:ntime
  for j=1:ntime
    z= abs(time[i]-time[j])*beta
#    kex[i,j]=(1.0+z)*exp(-z)
    kex[i,j]=exp(-z^2/2.0)
#    kap[i,j]=2.670143*exp(-1.500258*z)*cos(0.346731*z)-
#             0.458705*exp(-4.125588*z)*cos(1.390912*z)-
#             1.211438*exp(-1.744602*z)*cos(1.798436*z)
     for k=1:nc
       kap[i,j] +=amp[k]*exp(-b[k]*z)*cos(omega[k]*z)
     end
  end
  kex[i,i]+=w
  kap[i,i]+=w
end
println("Determinant: ",det(kex))
eigex,eigvec = eig(kex)
println("Minimum eigenvalue: ",minimum(eigex))
eigap,eigvec = eig(kap)
println("Minimum eigenvalue: ",minimum(eigap))
# Cholesky decomposition (square root) of kernels:
sqrt_kex = chol(kex)
sqrt_kap = chol(kap)

# Random normal deviates:
y=randn(ntime)
# Compute correlated noise component:
corrnoise_ex=*(transpose(sqrt_kex),y);
corrnoise_ap=*(transpose(sqrt_kap),y);

plot(time,corrnoise_ex)
plot(time,corrnoise_ap,"r.")


text:
@testset "MutableLinkedList" begin

    @testset "empty list" begin
        l1 = MutableLinkedList{Int}()
        @test MutableLinkedList() == MutableLinkedList{Any}()
        @test iterate(l1) === nothing
        @test isempty(l1)
        @test length(l1) == 0
        @test lastindex(l1) == 0
        @test collect(l1) == Int[]
        @test eltype(l1) == Int
        @test eltype(typeof(l1)) == Int
        @test_throws ArgumentError pop!(l1)
        @test_throws ArgumentError popfirst!(l1)
    end

    @testset "core functionality" begin
        n = 10

        @testset "push back / pop back" begin
            l = MutableLinkedList{Int}()

            @testset "push back" begin
                for i = 1:n
                    push!(l, i)
                    @test last(l) == i
                    if i > 4
                        @test getindex(l, i) == i
                        @test getindex(l, 1:floor(Int, i/2)) == MutableLinkedList{Int}(1:floor(Int, i/2)...)
                        @test l[1:floor(Int, i/2)] == MutableLinkedList{Int}(1:floor(Int, i/2)...)
                        setindex!(l, 0, i - 2)
                        @test l == MutableLinkedList{Int}(1:i-3..., 0, i-1:i...)
                        setindex!(l, i - 2, i - 2)
                    end
                    @test lastindex(l) == i
                    @test length(l) == i
                    @test isempty(l) == false
                    for (j, k) in enumerate(l)
                        @test j == k
                    end
                    if i > 3
                        l1 = MutableLinkedList{Int32}(1:i...)
                        io = IOBuffer()
                        @test sprint(io -> show(io, iterate(l1))) == "(1, DataStructures.ListNode{Int32}(2))"
                        @test sprint(io -> show(io, iterate(l1, l1.node.next.next))) == "(2, DataStructures.ListNode{Int32}(3))"
                    end
                    cl = collect(l)
                    @test isa(cl, Vector{Int})
                    @test cl == collect(1:i)
                end
            end

            @testset "pop back" begin
                for i = 1:n
                    x = pop!(l)
                    @test length(l) == n - i
                    @test isempty(l) == (i == n)
                    @test x == n - i + 1
                    cl = collect(l)
                    @test cl == collect(1:n-i)
                end
            end
        end

        @testset "push front / pop front" begin
            l = MutableLinkedList{Int}()

            @testset "push front" begin
                for i = 1:n
                    pushfirst!(l, i)
                    @test first(l) == i
                    @test length(l) == i
                    @test isempty(l) == false
                    cl = collect(l)
                    @test isa(cl, Vector{Int})
                    @test cl == collect(i:-1:1)
                end
            end

            @testset "pop front" begin
                for i = 1:n
                    x = popfirst!(l)
                    @test length(l) == n - i
                    @test isempty(l) == (i == n)
                    @test x == n - i + 1
                    cl = collect(l)
                    @test cl == collect(n-i:-1:1)
                end
            end

        end

        @testset "append / delete / copy / reverse" begin
            for i = 1:n
                l = MutableLinkedList{Int}(1:n...)

                @testset "append" begin
                    l2 = MutableLinkedList{Int}(n+1:2n...)
                    append!(l, l2)
                    @test l == MutableLinkedList{Int}(1:2n...)
                    l3 = MutableLinkedList{Int}(1:n...)
                    append!(l3, n+1:2n...)
                    @test l3 == MutableLinkedList{Int}(1:2n...)
                end

                @testset "delete" begin
                    delete!(l, n+1:2n)
                    @test l == MutableLinkedList{Int}(1:n...)
                    for i = n:-1:1
                        delete!(l, i)
                    end
                    @test l == MutableLinkedList{Int}()
                    l = MutableLinkedList{Int}(1:n...)
                    @test_throws BoundsError delete!(l, n-1:2n)
                    @test_throws BoundsError delete!(l, 2n)
                end

                @testset "copy" begin
                    l2 = copy(l)
                    @test l == l2
                end

                @testset "reverse" begin
                    l2 = MutableLinkedList{Int}(n:-1:1...)
                    @test l == reverse(l2)
                end
            end
        end

        @testset "map / filter" begin
            for i = 1:n
                @testset "map" begin
                    l = MutableLinkedList{Int}(1:n...)
                    @test map(x -> 2x, l) == MutableLinkedList{Int}(2:2:2n...)
                    l2 = MutableLinkedList{Float64}()
                    @test map(x -> x*im, l2) == MutableLinkedList{Complex{Float64}}()
                    @test map(Int32, l2) == MutableLinkedList{Int32}()
                    f(x) = x % 2 == 0 ? convert(Int8, x) : convert(Float16, x)
                    @test typeof(map(f, l)) == MutableLinkedList{Real}
                end

                @testset "filter" begin
                    l = MutableLinkedList{Int}(1:n...)
                    @test filter(x -> x % 2 == 0, l) == MutableLinkedList{Int}(2:2:n...)
                end

                @testset "show" begin
                    l = MutableLinkedList{Int32}(1:n...)
                    io = IOBuffer()
                    @test sprint(io -> show(io, l.node.next)) == "$(typeof(l.node.next))($(l.node.next.data))"
                    io1 = IOBuffer()
                    write(io1, "MutableLinkedList{Int32}(");
                    write(io1, join(l, ", "));
                    write(io1, ")")
                    seekstart(io1)
                    @test sprint(io -> show(io, l)) == read(io1, String)
                end
            end
        end
    end

    @testset "random operations" begin
        l = MutableLinkedList{Int}()
        r = Int[]
        m = 100

        for k = 1 : m
            la = rand(1:20)
            x = rand(1:1000, la)

            for i = 1 : la
                if rand(Bool)
                    push!(r, x[i])
                    push!(l, x[i])
                else
                    pushfirst!(r, x[i])
                    pushfirst!(l, x[i])
                end
            end

            @test length(l) == length(r)
            @test collect(l) == r

            lr = rand(1:length(r))
            for i = 1 : lr
                if rand(Bool)
                    pop!(r)
                    pop!(l)
                else
                    popfirst!(r)
                    popfirst!(l)
                end
            end

            @test length(l) == length(r)
            @test collect(l) == r
        end
    end
end


text:
export SliceIterator

using QXContexts.ComputeGraphs
"""
Data structure that implements iterator interface for iterating over multi-dimensional
objects with configurable start and end points.
"""
struct SliceIterator
    iter::CartesianIndices
    start::Int
    stop::Int
end

"""
    SliceIterator(dims::Vector{Int}, start::Int=1, stop::Int=-1)

Constructor for slice iterator which takes dimensions as argument
"""
function SliceIterator(dims::Vector{Int}, start::Int=1, stop::Int=-1)
    iter = CartesianIndices(Tuple(dims))
    if stop == -1 stop = length(iter) end
    SliceIterator(iter, start, stop)
end

"""
    SliceIterator(cg::ComputeGraph, args..., num_bonds=nothing)

Constructor to initialise instance from compute graph object. Optional num_bonds argument
allows number of bonds used to be limited.
"""
function SliceIterator(cg::ComputeGraph, args...; max_slices=nothing)
    slice_params = params(cg, ViewCommand)
    num_slices = (max_slices === nothing) ? length(slice_params) : min(max_slices, length(slice_params))
    dims = map(x -> slice_params[Symbol("v$(x)")], 1:num_slices)
    SliceIterator(dims, args...)
end

"""
    SliceIterator(si::SliceIterator, start::Int=1, end::Int=-1)

Constructor to initialise instance from an existing instance.
"""
function SliceIterator(si::SliceIterator, start::Int=1, stop::Int=-1)
    new_start = si.start + start - 1
    new_stop = si.start + stop - 1
    @assert new_stop <= length(si.iter) "Stop index out of range"
    SliceIterator(si.iter, new_start, new_stop)
end

"""Implement required iterator interface functions"""

Base.iterate(a::SliceIterator) = length(a) == 0 ? nothing : (Int[Tuple(a.iter[a.start])...], a.start)
Base.iterate(a::SliceIterator, state) = length(a) <= (state + 1 - a.start) ? nothing : (Int[Tuple(a.iter[state + 1])...], state + 1)
Base.length(a::SliceIterator) = a.stop - a.start + 1
Base.eltype(::SliceIterator) = Vector{Int}
Base.getindex(a::SliceIterator, i...) = collect(Tuple(a.iter[i...]))


text:
module MetaImageFormat

using ImageCore
using ImageCore.FixedPointNumbers
using ImageCore.Colors
using ImageCore.MappedArrays
using ColorVectorSpace
# Other packages
using AxisArrays
using AxisArrays: HasAxes
using ImageAxes
using FileIO
# import Libz   # what compression does MetaImageFormat use? seems undocumented
using LinearAlgebra

string2type = Dict(
    "MET_UCHAR" => UInt8,
    "MET_SHORT" => Int16,
    "MET_USHORT" => UInt16,
    "MET_INT" => Int32,
    "MET_UINT" => UInt32,
    "MET_FLOAT" => Float32,
    "MET_UCHAR_ARRAY" => RGB{N0f8},
    "MET_USHORT_ARRAY" => RGB{N0f16},
)

# The opposite of string2type
type2string(::Type{UInt8})  = "MET_UCHAR"
type2string(::Type{Int16})  = "MET_SHORT"
type2string(::Type{UInt16}) = "MET_USHORT"
type2string(::Type{Int32})  = "MET_INT"
type2string(::Type{UInt32}) = "MET_UINT"
type2string(::Type{Float32}) = "MET_FLOAT"
type2string(::Type{RGB{T}}) where T = type2string(T)*"_ARRAY"

function space2axes(str)
    syms = Symbol.(tuple(str...))
    if length(syms) != 3 || !allunique(syms)
        @warn "AnatomicalOrientation string unrecognized, got $str. Using xyz instead."
        (:x, :y, :z)
    else
        syms
    end
end

axes2space(syms::NTuple{3,Symbol}) = join(syms)

# Parsing
VTuple{T} = Tuple{Vararg{T}}          # space-delimited tuple: 80 150

# This should list anything that DOESN'T parse to a string
# We replace True/False with the corresponding Bool
const parse_type = Dict(
    "NDims"=>Int,
    "ID"=>Int,
    "ParentID"=>Int,
    "BinaryData"=>Bool,
    "CompressedData"=>Bool,
    "BinaryDataByteOrderMSB"=>Bool,
    "ElementByteOrderMSB"=>Bool,
    "Color"=>RGBA{Float32},
    "Origin"=>VTuple{Float32},
    "Position"=>VTuple{Float32},
    "Offset"=>VTuple{Float32},
    "Orientation"=>Matrix{Float32},
    "Rotation"=>Matrix{Float32},
    "TransformMatrix"=>Matrix{Float32},
    "ElementSpacing"=>VTuple{Float32},
    "DimSize"=>VTuple{Int},
    "HeaderSize"=>Int,
    "SequenceID"=>VTuple{Int},
    "ElementMin"=>Float32,
    "ElementMax"=>Float32,
    "ElementNumberOfChannels"=>Int,
    "ElementSize"=>VTuple{Float32},
)

# const fieldorder = ["ObjectType", "NDims", "BinaryData",
#                     "BinaryDataByteOrderMSB", "ElementByteOrderMSB",
#                     "CompressedData",
#                     "TransformMatrix", "Offset", "CenterOfRotation",
#                     "sizes", "spacings", "space directions", "kinds",
#                     "centers", "centerings", "thickness",
#                     "axis mins", "axismins", "axis maxs", "axismaxs",
#                     "labels", "units",
#                     "min", "max", "old min", "oldmin", "old max", "oldmax",
#                     "block size", "blocksize", "endian", "encoding",
#                     "space units", "space origin", "measurement frame",
#                     "line skip", "lineskip", "byte skip", "byteskip",
#                     "sample units", "sampleunits",
#                     "data file", "datafile"]

const mymsb = ENDIAN_BOM == 0x01020304

# Don't extend FileIO.load
# Set mode to "r+" if you want to be able to modify values in the
# image and have them update in the disk file
function load(f::File{format"MetaImage"}, args...; mode="r", mmap=:auto)
    open(f, mode) do io
        load(io, args...; mode=mode, mmap=mmap)
    end
end

function load(io::Stream{format"MetaImage"}, Tuser::Type=Any; mode="r", mmap=:auto)
    # Assemble all the information about the array we're about to
    # read: element type, size, and the "meaning" of axes
    header = parse_header(io)
    Traw, need_bswap = raw_eltype(header)
    szraw = (header["DimSize"]...,)  # "sizes" may change in outer_eltype!, grab it now
    T, nd, perm = outer_eltype!(header, Traw, Tuser)
    axs = get_axes(header, nd)
    sz = get_size(axs)

    # Read the data
    iodata = find_datafile(io, header; mode=mode)
    compressed = get(header, "CompressedData", false)
    if compressed
        error("The MetaImage format doesn't appear to specifiy the compression algorithm, so not supported")
        # iodata = Libz.ZlibInflateInputStream(iodata)
    end

    can_mmap = !compressed

    if mmap == true && (!can_mmap)
        error("Cannot use memory-mapped for reading a compressed file")
    end

    # Use memory-mapping for large files
    do_mmap = can_mmap && (prod(szraw) > 10^8) && (mmap == :auto)
    do_mmap |= can_mmap && (mmap == true)

    if !compressed
        szraw = checked_size(Traw, szraw, sz, iodata)
    end

    if do_mmap
        A = Mmap.mmap(iodata, Array{Traw,length(szraw)}, szraw, position(iodata);
                      grow=false)
        if need_bswap
            f = mode == "r+" ? (bswap, bswap) : bswap
            A = mappedarray(f, A)
        end
    else
        A = Array{Traw}(undef, szraw...)
        read!(iodata, A)

        if need_bswap
            A = [bswap(a) for a in A]
        end
    end

    if perm == ()
        if T != eltype(A)
            A = need_bswap ? A = mappedarray(x->T(x), A) : reshape(reinterpret(T, vec(A)), sz)
        end
    else
        A = PermutedDimsArray(A, perm)
        if T<:Color
            A = colorview(T, A)
        end
    end

    isa(axs, Dims) ? A : AxisArray(A, axs)
end

# function save(f::File{format"MetaImage"}, img::AbstractArray; props::Dict = Dict{String,Any}(), keyvals=nothing, comments=nothing)
#     open(f, "w") do io
#         write(io, magic(format"NRRD"))
#         save(io, img; props=props, keyvals=keyvals, comments=comments)
#     end
# end

# function save{T}(io::Stream{format"NRRD"}, img::AbstractArray{T}; props::Dict = Dict{String,Any}(), keyvals=nothing, comments=nothing)
#     axs = axisinfo(img)
#     header = headerinfo(T, axs)
#     header_eltype!(header, T)
#     # copy fields from props to override those in header
#     for (k, v) in props
#         header[k] = v
#     end
#     v = version(header, !(keyvals==nothing || isempty(keyvals)))
#     write_header(io, v, header, keyvals, comments)
#     datafilename = get(props, "datafile", "")
#     if isempty(datafilename)
#         datafilename = get(props, "data file", "")
#     end
#     if isempty(datafilename)
#         nrrd_write(io, img)
#     else
#         println(io, "data file: ", datafilename)
#         if !get(props, "headeronly", false)
#             open(datafilename, "w") do file
#                 nrrd_write(file, img)
#             end
#         end
#     end
# end

# axisinfo(img) = axisinfo(HasAxes(img), img)
# axisinfo(::HasAxes{true}, img) = axes(img)
# axisinfo(::HasAxes, img) = size(img)

### Interpreting header settings

"""
    arraytype!(header) -> T, axs, perm, need_bswap

Analyze the `header` dictionary to extract the element-type `T`, size
or axes information `axs`, the permutation `perm` (if any) that julia
should use for "wrapping" the read data, and a boolean `need_bswap`
indicating whether the data need to be byte-swapped (to account for
differences in endianness). `T` includes any color information (in
which case a dimension of the array will be "consumed"). `axs` will be
a Dims-tuple in simple cases, or an `Axes` tuple (from AxisArrays.jl)
if dimensions are labeled or have their spatial information
(pixelspacing, spacedirections, etc) specified. `perm` is the
permutation needed to move the color data to the first dimension, or
an empty tuple if no permutation is required.

This function may modify the `header` dictionary (the reason for the !
in the name), so make a copy first if necessary.
"""
function arraytype!(header)
    Traw, need_bswap = raw_eltype(header)
    T, nd, perm = outer_eltype!(header, Traw)
    axs = get_axes(header, nd)
    T, axs, perm, need_bswap
end

"""
    arraytype(filename)

Parse MetaImage header and call `arraytype!(header)`. See
`arraytype!` for information about the return values.
"""
function arraytype(filename)
    header = parse_header(filename)
    arraytype!(header)
end

# function headerinfo(T, axs)
#     header = Dict{String,Any}()
#     Traw = raw_eltype(T)
#     header["type"] = type2string(Traw)
#     header["endian"] = myendian()
#     header["encoding"] = "raw"
#     if T <: Gray
#         val = gray(one(T))
#         val = isa(val, FixedPoint) ? reinterpret(val) : val
#         header["sample units"] = string("gray ", val)
#     elseif T <: Union{RGB,RGBA}
#         val = red(one(T))
#         val = isa(val, FixedPoint) ? reinterpret(val) : val
#         valfmt = isa(val, Integer) ? 'd' : 'f'
#         if T <: RGB
#             colstr = "rgb"
#             valfmtstr = "(%$valfmt,%$valfmt,%$valfmt)"
#             vals = (val,val,val)
#         else
#             colstr = "rgba"
#             valfmtstr = "(%$valfmt,%$valfmt,%$valfmt,%$valfmt)"
#             vals = (val,val,val,val)
#         end
#         fmtstr = "%s $valfmtstr"
#         x = (colstr, vals...)
#         header["sample units"] = @eval @sprintf($fmtstr, $(x...))
#     end
#     # Do the axes information
#     header["dimension"] = length(axs)
#     if isa(axs, Base.Indices)
#         axs = map(length, axs)
#     end
#     if isa(axs, Dims)
#         header["sizes"] = [axs...]
#     else
#         # axs is an Axis-tuple
#         header["sizes"] = [map(length, axs)...]
#         axnames = map(ax->axisnames(ax)[1], axs)
#         isspace = map(s->!startswith(string(s), "time"), axnames)
#         if haskey(axes2space, axnames)
#             header["space"] = axes2space[axnames]
#         else
#             header["space dimension"] = sum(isspace)
#         end
#         header["kinds"] = [isspc ? "domain" : "time" for isspc in isspace]
#         if !all(isdefaultname, axnames)
#             header["labels"] = [string(s) for s in axnames]
#         end
#         rng = map(ax->axisvalues(ax)[1], axs)
#         stepval = map(step, rng)
#         unitstr = map(x->isa(x, Quantity) ? string(unit(x)) : "", stepval)
#         spacing = map(x->isa(x, Quantity) ? ustrip(x) : x,        stepval)
#         if !all(x->x=="", unitstr)
#             header["units"] = [unitstr...]
#         end
#         if !all(x->x==1, spacing)
#             header["spacings"] = [spacing...]
#         end
#         origin = map(x->isa(x, Quantity) ? ustrip(x) : x, map(first, rng))
#         if any(x->x!=0, origin)
#             header["space origin"] = [origin[[isspace...]]...]
#         end
#     end
#     # Adjust the axes for color
#     if T <: Colorant && !(T <: AbstractGray)
#         header["dimension"] = length(axs)+1
#         unshift!(header["sizes"], length(T))
#         if haskey(header, "spacings")
#             unshift!(header["spacings"], NaN)
#         end
#         if haskey(header, "labels")
#             unshift!(header["labels"], lowercase(string(T.name.name)))
#         end
#         if haskey(header, "units")
#             unshift!(header["units"], "")
#         end
#         if !haskey(header, "kinds")
#             header["kinds"] = ["domain" for d = 1:length(axs)]
#         end
#         if T <: RGB
#             colkind = "RGB-color"
#         elseif T <: HSV
#             colkind = "HSV-color"
#         elseif T <: XYZ
#             colkind = "XYZ-color"
#         elseif T <: RGBA
#             colkind = "RGBA-color"
#         else
#             colkind = string(length(T), "-color")
#         end
#         unshift!(header["kinds"], colkind)
#     end
#     header
# end

# function version(header, has_keyvalue::Bool=false)
#     for n = length(version_reqs):-1:1
#         vr = version_reqs[n]
#         for f in vr
#             if haskey(header, f)
#                 return n
#             end
#         end
#     end
#     has_keyvalue ? 2 : 1
# end

"""
    raw_eltype(header) -> Traw, need_bswap
    raw_eltype(::Type{T}) -> Traw

Get the "basic" element type of the data, e.g., `UInt16` or
`Float32`.

This function does not try to determine whether the image is color
(`Traw` does not contain any color information), nor does it try to
interpret `Traw` as a `Normed` type.

See also: outer_eltype!, fixedtype.
"""
function raw_eltype(header)
    Traw = string2type[header["ElementType"]]
    msb = haskey(header, "BinaryDataByteOrderMSB") ? header["BinaryDataByteOrderMSB"] :
          haskey(header, "ElementByteOrderMSB") ? header["ElementByteOrderMSB"] : false
    need_bswap = msb != mymsb && sizeof(Traw) > 1
    Traw, need_bswap
end

#raw_eltype{C<:Colorant}(::Type{C}) = raw_eltype(eltype(C))
raw_eltype(::Type{C}) where C <: Colorant = raw_eltype(eltype(C))
raw_eltype(::Type{T}) where T <:FixedPoint = FixedPointNumbers.rawtype(T)
raw_eltype(::Type{T}) where T = T

"""
    fixedtype(Traw, header) -> Tu

Attempt to interpret type `Traw` in terms of FixedPoint numbers.
If `Traw` cannot be interpreted as `Normed`, `Tu = Traw`.
"""
function fixedtype(::Type{Traw}, header) where Traw <: Unsigned
    Traw == UInt8 && return N0f8
    if Traw == UInt16
        maxval = UInt16(get(header, "ElementMax", typemax(UInt16)))
        maxval == 0xffff && return N0f16
        maxval == 0x3fff && return N2f14
        maxval == 0x3fff && return N2f14
        maxval == 0x0fff && return N4f12
        maxval == 0x03ff && return N6f10
        maxval == 0x00ff && return N8f8
        error("unknown maxval $maxval")
    end
    Traw
end
fixedtype(::Type{Traw}, header) where Traw = Traw

function fixedtype_max(::Type{Traw}, mx) where Traw <: Unsigned
    fmx = log2(mx+1)
    if round(fmx) == fmx
        return Normed{Traw,round(Int,fmx)}
    end
    Traw
end

# header_eltype!(header, ::Type) = header
# function header_eltype!{T<:FixedPoint}(header, ::Type{T})
#     header["sample units"] = string("gray ", reinterpret(one(T)))
#     header
# end
# function header_eltype!{C<:Colorant}(header, ::Type{C})
#     _header_eltype!(header, C, eltype(C))
#     header
# end
# function _header_eltype!{C<:AbstractGray,T<:FixedPoint}(header, ::Type{C}, ::Type{T})
#     header["sample units"] = string("gray ", reinterpret(one(T)))
# end
# function _header_eltype!{C<:AbstractGray,T}(header, ::Type{C}, ::Type{T})
#     header["sample units"] = string("gray ", one(T))
# end
# function _header_eltype!{C<:AbstractRGB,T<:FixedPoint}(header, ::Type{C}, ::Type{T})
#     o = reinterpret(one(T))
#     header["sample units"] = "rgb ($o,$o,$o)"
# end
# function _header_eltype!{C<:AbstractRGB,T}(header, ::Type{C}, ::Type{T})
#     o = one(T)
#     header["sample units"] = "rgb ($o,$o,$o)"
# end
# function _header_eltype!{C<:XYZ}(header, ::Type{C}, ::Type)
#     header["sample units"] = "xyz (95.047,100.000,108.883)"
# end
# function _header_eltype!{C<:HSV}(header, ::Type{C}, ::Type)
#     header["sample units"] = "hsv (360, 0, 1)"
# end

"""
    colorant_eltype(C, T) -> Tc

Return a valid "inner" element type `Tc` for colorant type `C`. When
`T` != `Tc`, values must be "converted" before they can be interpreted
as type `C`.
"""
colorant_eltype(::Type{C}, ::Type{T}) where {C <: Colorant, T <: AbstractFloat} = C{T}
colorant_eltype(::Type{C}, ::Type{T}) where {C <: Colorant, T} = C{Float32}

"""
    UnknownColor{T,N}

An unknown Color. This type gets returned when one of the "kind"
settings is "3-color" or "4-color".
"""
struct UnknownColor{T,N} <: Color{T,N}
    col::NTuple{N,T}
end

"""
    outer_eltype!(header, Traw) -> T, nd, perm

Extract the julia array `eltype` `T`, the number of dimensions `nd`
**excluding** color/complex/vector/matrix element data, and any
permutation needed to put the eltype dimension first. Any dimensions
in the header corresponding to color (or if "kind" is set to one of
the vector types) will be "consumed" upon exit. `Traw` is the
element type as determined by `raw_eltype`.

See also: raw_eltype.
"""
function outer_eltype!(header, Traw, Tuser=Any)
    nd = header["NDims"]
    sz = header["DimSize"]
    length(sz) == nd || error("parsing of sizes: $sz is inconsistent with $nd dimensions")
    perm = ()
    T = fixedtype(Traw, header)
    T, nd, perm
end

function get_axes(header, nd)
    if haskey(header, "AnatomicalOrientation")
        axnames = space2axes(header["AnatomicalOrientation"])
    else
        axnames = nd == 2 ? (:x, :y) :
                  nd == 3 ? (:x, :y, :z) : error("$nd dimension not handled")
    end
    tform = haskey(header, "TransformMatrix") ? header["TransformMatrix"] :
            haskey(header, "Rotation") ? header["Rotation"] : Matrix(I, nd, nd)
    if tform != Matrix(I, nd, nd)
        error("rotations aren't yet handled, got $tform")
    end
    offset = get(header, "Offset", zeros(nd))
    spacing = get(header, "ElementSpacing", ones(nd))
    sz = header["DimSize"]
    rng = map((o, s, l) -> range(o, step=s, length=l), offset, spacing, sz)
    (map((name, r) -> Axis{name}(r), axnames, rng)...,)
end

get_size(sz::Dims) = sz
get_size(axs) = map(length, axs)

### Parsing

"""
    header = parse_header(io)

Parse the MetaImage header (an .mhd file). `io` should be positioned
at the beginning of the header file.

See also: write_header.
"""
function parse_header(io)
    header = Dict{String, Any}()
    # Read until we encounter a blank line, which is the separator
    # between the header and data
    line = strip(readline(io))
    while !isempty(line)
        idx = findfirst(isequal('='), line)
        idx == 0 && error("no colon found in $line")
        key, value = strip(line[1:idx-1]), strip(line[idx+1:end])
        T = get(parse_type, key, String)
        header[key] = meta_parse(T, value)
        line = strip(readline(io))
    end
    header
end

parse_header(s::Stream{format"MetaImage"}) = parse_header(stream(s))

function parse_header(filename::AbstractString)
    f = File{format"MetaImage"}(filename)
    open(f) do io
        parse_header(io)
    end
end

# """
#     write_header(io, version, header, [keyvals, [comments]])

# Write an NRRD header, as the top of the .nrrd or the separate .nhdr
# file. `io` should be positioned just after the initial "NRRD" in the
# file. This writes the header and a blank line, so that at the end `io`
# is positioned at the first byte of the data (if present).

# Note that if you're writing a header for a "detached" data file
# (separate .nhdr and .raw files), `header` should contain a "data file"
# or "datafile" field storing the name of the .raw file.

# Inputs:
# - `version` is a 4-character string, e.g., "0002", giving the NRRD version of the header, or a integer corresponding to a recognized NRRD header version number
# - `header` is a `Dict{String,Any}` of `field=>setting` pairs (as returned by `parse_header`)
# - `keyvals` is a `Dict{String,String}` containing `key=>value` pairs (NRRD0002 or higher, lines like key:=value; many NRRD files do not contain any of these)
# - `comments` is an array containing lines of the header that began with `#` (but with the `#` and leading whitespace stripped out)

# See also: parse_header.
# """
# function write_header(io::IO, version, header, keyvals=nothing, comments=nothing)
#     writeversionstr(io, version)

#     for fn in fieldorder
#         if haskey(header, fn)
#             print(io, fn, ": ")
#             T = get(parse_type, fn, String)
#             nrrd_format(io, T, header[fn])
#             println(io, "")
#         end
#     end
#     if keyvals != nothing
#         for (k,v) in keyvals
#             println(io, k, ":=", v)
#         end
#     end
#     if comments != nothing
#         for c in comments
#             println(io, "# ", c)
#         end
#     end
#     println(io)
#     nothing
# end
# write_header(s::Stream{format"NRRD"}, args...) = write_header(stream(s), args...)

meta_parse(::Type{T}, str) where T = parse(T, str)  # fallback
meta_parse(::Type{String}, str) = str

function meta_parse(::Type{Bool}, str)
    str == "True" && return true
    str == "False" && return false
    error("unrecognized boolean string $str")
end

function meta_parse(::Type{VTuple{T}}, s::AbstractString) where T
    ss = split(s)
    v = Vector{T}(undef, length(ss))
    for i = 1:length(ss)
        v[i] = meta_parse(T, ss[i])
    end
    return v
end

function meta_parse(::Type{Matrix{T}}, s::AbstractString) where T
    v = meta_parse(VTuple{T}, s)
    l = length(v)
    sqrtl = round(Int, sqrt(l))
    sqrtl^2 == l || error("cannot interpret transformation matrix $v")
    reshape(v, sqrtl, sqrtl)
end

meta_format(io, ::Type{T}, x) where T = print(io, x)

function meta_format(io, ::Type{VTuple{T}}, container) where T
    isfirst = true
    for x in container
        if isfirst
            isfirst = false
        else
            print(io, ' ')
        end
        nrrd_format(io, T, x)
    end
    nothing
end

### Utilities

function chksize(sz, sztarget)
    sz == sztarget || error("dimension size should be $sztarget, got $sz")
    nothing
end

function spacings(spacedirections)
    spacings = zeros(length(spacedirections))
    can_convert = true
    d = 1
    for (i, v) in enumerate(spacedirections)
        if v == "none"
            spacings[i] = NaN
        else
            n = 0
            for (j,x) in enumerate(v)
                if x != 0
                    if j == d   # make certain they are in order
                        spacings[i] = x
                        d += 1
                    else
                        can_convert = false
                    end
                    n += 1
                end
            end
            can_convert &= n < 2
        end
    end
    can_convert, spacings
end

### File-related functions

# Adjust the array size if the file is not big enough for reading to succeed
function checked_size(Traw, szraw, sz, iodata)
    cpos = position(iodata)
    datalen = div(filesize(stat(iodata)) - cpos, sizeof(Traw))
    if datalen < prod(szraw)
        if szraw != sz
            # If we've dropped a dimension due to "absorbing"
            # color, etc, don't try to figure out how to correct
            # it, just punt to the user
            error("data length $datalen too small for size $sz array of $T")
        end
        # If the data are smaller than the header suggests, read in as
        # much "complete" data as are available
        sznew = [szraw...]
        strds = [1;cumprod(sznew)]
        k = length(sznew)
        sznew[k] = div(datalen, strds[k])
        while sznew[k] == 0 && k > 1
            pop!(sznew)
            k -= 1
            sznew[k] = div(datalen, strds[k])
        end
        tsznew = (sznew...,)
        warn("header indicates an array size $szraw, but the file size is consistent with at most $tsznew")
        szraw = tsznew
    end
    szraw
end

function stream2name(s::IO)
    name = s.name
    if !startswith(name, "<file ")
        error("io name ", name, " doesn't fit expected pattern")
    end
    name[7:end-1]
end

function find_datafile(iodata::IOStream, header; mode="r", path=nothing)
    fdata = header["ElementDataFile"]
    if path == nothing
        # Check current directory first
        if isfile(fdata)
            iodata = open(fdata)
        else
            path = dirname(stream2name(iodata))
            iodata = open(joinpath(path, fdata), mode)
        end
    else
        iodata = open(joinpath(path, fdata), mode)
    end
    iodata
end

function find_datafile(s::Stream{format"MetaImage"}, header; mode="r")
    find_datafile(stream(s), header; mode=mode)
end

# nrrd_write(io, A::AxisArray) = nrrd_write(io, A.data)
# nrrd_write(io, A::AbstractArray) = nrrd_write_elty(io, A)
# nrrd_write_elty{C<:Colorant}(io, A::AbstractArray{C}) = nrrd_write_elty(io, channelview(A))
# nrrd_write_elty{T<:Fixed}(io, A::AbstractArray{T}) = nrrd_write_elty(io, rawview(A))
# nrrd_write_elty(io, A::AbstractArray) = write(io, A)

end


text:
using Distributions

abstract DifferentialEvolutionOpt <: PopulationOptimizer

name(de::DifferentialEvolutionOpt) = de.name

DE_DefaultOptions = {
  "f" => 0.6,
  "cr" => 0.7,
  "NumParents" => 3,
  "SamplerRadius" => 8,
}

type DiffEvoOpt <: DifferentialEvolutionOpt
  name::ASCIIString

  # A population is a matrix of floats.
  population::Array{Float64, 2}

  # A search space is defined by the min and max values (in tuples) for each
  # of its dimenions. The dimension is the length of an individual, i.e. the
  # number of Float64 values in it.
  search_space::SearchSpace

  # Options
  options

  # Set of functions that together define a specific DE strategy.
  sample::Function
  mutate::Function
  crossover::Function
  bound::Function

  function DiffEvoOpt(name, pop, ss, options, sample, mutate, crossover, bound)
    new(name, pop, ss, merge(DE_DefaultOptions, options), sample, mutate, crossover, bound)
  end
end

popsize(opt::DifferentialEvolutionOpt) = Base.size(population(opt),1)
fconst(de::DifferentialEvolutionOpt, i) = de.options["f"]
crconst(de::DifferentialEvolutionOpt, i) = de.options["cr"]

# Ask for a new candidate object to be evaluated, and a list of individuals
# it should be ranked with. The individuals are supplied as an array of tuples
# with the individual and its index.
function ask(de::DifferentialEvolutionOpt)
  # Sample parents and target
  numparents = de.options["NumParents"]
  indices = de.sample(de, 1 + numparents)
  #print("indices = "); show(indices); println("")
  parent_indices = indices[1:numparents]
  #print("parent_indices = "); show(parent_indices); println("")
  target_index = indices[end]
  #print("target_index = "); show(target_index); println("")
  target = de.population[target_index,:]
  #print("target = "); show(target); println("")

  donor = de.mutate(de, target_index, parent_indices)
  #print("donor = "); show(donor); println("")

  # Crossover donor and target
  trial = de.crossover(de, target, target_index, donor)

  # Bound the trial vector according to search space bounds
  trial = de.bound(trial, target, de.search_space)

  # Return the candidates that should be ranked as tuples including their
  # population indices.
  return [(trial, target_index), (target, target_index)]
end

function random_sampler(de::DifferentialEvolutionOpt, numSamples)
  sample(1:popsize(de), numSamples; replace = false)
end

# This implements a "trivial geography" similar to Spector and Kline (2006)
# by first sampling an individual randomly and then selecting additional
# individuals for the same tournament within a certain deme of limited size
# for the sub-sequent individuals in the population. The version we implement
# here is from:
#  I. Harvey, "The Microbial Genetic Algorithm", in Advances in Artificial Life
#  Darwin Meets von Neumann, Springer, 2011.
function radius_limited_sampler(de::DifferentialEvolutionOpt, numSamples)
  # The radius must be at least as big as the number of samples + 2 so that
  # there is something to sample from.
  radius = max(de.options["SamplerRadius"], numSamples+2)
  psize = popsize(de)
  deme_start = rand(1:psize)
  indices = sample(deme_start:(deme_start+radius-1), numSamples; replace = false)
  # Ensure they are not out of bounds by wrapping over at the end.
  map(indices) do index
    if index > psize
      mod(index, psize) + 1 # We have to increase by 1 since Julia arrays start indices at 1
    else
      index
    end
  end
end

# DE/rand/1 mutation strategy
function de_mutation_rand_1(de::DifferentialEvolutionOpt, targetIndex, parentIndices)
  f = fconst(de, targetIndex)
  p = de.population[parentIndices,:]
  return p[3,:] + (f * (p[1,:] - p[2,:]))
end

# DE/rand/2 mutation strategy
function de_mutation_rand_2(de::DifferentialEvolutionOpt, targetIndex, parentIndices)
  f = fconst(de, targetIndex)
  p = de.population[parentIndices,:]
  return p[3,:] + (f * (p[1,:] - p[2,:])) + (f * (p[4,:] - p[5,:]))
end

# Binomial crossover for DE, i.e. DE/*/*/bin.
function de_crossover_binomial(de::DifferentialEvolutionOpt, target, targetIndex, donor)
  trial = copy(target)

  # Always ensure at least one value from donor is copied to trial vector
  jrand = rand(1:length(trial))
  trial[jrand] = donor[jrand]

  # Now crossover randomly for the rest of the indices
  switch = rand(length(trial)) .<= crconst(de, targetIndex)
  #print("switch = "); show(switch); println("")
  #print("trial = "); show(trial); println("")
  #print("donor = "); show(donor); println("")
  trial[:,switch] = donor[:,switch]
  #print("trial = "); show(trial); println("")

  return trial
end

# If we come out-of-bounds we randomly sample between the target value
# and the bound.
function rand_bound_from_target!(individual, target, searchSpace::SearchSpace)
  ssmins, ssmaxs = mins(searchSpace), maxs(searchSpace)
  for i in 1:length(individual)
    min, max = ssmins[i], ssmaxs[i]

    if individual[i] < min
      individual[i] = min + rand() * (target[i] - min)
    elseif individual[i] > ssmaxs[i]
      individual[i] = target[i] + rand() * (max - target[i])
    end
  end
  individual
end

function rand_bound_from_target!(individual, target, searchSpace::Array{(Float64,Float64),1})
  rand_bound_from_target!(individual, target, RangePerDimSearchSpace(searchSpace))
end

# Tell the optimizer about the ranking of candidates. Returns the number of
# better candidates that were inserted into the population.
function tell!(de::DiffEvoOpt,
  # archive::Archive, # Skip for now
  rankedCandidates)
  num_candidates = length(rankedCandidates)
  num_better = 0
  for i in 1:div(num_candidates, 2)
    candidate, index = rankedCandidates[i]
    if candidate != de.population[index, :]
      num_better += 1
      #print("candidate = "); show(candidate); println("")
      #print("index = "); show(index); println("")
      #print("target = "); show(de.population[index,:]); println("")
      old = de.population[index,:]
      de.population[index,:] = candidate
    end
  end
  num_better
end

# Now we can create specific DE optimizers that are commonly used in the
# literature.

# The most used DE/rand/1/bin.
function de_rand_1_bin(parameters = Dict(); sampler = random_sampler, name = "DE/rand/1/bin")
  params = Parameters(parameters, DE_DefaultOptions)
  # Ensure NumParents is 3 since de_mutation_rand_1 requires it.
  params["NumParents"] = 3
  DiffEvoOpt(name, params[:Population], params[:SearchSpace], params,
    sampler,
    de_mutation_rand_1,
    de_crossover_binomial,
    rand_bound_from_target!)
end

function de_rand_2_bin(parameters = Dict(); sampler = random_sampler, name = "DE/rand/2/bin")
  params = Parameters(parameters, DE_DefaultOptions)
  # Ensure NumParents is 5 since de_mutation_rand_2 requires it.
  params["NumParents"] = 5
  DiffEvoOpt(name, params[:Population], params[:SearchSpace], params,
    sampler,
    de_mutation_rand_2,
    de_crossover_binomial,
    rand_bound_from_target!)
end

# The most used DE/rand/1/bin with "local geography" via radius limited sampling.
function de_rand_1_bin_radiuslimited(parameters = Dict())
  de_rand_1_bin(parameters; sampler = radius_limited_sampler, name = "DE/rand/1/bin/radiuslimited")
end

function de_rand_2_bin_radiuslimited(parameters = Dict())
  de_rand_2_bin(parameters; sampler = radius_limited_sampler, name = "DE/rand/2/bin/radiuslimited")
end