text:
using Pkg
#Pkg.add("StatsBase")
using StatsBase

include("ant.jl")
include("graph.jl")

function algorithmAntColonyOptimization(graph::Graph, maxIteration:: Int64,antsCount:: Int64, pheromoneExponent::Int64, desirableExponent::Float64, pheromoneDepositFactor::Int64, evaporationRate:: Float64)

    bestFitness = Inf
    bestPath = []

    for i in 1:maxIteration
        ants = generateAntPopulation(antsCount,1)
        for j in 1:antsCount
            antReachDestination = false
            while (!antReachDestination)
                # trazimo mogucu putanju
                possiblePaths = findPossiblePaths(ants[j],graph)
                # racunamo verovatnocu
                probabilities = calculateProbabilities(graph, ants[j], possiblePaths, pheromoneExponent, desirableExponent)
                # biramo putanju
                chosenNode = sample(possiblePaths,Weights(probabilities))
                # pomeramo mrava
                push!(ants[j].route,ants[j].currentIndexNode)
                ants[j].currentIndexNode = chosenNode
                if chosenNode == graph.foodIndex
                    antReachDestination = true
                    push!(ants[j].route, ants[j].currentIndexNode)
                end
            end
            # trazenje najboljeg resenja za iteraciju
            fitness = calculateFitness(ants[j],graph)
            ants[j].fitness = fitness
            if fitness < bestFitness
                bestFitness = fitness
                bestPath = deepcopy(ants[j].route)
            end
        end
        # update feromona
        updatePheromones!(ants,graph, pheromoneDepositFactor, evaporationRate)
    end
    return bestPath, bestFitness, graph
end

function findPossiblePaths(ant:: Ant,graph:: Graph)
    possiblePaths = findall(graph.adjMatrix[ant.currentIndexNode,:].!= -1)
    if !isempty(ant.route)
        possiblePaths = deleteat!(possiblePaths,findall(possiblePaths.== ant.route[end]))
    end
    return possiblePaths
end

function calculateProbabilities(graph::Graph, ant::Ant, possiblePaths::Array{Int64,1}, pheromoneExponent, desirableExponent)
    pathSum = 0
    for possibleNode in possiblePaths
        pathSum += ((graph.pheromoneMatrix[ant.currentIndexNode, possibleNode])^pheromoneExponent) *
                    ((graph.adjMatrix[ant.currentIndexNode,possibleNode])^desirableExponent)
    end

    probabilities = Float64[]

    for possibleNode in possiblePaths
        probability = ((graph.pheromoneMatrix[ant.currentIndexNode, possibleNode]) ^ pheromoneExponent) *
        ((graph.adjMatrix[ant.currentIndexNode,possibleNode]) ^ desirableExponent) / pathSum

        push!(probabilities,probability)
    end
    return probabilities
end

function calculateFitness(ant::Ant,graph::Graph)
    pathLength = 0
    for i in 1:(length(ant.route)-1)
        pathLength += graph.adjMatrix[ant.route[i],ant.route[i+1]]
    end
    return pathLength
end

function updatePheromones!(ants, graph, pheromoneDepositFactor, evaporationRate)
    pheromoneDepositMatrix = zeros(Float64, graph.nodesCount, graph.nodesCount)

    for currentAnt in ants
        route = currentAnt.route
        for j in length(route)-1
            pheromoneDepositMatrix[route[j],route[j+1]] += pheromoneDepositFactor / currentAnt.fitness
        end
    end

    for i in 1:graph.nodesCount
        for j in 1:graph.nodesCount
           pheromone = graph.pheromoneMatrix[i,j]
           graph.pheromoneMatrix[i,j] = (1 - evaporationRate) * pheromone + pheromoneDepositMatrix[i,j]
           pheromoneDepositMatrix[j,i] = graph.pheromoneMatrix[i,j]
        end
    end
end

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
using SubmodularMaximization
using POMDPs
using MCTS

using Statistics
using Profile

default_steps = 50
grid_size = 10
horizon = SubmodularMaximization.default_horizon
iterations = 1000
information_samples = 1

function run_test(steps; print=true)
  ifprint(xs...) = if print println(xs...) end

  ifprint("Running test: ", steps, " steps")
  grid = Grid(grid_size, grid_size)
  sensor = RangingSensor()

  robot_state = random_state(grid)
  target_state = random_state(grid)

  target_states = Array{State}(undef, steps)
  target_states[1] = target_state

  robot_states = Array{State}(undef, steps)
  robot_states[1] = robot_state

  histogram_filter = Filter(grid, target_state)

  problem = SingleRobotTargetTrackingProblem(grid, sensor, horizon,
                                             [histogram_filter],
                                             num_information_samples =
                                               information_samples)

  time = @elapsed for ii = 2:steps
    # Before the target moves and the robot receives a measurement, execute robot
    # dynamics
    robot_state = solve_single_robot(problem, robot_state,
                                     n_iterations = iterations).action
    robot_states[ii] = robot_state

    # Then update the target and sample the observation
    target_state = target_dynamics(grid, target_states[ii-1])
    target_states[ii] = target_state

    range_observation = generate_observation(grid, sensor, robot_state,
                                             target_state)

    # After updating states, update the filter (in place for now)
    process_update!(histogram_filter, transition_matrix(grid))
    measurement_update!(histogram_filter, robot_state, get_states(grid), sensor,
                        grid, range_observation)
  end

  ifprint("Test duration: ", time, " seconds")
  ifprint("  ", time / (steps - 1), " seconds per step")
  ifprint("  ", time / ((steps - 1) * iterations), " seconds per MCTS sample")
  ifprint("  ", time / ((steps - 1) * iterations * information_samples),
          " seconds per information sample")
  ifprint("Note: horizon length is ", horizon)
  ifprint("Note: total number of information samples is ",
          (steps-1) * iterations * information_samples)
end

# Show timing without profiling
run_test(2)
# Timing should produce useful memory data
@time run_test(default_steps)

# Run the profiler
Profile.init(n=1000000)
@profile run_test(2, print=false)

Profile.clear()
Profile.init(n=1000000)
@profile run_test(default_steps, print=false)

println()
Profile.print(mincount=100)

nothing

ext: jl

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
using ForwardDiff, NLopt

l = x->loss(x,full=false) #from uncertainty.jl

function loss_nlopt(x,∇)
    length(∇) > 0 ? ForwardDiff.gradient!(∇, l,x) : nothing
    l(x)
end

function gimbal_optimize(x₀)
    #x₀ = [0.01,0.029,0.16,0.06,0.009,0.01,7000.0]
    #x₀ = [0.011,0.0319,0.0319,0.066,0.0081,0.011,7000.0]
    opt = Opt(:LD_MMA, 7)
    opt.lower_bounds = x₀*0.8
    opt.upper_bounds = x₀*1.2
    opt.xtol_rel = 1e-3
    #opt.maxtime = 10 #time limit in seconds
    opt.min_objective = loss_nlopt
    (minf,minx,ret) = NLopt.optimize(opt, x₀)
end

export gimbal_optimize


text:
using BenchmarkTools
using Ramnet
using MLDatasets

# Should I make this a standard matrix? I may be benchmarking the transpose operation if it's lazy
X_train = permutedims(reshape(MNIST.traintensor(), 784, :)) .> 0.5
y_train = MNIST.trainlabels()

# b_train = @benchmarkable train!(model, $X_train, $y_train) setup = (model = MultiDiscriminatorClassifier{Int64}(784, 28; seed=1))

# display(run(b_train))

# println("\n2) MultiDiscriminatorClassifier inference with the MNIST dataset(10.000 test images):")

# model = MultiDiscriminatorClassifier{Int64}(784, 28; seed=1)
# train!(model, X_train, y_train)

X_test = permutedims(reshape(MNIST.testtensor(), 784, :)) .> 0.5

# b_predict = @benchmarkable predict($model, $X_test)

# display(run(b_predict))

# println()

########

const suite = BenchmarkGroup()

# Subgroups
suite["mapper"] = BenchmarkGroup(["mappers", "mapping"])
suite["mapper"]["random mapper"] = BenchmarkGroup(["random mapping"])
suite["node"] = BenchmarkGroup()
suite["node"]["dict node"] = BenchmarkGroup()
suite["discriminator"] = BenchmarkGroup(["discriminators"])

# Random mapping tests
suite["mapper"]["random mapper"]["Instantiation"] = @benchmarkable RandomMapper(784, 28; seed=1)

mapper = RandomMapper(784, 28; seed=1)
mapper_itr = map(mapper, X_train)

suite["mapper"]["random mapper"]["Iteration"] = @benchmarkable collect(map($mapper, $X_train))
# suite["mapper"]["random mapper"]["Iteration"] = @benchmarkable begin
#     for (i, tuple) in Iterators.enumerate(map($mapper, $X_train))
#         tuples[i] = tuple
#     end
# end setup = (tuples = Array{eltype(mapper_itr)}(undef, length(mapper_itr)))

# Node tests
zero_digits = X_train[findall(y_train .== 0), :]
tuple = first(map(RandomMapper(784, 28; seed=1), zero_digits))

suite["node"]["dict node"]["train (Vector{Bool})"] = @benchmarkable Ramnet.Nodes.train!(node, tuple) setup = (node = Ramnet.Nodes.DictNode{Vector{Bool}}())
suite["node"]["dict node"]["train (BitVector)"] = @benchmarkable Ramnet.Nodes.train!(node, tuple) setup = (node = Ramnet.Nodes.DictNode{BitVector}())

paramspath = joinpath(dirname(@__FILE__), "params.json")

if isfile(paramspath)
    loadparams!(suite, BenchmarkTools.load(paramspath)[1], :evals);
else
    tune!(suite)
    BenchmarkTools.save(paramspath, params(suite));
end

max_stars_repo_name: RafaelFK/ramnet

max_issues_repo_name: RafaelFK/ramnet

max_forks_repo_name: RafaelFK/ramnet


text:
# File: barchart.jl
# Author:
# Date: 29-June-2020

function barchart()
    # import XL data
    [~,dates,returns] = getXLData('hfma_matlab_data.xlsx','CTA Index');

    n = size(returns,1); # # data points

    figure; # create figure
    sDate = datenum(dates(1)); # set start date for x-axis
    eDate = datenum(dates(end)); # set end date for x-axis
    xData = linspace(sDate,eDate,n);

    bar(xData, returns); # plot bar chart
    xlabel('Date'); # add x label
    ylabel('RoR (# )'); # add y label
    box off;
    ytick = get(gca,'YTick'); # format axes
    set(gca,'YTickLabel', sprintf('# .2f|',ytick))
    datetick('x','yyyy','keeplimits');
    return bar
end


text:
# test graph wrapper for testing @wrap_graph!,
# should forward everything to the wrapped graph
struct TestGraphWrapper{V <: Integer, V_VALS, E_VALS, G_VALS, GI <: AbstractValGraph{V, V_VALS, E_VALS, G_VALS}} <: AbstractValGraph{V, V_VALS, E_VALS, G_VALS}

    inner::GI
end

SimpleValueGraphs.wrapped_graph(gw::TestGraphWrapper) = gw.inner

SimpleValueGraphs.wrapped_graph_type(::Type{<:TestGraphWrapper{V, V_VALS, E_VALS, G_VALS, GI}}) where {V, V_VALS, E_VALS, G_VALS, GI} = GI

SimpleValueGraphs.@wrap_graph! TestGraphWrapper include=[ne, outneighbors, inneighbors, outedgevals, inedgevals, edges]

@testset "graph wrapper" begin

gi = ValDiGraph{Int8}(
        4,
        vertexval_types=(a=Int, b=String),
        vertexval_init=v -> (a=v, b="$v"),
        edgeval_types=(c=String, d=Float64),
        graphvals=(e='x', f="y")
)
add_edge!(gi, 1, 2, ("1-2", 1.2))
add_edge!(gi, 2, 1, ("2-1", 2.1))
add_edge!(gi, 1, 4, ("1-4", 1.4))
add_edge!(gi, 4, 4, ("4-4", 4.4))

gw = TestGraphWrapper(gi)

@test nv(gw) == nv(gi)

@test is_directed(typeof(gw)) == is_directed(typeof(gi))

@test has_edge(gw, 1, 2) == has_edge(gi, 1, 2)
@test has_edge(gw, 1, 3) == has_edge(gi, 1, 3)

@test get_graphval(gw, 1) == get_graphval(gi, 1)
@test get_graphval(gw, :f) == get_graphval(gi, :f)
@test get_graphval(gw, :) == get_graphval(gi, :)

@test get_vertexval(gw, 1, 1) == get_vertexval(gi, 1, 1)
@test get_vertexval(gw, 2, :a) == get_vertexval(gi, 2, :a)
@test get_vertexval(gw, 3, :) == get_vertexval(gi, 3, :)

@test get_edgeval(gw, 1, 2, 1) == get_edgeval(gi, 1, 2, 1)
@test get_edgeval(gw, 2, 1, :d) == get_edgeval(gi, 2, 1, :d)
@test get_edgeval(gw, 4, 4, :) == get_edgeval(gi, 4, 4, :)

# TODO tests for modifiers

@test ne(gw) == ne(gi)

@test outneighbors(gw, 1) == outneighbors(gi, 1)
@test outneighbors(gw, 3) == outneighbors(gi, 3)

@test inneighbors(gw, 1) == inneighbors(gi, 1)
@test inneighbors(gw, 3) == inneighbors(gi, 3)

@test outedgevals(gw, 1, 2) == outedgevals(gi, 1, 2)
@test outedgevals(gw, 2, :c) == outedgevals(gi, 2, :c)
@test outedgevals(gw, 4, :) == outedgevals(gi, 4, :)

@test inedgevals(gw, 1, 2) == inedgevals(gi, 1, 2)
@test inedgevals(gw, 2, :c) == inedgevals(gi, 2, :c)
@test inedgevals(gw, 4, :) == inedgevals(gi, 4, :)

@test edges(gw) == edges(gi)
@test edges(gw, :) == edges(gi, :)

end # testset


text:
const libdex = joinpath(dirname(@__DIR__), "deps", "libDex.so")
isfile(libdex) || error("libDex not found in $libdex. Please run `Pkg.build()`")

##########################################################################################
# Global State

# These can only be called once in life-time of program, can not re-init after fini etc
init() = @ccall libdex.dexInit()::Nothing
fini() = @ccall libdex.dexFini()::Nothing

# No reason to call these more than once in life-time of program
create_JIT() = @ccall libdex.dexCreateJIT()::Ptr{HsJIT}
destroy_JIT(jit) = NO_FREE[] || @ccall libdex.dexDestroyJIT(jit::Ptr{HsJIT})::Nothing

##########################################################################################
"An error thrown from with-in Dex"
struct DexError <: Exception
    msg::String
end
function Base.showerror(io::IO, err::DexError)
    if '\n' ∈ err.msg
        # If message is multiline then it may dend on exact alignment
        println(io, "DexError:\n", err.msg)
    else
        # If one line then short enough to happen on same line as everuthing else
        println(io, "DexError: ", err.msg)
    end
end

get_error_msg() = unsafe_string(@ccall libdex.dexGetError()::Cstring)
throw_from_dex() = throw(DexError(get_error_msg()))

create_context() = @ccall libdex.dexCreateContext()::Ptr{HsContext}
destroy_context(ctx) = NO_FREE[] || @ccall libdex.dexDestroyContext(ctx::Ptr{HsContext})::Nothing

function context(f)
    ctx = create_context()
    try
        f(ctx)
    finally
        destroy_context(ctx)
    end
end

dex_eval(ctx, str) = @ccall libdex.dexEval(ctx::Ptr{HsContext}, str::Cstring)::Ptr{HsContext}

insert(ctx, str, atm) = @ccall libdex.dexInsert(ctx::Ptr{HsContext}, str::Cstring, atm::Ptr{HsAtom})::Ptr{HsContext}
eval_expr(ctx, str) = @ccall libdex.dexEvalExpr(ctx::Ptr{HsContext}, str::Cstring)::Ptr{HsAtom}
lookup(ctx, str) = @ccall libdex.dexLookup(ctx::Ptr{HsContext}, str::Cstring)::Ptr{HsAtom}

print(atm) = unsafe_string(@ccall libdex.dexPrint(atm::Ptr{HsAtom})::Cstring)

compile(ctx, atm, jit=JIT) = @ccall libdex.dexCompile(jit::Ptr{HsJIT}, ctx::Ptr{HsContext}, atm::Ptr{HsAtom})::Ptr{NativeFunctionObj}
unload(f, jit=JIT) = NO_FREE[] || @ccall libdex.dexUnload(jit::Ptr{HsJIT}, f::Ptr{NativeFunctionObj})::Nothing

get_function_signature(f, jit=JIT) = @ccall libdex.dexGetFunctionSignature(jit::Ptr{HsJIT}, f::Ptr{NativeFunctionObj})::Ptr{NativeFunctionSignature}
free_function_signature(s) = NO_FREE[] || @ccall libdex.dexFreeFunctionSignature(s::Ptr{NativeFunctionSignature})::Nothing

to_CAtom(src, dest) = @ccall libdex.dexToCAtom(src::Ptr{HsAtom}, dest::Ptr{CAtom})::Int32
from_CAtom(src) = @ccall libdex.dexFromCAtom(src::Ptr{CAtom})::Ptr{HsAtom}

max_stars_repo_name: justinchiu/dex-lang

max_issues_repo_name: justinchiu/dex-lang

max_forks_repo_name: justinchiu/dex-lang


text:
###############################################################################
#
# Example of how to add conditional formatting to an XlsxWriter file.
#
# Conditional formatting allows you to apply a format to a cell or a
# range of cells based on certain criteria.

# Original Python Copyright 2013-2016, John McNamara, jmcnamara@cpan.org
# https://github.com/jmcnamara/XlsxWriter

using XlsxWriter

function test()
	wb = Workbook("conditional_format.xlsx")

	ws1 = add_worksheet!(wb)
	ws2 = add_worksheet!(wb)
	ws3 = add_worksheet!(wb)
	ws4 = add_worksheet!(wb)
	ws5 = add_worksheet!(wb)
	ws6 = add_worksheet!(wb)
	ws7 = add_worksheet!(wb)
	ws8 = add_worksheet!(wb)

	# Add a format. Light red fill with dark red text.
	format1 = add_format!(wb, Dict("bg_color"=> "#FFC7CE", "font_color"=> "#9C0006"))

	# Add a format. Green fill with dark green text.
	format2 = add_format!(wb, Dict("bg_color"=> "#C6EFCE", "font_color"=> "#006100"))

	# Some sample data to run the conditional formatting against.
	data = [
		[34 72 38 30 75 48 75 66 84 86];
		[6 24 1 84 54 62 60 3 26 59];
		[28 79 97 13 85 93 93 22 5 14];
		[27 71 40 17 18 79 90 93 29 47];
		[88 25 33 23 67 1 59 79 47 36];
		[24 100 20 88 29 33 38 54 54 88];
		[6 57 88 28 10 26 37 7 41 48];
		[52 78 1 96 26 45 47 33 96 36];
		[60 54 81 66 81 90 80 93 12 55];
		[70 5 46 14 71 19 66 36 41 21];
	]
	###############################################################################
	#
	# Example 1.
	#
	caption = "Cells with values >= 50 are in light red. Values < 50 are in light green."

	# Write the data.
	write!(ws1, "A1", caption)

	for row in size(data, 2)
		write_row!(ws1, row + 2, 1, data[row,:])
	end

	# Write a conditional format over a range.
	conditional_format!(ws1, "B3:K12", Dict("type"=> "cell",
											 "criteria"=> ">=",
											 "value"=> 50,
											 "format"=> format1))

	# Write another conditional format over the same range.
	conditional_format!(ws1, "B3:K12", Dict("type"=> "cell",
											 "criteria"=> "<",
											 "value"=> 50,
											 "format"=> format2))
	###############################################################################
	#
	# Example 2.
	#
	caption = "Values between 30 and 70 are in light red. Values outside that range are in light green."

	write!(ws2, "A1", caption)

	for row in size(data, 2)
		write_row!(ws2, row + 2, 1, data[row,:])
	end

	conditional_format!(ws2, "B3:K12", Dict("type"=> "cell",
											 "criteria"=> "between",
											 "minimum"=> 30,
											 "maximum"=> 70,
											 "format"=> format1))

	conditional_format!(ws2, "B3:K12", Dict("type"=> "cell",
											 "criteria"=> "not between",
											 "minimum"=> 30,
											 "maximum"=> 70,
											 "format"=> format2))
	###############################################################################
	#
	# Example 3.
	#
	caption = "Duplicate values are in light red. Unique values are in light green."

	write!(ws3, "A1", caption)

	for row in size(data, 2)
		write_row!(ws3, row + 2, 1, data[row,:])
	end
	conditional_format!(ws3, "B3:K12", Dict("type"=> "duplicate",
											 "format"=> format1))

	conditional_format!(ws3, "B3:K12", Dict("type"=> "unique",
											 "format"=> format2))
	###############################################################################
	#
	# Example 4.
	#
	caption = "Above average values are in light red. Below average values are in light green."

	write!(ws4, "A1", caption)

	for row in size(data, 2)
		write_row!(ws4, row + 2, 1, data[row,:])
	end
	conditional_format!(ws4, "B3:K12", Dict("type"=> "average",
											 "criteria"=> "above",
											 "format"=> format1))

	conditional_format!(ws4, "B3:K12", Dict("type"=> "average",
											 "criteria"=> "below",
											 "format"=> format2))

	###############################################################################
	#
	# Example 5.
	#
	caption = "Top 10 values are in light red. Bottom 10 values are in light green."

	write!(ws5, "A1", caption)

	for row in size(data, 2)
		write_row!(ws5, row + 2, 1, data[row,:])
	end
	conditional_format!(ws5, "B3:K12", Dict("type"=> "top",
											 "value"=> "10",
											 "format"=> format1))

	conditional_format!(ws5, "B3:K12", Dict("type"=> "bottom",
											 "value"=> "10",
											 "format"=> format2))
	###############################################################################
	#
	# Example 6.
	#
	caption = "Cells with values >= 50 are in light red. Values < 50 are in light green. Non-contiguous ranges."

	# Write the data.
	write!(ws6, "A1", caption)

	for row in size(data, 2)
		write_row!(ws6, row + 2, 1, data[row,:])
	end
	# Write a conditional format over a range.
	conditional_format!(ws6, "B3:K6", Dict("type"=> "cell",
											"criteria"=> ">=",
											"value"=> 50,
											"format"=> format1,
											"multi_range"=> "B3:K6 B9:K12"))

	# Write another conditional format over the same range.
	conditional_format!(ws6, "B3:K6", Dict("type"=> "cell",
											"criteria"=> "<",
											"value"=> 50,
											"format"=> format2,
											"multi_range"=> "B3:K6 B9:K12"))

	###############################################################################
	#
	# Example 7.
	#
	caption = "Examples of color scales and data bars. Default colors."

	data = collect(1:13)

	write!(ws7, "A1", caption)

	write!(ws7, "B2", "2 Color Scale")
	write!(ws7, "D2", "3 Color Scale")
	write!(ws7, "F2", "Data Bars")

	for row in 1:length(data)
		write!(ws7, row + 2, 1, data[row])
		write!(ws7, row + 2, 3, data[row])
		write!(ws7, row + 2, 5, data[row])
	end

	conditional_format!(ws7, "B3:B14", Dict("type"=> "2_color_scale"))
	conditional_format!(ws7, "D3:D14", Dict("type"=> "3_color_scale"))
	conditional_format!(ws7, "F3:F14", Dict("type"=> "data_bar"))

	###############################################################################
	#
	# Example 8.
	#
	caption = "Examples of color scales and data bars. Modified colors."

	data = collect(1:13)

	write!(ws8, "A1", caption)

	write!(ws8, "B2", "2 Color Scale")
	write!(ws8, "D2", "3 Color Scale")
	write!(ws8, "F2", "Data Bars")

	for row in 1:length(data)
		write!(ws8, row + 2, 1, data[row])
		write!(ws8, row + 2, 3, data[row])
		write!(ws8, row + 2, 5, data[row])
	end
	conditional_format!(ws8, "B3:B14", Dict("type"=> "2_color_scale",
											 "min_color"=> "#FF0000",
											 "max_color"=> "#00FF00"))

	conditional_format!(ws8, "D3:D14", Dict("type"=> "3_color_scale",
											 "min_color"=> "#C5D9F1",
											 "mid_color"=> "#8DB4E3",
											 "max_color"=> "#538ED5"))

	conditional_format!(ws8, "F3:F14", Dict("type"=> "data_bar",
											 "bar_color"=> "#63C384"))

	close(wb)
	isfile("conditional_format.xlsx")

end

test()


text:
###############################
#
# Correctness testing
#
################################

# Standard POMDP
# All values calculated using POMDPs on Julia 1.0
#
# Tiger Problem -------------------------------------------------------
# TigerPOMDP(-1.0, -100.0, 10.0, 0.85, 0.95)
# tolerance = 0.001
#
# SARSOP
# Expected: 19.372
# Simulator: Sim ETR 19.3611, Sim ETR 95% CI (17.48, 21.241)
# Evaluator: Eval ETR 19.4831, Sim ETR 95% CI (19.2147, 19.7516)
# Policy: POMDPXFiles.POMDPAlphas([-81.5975 3.01448 … 28.4025 19.3711; 28.4025 24.6954 … -81.5975 19.3711], [2, 0, 0, 1, 0])
# alphasref_tiger = [[-81.5975, 28.4025],
#         [3.01448, 24.6954],
#         [24.6954, 3.01452],
#         [28.4025, -81.5975],
#         [19.3711, 19.3711]]
#
# Baby Problem -------------------------------------------------------
# BabyPOMDP(-5.0, -10.0, 0.1, 0.8, 0.1, 0.9)
# tolerance = 0.001
#
# SARSOP:
# Expected: -16.305
# Simulator: Sim ETR -16.505, Sim ETR 95% CI (-17.075, -15.935)
# Evaluator: Eval ETR -16.384, Sim ETR 95% CI (-16.749, -16.019)
# Policy: POMDPXFiles.POMDPAlphas([-19.6749 -16.3055; -29.6749 -38.2512], [1, 0])
# alphasref_baby = [[-19.6749, -29.6749], [-16.3055, -38.2512]]

# compare RPOMDP to POMDP
using RPOMDPs, RPOMDPModels, RPOMDPToolbox
using RobustValueIteration
using Plots; gr()
const RPBVI = RobustValueIteration

# PBVI POMDPs
srand(8473272)
bs = [[b, 1-b] for b in 0.0:0.01:1.0]
prob = TigerPOMDP(0.95)
prob2 = Baby2POMDP(-5.0, -10.0, 0.9)
sr = RPBVISolver(beliefpoints = bs, max_iterations = 200)
sr2 = RPBVISolver(beliefpoints = bs, max_iterations = 200)
polr = RobustValueIteration.solve(sr, prob)
polr2 = RobustValueIteration.solve(sr2, prob2)
bur = updater(polr)
bur2 = updater(polr2)
value(polr, [0.5,0.5])
value(polr2, [0.0, 1.0])
N = 1_000
N2 = 1_000
simvals = [discounted_reward(simulate(HistoryRecorder(max_steps=1_000),
                                                prob, polr, bur)) for i = 1:N]
m = mean(simvals)
s = std(simvals)
tstar = 1.962
ci = (m - tstar * s / sqrt(N), m + tstar * s / sqrt(N))

simvals2 = [discounted_reward(simulate(HistoryRecorder(max_steps=100),
                                                prob2, polr2, bur2)) for i = 1:N2]
m2 = mean(simvals2)
s2 = std(simvals2)
tstar = 1.962
ci2 = (m2 - tstar * s2 / sqrt(N2), m2 + tstar * s2 / sqrt(N2))

using Plots; gr()
plotvals0_tiger = [maximum(dot(polr.alphas[i], [1-b, b]) for i = 1:length(polr.alphas)) for b = 0.0:0.01:1.0]
alphasref_tiger = [[-81.5975, 28.4025],
        [3.01448, 24.6954],
        [24.6954, 3.01452],
        [28.4025, -81.5975],
        [19.3711, 19.3711]]

valsref_tiger = [maximum(dot(alphasref_tiger[i], [b, 1-b]) for i = 1:length(alphasref_tiger))
        for b = 0.0:0.01:1.0]

x = collect(0:0.01:1.0)
p_tiger = Plots.plot(x, plotvals0_tiger, xticks = 0:0.1:1,
        legend = :top, label = "PBVI",
        # title = "Tiger POMDP Value Function",
        xlab = "Belief, P(State = Tiger Left)",
        ylab = "Expected Total Discounted Reward",
        line = (4, :solid, :red), alpha = 0.5)
plot!(x, valsref_tiger,
    label = "SARSOP", line = :dash, alpha = 1, linewidth = 2,
    linecolor = :blue)
fn = "value_function_tiger_2.pdf"
path = joinpath(homedir(),".julia\\v0.6\\RobustInfoPOMDP\\data\\figures\\",fn)
savefig(p_tiger, path)

plotvals0_baby = [maximum(dot(polr2.alphas[i], [b, 1-b]) for i = 1:length(polr2.alphas)) for b = 0.0:0.01:1.0]
alphasref_baby = [[-19.6749, -29.6749], [-16.3055, -38.2512]]
valsref_baby = [maximum(dot(alphasref_baby[i], [1-b, b]) for i = 1:length(alphasref_baby))
        for b = 0.0:0.01:1.0]

x = collect(0:0.01:1.0)
p_baby = Plots.plot(x, plotvals0_baby, xticks = 0:0.1:1,
        legend = :topright, label = "PBVI",
        # title = "Baby POMDP Value Function",
        xlab = "Belief, P(State = Hungry)",
        ylab = "Expected Total Discounted Reward",
        line = (4, :solid, :red, 0.5))
plot!(x, valsref_baby, label = "SARSOP",
    line = (2, :dash, :blue, 1.0))
fn = "value_function_baby_2.pdf"
path = joinpath(homedir(),".julia\\v0.6\\RobustInfoPOMDP\\data\\figures\\",fn)
savefig(p_baby, path)

#########################################
# RPBVI Robust POMDPs Tiger
# unc size = 0.001, 0.01, 0.05
srand(8473272)
nr = 5
bs = fill([[b, 1-b] for b in 0.0:0.05:1.0], nr)
ambiguity = [0.001, 0.01, 0.05, 0.1, 0.2]
maxiter = [100, 100, 100, 100, 100]
nreps = [1, 1, 500, 1, 1]
maxstep = [100, 100, 100, 100, 100]
policyvalue_tiger = Vector{Float64}(nr)
simvals_tiger = Vector{Vector{Float64}}(nr)
simps_tiger = Vector{Vector{Float64}}(nr)
simvals_tiger_nom = Vector{Vector{Float64}}(nr)
simps_tiger_nom = Vector{Vector{Float64}}(nr)
m_tiger = Vector{Float64}(nr)
m_nom_tiger = Vector{Float64}(nr)
s_tiger = Vector{Float64}(nr)
ci_tiger = Vector{Tuple{Float64, Float64}}(nr)
ci_nom_tiger = Vector{Tuple{Float64, Float64}}(nr)
plotvals_tiger = Vector{Vector{Float64}}(nr)

# plot value
for i = 1:nr
    println("Starting loop $i...")
    prob = TigerRPOMDP(0.95, ambiguity[i])
    sr = RPBVISolver(beliefpoints = bs[i], max_iterations = maxiter[i])
    println("Calculating policy...")
    polr = RobustValueIteration.solve(sr, prob)
    plotvals_tiger[i] = [maximum(dot(polr.alphas[j], [1-b, b]) for j = 1:length(polr.alphas))
        for b = 0.0:0.01:1.0]
end

# reference values
alphasref_tiger = [[-81.5975, 28.4025],
        [3.01448, 24.6954],
        [24.6954, 3.01452],
        [28.4025, -81.5975],
        [19.3711, 19.3711]]
valsref_tiger = [maximum(dot(alphasref_tiger[i], [b, 1-b]) for i = 1:length(alphasref_tiger))
        for b = 0.0:0.01:1.0]

x = collect(0:0.01:1.0)
p_tiger_robust = plot(x, [valsref_tiger, plotvals_tiger],
        xticks = 0:0.1:1, label = "SARSOP",
        xlab = "Belief, P(State = Tiger Left)",
        ylab = "Expected Total Discounted Reward",
        legend = :none,
        linestyle = [:dot :solid :solid :solid :solid :solid],
        linewidth = [4 2 2 2 2 2],
        linealpha = [1 0.9 0.7 0.5 0.3 0.1],
        linecolor = [:blue :red :red :red :red :red])
# plot!(x, plotvals0_tiger, line = (2, :red, 1.0), label = "PBVI")
tigerrobustlabels = ["SARSOP", "RPBVI: 0.001", "RPBVI: 0.01", "RPBVI: 0.05", "RPBVI: 0.1", "RPBVI: 0.2"]
l = @layout [a b{0.2w}]
p1 = p_tiger_robust
p2 = plot(x, [valsref_tiger, plotvals_tiger],
        linestyle = [:dot :solid :solid :solid :solid :solid],
        linewidth = [4 2 2 2 2 2],
        linealpha = [1 0.9 0.7 0.5 0.3 0.1],
        linecolor = [:blue :red :red :red :red :red],
        label=tigerrobustlabels, grid=false, xlims=(20,3), showaxis=false)
p_tiger_robust_out = plot(p1,p2,layout=l)
fn = "value_function_tiger_robust_4.pdf"
path = joinpath(homedir(),".julia\\v0.6\\RobustInfoPOMDP\\data\\figures\\",fn)
savefig(p_tiger_robust_out, path)

# RPBVI Robust POMDPs Baby
# unc size = 0.001, 0.1
srand(8473272)
nr = 5
bs = [[[b, 1-b] for b in 0.0:0.05:1.0],
        [[b, 1-b] for b in 0.0:0.05:1.0],
        [[b, 1-b] for b in 0.0:0.05:1.0],
        [[b, 1-b] for b in 0.0:0.05:1.0],
        [[b, 1-b] for b in 0.0:0.05:1.0]]
ambiguity = [0.001, 0.01, 0.05, 0.1, 0.2]
maxiter = [100, 100, 100, 100, 100]
nreps = [1, 1, 500, 1, 1]
maxstep = [100, 100, 100, 100, 100]
policy_baby = Vector{}(nr)
policyvalue_baby = Vector{}(nr)
simvals_baby = Vector{}(nr)
m_baby = Vector{}(nr)
s_baby = Vector{}(nr)
ci_baby = Vector{}(nr)
plotvals_baby = Vector{}(nr)

# plot value
for i = 1:nr
    println("Starting loop $i...")
    prob = Baby2RPOMDP(-5.0, -10.0, 0.9, ambiguity[i])
    sr = RPBVISolver(beliefpoints = bs[i], max_iterations = maxiter[i])
    println("Calculating policy...")
    policy_baby[i] = RobustValueIteration.solve(sr, prob)
    policyvalue_baby[i] = value(policy_baby[i], [0.0,1.0])
    plotvals_baby[i] = [maximum(dot(policy_baby[i].alphas[j], [b, 1-b]) for j = 1:length(policy_baby[i].alphas))
        for b = 0.0:0.01:1.0]
end

# reference values
alphasref_baby = [[-19.6749, -29.6749], [-16.3055, -38.2512]]
valsref_baby = [maximum(dot(alphasref_baby[i], [1-b, b]) for i = 1:length(alphasref_baby))
        for b = 0.0:0.01:1.0]

x = collect(0:0.01:1.0)
p_baby_robust = plot(x, valsref_baby, xticks = 0:0.1:1, label = "SARSOP",
        # title = "Baby POMDP Value Function",
        xlab = "Belief, P(State = Hungry)",
        ylab = "Expected Total Discounted Reward",
        legend = :bottomright,
        color = :blue,
        line = :dot,
        linealpha = 1.0,
        linewidth = 4)
# plot!(x, plotvals0_baby, line = (2, :red, 1.0), label = "PBVI")
plot!(x, plotvals_baby[1], line = (2, :red, 0.9), label = "RPBVI: 0.001")
plot!(x, plotvals_baby[2], line = (2, :red, 0.7), label = "RPBVI: 0.01")
plot!(x, plotvals_baby[3], line = (2, :red, 0.5), label = "RPBVI: 0.05")
plot!(x, plotvals_baby[4], line = (2, :red, 0.3), label = "RPBVI: 0.1")
plot!(x, plotvals_baby[5], line = (2, :red, 0.1), label = "RPBVI: 0.2")
fn = "value_function_baby_robust_3.pdf"
path = joinpath(homedir(),".julia\\v0.6\\RobustInfoPOMDP\\data\\figures\\",fn)
savefig(p_baby_robust, path)

p_baby_robust = plot(x, [valsref_baby, plotvals_baby],
        xticks = 0:0.1:1, label = "SARSOP",
        xlab = "Belief, P(State = Hungry)",
        ylab = "Expected Total Discounted Reward",
        legend = :none,
        linestyle = [:dot :solid :solid :solid :solid :solid],
        linewidth = [4 2 2 2 2 2],
        linealpha = [1 0.9 0.7 0.5 0.3 0.1],
        linecolor = [:blue :red :red :red :red :red])
# plot!(x, plotvals0_tiger, line = (2, :red, 1.0), label = "PBVI")
babyrobustlabels = ["SARSOP", "RPBVI: 0.001", "RPBVI: 0.01", "RPBVI: 0.05", "RPBVI: 0.1", "RPBVI: 0.2"]
l = @layout [a b{0.2w}]
p1 = p_baby_robust
p2 = plot(x, [valsref_baby, plotvals_baby],
        linestyle = [:dot :solid :solid :solid :solid :solid],
        linewidth = [4 2 2 2 2 2],
        linealpha = [1 0.9 0.7 0.5 0.3 0.1],
        linecolor = [:blue :red :red :red :red :red],
        label=babyrobustlabels, grid=false, xlims=(20,3), showaxis=false)
p_baby_robust_out = plot(p1,p2,layout=l)
fn = "value_function_baby_robust_4.pdf"
path = joinpath(homedir(),".julia\\v0.6\\RobustInfoPOMDP\\data\\figures\\",fn)
savefig(p_baby_robust_out, path)

############################################
#
# Baby2RPOMDP Sim Value
#
############################################
using RPOMDPs, RPOMDPModels, RPOMDPToolbox
using RobustValueIteration
using Plots; gr()
const RPBVI = RobustValueIteration

# RPBVI Robust POMDPs Baby2
# unc size = 0.001, 0.01, 0.05
srand(8473272)
nr = 5
bs = fill([[b, 1-b] for b in 0.0:0.05:1.0], nr)
ambiguity = [0.001, 0.01, 0.05, 0.1, 0.2]
maxiter = [100, 100, 100, 100, 100]
nreps = [50, 50, 50, 50, 50]
maxstep = [100, 100, 100, 100, 100]

function meanci(data::Vector{Float64})
    n = length(data)
    m = mean(data)
    s = std(data)
    tstar = 1.962
    hw = tstar * s / sqrt(n)
    m, (m - hw, m + hw)
end

nreps[1] = 200
uncsize = 0.4
rfeed = -15.0
prob_nom = Baby2POMDP(rfeed, -10.0, 0.9)
# prob_nom = TigerPOMDP(0.95)
solver_nom = RPBVISolver(beliefpoints = bs[1], max_iterations = maxiter[1])
pol_nom = RPBVI.solve(solver_nom, prob_nom)
bu_nom = updater(pol_nom)
prob_r_nom = SimpleBaby2RPOMDP(rfeed, -10.0, 0.9, uncsize)
pol_r_nom = RPBVI.solve(solver_nom, prob_r_nom)
@show pol_nom.action_map
@show pol_r_nom.action_map
# prob_r_nom = SimpleTigerRPOMDP(0.95, uncsize)
# prob_r_nom = SimpleTigerRPOMDP(0.95, uncsize)
bu_r_nom = updater(pol_r_nom)
simulator = RolloutSimulator(max_steps = maxstep[1])
simvals_nn = Array{Float64}(nr, nreps[1],2)
simvals_nr = Array{Float64}(nr, nreps[1],2)
simvals_rn = Array{Float64}(nr, nreps[1],2)
simvals_rr = Array{Float64}(nr, nreps[1],2)
for j = 1:nreps[1]
    (j % 10 == 0) && print("\rRep $j")
    # nominal policy against nominal dynamics
    simvals_nn[1,j,1], simvals_nn[1,j,2] = simulate(simulator,
        prob_nom, pol_nom, bu_nom)
    # nominal policy against worst-case dynamics
    simvals_nr[1,j,1], simvals_nr[1,j,2]  = simulate_worst(simulator,
        prob_r_nom, pol_nom, bu_nom, pol_r_nom.alphas)
    # robust policy against nominal dynamics
    simvals_rn[1,j,1], simvals_rn[1,j,2] = simulate(simulator,
        prob_nom, pol_r_nom, bu_r_nom)
    # robust policy against worst-case dynamics
    simvals_rr[1,j,1], simvals_rr[1,j,2]  = simulate_worst(simulator,
        prob_r_nom, pol_r_nom, bu_r_nom, pol_r_nom.alphas)
end
value(pol_nom, [0.0, 1.0])
value(pol_r_nom, [0.0, 1.0])
value(pol_nom, [0.5, 0.5])
value(pol_r_nom, [0.5, 0.5])
m_nn, ci_nn = meanci(simvals_nn[1,:,1])
m_nr, ci_nr = meanci(simvals_nr[1,:,1])
m_rn, ci_rn = meanci(simvals_rn[1,:,1])
m_rr, ci_rr = meanci(simvals_rr[1,:,1])

mp_nn, cip_nn = meanci(simvals_nn[1,:,2])
mp_nr, cip_nr = meanci(simvals_nr[1,:,2])
mp_rn, cip_rn = meanci(simvals_rn[1,:,2])
mp_rr, cip_rr = meanci(simvals_rr[1,:,2])

############################################
#
# RockDiagnosis Sim Value
#
############################################
using RPOMDPs, RPOMDPModels, RPOMDPToolbox
using RobustValueIteration
using SimpleProbabilitySets
# using Plots; gr()
const RPBVI = RobustValueIteration

# RPBVI Robust POMDPs Baby2
# unc size = 0.001, 0.01, 0.05
srand(8473272)
nr = 5
# nb = 30
# bs1 = [vcat(psample(zeros(2), ones(2)),zeros(2)) for i = 1:nb]
# bs2 = [vcat(zeros(2), psample(zeros(2), ones(2))) for i = 1:nb]
# bs0 = vcat(bs1, bs2)
bdelt = 0.05
nb2 = 10
bs1 = [vcat([b, 1-b],zeros(2)) for b in 0.0:bdelt:1.0]
bs2 = [vcat(zeros(2), [b, 1-b]) for b in 0.0:bdelt:1.0]
bs3 = [psample(zeros(4), ones(4)) for i = 1:nb2]
bs0 = vcat(bs1, bs2, bs3)
bs = fill(bs0, nr)
ambiguity = [0.001, 0.1, 0.2, 0.3, 0.4]
maxiter = fill(150, nr)
nreps = fill(150, nr)
maxstep = fill(150, nr)

function meanci(data::Vector{Float64})
    n = length(data)
    m = mean(data)
    s = std(data)
    tstar = 1.962
    hw = tstar * s / sqrt(n)
    m, (m - hw, m + hw)
end

# simple info reward function
ralphas1 = [[1.0, -1/3, -1/3, -1/3],
         [-1/3, 1.0, -1/3, -1/3],
         [-1/3, -1/3, 1.0, -1/3],
         [-1/3, -1/3, -1/3, 1.0]]

# complicated info reward function
vhi = -10.0
vlo = 0.1
ralphas2 = [[1.0, vhi, vhi, vhi],
        [vhi, 1.0, vhi, vhi],
        [vhi, vhi, 1.0, vhi],
        [vhi, vhi, vhi, 1.0],
        [vlo, -vlo/3, -vlo/3, -vlo/3],
        [-vlo/3, vlo, -vlo/3, -vlo/3],
        [-vlo/3, -vlo/3, vlo, -vlo/3],
        [-vlo/3, -vlo/3, -vlo/3, vlo]]

ra = ralphas2
disc = 0.95
uncsize = 0.2
prob_nom = RockIPOMDP(ra, disc)
solver_nom = RPBVISolver(beliefpoints = bs[1],
    max_iterations = maxiter[1])
pol_nom = RPBVI.solve(solver_nom, prob_nom, verbose = true)
bu_nom = updater(pol_nom)
binit_nom = initial_belief_distribution(prob_nom)
prob_r_nom = RockRIPOMDP(ra, disc, uncsize)
pol_r_nom = RPBVI.solve(solver_nom, prob_r_nom, verbose = true)
@show pol_nom.action_map
@show pol_r_nom.action_map
bu_r_nom = updater(pol_r_nom)
binit_r_nom = initial_belief_distribution(prob_r_nom)
simulator = RolloutSimulator(max_steps = maxstep[1])

#########################################
# robust vs nom (exp and worst case)
simvals_nn = Array{Float64}(nr, nreps[1],2)
simvals_nr = Array{Float64}(nr, nreps[1],2)
simvals_rn = Array{Float64}(nr, nreps[1],2)
simvals_rr = Array{Float64}(nr, nreps[1],2)
for j = 1:nreps[1]
    print("\rRep $j")
    s_nom = rand(simulator.rng, initial_state_distribution(prob_nom))
    s_r_nom = rand(simulator.rng, initial_state_distribution(prob_r_nom))
    # nominal policy against nominal dynamics
    simvals_nn[1,j,1], simvals_nn[1,j,2] = simulate(simulator,
        prob_nom, pol_nom, bu_nom, binit_nom, s_nom)
    # nominal policy against worst-case dynamics
    simvals_nr[1,j,1], simvals_nr[1,j,2]  = simulate_worst(simulator,
        prob_r_nom, pol_nom, bu_nom, binit_nom, s_r_nom, pol_r_nom.alphas)
    # robust policy against nominal dynamics
    simvals_rn[1,j,1], simvals_rn[1,j,2] = simulate(simulator,
        prob_nom, pol_r_nom, bu_r_nom, binit_r_nom, s_nom)
    # robust policy against worst-case dynamics
    simvals_rr[1,j,1], simvals_rr[1,j,2]  = simulate_worst(simulator,
        prob_r_nom, pol_r_nom, bu_r_nom, binit_r_nom, s_r_nom, pol_r_nom.alphas)
end
value(pol_nom, [0.5, 0.5, 0.0, 0.0])
value(pol_r_nom, [0.5, 0.5, 0.0, 0.0])
m_nn, ci_nn = meanci(simvals_nn[1,:,1])
m_nr, ci_nr = meanci(simvals_nr[1,:,1])
m_rn, ci_rn = meanci(simvals_rn[1,:,1])
m_rr, ci_rr = meanci(simvals_rr[1,:,1])

mp_nn, cip_nn = meanci(simvals_nn[1,:,2])
mp_nr, cip_nr = meanci(simvals_nr[1,:,2])
mp_rn, cip_rn = meanci(simvals_rn[1,:,2])
mp_rr, cip_rr = meanci(simvals_rr[1,:,2])

#################################################
#################################################
using RPOMDPs, RPOMDPModels, RPOMDPToolbox
using RobustValueIteration
using SimpleProbabilitySets
# using Plots; gr()
const RPBVI = RobustValueIteration

function meanci(data::Vector{Float64})
    n = length(data)
    m = mean(data)
    s = std(data)
    tstar = 1.962
    hw = tstar * s / sqrt(n)
    m, (m - hw, m + hw)
end

# random beliefs
nr = 1
bdelt = 0.05
nb2 = 20
bs1 = [vcat([b, 1-b],zeros(2)) for b in 0.0:bdelt:1.0]
bs2 = [vcat(zeros(2), [b, 1-b]) for b in 0.0:bdelt:1.0]
bs3 = [psample(zeros(4), ones(4)) for i = 1:nb2]
bs = vcat(bs1, bs2, bs3)
maxiter = 200
maxstep = 200
nreps = 100

# simple info reward function
ralphas1 = [[1.0, -1/3, -1/3, -1/3],
         [-1/3, 1.0, -1/3, -1/3],
         [-1/3, -1/3, 1.0, -1/3],
         [-1/3, -1/3, -1/3, 1.0]]

# complicated info reward function
vhi = -10.0
vlo = 0.1
ralphas2 = [[1.0, vhi, vhi, vhi],
        [vhi, 1.0, vhi, vhi],
        [vhi, vhi, 1.0, vhi],
        [vhi, vhi, vhi, 1.0],
        [vlo, -vlo/3, -vlo/3, -vlo/3],
        [-vlo/3, vlo, -vlo/3, -vlo/3],
        [-vlo/3, -vlo/3, vlo, -vlo/3],
        [-vlo/3, -vlo/3, -vlo/3, vlo]]

prob_nom = RockIPOMDP(ralphas1, 0.95)
solver_nom = RPBVISolver(beliefpoints = bs, max_iterations = maxiter)
pol_nom = RPBVI.solve(solver_nom, prob_nom)
println(pol_nom.action_map)
bu_nom = updater(pol_nom)
binit_nom = initial_belief_distribution(prob_nom)

srand(11)
nreps = 100
maxstep = 200
simvals_nn = Array{Float64}(nr, nreps, 2)
simulator = RolloutSimulator(max_steps = maxstep)
for j = 1:nreps
    (j % 10 == 0) && print("\rRep $j")
    s_nom = rand(simulator.rng, initial_state_distribution(prob_nom))
    # nominal policy against nominal dynamics
    simvals_nn[1,j,1], simvals_nn[1,j,2] = simulate(simulator,
        prob_nom, pol_nom, bu_nom, binit_nom, s_nom)
end

value(pol_nom, [0.5, 0.5, 0.0, 0.0])
value(pol_nom, [1.0, 0.0, 0.0, 0.0])
value(pol_nom, [0.0, 1.0, 0.0, 0.0])
value(pol_nom, [0.0, 0.0, 0.5, 0.5])
m_nn, ci_nn = meanci(simvals_nn[1,:,1])
mp_nn, cip_nn = meanci(simvals_nn[1,:,2])

xb = 0.0:0.1:1.0
anom = [action(pol_nom, DiscreteBelief(prob_nom, [b, 1-b, 0.0, 0.0])) for b in xb]

using Plots; gr()
valfunc1 = [maximum(dot(pol_nom.alphas[i], [b, 1-b, 0.0, 0.0]) for i = 1:length(pol_nom.alphas))
        for b = 0.0:0.01:1.0]
x = collect(0:0.01:1.0)
vplot = plot(x, valfunc1, xticks = 0:0.1:1, label = "Nominal",
        # title = "Tiger POMDP Value Function",
        xlab = "Belief, P(State = [b, 1-b, 0 0])",
        # xlab = "Belief, P(State = Hungry)",
        ylab = "Expected Total Discounted Reward",
        legend = :bottomleft,
        # legend = :topright,
        line = :dash,
        linealpha = 0.9)
# plot!(x, valfunc13, color = :red, linealpha = 0.8, label = "Robust: 0.001")
# plot!(x, valfunc21, color = :red, linealpha = 0.4, label = "Robust: 0.4")
# fn = string("value_function_", sname, "_", sversion, ".pdf")
# path = joinpath(homedir(),".julia\\v0.6\\RobustInfoPOMDP\\data\\figures\\",fn)
# savefig(vplot, path)

psim2 = RolloutSimulator(max_steps = 10)
rnum = 28
prob = probs[rnum]
simprob = simprobs[rnum]
sol = policies[rnum]
bu = updater(sol)
binit = initial_belief_distribution(prob)
sinit = rand(psim2.rng, initial_state_distribution(simprob))
simulate_worst(psim2, simprob, sol, bu, binit, sinit, soldynamics[rnum].alphas)

value(policies[26], initial_belief(probs[26]))
value(policies[24], initial_belief(probs[24]))
initial_belief(probs[28]) == initial_belief(probs[26])
policies[26].action_map == policies[28].action_map

max_stars_repo_name: ajkeith/RobustInfoPOMDP

max_issues_repo_name: ajkeith/RobustInfoPOMDP

max_forks_repo_name: ajkeith/RobustInfoPOMDP


text:
## Common implementation of each block:
## returns a tuple of (output, internal_derivative)

# Low pass modified with denominator K_den instead of 1.0
"""
Low Pass Filter Modified
     ┌─────────────┐
     │      K      │
u -> │ ────────────│ -> y
     │ K_den + sT  │
     └─────────────┘
"""
function low_pass_modified_mass_matrix(
    u::Z,
    y::V,
    K::Float64,
    K_den::W,
    ::Float64,
) where {V <: ACCEPTED_REAL_TYPES, W <: ACCEPTED_REAL_TYPES, Z <: ACCEPTED_REAL_TYPES}
    return y, K * u - K_den * y
end

function low_pass_modified(
    u::Z,
    y::V,
    K::Float64,
    K_den::W,
    T::Float64,
) where {V <: ACCEPTED_REAL_TYPES, W <: ACCEPTED_REAL_TYPES, Z <: ACCEPTED_REAL_TYPES}
    return y, (1.0 / T) * low_pass_modified_mass_matrix(u, y, K, K_den, T)[2]
end

"""
Low Pass Filter
     ┌────────┐
     │    K   │
u -> │ ────── │ -> y
     │ 1 + sT │
     └────────┘
"""

# Use this one if T = 0 is allowed, and let the mass matrix take care of it.
function low_pass_mass_matrix(
    u::Z,
    y::V,
    K::Float64,
    T::Float64,
) where {V <: ACCEPTED_REAL_TYPES, Z <: ACCEPTED_REAL_TYPES}
    return low_pass_modified_mass_matrix(u, y, K, 1.0, T)
end

function low_pass(
    u::Z,
    y::V,
    K::Float64,
    T::Float64,
) where {V <: ACCEPTED_REAL_TYPES, Z <: ACCEPTED_REAL_TYPES}
    return low_pass_modified(u, y, K, 1.0, T)
end

"""
Low Pass Filter with non-windup
             y_max
           /¯¯¯¯¯¯
     ┌────────┐
     │    K   │
u -> │ ────── │ -> y
     │ 1 + sT │
     └────────┘
   ______/
   y_min
"""

function low_pass_nonwindup_mass_matrix(
    u::Z,
    y::V,
    K::Float64,
    ::Float64,
    y_min::Float64,
    y_max::Float64,
) where {V <: ACCEPTED_REAL_TYPES, Z <: ACCEPTED_REAL_TYPES}
    dydt = K * u - y
    y_sat = clamp(y, y_min, y_max)
    # Non Windup logic from IEEE Std 421.5
    binary_logic = ((y >= y_max) && (dydt > 0)) || ((y <= y_min) && (dydt < 0)) ? 0.0 : 1.0
    return y_sat, binary_logic * dydt
end

# Does not accept T = 0
function low_pass_nonwindup(
    u::Z,
    y::V,
    K::Float64,
    T::Float64,
    y_min::Float64,
    y_max::Float64,
) where {V <: ACCEPTED_REAL_TYPES, Z <: ACCEPTED_REAL_TYPES}
    y_sat, dydt_scaled = low_pass_nonwindup_mass_matrix(u, y, K, T, y_min, y_max)
    return y_sat, (1.0 / T) * dydt_scaled
end

# Does not accept T = 0
function low_pass_nonwindup_ramp_limits(
    u::Z,
    y::V,
    K::Float64,
    T::Float64,
    y_min::Float64,
    y_max::Float64,
    dy_min::Float64,
    dy_max::Float64,
) where {V <: ACCEPTED_REAL_TYPES, Z <: ACCEPTED_REAL_TYPES}
    y, dydt = low_pass(u, y, K, T)
    y_sat = clamp(y, y_min, y_max)
    dydt_sat = clamp(dydt, dy_min, dy_max)
    # Non Windup logic from IEEE Std 421.5
    binary_logic =
        ((y >= y_max) && (dydt_sat > 0)) || ((y <= y_min) && (dydt_sat < 0)) ? 0.0 : 1.0
    return y_sat, binary_logic * dydt_sat
end

"""
High Pass Filter
     ┌────────┐
     │   sK   │
u -> │ ────── │ -> y
     │ 1 + sT │
     └────────┘

Internal State: x
"""

function high_pass_mass_matrix(
    u::Z,
    x::Z,
    K::Float64,
    T::Float64,
) where {Z <: ACCEPTED_REAL_TYPES}
    K_T = T < eps() ? 0.0 : (K / T)
    return x + K_T * u, -(K_T * u + x)
end

function high_pass(u::Z, x::Z, K::Float64, T::Float64) where {Z <: ACCEPTED_REAL_TYPES}
    y, dxdt_scaled = high_pass_mass_matrix(u, x, K, T)
    return y, (1.0 / T) * dxdt_scaled
end

"""
Lead-Lag Block:
     ┌───────────┐
     │   1 + sT1 │
u -> │ K ─────── │ -> y
     │   1 + sT2 │
     └───────────┘

Internal State: x
"""

function lead_lag_mass_matrix(
    u::Z,
    x::Z,
    K::Float64,
    T1::Float64,
    T2::Float64,
) where {Z <: ACCEPTED_REAL_TYPES}
    T1_T2 = T2 < eps() ? 0.0 : (T1 / T2)
    return x + (K * T1_T2) * u, K * (1 - T1_T2) * u - x
end

function lead_lag(
    u::Z,
    x::Z,
    K::Float64,
    T1::Float64,
    T2::Float64,
) where {Z <: ACCEPTED_REAL_TYPES}
    y, dxdt_scaled = lead_lag_mass_matrix(u, x, K, T1, T2)
    return y, (1.0 / T2) * dxdt_scaled
end

"""
2nd Order Low Pass Filter
     ┌──────────────────┐
     │         K        │
u -> │ ──────────────── │ -> y
     │ 1 + sT1 + s^2 T2 │
     └──────────────────┘

Internal state: x
"""

function low_pass_2nd(
    u::Z,
    x::Z,
    y::Z,
    K::Float64,
    T1::Float64,
    T2::Float64,
) where {Z <: ACCEPTED_REAL_TYPES}
    dxdt = (1.0 / T2) * (K * u - T1 * x - y)
    dydt = x
    return y, dxdt, dydt
end

function low_pass_2nd_mass_matrix(
    u::Z,
    x::Z,
    y::Z,
    K::Float64,
    T1::Float64,
    T2::Float64,
) where {Z <: ACCEPTED_REAL_TYPES}
    dxdt = K * u - T1 * x - y
    dydt = x
    return y, dxdt, dydt
end

"""
2nd Order Lead-Lag Block
     ┌──────────────────┐
     │ 1 + sT3 + s^2 T4 │
u -> │ ──────────────── │ -> y
     │ 1 + sT1 + s^2 T2 │
     └──────────────────┘

Internal States: x1, x2
"""

function lead_lag_2nd_mass_matrix(
    u::Z,
    x1::Z,
    x2::Z,
    T1::Float64,
    T2::Float64,
    T3::Float64,
    T4::Float64,
) where {Z <: ACCEPTED_REAL_TYPES}
    dx1dt = u - T1 * x1 - x2
    dx2dt = x1
    T4_T2 = T2 < eps() ? 0.0 : (T4 / T2)
    y = T4_T2 * u + (T3 - T1 * T4_T2) * x1 + (1.0 - T4_T2) * x2
    return y, dx1dt, dx2dt
end

function lead_lag_2nd(
    u::Z,
    x1::Z,
    x2::Z,
    T1::Float64,
    T2::Float64,
    T3::Float64,
    T4::Float64,
) where {Z <: ACCEPTED_REAL_TYPES}
    y, dx1dt_scaled, dx2dt = lead_lag_2nd_mass_matrix(u, x1, x2, T1, T2, T3, T4)
    return y, (1.0 / T2) * dx1dt_scaled, dx2dt
end

"""
Proportional-Integral Block
             y_max
            /¯¯¯¯¯¯
     ┌──────────┐
     │      ki  │
u -> │kp + ───  │ -> y
     │      s   │
     └──────────┘
    ______/
     y_min

Internal State: x
"""

function pi_block(u::Z, x::Z, kp::Float64, ki::Float64) where {Z <: ACCEPTED_REAL_TYPES}
    return kp * u + ki * x, u
end

function pi_block_nonwindup(
    u::Z,
    x::Z,
    kp::Float64,
    ki::Float64,
    y_min::Float64,
    y_max::Float64,
) where {Z <: ACCEPTED_REAL_TYPES}
    y, _ = pi_block(u, x, kp, ki)
    y_sat = clamp(y, y_min, y_max)
    binary_logic = y_min < y < y_max ? 1.0 : 0.0
    return y_sat, binary_logic * u
end

"""
Integrator with windup limits
                                 y_max
                                _ _ _
        ┌────────┐             /
        │    K   │   y        /
 u - -->│ ────── │ - - - - - /- - - - - --> ysat
        │   sT   │          /
        └────────┘   _ _ _ /
                     y_min
 """

function integrator_windup_mass_matrix(
    u::Z,
    y::Z,
    K::Float64,
    ::Float64,
    y_min::Float64,
    y_max::Float64,
) where {Z <: ACCEPTED_REAL_TYPES}
    dydt_scaled = K * u
    y_sat = clamp(y, y_min, y_max)
    return y_sat, dydt_scaled
end

# Does not accept T = 0
function integrator_windup(
    u::Z,
    y::Z,
    K::Float64,
    T::Float64,
    y_min::Float64,
    y_max::Float64,
) where {Z <: ACCEPTED_REAL_TYPES}
    y_sat = clamp(y, y_min, y_max)
    return y_sat, (1.0 / T) * integrator_windup_mass_matrix(u, y, K, T, y_min, y_max)[2]
end

"""
Integrator with non-windup limits
             y_max
           /¯¯¯¯¯¯
     ┌────────┐
     │    K   │
u -> │ ────── │ -> y
     │   sT   │
     └────────┘
   ______/
   y_min
"""

function integrator_nonwindup_mass_matrix(
    u::Z,
    y::Z,
    K::Float64,
    ::Float64,
    y_min::Float64,
    y_max::Float64,
) where {Z <: ACCEPTED_REAL_TYPES}
    dydt_scaled = K * u
    y_sat = clamp(y, y_min, y_max)
    upper_lim = (y >= y_max) && (dydt_scaled > 0)
    lower_lim = (y <= y_min) && (dydt_scaled < 0)
    binary_logic = upper_lim || lower_lim ? 0.0 : 1.0
    return y_sat, binary_logic * dydt_scaled
end

# Does not accept T = 0
function integrator_nonwindup(
    u::Z,
    y::Z,
    K::Float64,
    T::Float64,
    y_min::Float64,
    y_max::Float64,
) where {Z <: ACCEPTED_REAL_TYPES}
    y_sat, dydt_scaled = integrator_nonwindup_mass_matrix(u, y, K, T, y_min, y_max)
    return y_sat, (1.0 / T) * dydt_scaled
end


text:
# Autogenerated wrapper script for ripgrep_jll for x86_64-linux-musl
export rg

JLLWrappers.@generate_wrapper_header("ripgrep")
JLLWrappers.@declare_executable_product(rg)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        rg,
        "bin/rg",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using CairoMakie
using CSV
using Colors

color_lst = [colorant"#332288" colorant"#555555" colorant"#117733" colorant"#882255"]
# Blue, Grey, Green, Red, Yellow
# Colorblind friendly color palette

f = CSV.File("HopfSampling/HopfSampling.csv")

Hopf_random_sampling_trace  = f.random_sampling_trace
Hopf_optimal_sampling_trace = f.optimal_sampling_trace
Hopf_optimal_sampling_trace_ekf = f.optimal_sampling_trace_ekf
Hopf_approx_sampling_trace  = f.approx_sampling_trace
Hopf_approx_sampling_trace_ekf = f.approx_sampling_trace_ekf
Hopf_samples = f.samples

f = CSV.File("LorenzSampling/LorenzSampling.csv")

Lorenz_random_sampling_trace = f.random_sampling_trace
Lorenz_optimal_sampling_trace = f.optimal_sampling_trace
Lorenz_optimal_sampling_trace_ekf = f.optimal_sampling_trace_ekf
Lorenz_approx_sampling_trace = f.approx_sampling_trace
Lorenz_approx_sampling_trace_ekf = f.approx_sampling_trace_ekf
Lorenz_samples = f.samples

max_steps = length(Lorenz_samples)

noto_sans = "./resources/NotoSans-Regular.ttf"
noto_sans_bold = "./resources/NotoSans-Bold.ttf"

tickfontsize    = 32
labelfontsize   = 38
basewidth = 2

fig = Figure(font=noto_sans, resolution=(1600,560), figure_padding=5)

ax1 = Axis( fig[1,2],
            xticklabelsize=tickfontsize,
            yticklabelsize=tickfontsize,
            yticklabelpad=2,
            xlabel="Sample", ylabel = "Tr(Σ)",
            xlabelsize=labelfontsize,
            ylabelsize=labelfontsize,
            yscale = log10,
            yminorticksvisible = true,
            yminorgridvisible = true,
            yminorticks = IntervalsBetween(10),
            xminorticksvisible = true,
            xminorgridvisible = true,
            xminorticks = IntervalsBetween(10),
            title="Lorenz System",
            titlesize=labelfontsize)

ylims!(ax1,(1e-1,1e1))
xlims!(ax1,(0,1000))

l2 = lines!(ax1, 1:max_steps, Lorenz_optimal_sampling_trace ./ Lorenz_optimal_sampling_trace,     color=color_lst[1], linewidth=basewidth)
l3 = lines!(ax1, 1:max_steps, Lorenz_optimal_sampling_trace_ekf ./ Lorenz_optimal_sampling_trace, color=color_lst[2], linewidth=basewidth)
l1 = lines!(ax1, 1:max_steps, Lorenz_approx_sampling_trace ./ Lorenz_optimal_sampling_trace,      color=color_lst[3], linewidth=basewidth)
l4 = lines!(ax1, 1:max_steps, Lorenz_approx_sampling_trace_ekf./ Lorenz_optimal_sampling_trace,  color=color_lst[4], linewidth=basewidth)
#l5 = lines!(ax1, 1:max_steps, Lorenz_random_sampling_trace./ Lorenz_optimal_sampling_trace, color=:black)

ax2= Axis(  fig[1,1],
            xticklabelsize=tickfontsize,
            yticklabelsize=tickfontsize,
            yticklabelpad=2,
            xlabel="Sample", ylabel = "Tr(Σ)",
            xlabelsize=labelfontsize,
            ylabelsize=labelfontsize,
            yscale = log10,
            yminorticksvisible = true,
            yminorgridvisible = true,
            yminorticks = IntervalsBetween(10),
            xminorticksvisible = true,
            xminorgridvisible = true,
            xminorticks = IntervalsBetween(10),
            title="Hopf Bifurcation",
            titlesize=labelfontsize)

ylims!(ax2,(1e-1,1e1))
xlims!(ax2,(0,1000))

l2 = lines!(ax2, 1:max_steps, Hopf_optimal_sampling_trace ./ Hopf_optimal_sampling_trace,     color=color_lst[1], linewidth=basewidth)
l3 = lines!(ax2, 1:max_steps, Hopf_optimal_sampling_trace_ekf./ Hopf_optimal_sampling_trace, color=color_lst[2], linewidth=basewidth)
l1 = lines!(ax2, 1:max_steps, Hopf_approx_sampling_trace./ Hopf_optimal_sampling_trace,      color=color_lst[3], linewidth=basewidth)
l4 = lines!(ax2, 1:max_steps, Hopf_approx_sampling_trace_ekf./ Hopf_optimal_sampling_trace,  color=color_lst[4], linewidth=basewidth)

l1elem = LineElement(color=color_lst[3], linewidth=10)
l2elem = LineElement(color=color_lst[1], linewidth=10)
l3elem = LineElement(color=color_lst[2], linewidth=10)
l4elem = LineElement(color=color_lst[4], linewidth=10)

Legend(fig[2,:],
       [l1elem, l4elem, l2elem, l3elem],
       ["1D Limit Set Oracle", "1D Limit Set EKF", "Dynamic Programming Oracle", "Dynamic Programming EKF"],
       orientation=:horizontal,
       nbanks=2,
       labelsize=labelfontsize,
       colgap=40)

colgap!(fig.layout,1,60)

Label(fig[1,1,TopLeft()], "A", font=noto_sans_bold, textsize = 40, halign=:left, valign=:bottom, padding=(25,0,10,0))
Label(fig[1,2,TopLeft()], "B", font=noto_sans_bold, textsize = 40, halign=:left, valign=:bottom, padding=(25,0,10,0))

# Generate Timing Plot

color_lst = [colorant"#332288" colorant"#555555" colorant"#117733" colorant"#882255"]

f = CSV.File("DimensionScaling/DimensionScaling.csv")

crlb_trace = f.crlb_trace
random_crlb_trace = f.random_crlb_trace

noto_sans = "./resources/NotoSans-Regular.ttf"

ax = Axis(  fig[1,3],
            xticklabelsize=tickfontsize,
            yticklabelsize=tickfontsize,
            yticklabelpad=2,
            xlabel="Ambient Dimension", ylabel = "Tr(Σ)",
            xlabelsize=labelfontsize,
            ylabelsize=labelfontsize,
#            yscale = log10,
#            xscale = log10,
            title = "Dimension Scaling",
            titlesize=labelfontsize,
            yminorticksvisible = true,
            yminorgridvisible = true,
            yminorticks = IntervalsBetween(10),
            xminorticksvisible = true,
            xminorgridvisible = true,
            xminorticks = IntervalsBetween(10))

scaling = sum(crlb_trace)/length(crlb_trace)
l1 = lines!(ax, 2 .+ (1:length(crlb_trace)), crlb_trace/scaling,        color = color_lst[4], linewidth=2)
l2 = lines!(ax, 2 .+ (1:length(crlb_trace)), random_crlb_trace/scaling, color = color_lst[1], linewidth=2)

xlims!(ax,(0,52))

l1elem = LineElement(color=color_lst[4], linewidth=10)
l2elem = LineElement(color=color_lst[1], linewidth=10)

Legend(fig[2,3], [l2elem,l1elem], ["Random Sampling","1D Limit Set"], labelsize=labelfontsize)
Label(fig[1,3,TopLeft()], "C", font=noto_sans_bold, textsize = 40, halign=:left, valign=:bottom, padding=(25,0,10,0))

colgap!(fig.layout,2,150)

save("out/FinalSim4.pdf",fig)


text:
get_sform(x::NIVolume) = get_sform(x.header)
function get_sform(hdr::NIfTI1Header)
    if hdr.sform_code > 0
        return @inbounds Float32[
            hdr.srow_x[1]  hdr.srow_x[2]  hdr.srow_x[3]  hdr.srow_x[4]
            hdr.srow_y[1]  hdr.srow_y[2]  hdr.srow_y[3]  hdr.srow_y[4]
            hdr.srow_z[1]  hdr.srow_z[2]  hdr.srow_z[3]  hdr.srow_z[4]
            0              0              0              1
        ]
    else
        # TODO there's no sform in this case, should we just throw an error?
        return nothing
    end
end

get_qform(x::NIVolume) = get_qform(x.header)
function get_qform(hdr::NIfTI1Header)
    if hdr.qform_code <= 0
        return @inbounds Float32[
            hdr.pixdim[2]   0              0              0
            0               hdr.pixdim[3]  0              0
            0               0              hdr.pixdim[4]  0
            0               0              0              1
        ]
    else
        dx = hdr.pixdim[2]
        dy = hdr.pixdim[3]
        # aka qfac left handedness
        if hdr.pixdim[1] < 0
            dz = -hdr.pixdim[4]
        else
            dz = hdr.pixdim[4]
        end
        b = hdr.quatern_b
        c = hdr.quatern_c
        d = hdr.quatern_d
        b2 = b*b
        c2 = c*c
        d2 = d*d
        a = 1 - b2 - c2 - d2
        if a < 1.e-7
            a = 1 / sqrt(b2 + c2 + d2)
            b *= a
            c *= a
            d *= a       # normalize (b,c,d) vector
            a = zero(a)  # a = 0 ==> 180 degree rotation
        else
            a = sqrt(a)   # angle = 2*arccos(a)
        end
        return Float32[
            (a*a+b*b-c*c-d*d)*dx   (2*b*c-2*a*d)*dy       (2*b*d+2*a*c)*dz      hdr.qoffset_x
            (2*b*c+2*a*d)*dx       (a*a+c*c-b*b-d*d)*dy   (2*c*d-2*a*b)*dz      hdr.qoffset_y
            (2*b*d-2*a*c)*dx       (2*c*d+2*a*b)*dy       (a*a+d*d-c*c-b*b)*dz  hdr.qoffset_z
            0                      0                      0                     1
        ]
    end
end

getaffine(x::NIVolume) = getaffine(x.header)

# Convert a NIfTI header to a 4x4 affine transformation matrix
function getaffine(hdr::NIfTI1Header)
    if hdr.sform_code > 0
        return get_sform(hdr)
    else
        return get_qform(hdr)
    end
end

# Set affine matrix of NIfTI header
function setaffine(h::NIfTI1Header, affine::Array{T,2}) where {T}
    size(affine, 1) == size(affine, 2) == 4 ||
        error("affine matrix must be 4x4")
    affine[4, 1] == affine[4, 2] == affine[4, 3] == 0 && affine[4, 4] == 1 ||
        error("last row of affine matrix must be [0 0 0 1]")
    h.qform_code = one(Int16)
    h.sform_code = one(Int16)
    h.pixdim = (zero(Float32), h.pixdim[2:end]...,)
    h.quatern_b = zero(Float32)
    h.quatern_c = zero(Float32)
    h.quatern_d = zero(Float32)
    h.qoffset_x = zero(Float32)
    h.qoffset_y = zero(Float32)
    h.qoffset_z = zero(Float32)
    h.srow_x = convert(Tuple{Vararg{Float32}}, (affine[1, :]...,))
    h.srow_y = convert(Tuple{Vararg{Float32}}, (affine[2, :]...,))
    h.srow_z = convert(Tuple{Vararg{Float32}}, (affine[3, :]...,))
    h
end


text:
using Distributions: Uniform, Normal, MultivariateNormal
using LazySets: DefaultUniform

for N in [Float64]
    P1 = BallInf([0.6, 0.1, -1.3, -0.4, 0.2], 0.6)
    A = [1.0  0.0;
         0.0  1.0;
        -1.0  0.0;
         0.0 -1.0]
    b = [1.0, 2.0, 3.0, 4.0]
    P2 = HPolyhedron(A, b)
    P3 = HPolyhedron([0.0 0.1], [3.0])
    P4 = Ball2([0.2, -0.3, -1.1, 0.6, -0.7], 0.4)

    # Test rand samples are contained in the set
    p1 = LazySets.sample(P1)
    @test p1 ∈ P1
    p1_samples = LazySets.sample(P1, 100)
    @test sum(p1_samples .∈ Ref(P1)) == length(p1_samples)

    # default distribution
    @test LazySets.RejectionSampler(P2).distribution == [DefaultUniform(-3.0,1.0), DefaultUniform(-4.0,2.0)]

    # polytope sampler
    p1_samples = LazySets.sample(P1, 100, sampler=LazySets.RandomWalkSampler())
    @test all(v ∈ P1 for v in p1_samples)

    # specifying a distribution from Distributions.jl
    @test LazySets.RejectionSampler(P2, Uniform).distribution == [Uniform(-3.0, 1.0), Uniform(-4.0,2.0)]
    @test LazySets.RejectionSampler(P2, Normal).distribution == [Normal(-3.0, 1.0), Normal(-4.0, 2.0)]

    # test univariate distributions
    X = Interval(1, 2.)
    v = sample(X, 5)
    @test all(vi ∈ X for vi in v)
    d = Normal(1.5, 0.5)

    v = sample(X, 5, sampler=LazySets.RejectionSampler(d))
    @test all(vi ∈ X for vi in v)

    v = sample(X, 5, sampler=LazySets.RejectionSampler([d]))
    @test all(vi ∈ X for vi in v)

    v = sample(X, 5, sampler=LazySets.RejectionSampler([DefaultUniform(1, 2)]))
    @test all(vi ∈ X for vi in v)

    v = sample(X, 5, sampler=LazySets.RejectionSampler(DefaultUniform(1, 2)))
    @test all(vi ∈ X for vi in v)

    # test multivariate distribution
    H = rand(Hyperrectangle, dim=2)
    sampler = LazySets.RejectionSampler(MultivariateNormal(N[1 0; 0 1.]), tight=true)
    v = sample(H, 5, sampler=sampler)

    # including vertices
    for k in 0:4
        p1 = sample(P1, 10; include_vertices=k)
        @test length(p1) == 10 + k
    end
    @test length(sample(P1, 10; include_vertices=false)) == 10
    @test length(sample(P1, 10; include_vertices=true)) == 42
end


text:
### A Pluto.jl notebook ###
# v0.14.8

using Markdown
using InteractiveUtils

# ╔═╡ 5168c034-1879-11eb-2df4-e74ba7880da1
using Pkg, DrWatson

# ╔═╡ 5168ed8e-1879-11eb-31cd-d3a727fff08d
begin
	@quickactivate "SR2TuringPluto"
	using Turing
	using StatisticalRethinking
	Turing.setprogress!(false);
end;

# ╔═╡ e6220984-1878-11eb-275c-db7fe829037a
md"## Intro clip 2_chains.jl"

# ╔═╡ 51695ecc-1879-11eb-2fd1-b1d0abfb4b8f
begin
	N = 200
	df = DataFrame()
	df.weight = rand(Uniform(20, 90), N)
	f(x) = mean(df.weight) + 1.6x + rand(Normal(0, 10))
	df.height = f.(df.weight)
	Text(precis(df; io=String))
end

# ╔═╡ 5174d414-1879-11eb-28ed-b3f36412a25f
@model function ppl4_3(weights, heights)
    a ~ Normal(178, 20)
    b ~ LogNormal(0, 1)
    σ ~ Uniform(0, 50)
    μ = a .+ b .* weights
    for i in eachindex(heights)
        heights[i] ~ Normal(μ[i], σ)
    end
end

# ╔═╡ 517544f8-1879-11eb-3559-81508191d427
begin
	m4_3t = ppl4_3(df.weight, df.height)
    sampler=NUTS(0.65)
	nsamples=2000
	nchains=4
    chns4_3t = mapreduce(c -> sample(m4_3t, sampler, nsamples), chainscat, 1:nchains)
	Text(sprint(show, "text/plain", chns4_3t))
end

# ╔═╡ 2c0b2e42-187b-11eb-0c00-43aeac8b5045
plot(chns4_3t; seriestype=:traceplot)

# ╔═╡ 9110be74-187b-11eb-34e8-c1d4da25b29e
plot(chns4_3t; seriestype=:density)

# ╔═╡ b1690f50-19fa-11eb-15c7-87b89cf8ffa5
md"##### For additional diagnostics, see the [MCMCChains README document](https://github.com/TuringLang/MCMCChains.jl/blob/master/README.md)"

# ╔═╡ 0365d298-187f-11eb-1dc1-b7f61a13c7b3
md"## End of Intro clip 2_chains.jl"

# ╔═╡ Cell order:
# ╟─e6220984-1878-11eb-275c-db7fe829037a
# ╠═5168c034-1879-11eb-2df4-e74ba7880da1
# ╠═5168ed8e-1879-11eb-31cd-d3a727fff08d
# ╠═51695ecc-1879-11eb-2fd1-b1d0abfb4b8f
# ╠═5174d414-1879-11eb-28ed-b3f36412a25f
# ╠═517544f8-1879-11eb-3559-81508191d427
# ╠═2c0b2e42-187b-11eb-0c00-43aeac8b5045
# ╠═9110be74-187b-11eb-34e8-c1d4da25b29e
# ╟─b1690f50-19fa-11eb-15c7-87b89cf8ffa5
# ╟─0365d298-187f-11eb-1dc1-b7f61a13c7b3


text:
module TestUtils

using Test
using SyncBarriers.Internal: ceillog2

function test_ceillog2()
    xs = 1:2^10
    @test ceillog2.(xs) == ceil.(Int, log2.(xs))
end

end  # module


text:
"""

```
mddf_linkedcells_parallel(trajectory::Trajectory, options::Options,samples::Samples, RNG, mddf_compute!)
```

Computes the MDDF using linked cells, in parallel.

"""
function mddf_linkedcells_parallel(
    trajectory::Trajectory,
    options::Options,
    samples::Samples,
    RNG,
    mddf_compute!,
)

    # Simplify code by assigning some shortened names
    solute = trajectory.solute
    solvent = trajectory.solvent
    x_solute = trajectory.x_solute
    x_solvent = trajectory.x_solvent

    # Number of threads
    if options.nthreads < 0
        nspawn = Threads.nthreads() - 1
    else
        nspawn = options.nthreads - 1
    end
    if nspawn == 0
        error(" Parallel version must be executed only with more than 1 thread. ")
    end

    # Initializing the structure that carries the result per thread
    R0 = Result(trajectory, options)
    R = [Result(trajectory, options, irefatom = R0.irefatom) for i = 1:nspawn]

    # Safe passing of frame counter to the threads
    tframe = zeros(Int, nspawn)

    # Create data structures required for multithreading
    framedata = Vector{FrameData}(undef, nspawn)
    for ispawn = 1:nspawn
        framedata[ispawn] = FrameData(deepcopy(trajectory), R[1])
    end

    # Tasks
    t = Vector{Task}(undef, nspawn)
    free = ones(Bool, nspawn)

    # Skip initial frames if desired
    iframe = 0 # Counter for all frames of the input file
    while iframe < options.firstframe - 1
        nextframe!(trajectory)
        iframe += 1
    end

    # Print some information about this run
    options.silent || title(R[1], solute, solvent, nspawn)

    ndone = 0
    iframe_read = 0 # Counter for the frames that are actually being considered
    if !options.silent
        progress = Progress(R[1].nframes_read, 1)
    end
    while ndone < R[1].nframes_read

        # Launch for each free thread the computation of a frame
        while iframe_read < R[1].nframes_read && count(free) > 0
            ifree = findfirst(x -> x == true, free)
            # Skip frames that are not going to be considered
            while (iframe + 1) % options.stride != 0
                nextframe!(trajectory)
                iframe += 1
            end
            # Actually read the data of the frame that will be considered
            nextframe!(trajectory)
            iframe += 1
            iframe_read += 1
            tframe[ifree] = iframe
            @. framedata[ifree].trajectory.x_solute = trajectory.x_solute
            @. framedata[ifree].trajectory.x_solvent = trajectory.x_solvent
            @. framedata[ifree].trajectory.sides = trajectory.sides
            # Spawn the calculations for this frame
            t[ifree] = ThreadPools.@tspawnat ifree + 1 mddf_compute!(
                tframe[ifree],
                framedata[ifree],
                options,
                RNG,
                R[ifree],
            )
            free[ifree] = false
        end

        # Wait a little bit before checking
        sleep(options.sleep)

        # Check thread status
        for ispawn = 1:nspawn
            if !free[ispawn]
                if istaskfailed(t[ispawn])
                    error(
                        " Computation of MDDF failed in thread: $ispawn",
                        fetch(t[ispawn]),
                    )
                end
                if istaskdone(t[ispawn])
                    ndone += 1
                    free[ispawn] = true
                    options.silent || next!(progress)
                    if options.GC &&
                       (Sys.free_memory() / Sys.total_memory() < options.GC_threshold)
                        GC.gc() # why we need this anyway??? There should not be so much garbage.
                    end
                end
            end
        end

    end
    closetraj(trajectory)

    # Sum up the results of all threads into the data of thread one (R1<-R1+R2)
    for ispawn = 2:nspawn
        sum!(R[1], R[ispawn])
    end

    # Setup the final data structure with final values averaged over the number of frames,
    # sampling, etc, and computes final distributions and integrals
    finalresults!(R[1], options, trajectory, samples)
    options.silent || println(bars)

    return R[1]

end


text:
#   Unit "plots" of the FourierAnalysis Package for julia language
#
#   MIT License
#   Copyright (c) 2019-2021,
#   Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This unit implements recipes for plotting data generated
#   by FourierAnalysis without having any dependency to plot
#   packages

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

@recipe function f(::Type{Taper}, taper::Taper)
    title               --> taperinfo(taper)
    legend              --> false
    dpi                 --> 300
    background_color    --> :white
    gridcolor           --> :white
    taper.y
end

@recipe function f(::Type{Spectra}, S::Spectra;
            fmax   = S.sr/2,
            xspace = 4,
            ytitle = "Power (\\muV²)")

    b=f2b(fmax, S.sr, S.wl)
    xticstring=["$(xspace*i)" for i=1:fmax÷xspace]
    pushfirst!(xticstring, "$(fres(S.sr, S.wl))")
    xticrange=[1:S.wl/S.sr*xspace:b+1;]
    for i=2:length(xticrange) xticrange[i]=xticrange[i]-1 end

    legend         --> false
    dpi            --> 300
    xlabel         --> "Frequency (Hz)"
    xticks         --> (xticrange, xticstring)
    ylabel         --> ytitle
    delete!(plotattributes, :fmax)
    delete!(plotattributes, :xspace)
    delete!(plotattributes, :ytitle)
    S.y[1:(f2b(fmax, S.sr, S.wl; DC=S.DC)), :]
end

"""
Generate labels for the axes to be used in heatmap plots of [TFobjects](@ref).
On the y-axis the labels correspond to the center frequencies of the band-pass
filter bank used to obtain the time-frequency representation, that is,
to the `.flabels` field of the object.

See [plot time-frequency objects](@ref) for examples.
"""
function tfAxes(Y::TFobjects)
    function formatf(r::Real)
        s="$r"
        ss=split(s, ".")
        length(ss)>1 && parse(Float64, ss[2])==0.0 ? (return ss[1]) : (return s)
    end
    xs = [string(i) for i = 1:size(Y.y, 2)]
    ys = Vector{String}([formatf(f) for f in Y.flabels])
    return xs, ys
end


text:
g = CompleteDiGraph(5)
@test nv(g) == 5 && ne(g) == 20
g = CompleteGraph(5)
@test nv(g) == 5 && ne(g) == 10

g = CompleteBipartiteGraph(5, 8)
@test nv(g) == 13 && ne(g) == 40

g = StarDiGraph(5)
@test nv(g) == 5 && ne(g) == 4
g = StarGraph(5)
@test nv(g) == 5 && ne(g) == 4
g = StarGraph(1)
@test nv(g) == 1 && ne(g) == 0

g = PathDiGraph(5)
@test nv(g) == 5 && ne(g) == 4
g = PathGraph(5)
@test nv(g) == 5 && ne(g) == 4

g = CycleDiGraph(5)
@test nv(g) == 5 && ne(g) == 5
g = CycleGraph(5)
@test nv(g) == 5 && ne(g) == 5

g = WheelDiGraph(5)
@test nv(g) == 5 && ne(g) == 8
g = WheelGraph(5)
@test nv(g) == 5 && ne(g) == 8

g = Grid([3,3,4])
@test nv(g) == 3*3*4
@test ne(g) == 75
@test maximum(degree(g)) == 6
@test minimum(degree(g)) == 3

g = CliqueGraph(3,5)
@test nv(g) == 15 && ne(g) == 20
@test induced_subgraph(g,[1 2 3]) == CompleteGraph(3)

g = crosspath(3, BinaryTree(2))
# f = Vector{Vector{Int}}[[2 3 4];
# [1 5];
# [1 6];
# [1 5 6 7];
# [2 4 8];
# [3 4 9];
# [4 8 9];
# [5 7];
# [6 7]
# ]
I = [1,1,1,2,2,3,3,4,4,4,4,5,5,5,6,6,6,7,7,7,8,8,9,9]
J = [2,3,4,1,5,1,6,1,5,6,7,2,4,8,3,4,9,4,8,9,5,7,6,7]
V = ones(Int, length(I))
Adj = sparse(I,J,V)
@test Adj == sparse(g)

g = DoubleBinaryTree(3)
# [[3,2,8]
# [4,1,5]
# [1,6,7]
# [2]
# [2]
# [3]
# [3]
# [10,9,1]
# [11,8,12]
# [8,13,14]
# [9]
# [9]
# [10]
# [10]]
I = [1,1,1,2,2,2,3,3,3,4,5,6,7,8,8,8,9,9,9,10,10,10,11,12,13,14]
J = [3,2,8,4,1,5,1,6,7,2,2,3,3,10,9,1,11,8,12,8,13,14,9,9,10,10]
V = ones(Int, length(I))
Adj = sparse(I,J,V)
@test Adj == sparse(g)

rg3 = RoachGraph(3)
# [3]
# [4]
# [1,5]
# [2,6]
# [3,7]
# [4,8]
# [9,8,5]
# [10,7,6]
# [11,10,7]
# [8,9,12]
# [9,12]
# [10,11]
I = [ 1,2,3,3,4,4,5,5,6,6,7,7,7,8,8,8,9,9,9,10,10,10,11,11,12,12 ]
J = [ 3,4,1,5,2,6,3,7,4,8,9,8,5,10,7,6,11,10,7,8,9,12,9,12,10,11 ]
V = ones(Int, length(I))
Adj = sparse(I,J,V)
@test Adj == sparse(rg3)


text:
module Calc

using ..Mattermost

reduce_initialized = false

function calc(req)
    startswith(req.text, "!calc --wolfram ") && return req.text |> wcalc_lazy
    startswith(req.text, "!calc --reduce ")  && return req.text |> rcalc_lazy
    startswith(req.text, "!calc --maxima ")  && return req.text |> mcalc_lazy
    req.text |> wcalc_lazy
end

function wcalc_lazy(question)
    @eval begin
        using MathLink
        wcalc($question)
    end
end

function wcalc(question)
    question = replace(question, "!calc --wolfram " => "")
    question = replace(question, "!calc " => "")
    question = strip(question, ['`', ' '])

    question = question * " // TeXForm // ToString"

    try "Remove[\"Global`*\"]" |> MathLink.parseexpr |> weval catch; end

    local answer
    try
        answer = question |> MathLink.parseexpr |> weval |> string
    catch ex
        if isa(ex, MathLink.MathLinkError)
            return OutgoingWebhookResponse("```\n" * string(typeof(ex)) * "\n" * ex.msg * "```")
        else
            rethrow
        end
    end

    if occursin('^', answer) || occursin('_', answer) || occursin('\\', answer) || occursin('\n', answer)
        return OutgoingWebhookResponse("```latex\n" * answer * "\n```")
    else
        return OutgoingWebhookResponse("`" * answer * "`")
    end
end

function rcalc_lazy(question)
    @eval begin
        using Reduce
        rcalc($question)
    end
end

function rcalc(question)
    question = replace(question, "!calc --reduce " => "")
    question = replace(question, "!calc " => "")
    question = strip(question, ['`', ' '])

    if occursin(r"(in|out)\s+"i, question)
        throw(ReduceError("Forbidden"))
    end

    global reduce_initialized
    if !reduce_initialized
        rcall("load_package RESET")
        rcall("load_package RLFI")
        reduce_initialized = true
    end

#   try rcall("RESETREDUCE") catch; end

    local answer
    try
        answer = rcall(String(question), :latex) |> string
    catch ex
        if isa(ex, ReduceError)
            return OutgoingWebhookResponse("```\n" * string(typeof(ex)) * ex.errstr * "```")
        else
            rethrow
        end
    end

    answer = replace(answer, "\\begin{displaymath}" => "")
    answer = replace(answer, "\\end{displaymath}" => "")
    answer = strip(answer)

    if occursin('^', answer) || occursin('_', answer) || occursin('\\', answer) || occursin('\n', answer)
        return OutgoingWebhookResponse("```latex\n" * answer * "\n```")
    else
        return OutgoingWebhookResponse("`" * answer * "`")
    end
end

function mcalc_lazy(question)
    @eval begin
        using Maxima
        mcalc($question)
    end
end

function mcalc(question)
    question = replace(question, "!calc --maxima " => "")
    question = replace(question, "!calc " => "")
    question = strip(question, ['`', ' '])

    if occursin(r"(batch[a-z_]*|file[a-z_]*|load[a-z_]+|pathname[a-z_]*|save|stringout|with_stdout|[a-z_]*file)\s*\("i, question)
        throw(MaximaError("Forbidden"))
    end

    try
        mcall("kill (all)")
        mcall("reset ()")
        mcall("display2d: false")
    catch
    end

    local answer
    try
        answer = mcall(String(question)) |> string
    catch ex
        if isa(ex, MaximaError) || isa(ex, MaximaSyntaxError)
            return OutgoingWebhookResponse("```\n" * string(typeof(ex)) * ex.errstr * "```")
        else
            rethrow
        end
    end

    if startswith(answer, "\$\$")
        answer = replace(answer, r"^[^$]*[$]{2}" => "")
        answer = replace(answer, r"[$]{2}[^$]*$" => "")
        return OutgoingWebhookResponse("```latex\n" * answer * "\n```")
    else
        return OutgoingWebhookResponse("`" * answer * "`")
    end
end

end


text:
"""
      linear_mapping_approximation(rn_nonlinear::T, rn_linear::T, binary_vars::Array{Int,1}=Int[], m_order::Int=0;
                                   combinatoric_ratelaw = true) where T <: ReactionSystem

Given a *nonlinear* [`ReactionSystem`](https://catalyst.sciml.ai/stable/api/catalyst_api/#ModelingToolkit.ReactionSystem)
and an equivalent *linear* `ReactionSystem`, perform the Linear Mapping Approximation (LMA)
and return the corresponding linear [`RawMomentEquations`](@ref) of the system as well as
a Dictionary of reaction parameter substitutions obtained using LMA that are used to generate
the moment equations. See the [LMA theory section](@ref linear_mapping_approximation) for more details.

Notes:
- `rn_nonlinear` and `rn_linear` must be *identical* in layout in order to be interpreted correctly, and the
  nonlinear reactions contained in `rn_nonlinear` must all be linearised in `rn_linear` with rate coefficients
  updated accordingly. Although this requires a lot of manual input, automating the linearisation further is
  difficult due to arbitrary choices that may be mane in constructing the reaction networks.
- `binary_vars` *must* be specified for conditional closures as an array of indices of all species
  (as in [`speciesmap`](@ref)) which molecule number is a Bernoulli variable. Note that `rn_nonlinear`
  and `rn_linear` may internally order the species differently: `binary_vars` must be consistent with
  the ordering in the *nonlinear* network.
- By default the moment equations will be generated up to the order determined by the degree of nonlinearity
  of the nonlinear system's reactions. However, if higher order moment information is required, the optional
  `m_order` argument may be provided to increase the expansion order manually.
- `combinatoric_ratelaw=true` uses binomials in calculating the propensity functions
  of a `ReactionSystem`, see the notes for [`ModelingToolkit.jumpratelaw`]
  (https://mtk.sciml.ai/stable/systems/ReactionSystem/#ModelingToolkit.jumpratelaw).
-  [`ReactionSystemMod`](@ref) is currently unsupported as by design it does not contain information of reaction
  substrates and hence requires an API change. To be updated in the future!
"""
function linear_mapping_approximation(rn_nonlinear::T, rn_linear::T, binary_vars::Array{Int,1}=Int[], m_order::Int=0;
                                      combinatoric_ratelaw = true) where T <: ReactionSystem

      # check that all necessary information is provided and that LMA is applicable on the given nonlinear system

      # TODO: add further checks that both rn_nonlinear and rn_linear have reactions ordered identically
      @assert !isempty(binary_vars) "LMA does not work if there are no binary species"
      submat = substoichmat(rn_nonlinear)'
      no_substrates = vcat(sum(submat, dims=2)...)
      no_binary = vcat(sum(submat[:, binary_vars], dims=2)...)
      nonlinear_rs_inds = findall(no_substrates .> 1)
      @assert all(x == 1 for x in no_binary[nonlinear_rs_inds]) "non-linear reactions must involve one type of binary species"
      @assert all(submat[nonlinear_rs_inds, binary_vars] .< 2) "cannot have more than 1 molecules of each binary species"
      @assert all(sum(substoichmat(rn_linear), dims=1) .<= 1) "the linear network cannot contain nonlinear reactions"

      linearised_rs = reactions(rn_linear)[nonlinear_rs_inds]
      @assert all(!reaction.only_use_rate for reaction in linearised_rs) "linearised nonlinear reactions must follow the law of mass action (defining reactions using Catalyst's →)"
      coeffs = (reaction.rate for reaction in linearised_rs)

      # apply LMA: substitute reaction parameters in the linear network moment equations to reflect
      # the nonlinear interactions according to the LMA methodology

      term_factors, term_powers, poly_order = polynomial_propensities(propensities(rn_nonlinear; combinatoric_ratelaw)[nonlinear_rs_inds],
                                                                      get_iv(rn_nonlinear), speciesmap(rn_nonlinear))

      order = iszero(m_order) ? poly_order : max(poly_order, m_order)
      try sys = generate_raw_moment_eqs(rn_linear, order; combinatoric_ratelaw, smap=speciesmap(rn_nonlinear))
      catch e
            error("LMA cannot handle reactions with non-polynomial rates\n $e")
      end
      sys = bernoulli_moment_eqs(sys, binary_vars)

      which_binary = [binary_vars[findfirst(submat[i, binary_vars] .> 0)] for i in nonlinear_rs_inds]
      sub_params = OrderedDict()
      μ = sys.μ
      for (coeff, factors, powers, binary_ind) in zip(coeffs, term_factors, term_powers, which_binary)
            sub_params[coeff] = sum(factor*μ[Tuple(power)] for (factor, power) in zip(factors, powers))
            #sub_params[coeff] /= μ[sys.iter_1[binary_ind]]
            sub_params[coeff] *= μ[sys.iter_1[binary_ind]]^-1
            sub_params[coeff] = simplify(sub_params[coeff])
      end

      LMA_eqs = Equation[]
      for eq in get_eqs(sys.odes)
            rhs = substitute(eq.rhs, sub_params)
            rhs = expand(rhs)
            push!(LMA_eqs, Equation(eq.lhs, rhs))
      end

      field_values = [getfield(sys, field) for field in fieldnames(typeof(sys))]

      iv = get_iv(sys.odes)
      ps = reactionparams(rn_nonlinear)
      vars = states(sys.odes)
      odename = Symbol(nameof(sys), "_LMA")
      odes = ODESystem(LMA_eqs, iv, vars, ps; name=odename)
      new_system = typeof(sys)(odes, field_values[2:end]...)

      new_system, sub_params

end

#= TODO: consider adding
  • steady-state stuff
  • identification of whether a linear system has a time-dependendent or steady-state solution
  • construction of probability distributions if a solution exists
=#


text:
""" [Stevens2005](@cite) """
Dycoms_RF01_q_tot(::Type{FT}) where {FT} = z -> if z <= 840.0
        FT(9.0) / FT(1000.0)
    else
        FT(1.5) / FT(1000.0)
    end

""" [Stevens2005](@cite) """
Dycoms_RF01_θ_liq_ice(::Type{FT}) where {FT} = z -> if z <= 840.0
        FT(289)
    else
        FT(297.5) + (z - FT(840))^FT(1.0 / 3.0)
    end

""" [Stevens2005](@cite) """
Dycoms_RF01_u0(::Type{FT}) where {FT} = z -> FT(7)

""" [Stevens2005](@cite) """
Dycoms_RF01_v0(::Type{FT}) where {FT} = z -> FT(-5.5)

""" [Stevens2005](@cite) """
Dycoms_RF01_tke(::Type{FT}) where {FT} = z -> if z <= 800.0
        FT(1) - z / FT(1000)
    else
        FT(0)
    end


text:
@testset "string" begin

  @test string(Double64(pi)) == "3.1415926535897932384626433832795058"
  @test string(Double32(pi)) == "3.141592653589796"
  @test string(Double16(pi)) == "3.1415925"

  @test string(ComplexD64(1.0,1.0)) == "1.0 + 1.0im"
  @test string(ComplexD32(1.0,1.0)) == "1.0 + 1.0im"
  @test string(ComplexD16(1.0,1.0)) == "1.0 + 1.0im"

  @test stringtyped(Double64(pi)) == "Double64(3.141592653589793, 1.2246467991473532e-16)"
  @test stringtyped(Double32(pi)) == "Double32(3.1415927, -8.742278e-8)"
  @test stringtyped(Double16(pi)) == "Double16(3.14, 0.0009675)"

  @test stringtyped(ComplexD64(1.0,1.0)) == "ComplexD64(Double64(1.0, 0.0), Double64(1.0, 0.0))"
  @test stringtyped(ComplexD32(1.0,1.0)) == "ComplexD32(Double32(1.0, 0.0), Double32(1.0, 0.0))"
  @test stringtyped(ComplexD16(1.0,1.0)) == "ComplexD16(Double16(1.0, 0.0), Double16(1.0, 0.0))"

end


text:
using Distributed

mutable struct DataFrameWriter{C,T} <: Writer{C,T}
    results::DataFrame
    collector::C # The collector used by the Writer.
    localbuffer::TypedBuffer{T} # Stores data entries until the buffer is full.
    remotequeue::RemoteChannel{Channel{TypedBuffer{T}}} # References the channel on the master
                                                        # to which the buffer should be flushed.
end

"""
    DataFrameWriter(collector; chunksz = 15000, channels_per_worker = 2)

Return a `Writer` that returns the `collector`'s data as a DataFrame.

# Arguments
* `collector::Collector`: the Collector to use to collect data.
* `chunksz::Integer`: specifies the number of sets of data that the collector's buffer can store
    before sending it off to the master process to write it to disk.
* `channels_per_worker::Integer`: specifies the number of channels that the master process should
    allocate for each worker.
"""
function DataFrameWriter(collector::Collector{C};
                         chunksz::Integer = 15000, channels_per_worker::Integer = 2) where C
    localbuffer = TypedBuffer(collector, chunksz)
    buffertype = typeof(localbuffer)
    remotequeue = RemoteChannel(() -> Channel{buffertype}(channels_per_worker * nworkers()), 1)

    df_columns = [Array{coltype}(undef, 0) for coltype in C.types]
    df_names = [Symbol(colname) for colname in collector.names]
    results = DataFrame(df_columns, df_names)

    return DataFrameWriter(results, collector, localbuffer, remotequeue)
end

function open_sink(writer::DataFrameWriter)
    return nothing
end

function close_sink(sink::Nothing, writer::DataFrameWriter)
    return nothing
end

function get_return(writer::DataFrameWriter)
    return writer.results
end

function insert!(sink::Nothing, writer::DataFrameWriter, buffer::TypedBuffer{T}) where T
    n_datavectors = length(T.parameters)

    # Hoist from the loop.
    results = writer.results
    data = buffer.data
    size = buffer.size
    for datavector_idx in 1:n_datavectors
        @inbounds append!(results[datavector_idx], view(data[datavector_idx], 1:size))
    end

    return nothing
end


text:
import Base.Meta: quot

export @resource, @mixin

macro codegen(ex)
    :( push!(result.args, $(quot(ex))) )
end

macro resource(declaration, content)
    result = Expr(:block)
    if isa(declaration, Expr) && declaration.head == :<:
        this, super = declaration.args
        @codegen $this = Restful.Resource($(string(this)))
        @codegen Restful.addsubresource!($super, $this)
    elseif isa(declaration, Symbol)
        this = declaration
        @codegen $this = Restful.Resource($(string(this)))
    else
        error("unexpected $(declaration)")
    end

    definations = if content.head == :let
        content.args[1].args
    elseif content.head == :block
        content.args
    else
        error("unrecognized content type")
    end

    local description = ""
    for i in definations
        if isa(i, AbstractString) && isempty(description)
            description = i
        elseif i.head == :line
            continue
        elseif i.head == :(=>)
            key = i.args[1].args[1]
            if key == :name
                @codegen $(this).name = $(i.args[2])
            elseif key == :route
                @codegen $(this).route = $(i.args[2])
            elseif key == :mixin
                @codegen Restful.addmixin!($this, $(i.args[2]))
            elseif key in (:children, :subresources)
                @codegen Restful.addsubresource!($this, $(i.args[2]))
            elseif key in setdiff(HOOKS, METHODS)
                @codegen Restful.hook!($this, $(quot(key)), $(i.args[2]))
            else
                (hooks, method) = parsepipe(i.args[1])
                @codegen Restful.hook!($this, $method, $(Expr(:ref, :Function, hooks...)))
                if description==""
                    @codegen Restful.addmethod!($this, $method) do req, id $(i.args[2]) end
                else
                    @codegen Restful.addmethod!($this, $method, $description) do req, id $(i.args[2]) end
                    description = ""
                end
            end
        else
            error("unexpected $(i.head)")
        end
    end
    esc(result)
end

macro mixin(this, content)
    result = Expr(:block)
    @codegen $this = Restful.Mixin()
    for i in content.args
        if i.head == :line
            continue
        elseif i.head == :(=>)
            key = i.args[1].args[1]
            if key in HOOKS
                @codegen Restful.hook!($this, $(quot(key)), $(i.args[2]))
            else
                error("no hooks called $key")
            end
        else
            error("unexpected $(i.head)")
        end
    end
    esc(result)
end

function parsepipe(x::Expr)
    if x.head == :call && x.args[1] == :|
        (hooks, method) = parsepipe(x.args[2])
        (push!(hooks, x.args[3]), method)
    elseif x.head == :quote
        (Union{Symbol, Expr}[], x)
    else
        error("unexpected $x")
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiExtensionsV1beta1HTTPIngressRuleValue <: SwaggerModel
    paths::Any # spec type: Union{ Nothing, Vector{IoK8sApiExtensionsV1beta1HTTPIngressPath} } # spec name: paths

    function IoK8sApiExtensionsV1beta1HTTPIngressRuleValue(;paths=nothing)
        o = new()
        validate_property(IoK8sApiExtensionsV1beta1HTTPIngressRuleValue, Symbol("paths"), paths)
        setfield!(o, Symbol("paths"), paths)
        o
    end
end # type IoK8sApiExtensionsV1beta1HTTPIngressRuleValue

const _property_map_IoK8sApiExtensionsV1beta1HTTPIngressRuleValue = Dict{Symbol,Symbol}(Symbol("paths")=>Symbol("paths"))
const _property_types_IoK8sApiExtensionsV1beta1HTTPIngressRuleValue = Dict{Symbol,String}(Symbol("paths")=>"Vector{IoK8sApiExtensionsV1beta1HTTPIngressPath}")
Base.propertynames(::Type{ IoK8sApiExtensionsV1beta1HTTPIngressRuleValue }) = collect(keys(_property_map_IoK8sApiExtensionsV1beta1HTTPIngressRuleValue))
Swagger.property_type(::Type{ IoK8sApiExtensionsV1beta1HTTPIngressRuleValue }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiExtensionsV1beta1HTTPIngressRuleValue[name]))}
Swagger.field_name(::Type{ IoK8sApiExtensionsV1beta1HTTPIngressRuleValue }, property_name::Symbol) =  _property_map_IoK8sApiExtensionsV1beta1HTTPIngressRuleValue[property_name]

function check_required(o::IoK8sApiExtensionsV1beta1HTTPIngressRuleValue)
    (getproperty(o, Symbol("paths")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiExtensionsV1beta1HTTPIngressRuleValue }, name::Symbol, val)
end


text:
# ---
# title: 1309. Decrypt String from Alphabet to Integer Mapping
# id: problem1309
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/decrypt-string-from-alphabet-to-integer-mapping/description/>
# hidden: true
# ---
#
# Given a string `s` formed by digits (`'0'` \- `'9'`) and `'#'` . We want to
# map `s` to English lowercase characters as follows:
#
#   * Characters (`'a'` to `'i')` are represented by (`'1'` to `'9'`) respectively.
#   * Characters (`'j'` to `'z')` are represented by (`'10#'` to `'26#'`) respectively.
#
# Return the string formed after mapping.
#
# It's guaranteed that a unique mapping will always exist.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "10#11#12"
#     Output: "jkab"
#     Explanation: "j" -> "10#" , "k" -> "11#" , "a" -> "1" , "b" -> "2".
#
#
# **Example 2:**
#
#
#
#     Input: s = "1326#"
#     Output: "acz"
#
#
# **Example 3:**
#
#
#
#     Input: s = "25#"
#     Output: "y"
#
#
# **Example 4:**
#
#
#
#     Input: s = "12345678910#11#12#13#14#15#16#17#18#19#20#21#22#23#24#25#26#"
#     Output: "abcdefghijklmnopqrstuvwxyz"
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 1000`
#   * `s[i]` only contains digits letters (`'0'`-`'9'`) and `'#'` letter.
#   * `s` will be valid string such that mapping is always possible.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end