text:
using Test
using ZXCalculus
using YaoHIR, YaoLocations
using YaoHIR.IntrinsicOperation
using CompilerPluginTools

chain = Chain()
push_gate!(chain, Val(:Sdag), 1)
push_gate!(chain, Val(:Z), 1)
push_gate!(chain, Val(:H), 1)
push_gate!(chain, Val(:S), 1)
push_gate!(chain, Val(:S), 2)
push_gate!(chain, Val(:H), 4)
push_gate!(chain, Val(:CNOT), 3, 2)
push_gate!(chain, Val(:CZ), 4, 1)
push_gate!(chain, Val(:H), 2)
push_gate!(chain, Val(:T), 2)
push_gate!(chain, Val(:CNOT), 3, 2)
push_gate!(chain, Val(:Tdag), 2)
push_gate!(chain, Val(:CNOT), 1, 4)
push_gate!(chain, Val(:H), 1)
push_gate!(chain, Val(:T), 2)
push_gate!(chain, Val(:S), 3)
push_gate!(chain, Val(:H), 4)
push_gate!(chain, Val(:T), 1)
push_gate!(chain, Val(:T), 2)
push_gate!(chain, Val(:T), 3)
push_gate!(chain, Val(:X), 3)
push_gate!(chain, Val(:H), 2)
push_gate!(chain, Val(:H), 3)
push_gate!(chain, Val(:Sdag), 4)
push_gate!(chain, Val(:S), 3)
push_gate!(chain, Val(:X), 4)
push_gate!(chain, Val(:CNOT), 3, 2)
push_gate!(chain, Val(:H), 1)
push_gate!(chain, Val(:shift), 4, ZXCalculus.Phase(1//2))
push_gate!(chain, Val(:Rx), 4, ZXCalculus.Phase(1//1))
push_gate!(chain, Val(:Rx), 3, ZXCalculus.Phase(1//4))
push_gate!(chain, Val(:Rx), 2, ZXCalculus.Phase(1//4))
push_gate!(chain, Val(:S), 3)

ir = @make_ircode begin
end
bir = BlockIR(ir, 4, chain)
zxd = convert_to_zxd(bir)
convert_to_chain(zxd)
pt_zxd = phase_teleportation(zxd)
@test tcount(pt_zxd) <= tcount(zxd)
pt_chain = convert_to_chain(pt_zxd)
@test length(pt_chain) <= length(chain)

zxg = clifford_simplification(zxd)
cl_chain = circuit_extraction(zxg)

zxg = full_reduction(zxd)
fl_chain = circuit_extraction(zxg)
ZXCalculus.generate_layout!(zxg)
@test ZXCalculus.qubit_loc(zxg, 40) == 0//1
ZXCalculus.spider_sequence(zxg)

pt_bir = phase_teleportation(bir)
cl_bir = clifford_simplification(bir)
fl_bir = full_reduction(bir)

@test length(pt_chain) == length(pt_bir.circuit)
@test length(cl_chain) == length(cl_bir.circuit)
@test length(fl_chain) == length(fl_bir.circuit)

@testset "issue#80" begin
    ir = @make_ircode begin
    end

    circuit = Chain(Gate(X, Locations(1)), Gate(X, Locations(1)))
    bir = BlockIR(ir, 1, circuit)
    bir = clifford_simplification(bir)
    bir = clifford_simplification(bir)
    @test bir.circuit == Chain(Gate(H, Locations((1, ))), Gate(H, Locations((1, ))))
end

@testset "generate_layout!" begin
    circ = Chain(Gate(H, Locations(2)), Gate(T, Locations(4)), Gate(H, Locations(1)), Gate(AdjointOperation{SGate}(S), Locations(2)), Gate(H, Locations(2)), Gate(X, Locations(3)), Gate(AdjointOperation{SGate}(S), Locations(1)), Gate(Z, Locations(1)), Gate(H, Locations(2)), Gate(X, Locations(1)), Gate(Z, Locations(1)), Gate(T, Locations(5)), Ctrl(Gate(X, Locations(5)), CtrlLocations(1)), Gate(H, Locations(1)), Gate(T, Locations(1)), Ctrl(Gate(X, Locations(3)), CtrlLocations(5)), Gate(H, Locations(1)), Gate(X, Locations(4)), Ctrl(Gate(X, Locations(5)), CtrlLocations(4)), Gate(S, Locations(2)), Ctrl(Gate(X, Locations(1)), CtrlLocations(3)), Gate(X, Locations(2)), Ctrl(Gate(X, Locations(3)), CtrlLocations(1)), Gate(X, Locations(2)), Gate(S, Locations(3)), Gate(Z, Locations(2)), Gate(Z, Locations(5)), Gate(X, Locations(2)), Gate(X, Locations(1)), Ctrl(Gate(X, Locations(3)), CtrlLocations(5)), Gate(S, Locations(4)), Gate(X, Locations(3)), Ctrl(Gate(X, Locations(1)), CtrlLocations(2)), Gate(AdjointOperation{SGate}(S), Locations(4)), Gate(Z, Locations(2)), Gate(AdjointOperation{SGate}(S), Locations(5)), Ctrl(Gate(X, Locations(1)), CtrlLocations(2)), Gate(Z, Locations(5)), Ctrl(Gate(X, Locations(1)), CtrlLocations(4)), Ctrl(Gate(X, Locations(3)), CtrlLocations(4)), Gate(H, Locations(4)), Gate(Z, Locations(1)), Gate(X, Locations(4)), Gate(Z, Locations(3)), Gate(H, Locations(4)), Ctrl(Gate(X, Locations(1)), CtrlLocations(4)), Ctrl(Gate(X, Locations(4)), CtrlLocations(1)), Gate(X, Locations(4)), Gate(S, Locations(3)), Gate(AdjointOperation{SGate}(S), Locations(2)), Gate(Z, Locations(3)), Gate(S, Locations(5)), Ctrl(Gate(X, Locations(3)), CtrlLocations(5)), Gate(H, Locations(2)), Gate(Z, Locations(4)), Gate(H, Locations(1)), Gate(X, Locations(1)), Gate(X, Locations(2)), Ctrl(Gate(X, Locations(5)), CtrlLocations(3)), Ctrl(Gate(X, Locations(1)), CtrlLocations(3)), Gate(Z, Locations(4)), Gate(S, Locations(5)), Gate(S, Locations(5)), Ctrl(Gate(X, Locations(5)), CtrlLocations(1)), Gate(T, Locations(4)), Gate(Z, Locations(2)), Gate(X, Locations(4)), Gate(H, Locations(2)), Gate(AdjointOperation{SGate}(S), Locations(3)), Gate(H, Locations(5)), Gate(T, Locations(2)), Gate(AdjointOperation{SGate}(S), Locations(5)), Gate(AdjointOperation{SGate}(S), Locations(4)), Gate(H, Locations(2)), Gate(S, Locations(2)), Gate(AdjointOperation{SGate}(S), Locations(3)), Ctrl(Gate(X, Locations(5)), CtrlLocations(1)), Gate(H, Locations(5)), Ctrl(Gate(X, Locations(4)), CtrlLocations(3)), Gate(H, Locations(4)), Gate(Z, Locations(2)), Gate(H, Locations(2)), Gate(S, Locations(2)), Gate(H, Locations(4)), Ctrl(Gate(X, Locations(1)), CtrlLocations(2)), Gate(X, Locations(4)), Gate(H, Locations(1)), Gate(AdjointOperation{SGate}(S), Locations(5)), Gate(X, Locations(1)), Gate(X, Locations(3)), Gate(Z, Locations(1)), Ctrl(Gate(X, Locations(1)), CtrlLocations(2)), Gate(T, Locations(3)), Gate(X, Locations(4)), Gate(X, Locations(3)), Gate(AdjointOperation{SGate}(S), Locations(1)), Gate(AdjointOperation{SGate}(S), Locations(3)), Gate(H, Locations(5)), Gate(S, Locations(4)), Gate(Z, Locations(4)))
    ir = @make_ircode begin end
    bir = BlockIR(ir, 5, circ)
    zxd = ZXDiagram(bir)
    zxg = ZXGraph(zxd)
    full_reduction(zxg)
    ZXCalculus.generate_layout!(zxg)
    @test ZXCalculus.generate_layout!(zxg) !== nothing
end

function random_circuit(nbits, ngates; T = 0.1, CZ = 0.0, CNOT = 0.1)
    ir = @make_ircode begin
    end
    CLIFF = 1 - T - CZ - CNOT
    circ = Chain()
    for _ = 1:ngates
        x = rand()
        nbits == 1 && (x = x*(CLIFF+T))
        if x <= CLIFF
            g = rand([:X, :X, :Z, :Z, :S, :Sdag, :H, :H])
            push_gate!(circ, Val(g), rand(1:nbits))
        elseif x - CLIFF <= T
            g = rand([:T, :Tdag])
            push_gate!(circ, Val(g), rand(1:nbits))
        elseif x - CLIFF - T <= CZ
            loc = rand(1:nbits)
            ctrl = loc
            while ctrl == loc
                ctrl = rand(1:nbits)
            end
            push_gate!(circ, Val(:CZ), loc, ctrl)
        else
            loc = rand(1:nbits)
            ctrl = loc
            while ctrl == loc
                ctrl = rand(1:nbits)
            end
            push_gate!(circ, Val(:CNOT), loc, ctrl)
        end
    end
    return BlockIR(ir, nbits, circ)
end


text:
import Test: AbstractTestSet, DefaultTestSet, record, finish, Result, Pass
import Actors: Actor

"Allows the Test macros to be safely used in any actor"
mutable struct LuvvyTestSet <: AbstractTestSet
    desc::String
    results::Vector{Result}
    expect::Int

    myself::Id{LuvvyTestSet}

    LuvvyTestSet(desc; expect=0) = new(desc, [], expect)
end

hear(s::Scene{LuvvyTestSet}, res::Result) = push!(my(s).results, res)

record(ts::LuvvyTestSet, res::Result) =
    say(task_local_storage(:scene), ts.myself, res)

finish(ts::LuvvyTestSet) = let dts = DefaultTestSet(ts.desc)
    wrong_length = ts.expect > 0 && ts.expect != length(ts.results)

    if any(res -> !(res isa Pass), ts.results) || wrong_length
        sleep(0.5)
    end

    foreach(res -> record(dts, res), ts.results)
    finish(dts)

    if wrong_length
        printstyled("Test Set Error: "; bold=true, color=Base.error_color())
        printstyled("Expected $(ts.expect) '$(ts.desc)' tests in total!\n";
                    color=Base.error_color())
    end
end

# Note that overriding the following methods is considered a last resort for
# integrating with another library.

struct TestEnvironment
    ts::LuvvyTestSet
end

# Get our test set from the main task's (The Stage's task) local storage
Actors.capture_environment(::Id) = let ts = Test.get_testset()
    ts isa LuvvyTestSet ? TestEnvironment(ts) : nothing
end

# Inject the test set actor into our play when the Stage starts
Actors.prologue!(s::Scene{Stage}, env::TestEnvironment) =
    env.ts.myself = enter!(s, env.ts)

# Inject the test set into a new actor's Task local storage
function Actors.prologue!(s::Actors.AbsScene{S}, env::TestEnvironment) where S
    task_local_storage(:scene, s)
    Test.push_testset(env.ts)
end

ext: jl


text:
export @pretty

"""
    @pretty @shape A[...] := B[...]

Prints an approximately equivalent expression with the macro expanded.
Compared to `@macroexpand1`, generated symbols are replaced with animal names (from MacroTools),
comments are deleted, module names are removed from functions,
and the final expression is fed to `println()`.

To copy and run the printed expression, you may need various functions which aren't exported.
Try something like `using TensorCast: orient, star, rview, @assert_, red_glue, sliceview, batchmul`
"""
macro pretty(ex)
    if @capture(ex, @cast_str str_)
        full = "@cast " * cast_string(str)
        println("# " * full)
        return :(@pretty $(Meta.parse(full)) )

    elseif @capture(ex, @reduce_str str_)
        full = reduce_string(str)
        println("# " * full)
        return :(@pretty $(Meta.parse(full)) )

    else
        :( macroexpand($(__module__), $(ex,)[1], recursive=false) |> pretty |> println )
    end
end

function pretty(ex::Union{Expr,Symbol})
    ex = MacroTools.alias_gensyms(ex)
    ex = MacroTools.striplines(ex) # remove line references

    str = string(ex)

    # str = replace(str, r"(#=.+=#)\s" => "")  # remove line references
    str = replace(str, r"\n(\s+\n)" => "\n") # and empty lines

    str = replace(str, r"\w+\.(\w+)\(" => s"\1(") # remove module names on functions
    str = replace(str, r"\w+\.(\w+)!\(" => s"\1!(")
    str = replace(str, r"\w+\.(@\w+)" => s"\1")   # and on macros
    str = replace(str, r"\w+\.(\w+){" => s"\1{")  # and on structs...

    str = replace(str, "Colon()" => ":")
end

pretty(tup::Tuple) = replace(string(tup), "Colon()" => ":")

pretty(vec::Vector) = join(something.(vec, "nothing"), ", ")

pretty(x) = string(x)


text:
using HDF5

file_id = h5open("test_unlimited.h5", "cw")

a = rand(2, 3, 4)

a_id = create_dataset(file_id, "a", Float64, ((size(a)..., 1), (size(a)..., -1)), chunk=(size(a)..., 1))
a_id[:, :, :, 1] .= a[:, :, :]

for i in 1:10
    a[:, :, :] .+= 1
    HDF5.set_extent_dims(a_id, (size(a)..., i))
    a_id[:, :, :, i] = a[:, :, :]
end

close(file_id)


text:
export PassiveADPAgentMDP, PassiveADPAgentProgram,
        PassiveTDAgentProgram, QLearningAgentProgram;

#=

	PassiveADPAgentMDP is a MDP implementation of AbstractMarkovDecisionProcess

    that consists of a MarkovDecisionProcess 'mdp'.

=#
struct PassiveADPAgentMDP{T} <: AbstractMarkovDecisionProcess
	mdp::MarkovDecisionProcess{T}

    function PassiveADPAgentMDP{T}(initial::T, actions_list::Set{T}, terminal_states::Set{T}, gamma::Float64) where T
		return new(MarkovDecisionProcess(initial, actions_list, terminal_states, Dict(), gamma=gamma));
	end
end

PassiveADPAgentMDP(initial, actions_list::Set, terminal_states::Set, gamma::Float64) = PassiveADPAgentMDP{typeof(initial)}(initial, actions_list, terminal_states, gamma);

"""
    reward(mdp::PassiveADPAgentMDP, state)

Return a reward based on the given 'state'.
"""
function reward(mdp::PassiveADPAgentMDP, state)
    return mdp.mdp.reward[state];
end

"""
    transition_model(mdp::PassiveADPAgentMDP, state, action)

Return a list of (P(s'|s, a), s') pairs given the state 's' and action 'a'.
"""
function transition_model(mdp::PassiveADPAgentMDP, state, action)
    return collect((v, k) for (k, v) in get!(mdp.mdp.transitions, (state, action), Dict()));
end

"""
    actions(mdp::PassiveADPAgentMDP, state)

Return a set of actions that are possible in the given state.
"""
function actions(mdp::PassiveADPAgentMDP, state)
    if (state in mdp.mdp.terminal_states)
        return Set{Nothing}([nothing]);
    else
        return mdp.mdp.actions;
    end
end

"""
    policy_evaluation(pi::Dict, U::Dict, mdp::PassiveADPAgentMDP; k::Int64=20)

Return the updated utilities of the MDP's states by applying the modified policy iteration
algorithm on the given Markov decision process 'mdp', utility function 'U', policy 'pi',
and number of Bellman updates to use 'k'.
"""
function policy_evaluation(pi::Dict, U::Dict, mdp::PassiveADPAgentMDP; k::Int64=20)
    for i in 1:k
        for state in mdp.mdp.states
            if (length(transition_model(mdp, state, pi[state])) != 0)
                U[state] = (reward(mdp, state)
                            + (mdp.mdp.gamma
                            * sum((p * U[state_prime] for (p, state_prime) in transition_model(mdp, state, pi[state])))));
            else
                U[state] = (reward(mdp, state) + (mdp.mdp.gamma * 0));
            end
        end
    end
    return U;
end

#=

    PassiveADPAgentProgram is a passive reinforcement learning agent based on

    adaptive dynamic programming (Fig. 21.2).

=#
mutable struct PassiveADPAgentProgram <: AgentProgram
    state   # can be any DataType, but check for Nothing DataType later
    action  # can be any DataType, but check for Nothing DataType later
    U::Dict
    pi::Dict
    mdp::PassiveADPAgentMDP
    N_sa::Dict
    N_s_prime_sa::Dict

    function PassiveADPAgentProgram(pi::Dict, mdp::T) where {T <: AbstractMarkovDecisionProcess}
        return new(nothing,
                    nothing,
                    Dict(),
                    pi,
                    PassiveADPAgentMDP(mdp.initial, mdp.actions, mdp.terminal_states, mdp.gamma),
                    Dict(),
                    Dict());
    end
end

function execute(padpap::PassiveADPAgentProgram, percept::Tuple{Any, Any})
    local r_prime::Float64;
    s_prime, r_prime = percept;

    push!(padpap.mdp.mdp.states, s_prime);
    if (!haskey(padpap.mdp.mdp.reward, s_prime))
        padpap.U[s_prime] = r_prime;
        padpap.mdp.mdp.reward[s_prime] = r_prime;
    end
    if (!(padpap.state === nothing))
        padpap.N_sa[(padpap.state, padpap.action)] = get!(padpap.N_sa, (padpap.state, padpap.action), 0) + 1;
        padpap.N_s_prime_sa[(s_prime, padpap.state, padpap.action)] = get!(padpap.N_s_prime_sa, (s_prime, padpap.state, padpap.action), 0) + 1;
        for t in collect(result_state
                        for ((result_state, state, action), occurrences) in padpap.N_s_prime_sa
                        if (((state, action) == (padpap.state, padpap.action)) && (occurrences != 0)))
            get!(padpap.mdp.mdp.transitions, (padpap.state, padpap.action), Dict())[t] = padpap.N_s_prime_sa[(t, padpap.state, padpap.action)] / padpap.N_sa[(padpap.state, padpap.action)];
        end
    end
    local U::Dict = policy_evaluation(padpap.pi, padpap.U, padpap.mdp);
    if (s_prime in padpap.mdp.mdp.terminal_states)
        padpap.state = nothing;
        padpap.action = nothing;
    else
        padpap.state = s_prime;
        padpap.action = padpap.pi[s_prime];
    end
    return padpap.action;
end

function update_state(padpap::PassiveADPAgentProgram, percept::Tuple{Any, Any})
    return percept;
end

#=

    PassiveTDAgentProgram is a passive reinforcement learning agent that learns

    utility estimates by using temporal differences (Fig. 21.4).

=#
mutable struct PassiveTDAgentProgram <: AgentProgram
    state   # can be any DataType, but check for Nothing DataType later
    action  # can be any DataType, but check for Nothing DataType later
    reward  # can be any DataType, but check for Nothing DataType later
    gamma::Float64
    U::Dict
    pi::Dict
    N_s::Dict
    terminal_states::Set
    alpha::Function

    function PassiveTDAgentProgram(pi::Dict, mdp::T; alpha::Union{Nothing, Function}=nothing) where {T <: AbstractMarkovDecisionProcess}
        local gamma::Float64;
        local terminal_states::Set;
        local new_alpha::Function;
        if (typeof(mdp) <: PassiveADPAgentMDP)
            gamma = mdp.mdp.gamma;
            terminal_states = mdp.mdp.terminal_states;
        else
            gamma = mdp.gamma;
            terminal_states = mdp.terminal_states;
        end
        if (typeof(alpha) <: Nothing)
            new_alpha = (function(n::Number)
                            return (1/(n + 1));
                        end);
        else
            new_alpha = alpha;
        end
        return new(nothing,
                    nothing,
                    nothing,
                    gamma,
                    Dict(),
                    pi,
                    Dict(),
                    terminal_states,
                    new_alpha);
    end
end

function execute(ptdap::PassiveTDAgentProgram, percept::Tuple{Any, Any})
    local r_prime::Float64;
    s_prime, r_prime = update_state(ptdap, percept);
    if (!haskey(ptdap.N_s, s_prime))
        ptdap.U[s_prime] = r_prime;
    end
    if (!(ptdap.state === nothing))
        ptdap.N_s[ptdap.state] = get!(ptdap.N_s, ptdap.state, 0) + 1;
        ptdap.U[ptdap.state] = (get!(ptdap.U, ptdap.state, 0.0)
                                    + ptdap.alpha(get!(ptdap.N_s, ptdap.state, 0))
                                    * (ptdap.reward
                                    + (ptdap.gamma * get!(ptdap.U, s_prime, 0.0))
                                    - get!(ptdap.U, ptdap.state, 0.0)));
    end
    if (s_prime in ptdap.terminal_states)
        ptdap.state = nothing;
        ptdap.action = nothing;
        ptdap.reward = nothing;
    else
        ptdap.state = s_prime;
        ptdap.action = ptdap.pi[s_prime];
        ptdap.reward = r_prime;
    end
    return ptdap.action;
end

function update_state(ptdap::PassiveTDAgentProgram, percept::Tuple{Any, Any})
    return percept;
end

#=

    QLearningAgentProgram is an exploratory Q-learning agent that learns the value

    Q(state, action) for each action in each situation (Fig. 21.8). The agent uses the

    same exploration function as the exploratory ADP agent, but avoid learning the

    transition model because the Q-value of a state can be related directly to those of

    its neighbor.

=#
mutable struct QLearningAgentProgram <: AgentProgram
    state   # can be any DataType, but check for Nothing DataType later
    action  # can be any DataType, but check for Nothing DataType later
    reward  # can be any DataType, but check for Nothing DataType later
    gamma::Float64
    Q::Dict
    N_sa::Dict
    actions::Set
    terminal_states::Set
    R_plus::Float64 # optimistic estimate of the best possible reward obtainable
    N_e::Int64  # try action-state pair at least N_e times
    f::Function
    alpha::Function

    function QLearningAgentProgram(mdp::T, N_e::Int64, R_plus::Number; alpha::Union{Nothing, Function}=nothing) where {T <: AbstractMarkovDecisionProcess}
        local new_alpha::Function;
        local gamma::Float64;
        local actions::Set;
        local terminal_states::Set;
        if (typeof(mdp) <: PassiveADPAgentMDP)
            gamma = mdp.mdp.gamma;
            actions = mdp.mdp.actions;
            terminal_states = mdp.mdp.terminal_states;
        else
            gamma = mdp.gamma;
            actions = mdp.actions;
            terminal_states = mdp.terminal_states;
        end
        if (typeof(alpha) <: Nothing)
            new_alpha = (function(n::Number)
                            return (1/(n + 1));
                        end);
        else
            new_alpha = alpha;
        end
        return new(nothing,
                    nothing,
                    nothing,
                    gamma,
                    Dict(),
                    Dict(),
                    actions,
                    terminal_states,
                    R_plus,
                    N_e,
                    exploration_function,
                    new_alpha);
    end
end

function exploration_function(qlap::QLearningAgentProgram, u::Number, n::Number)
    if (n < qlap.N_e)
        return qlap.R_plus;
    else
        return u;
    end
end

function actions(qlap::QLearningAgentProgram, state)
    if (state in qlap.terminal_states)
        return Set([nothing]);
    else
        return qlap.actions;
    end
end

function execute(qlap::QLearningAgentProgram, percept::Tuple{Any, Any})
    local r_prime::Float64;
    s_prime, r_prime = update_state(qlap, percept);
    if (!(qlap.state === nothing))
        if (qlap.state in qlap.terminal_states)
            qlap.Q[(qlap.state, nothing)] = r_prime;
        end
        qlap.N_sa[(qlap.state, qlap.action)] = get!(qlap.N_sa, (qlap.state, qlap.action), 0) + 1;
        # Default value for Q keys is 0.0.
        get!(qlap.Q, (qlap.state, qlap.action), 0.0);
        qlap.Q[(qlap.state, qlap.action)] = (qlap.Q[(qlap.state, qlap.action)]
                                                    + (qlap.alpha(qlap.N_sa[(qlap.state, qlap.action)]) *
                                                    (qlap.reward +
                                                    (qlap.gamma * reduce(max, collect(get!(qlap.Q, (s_prime, a_prime), 0.0)
                                                                                for a_prime in actions(qlap, s_prime))))
                                                    - qlap.Q[(qlap.state, qlap.action)])));
    end
    if (!(qlap.state === nothing) && (qlap.state in qlap.terminal_states))
        qlap.state = nothing;
        qlap.action = nothing;
        qlap.reward = nothing;
    else
        qlap.state = s_prime;
        qlap.action = argmax(collect(actions(qlap, s_prime)),
                            (function(a_prime)
                                return qlap.f(qlap, get!(qlap.Q, (s_prime, a_prime), 0.0), get!(qlap.N_sa, (s_prime, a_prime), 0));
                            end));
        qlap.reward = r_prime;
    end
    return qlap.action;
end

function update_state(qlap::QLearningAgentProgram, percept::Tuple{Any, Any})
    return percept;
end

"""
    take_single_action(mdp::T, state, action) where {T <: AbstractMarkovDecisionProcess}

Return the next state by choosing a weighted sample of the resulting states for
taking the action 'action' in state 'state'.
"""
function take_single_action(mdp::T, state, action) where {T <: AbstractMarkovDecisionProcess}
    local x::Float64 = rand(RandomDeviceInstance);
    local cumulative_probability::Float64 = 0.0;
    for (p, state_p) in transition_model(mdp, state, action)
        cumulative_probability = cumulative_probability + p;
        if (x < cumulative_probability)
            return state_p;
        end
    end
    error("take_single_action(): Could not find a valid resulting state for the state ", state,
        " and action ", action, "!");
end

"""
    run_single_trial(ap::T1, mdp::T2) where {T1 <: AgentProgram, T2 <: AbstractMarkovDecisionProcess}

The agent program 'ap' executes a trial in the environment represented by the MDP 'mdp'.
"""
function run_single_trial(ap::T1, mdp::T2) where {T1 <: AgentProgram, T2 <: AbstractMarkovDecisionProcess}
    current_state = mdp.initial;
    while (true)
        local current_reward::Float64 = reward(mdp, current_state);
        local percept::Tuple = (current_state, current_reward);
        next_action = execute(ap, percept);
        if (typeof(next_action) <: Nothing)
            break;
        end
        current_state = take_single_action(mdp, current_state, next_action);
    end
    return nothing;
end

max_stars_repo_name: mikhail-j/aima-julia

max_issues_repo_name: mikhail-j/aima-julia

max_forks_repo_name: mikhail-j/aima-julia


text:
module TestRVineMatrix

using Copulas
using Base.Test

#######################
## constructor tests ##
#######################

kk = [1 0 0; 2 3 0; 3 1 2]
@test_throws Exception Copulas.RVMatrix(kk)

kk = [2 0 0; 2 3 0; 3 2 2]
@test_throws Exception Copulas.RVMatrix(kk)

kk = [2 0 4; 2 3 0; 3 2 2]
@test_throws Exception Copulas.RVMatrix(kk)

kk = [1 0 0; 2 3 0; 3 2 2]
Copulas.RVMatrix(kk)

#########################
## columnSubsets tests ##
#########################

kk = [1 0 0; 2 3 0; 3 1 2]
@test_throws Exception Copulas.chkColumnSubsets(kk)

kk = [2 0 0; 2 3 0; 3 2 2]
@test_throws Exception Copulas.chkColumnSubsets(kk)

kk = [1 0 0; 2 3 0; 3 2 2]
@test Copulas.chkColumnSubsets(kk)

########
## dg ##
########

kk = [1 0 0; 2 3 0; 3 2 2]
rvm = Copulas.RVMatrix(kk)
@test Copulas.dg(rvm) == kk

###############
## allExcept ##
###############

vn = Copulas.testvine(2)
paths = convert(Copulas.CTreePaths, vn)

actOut = Copulas.allExcept(paths, 4)
expOut = [paths[1], paths[2], paths[3], paths[5], paths[6]]

@test actOut == expOut

##################
## findMaxDepth ##
##################

vn = Copulas.testvine(2)
paths = convert(Copulas.CTreePaths, vn)

expOut = 1
@test expOut == Copulas.findMaxDepth(paths)

## next example
vn = Copulas.testvine(6)
paths = convert(Copulas.CTreePaths, vn)

expOut = 10
@test expOut == Copulas.findMaxDepth(paths)

####################
## getLongestPath ##
####################

expOut = [5,6,3,1,4,2,14,13,12,11]
actOut = Copulas.getLongestPath(paths[7])
@test expOut == actOut

expOut = [6,4,2]
actOut = Copulas.getLongestPath(paths[5])
@test expOut == actOut

#############
## convert ##
#############

for ii=1:6
    vn = Copulas.testvine(ii)
    rmatr = convert(Copulas.RVMatrix, vn)
    vn2 = convert(Copulas.Vine, rmatr)
    @test vn.trees == vn2.trees
end

## convert two different RVines representing same vine
rmatr1 = Copulas.RVMatrix([4 0 0 0 0 0 0;
                           7 5 0 0 0 0 0;
                           6 7 1 0 0 0 0;
                           5 6 7 7 0 0 0;
                           1 1 6 2 6 0 0;
                           2 3 3 3 2 2 0;
                           3 2 2 6 3 3 3])

rmatr2 = Copulas.RVMatrix([7 0 0 0 0 0 0;
                           4 4 0 0 0 0 0;
                           5 6 6 0 0 0 0;
                           1 5 5 5 0 0 0;
                           2 1 1 1 1 0 0;
                           3 2 2 3 3 3 0;
                           6 3 3 2 2 2 2])

vn1 = convert(Copulas.Vine, rmatr1)
vn2 = convert(Copulas.Vine, rmatr2)

@test Copulas.dg(vn1) == Copulas.dg(vn2)

end


text:
using Test, OrdinaryDiffEq, HyperbolicDiffEq, Roots

function calc_error(balance_law, uₐₙₐ, tmin, tmax, fnum, reconstruction, N, parallel)
    u₀ = x -> uₐₙₐ(tspan[1], x)

    tspan = (tmin, tmax)
    mesh = UniformPeriodicMesh1D(0., 2., N)
    fv = UniformPeriodicReconstructedFV1D(balance_law, mesh, fnum, reconstruction, parallel)
    ode = semidiscretise(fv, u₀, tspan)
    if order(reconstruction) <= 4
        sol = solve(ode, SSPRK104(), dt=0.5*max_dt(tspan[1], ode.u0, fv),
                    save_everystep=false, dense=false)
    else
        sol = solve(ode, SSPRK104(), dt=0.01*max_dt(tspan[1], ode.u0, fv),
                    save_everystep=false, dense=false)
    end
    uana = compute_coefficients(x->uₐₙₐ(tspan[end],x), mesh, order(reconstruction))

    N, integrate(u->sum(abs.(u)), sol[end]-uana, mesh)
end

function calc_order_estimate(Ns_and_errors)
    orders = zeros(length(Ns_and_errors)-1)
    for i in 1:length(orders)
        orders[i] = -log(Ns_and_errors[i+1][2]/Ns_and_errors[i][2]) /
                     log(Ns_and_errors[i+1][1]/Ns_and_errors[i][1])
    end
    mean(orders)
end

function calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, reconstruction, Ns, parallel=Val{:serial}())
    calc_order_estimate(calc_error.(balance_law, uₐₙₐ, tspan..., fnum, reconstruction, Ns, parallel))
end

# test convergence orders
Ns_lo = 10 .* 2 .^ (2:6)
Ns_hi =  4 .* 2 .^ (2:5)

balance_law = ConstantLinearAdvection()
tspan = (0., 1.)
fnum = LocalLaxFriedrichsFlux()
uₐₙₐ = (t,x) -> sinpi(x-t)

println("CentralReconstruction")
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, CentralReconstruction(1), Ns_lo, Val{:serial}()) > 0.8
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, CentralReconstruction(3), Ns_lo, Val{:serial}()) > 2.8
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, CentralReconstruction(5), Ns_lo, Val{:serial}()) > 4.8
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, CentralReconstruction(7), Ns_hi, Val{:serial}()) > 6.8

@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, CentralReconstruction(1), Ns_lo, Val{:threads}()) > 0.8
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, CentralReconstruction(3), Ns_lo, Val{:threads}()) > 2.8
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, CentralReconstruction(5), Ns_lo, Val{:threads}()) > 4.8
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, CentralReconstruction(7), Ns_hi, Val{:threads}()) > 6.8

println("ModifiedENO")
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(1), Ns_lo, Val{:serial}()) > 0.8
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(2), Ns_lo, Val{:serial}()) > 1.6
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(3), Ns_lo, Val{:serial}()) > 2.8
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(4), Ns_lo, Val{:serial}()) > 3.6
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(5), Ns_lo, Val{:serial}()) > 4.8
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(6), Ns_hi, Val{:serial}()) > 5.6
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(7), Ns_hi, Val{:serial}()) > 6.8
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(8), Ns_hi, Val{:serial}()) > 7.5

@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(1), Ns_lo, Val{:threads}()) > 0.8
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(2), Ns_lo, Val{:threads}()) > 1.6
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(3), Ns_lo, Val{:threads}()) > 2.8
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(4), Ns_lo, Val{:threads}()) > 3.6
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(5), Ns_lo, Val{:threads}()) > 4.8
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(6), Ns_hi, Val{:threads}()) > 5.6
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(7), Ns_hi, Val{:threads}()) > 6.8
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(8), Ns_hi, Val{:threads}()) > 7.5

println("MinL2Choice")
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(Val{1}(), MinL2Choice()), Ns_lo, Val{:serial}()) > 0.8
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(Val{2}(), MinL2Choice()), Ns_lo, Val{:serial}()) > 1.6
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(Val{3}(), MinL2Choice()), Ns_lo, Val{:serial}()) > 2.8
@test 0.2 < calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(Val{4}(), MinL2Choice()), Ns_lo, Val{:serial}()) < 0.4
@test 0.2 < calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(Val{5}(), MinL2Choice()), Ns_lo, Val{:serial}()) < 0.4
@test 0.5 < calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(Val{6}(), MinL2Choice()), Ns_hi, Val{:serial}()) < 1.1
@test 0.5 < calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(Val{7}(), MinL2Choice()), Ns_hi, Val{:serial}()) < 1.1
@test 0.5 < calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(Val{8}(), MinL2Choice()), Ns_hi, Val{:serial}()) < 1.1

println("BiasedMinL2Choice")
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(Val{1}(), BiasedMinL2Choice()), Ns_lo, Val{:serial}()) > 0.8
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(Val{2}(), BiasedMinL2Choice()), Ns_lo, Val{:serial}()) > 1.6
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(Val{3}(), BiasedMinL2Choice()), Ns_lo, Val{:serial}()) > 2.8
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(Val{4}(), BiasedMinL2Choice()), Ns_lo, Val{:serial}()) > 3.6
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(Val{5}(), BiasedMinL2Choice()), Ns_lo, Val{:serial}()) > 4.8

# Order reduction observed by Rogerson & Meiburg and Shu (1990).
balance_law = ConstantLinearAdvection()
tspan = (0., 1.)
fnum = LocalLaxFriedrichsFlux()
uₐₙₐ = (t,x) -> (sinpi(x-t)^2)^2
println("ClassicalChoiceENO")
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(Val{2}(), ClassicalChoiceENO()), Ns_lo, Val{:serial}()) > 1.6
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(Val{3}(), ClassicalChoiceENO()), Ns_lo, Val{:serial}()) > 2.2
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(Val{4}(), ClassicalChoiceENO()), Ns_lo, Val{:serial}()) < 2.2
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(Val{5}(), ClassicalChoiceENO()), Ns_lo, Val{:serial}()) < 1.9

println("BiasedENOChoice")
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(Val{2}(), BiasedENOChoice()), Ns_lo, Val{:serial}()) > 1.6
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(Val{3}(), BiasedENOChoice()), Ns_lo, Val{:serial}()) > 2.6
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(Val{4}(), BiasedENOChoice()), Ns_lo, Val{:serial}()) > 3.4
@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(Val{5}(), BiasedENOChoice()), Ns_lo, Val{:serial}()) > 4.6

# test evaluation
mesh = UniformPeriodicMesh1D(0., 2., 50)

K = 1
fv = UniformPeriodicReconstructedFV1D(balance_law, mesh, fnum, ModifiedENO(K), Val{:serial})
u = compute_coefficients(sinpi, mesh, K)
xplot, uplot = evaluate_coefficients(u, fv)
@test norm(sinpi.(xplot) - uplot) < 5.e-1
K = 2
fv = UniformPeriodicReconstructedFV1D(balance_law, mesh, fnum, ModifiedENO(K), Val{:serial})
u = compute_coefficients(sinpi, mesh, K)
xplot, uplot = evaluate_coefficients(u, fv)
@test norm(sinpi.(xplot) - uplot) < 4.e-2
K = 3
fv = UniformPeriodicReconstructedFV1D(balance_law, mesh, fnum, ModifiedENO(K), Val{:serial})
u = compute_coefficients(sinpi, mesh, K)
xplot, uplot = evaluate_coefficients(u, fv)
@test norm(sinpi.(xplot) - uplot) < 2.e-3
K = 4
fv = UniformPeriodicReconstructedFV1D(balance_law, mesh, fnum, ModifiedENO(K), Val{:serial})
u = compute_coefficients(sinpi, mesh, K)
xplot, uplot = evaluate_coefficients(u, fv)
@test norm(sinpi.(xplot) - uplot) < 2.e-4
K = 5
fv = UniformPeriodicReconstructedFV1D(balance_law, mesh, fnum, ModifiedENO(K), Val{:serial})
u = compute_coefficients(sinpi, mesh, K)
xplot, uplot = evaluate_coefficients(u, fv)
@test norm(sinpi.(xplot) - uplot) < 7.e-6
K = 6
fv = UniformPeriodicReconstructedFV1D(balance_law, mesh, fnum, ModifiedENO(K), Val{:serial})
u = compute_coefficients(sinpi, mesh, K)
xplot, uplot = evaluate_coefficients(u, fv)
@test norm(sinpi.(xplot) - uplot) < 5.e-7
K = 7
fv = UniformPeriodicReconstructedFV1D(balance_law, mesh, fnum, ModifiedENO(K), Val{:serial})
u = compute_coefficients(sinpi, mesh, K)
xplot, uplot = evaluate_coefficients(u, fv)
@test norm(sinpi.(xplot) - uplot) < 3.e-8
K = 8
fv = UniformPeriodicReconstructedFV1D(balance_law, mesh, fnum, ModifiedENO(K), Val{:serial})
u = compute_coefficients(sinpi, mesh, K)
xplot, uplot = evaluate_coefficients(u, fv)
@test norm(sinpi.(xplot) - uplot) < 2.e-9


text:
struct DeviceDurationConstraintInfo <: PSI.AbstractStartConstraintInfo
    component_name::String
    duration_data::PSI.UpDown
    initial_duration::Tuple{PSI.InitialCondition, PSI.InitialCondition}
    multiplier::Float64
    timeseries::Vector{Float64}
end

struct DeviceOutageConstraintInfo <: PSI.AbstractStartConstraintInfo
    component_name::String
    initial_condition::Union{Nothing, PSI.InitialCondition}
    multiplier::Float64
    timeseries::Vector{Float64}
end


text:
export BrickletIndustrialAnalogOutV2Configuration
struct BrickletIndustrialAnalogOutV2Configuration
    voltage_range::Integer
    current_range::Integer
end

export BrickletIndustrialAnalogOutV2OutLEDStatusConfig
struct BrickletIndustrialAnalogOutV2OutLEDStatusConfig
    min::Integer
    max::Integer
    config::Integer
end

export BrickletIndustrialAnalogOutV2SPITFPErrorCount
struct BrickletIndustrialAnalogOutV2SPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletIndustrialAnalogOutV2Identity
struct BrickletIndustrialAnalogOutV2Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletIndustrialAnalogOutV2
"""
Generates configurable DC voltage and current, 0V to 10V and 4mA to 20mA
"""
mutable struct BrickletIndustrialAnalogOutV2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletIndustrialAnalogOutV2(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_industrial_analog_out_v2")
        deviceInternal = package.BrickletIndustrialAnalogOutV2(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export set_enabled
"""
    $(SIGNATURES)

Enables/disables the output of voltage and current.
"""
function set_enabled(device::BrickletIndustrialAnalogOutV2, enabled)
    device.deviceInternal.set_enabled(enabled)
end

export get_enabled
"""
    $(SIGNATURES)

Returns *true* if output of voltage and current is enabled, *false* otherwise.
"""
function get_enabled(device::BrickletIndustrialAnalogOutV2)
    return device.deviceInternal.get_enabled()
end

export set_voltage
"""
    $(SIGNATURES)

Sets the output voltage.

The output voltage and output current are linked. Changing the output voltage
also changes the output current.
"""
function set_voltage(device::BrickletIndustrialAnalogOutV2, voltage)
    device.deviceInternal.set_voltage(voltage)
end

export get_voltage
"""
    $(SIGNATURES)

Returns the voltage as set by :func:`Set Voltage`.
"""
function get_voltage(device::BrickletIndustrialAnalogOutV2)
    return device.deviceInternal.get_voltage()
end

export set_current
"""
    $(SIGNATURES)

Sets the output current.

The output current and output voltage are linked. Changing the output current
also changes the output voltage.
"""
function set_current(device::BrickletIndustrialAnalogOutV2, current)
    device.deviceInternal.set_current(current)
end

export get_current
"""
    $(SIGNATURES)

Returns the current as set by :func:`Set Current`.
"""
function get_current(device::BrickletIndustrialAnalogOutV2)
    return device.deviceInternal.get_current()
end

export set_configuration
"""
    $(SIGNATURES)

Configures the voltage and current range.

Possible voltage ranges are:

* 0V to 5V
* 0V to 10V

Possible current ranges are:

* 4mA to 20mA
* 0mA to 20mA
* 0mA to 24mA

The resolution will always be 12 bit. This means, that the
precision is higher with a smaller range.
"""
function set_configuration(device::BrickletIndustrialAnalogOutV2, voltage_range, current_range)
    device.deviceInternal.set_configuration(voltage_range, current_range)
end

export get_configuration
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Configuration`.
"""
function get_configuration(device::BrickletIndustrialAnalogOutV2)
    return device.deviceInternal.get_configuration()
end

export set_out_led_config
"""
    $(SIGNATURES)

You can turn the Out LED off, on or show a
heartbeat. You can also set the LED to "Out Status". In this mode the
LED can either be turned on with a pre-defined threshold or the intensity
of the LED can change with the output value (voltage or current).

You can configure the channel status behavior with :func:`Set Out LED Status Config`.
"""
function set_out_led_config(device::BrickletIndustrialAnalogOutV2, config)
    device.deviceInternal.set_out_led_config(config)
end

export get_out_led_config
"""
    $(SIGNATURES)

Returns the Out LED configuration as set by :func:`Set Out LED Config`
"""
function get_out_led_config(device::BrickletIndustrialAnalogOutV2)
    return device.deviceInternal.get_out_led_config()
end

export set_out_led_status_config
"""
    $(SIGNATURES)

Sets the Out LED status config. This config is used if the Out LED is
configured as "Out Status", see :func:`Set Out LED Config`.

For each channel you can choose between threshold and intensity mode.

In threshold mode you can define a positive or a negative threshold.
For a positive threshold set the "min" parameter to the threshold value in mV or
µA above which the LED should turn on and set the "max" parameter to 0. Example:
If you set a positive threshold of 5V, the LED will turn on as soon as the
output value exceeds 5V and turn off again if it goes below 5V.
For a negative threshold set the "max" parameter to the threshold value in mV or
µA below which the LED should turn on and set the "min" parameter to 0. Example:
If you set a negative threshold of 5V, the LED will turn on as soon as the
output value goes below 5V and the LED will turn off when the output value
exceeds 5V.

In intensity mode you can define a range mV or µA that is used to scale the brightness
of the LED. Example with min=2V, max=8V: The LED is off at 2V and below, on at
8V and above and the brightness is linearly scaled between the values 2V and 8V.
If the min value is greater than the max value, the LED brightness is scaled the
other way around.
"""
function set_out_led_status_config(device::BrickletIndustrialAnalogOutV2, min, max, config)
    device.deviceInternal.set_out_led_status_config(min, max, config)
end

export get_out_led_status_config
"""
    $(SIGNATURES)

Returns the Out LED status configuration as set by :func:`Set Out LED Status Config`.
"""
function get_out_led_status_config(device::BrickletIndustrialAnalogOutV2)
    return device.deviceInternal.get_out_led_status_config()
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletIndustrialAnalogOutV2)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletIndustrialAnalogOutV2, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletIndustrialAnalogOutV2)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletIndustrialAnalogOutV2, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletIndustrialAnalogOutV2, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletIndustrialAnalogOutV2, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletIndustrialAnalogOutV2)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletIndustrialAnalogOutV2)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletIndustrialAnalogOutV2)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletIndustrialAnalogOutV2, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletIndustrialAnalogOutV2)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletIndustrialAnalogOutV2)
    return device.deviceInternal.get_identity()
end


text:
"""Load from XML"""
function read_curve_xml(path::AbstractString)
    last(splitext(path)) == ".xml" || throw("$(path) is not XML file")
    # The structure is doc
    #                    mesh
    #                     nodes
    #                     cells
    xroot = root(parse_file(path))
    mesh = first(child_elements(xroot))

    celltype = Symbol(attribute(mesh, "celltype"))
    # We dispatch on the cell type: interval or else
    read_curve_xml(path, Val{celltype})
end

"""Load line mesh as points and segments(by indices)"""
function read_curve_xml(path::AbstractString, cell_type::Type{Val{:interval}})
    last(splitext(path)) == ".xml" || throw("$(path) is not XML file")
    # The structure is doc
    #                    mesh
    #                     nodes
    #                     cells
    xroot = root(parse_file(path))
    mesh = first(child_elements(xroot))

    dim = parse(Int, attribute(mesh, "dim"))

    nodes, cells = child_elements(mesh)

    npoints = parse(Int, attribute(nodes, "size"))
    points = Vector{SVector{dim, Float64}}(npoints)
    # Each node is an elemn with x, y [, z] attributes
    get_coords = (dim == 2) ? (
        e -> SVector(parse(Float64, attribute(e, "x")),
                     parse(Float64, attribute(e, "y")))
    ) : (
        e -> SVector(parse(Float64, attribute(e, "x")),
                     parse(Float64, attribute(e, "y")),
                     parse(Float64, attribute(e, "z"))))

    for (col, elm) in enumerate(child_elements(nodes))
        points[col] = get_coords(elm)
    end

    nsegments = parse(Int, attribute(cells, "size"))
    segments = Vector{Tuple{Int, Int}}(nsegments)

    for (i, cell) in enumerate(child_elements(cells))
        # One based indexing
        segments[i] = Tuple((1+parse(Int, attribute(cell, "v0")), 1+parse(Int, attribute(cell, "v1"))))
    end

    (points, segments)
end

"""Turn a triangle mesh into the edge graph mesh"""
function read_curve_xml(path::AbstractString, cell_type::Type{Val{:triangle}})
    last(splitext(path)) == ".xml" || throw("$(path) is not XML file")
    # The structure is doc
    #                    mesh
    #                     nodes
    #                     cells
    xroot = root(parse_file(path))
    mesh = first(child_elements(xroot))

    nodes, cells = child_elements(mesh)

    npoints = parse(Int, attribute(nodes, "size"))
    points = Vector{SVector{3, Float64}}(npoints)
    # Each node is an elemn with x, y z attributes
    for (col, elm) in enumerate(child_elements(nodes))
        points[col] = SVector(parse(Float64, attribute(elm, "x")),
                              parse(Float64, attribute(elm, "y")),
                              parse(Float64, attribute(elm, "z")))
    end

    # Local indides of triangle
    edge_indices = collect(Combinatorics.combinations((1, 2, 3), 2))
    # Now we break each cell into segments
    segments = Set{Tuple{Int, Int}}()
    for cell in imap(c -> Tuple((1+parse(Int, attribute(c, "v0")),
                                 1+parse(Int, attribute(c, "v1")),
                                 1+parse(Int, attribute(c, "v2")))),
                     child_elements(cells))
        # All the edges
        for (i0, i1) in edge_indices
            # One based indexing
            v0, v1 = cell[i0], cell[i1]   # Global
            push!(segments, (v0 < v1) ? (v0, v1) : (v1, v0))
        end
    end
    # Unique can be converted
    segments = collect(segments)
    # If we are not in the z plane
    norm(map(last, points), 1) > 0 && return (points, segments)

    # Project to 2d
    points2d = Vector{SVector{2, Float64}}(npoints)
    for (row, point) in enumerate(points)
        points2d[row] = SVector(point[1], point[2])
    end
    (points2d, segments)
end

"""Turn a tetrahedron mesh into the edge graph mesh"""
function read_curve_xml(path::AbstractString, cell_type::Type{Val{:tetrahedron}})
    last(splitext(path)) == ".xml" || throw("$(path) is not XML file")
    # The structure is doc
    #                    mesh
    #                     nodes
    #                     cells
    xroot = root(parse_file(path))
    mesh = first(child_elements(xroot))

    nodes, cells = child_elements(mesh)

    npoints = parse(Int, attribute(nodes, "size"))
    points = Vector{SVector{3, Float64}}(npoints)
    # Each node is an elemn with x, y z attributes
    for (col, elm) in enumerate(child_elements(nodes))
        points[col] = SVector(parse(Float64, attribute(elm, "x")),
                              parse(Float64, attribute(elm, "y")),
                              parse(Float64, attribute(elm, "z")))
    end

    # Cell local vertex indices
    edge_indices = collect(Combinatorics.combinations((1, 2, 3, 4), 2))
    # Now we break each cell into segments
    segments = Set{Tuple{Int, Int}}()
    for cell in imap(c -> Tuple((1+parse(Int, attribute(c, "v0")),
                                 1+parse(Int, attribute(c, "v1")),
                                 1+parse(Int, attribute(c, "v2")),
                                 1+parse(Int, attribute(c, "v3")))),
                     child_elements(cells))
        # Break cell to edges
        for (i0, i1) in edge_indices
            # One based indexing
            v0, v1 = cell[i0], cell[i1]
            push!(segments, (v0 < v1) ? (v0, v1) : (v1, v0))
        end
    end
    # Unique can be converted
    segments = collect(segments)

    (points, segments)
end


text:
all_F₁s = [25, 75, 400]#[2500, 1600, 900, 625, 400, 100, 25]#400:-25:25
all_hts = [1, 5, 10]#1:15#0#7#1:10

#ttype = "ridg"
ttype = "cosi"
#ttype = "rand"; fname_hrand = "hrand256Km2tk10filtnx32.mat"

SAVE_OUTPUT = true
PSI_HOVMOLLER = true

#-------------------
include("twolay.jl")

max_stars_repo_name: apaloczy/RoughOceanQGEnergetics

max_issues_repo_name: apaloczy/RoughOceanQGEnergetics

max_forks_repo_name: apaloczy/RoughOceanQGEnergetics


text:
@testset "1449.form-largest-integer-with-digits-that-add-up-to-target.jl" begin
    @test largest_number([4,3,2,5,6,7,2,5,5], 9) == big"7772"
    @test largest_number([7,6,5,5,5,6,8,7,8], 12) == big"85"
    @test largest_number([2,4,6,2,4,6,4,4,4], 5) == big"0"
    @test largest_number([6,10,15,40,40,40,40,40,40], 47) == big"32211"
end


text:
################################################################################
#  Copyright 2021, Marta Vanin                                                 #
################################################################################
function run_case_study_A(path_to_result_csv::String, nlsolver::Any; models::Array=[_PMDSE.ReducedIVRUPowerModel, _PMDSE.ReducedACRUPowerModel, _PMDSE.ReducedACPUPowerModel], abbreviation::Array = ["rIVR", "rACR", "rACP"], criteria::Array=["rwlav", "wls", "rwls"], set_rescaler::Int64 = 1000, power_base::Float64=1.0, start::Bool=true, vmin::Float64=0.7, vmax::Float64=1.3)

    # Input data
    rm_transfo = true
    rd_lines = true

    season = "summer"
    time_step = 144
    elm = ["load", "pv"]
    pfs = [0.95, 0.90]

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

    # Set path
    msr_path = joinpath(mktempdir(),"temp.csv")

    pf_solver = _PMD.optimizer_with_attributes(nlsolver...)

    se_solver = pf_solver

    df = _DF.DataFrame(ntw=Int64[], fdr=Int64[], solve_time=Float64[], n_bus=Int64[],
    termination_status=String[], objective=Float64[], criterion=String[], rescaler = Float64[], eq_model = String[],
            err_max_1=Float64[], err_max_2=Float64[], err_max_3=Float64[], err_avg_1 = Float64[], err_avg_2 = Float64[], err_avg_3 = Float64[], pbase=Float64[])

    for i in 1:length(models)
        mod = models[i]
        short = abbreviation[i]

        for criterion in criteria
            for ntw in 1:25 for fdr in 1:10
                data_path = _PMDSE.get_enwl_dss_path(ntw, fdr)
                if !isdir(dirname(data_path)) break end

                # Load the data
                data = _PMD.parse_file(_PMDSE.get_enwl_dss_path(ntw, fdr),data_model=_PMD.ENGINEERING);
                if rm_transfo _PMDSE.rm_enwl_transformer!(data) end
                if rd_lines _PMDSE.reduce_enwl_lines_eng!(data) end
                data["settings"]["sbase_default"] = power_base

                # Insert the load profiles
                _PMDSE.insert_profiles!(data, season, elm, pfs, t = time_step)

                # Transform data model
                data = _PMD.transform_data_model(data);

                # Solve the power flow
                pf_results = _PMD.solve_mc_pf(data, _PMD.ACPUPowerModel, se_solver)

                # Write measurements based on power flow
                v_pu = data["settings"]["vbases_default"]["1"]* data["settings"]["voltage_scale_factor"] # divider [V] to get the voltage in per units.
                v_max_err = 1.15 # maximum error of voltage measurement = 0.5% or 1.15 V
                σ_v = 1/3*v_max_err/v_pu

                p_pu = data["settings"]["sbase"] # divider [kW] to get the power in per units.
                p_max_err = 0.01 # maximum error of power measurement = 10W, or 0.01 kW
                σ_p = 1/3*p_max_err/p_pu

                # Write measurements based on power flow
                σ_dict = Dict("load" => Dict("load" => σ_p,
                                "bus"  => σ_v),
                                "gen"  => Dict("gen" => σ_p,
                                "bus"  => σ_v)
                                        )

                _PMDSE.write_measurements!(_PMD.ACPUPowerModel, data, pf_results, msr_path, σ= σ_dict)

                # Read-in measurement data and set initial values
                _PMDSE.add_measurements!(data, msr_path, actual_meas = true)
                _SEF.add_errors!(data)
                if start _PMDSE.assign_start_to_variables!(data) end

                _PMDSE.update_all_bounds!(data; v_min = vmin, v_max = vmax, pg_min=-Inf, pd_min=-Inf)

                # Set se settings
                data["se_settings"] = Dict{String,Any}("criterion" => criterion,
                                        "rescaler" => set_rescaler)

                # Solve the state estimation
                se_results = _PMDSE.solve_mc_se(data, _PMDSE.ReducedIVRUPowerModel, pf_solver)
                delta_1, delta_2, delta_3, max_1, max_2, max_3, mean_1, mean_2, mean_3 = _PMDSE.calculate_voltage_magnitude_error_perphase(se_results, pf_results)
                # store result
                push!(df, [ntw, fdr, se_results["solve_time"], length(data["bus"]),
                        string(se_results["termination_status"]),
                        se_results["objective"], criterion, set_rescaler, short, max_1, max_2, max_3, mean_1, mean_2, mean_3, power_base])
            end end #loop through feeder and network
        CSV.write(path_to_result_csv, df)
        end #criterion loop
    end #end models loop
end


text:
######################################################################
# lqr.jl
#
# Solves LQR problems
######################################################################

export LQR
export apply_LQR_gains

function LQR(A::MF, B::MF, Q::MF, R::MF, N::Int)
	P = Array(Matrix{Float64}, N+1)
	G = Array(Matrix{Float64}, N)
	K = Array(Matrix{Float64}, N)
	P[N+1] = Q

	for n = (N-1):-1:0
		G[n+1] = R + (B' * P[n+1+1] * B)
		K[n+1] = inv(G[n+1]) * B' * P[n+1+1] * A
		P[n+1] = Q + (A' * P[n+1+1] * A) - (K[n+1]' * G[n+1] * K[n+1])
	end

	return K
end

function apply_LQR_gains(A::MF, B::MF, K::Vector{MF}, x0::Vector{Float64})
	x = [x0]
	u = Array(Vector{Float64}, 0)

	N = length(K)

	for n = 1:N
		push!(u, -K[n]*x[n])
		push!(x, A*x[n] + B*u[n])
	end
	return x, u
end


text:
"""
Terrestrial water cycle: evapotranspiration, precipitation, and runoff.

* [`psychromet`](@ref)
* [`penman_monteith`](@ref)
* [`cond_sfc`](@ref)

# References

* [Shu12] Shuttleworth, W. J. (2012). *Terrestrial Hydrometeorology*. John
  Wiley & Sons. ISBN: 978-1119951933. <https://doi.org/10.1002/9781119951933>
* [MB14] Monson, R. and Baldocchi, D. (2014). *Terrestrial
  Biosphere--Atmosphere Fluxes*. Cambridge University Press, Cambridge, UK.
  ISBN: 978-1107040656.
"""
module WaterCycle

include("docstring_style.jl")

using Canopy.Constants: M_w, M_d
using Canopy.Water: latent_heat_vap, e_sat_deriv, vapor_deficit
using Canopy.Air: air_density
using Canopy.Transfer.Heat: heat_cap_mass_moistair

export psychromet,
    penman_monteith,
    cond_sfc

"""
Calculate the psychromet constant [Pa K^-1] from temperature [K], pressure
[Pa], and relative humidity [0--1].

# Examples

```jldoctest
julia> psychromet(298.15, 101325.0, 0.7)
67.91408269550081

julia> psychromet(283.15, 101325.0, 0.9)
66.5430908894645
```
"""
psychromet(temp, pressure, RH) =
    heat_cap_mass_moistair(temp, pressure, RH) * pressure *
    M_d / latent_heat_vap(temp)
# note: M_w cancels out because latent_heat_vap returns value in J mol^-1 K^-1

"""
Calculate latent heat flux [W m^-2] or evapotranspiration [mol m^-2 s^-1] from
the Penman--Monteith equation.

```math
LE = \\dfrac
    {\\Delta \\cdot (R_\\mathrm{net} - G) +
     G_\\mathrm{atm} \\rho_\\mathrm{a} c_p D}
    {\\Delta + \\gamma (1 + G_\\mathrm{a} / G_\\mathrm{c})}
```

where \$\\Delta := e_\\mathrm{sat}'(T)\$ and \$\\gamma\$ is the psychromet
constant.

# Arguments

* `temp`: Temperature [K].
* `pressure`: Ambient pressure [Pa].
* `RH`: Relative humidity [0--1].
* `heat_flux`: Net heat flux [W m^-2] which equals to net radiation
  (\$R_\\mathrm{net}\$) -- heat storage (\$G\$).
* `cond_atm`: Aerodynamic conductance [m s^-1].
* `cond_sfc`: Canopy surface conductance [m s^-1].
* `evap` (*kwarg*): Default is `false`. If `true`, return evapotranspiration
  flux instead of latent heat flux.
"""
function penman_monteith(temp, pressure, RH, heat_flux, cond_atm, cond_sfc;
                         evap::Bool=false)
    Delta = e_sat_deriv(temp)
    gamma = psychromet(temp, pressure, RH)
    rho_a = air_density(temp, pressure, RH)
    c_p = heat_cap_mass_moistair(temp, pressure, RH)
    D = vapor_deficit(temp, RH)
    LE = (Delta * heat_flux + rho_a * c_p * D * cond_atm) /
        (Delta + gamma * (1. + cond_atm / cond_sfc))
    return evap ? LE / latent_heat_vap(temp) : LE
end

"""
Calculate the canopy surface conductance [m s^-1] from the inverted
Penman--Monteith equation.

```math
G_\\mathrm{c} = \\gamma \\left[
    \\dfrac{G_\\mathrm{a}\\rho_\\mathrm{a} c_p D}{LE} +
    \\Delta \\cdot \\beta - \\gamma \\right]^{-1} G_\\mathrm{a}
```

where \$\\beta := \\dfrac{H}{LE}\$ is the Bowen ratio.

# Arguments

* `temp`: Temperature [K].
* `pressure`: Ambient pressure [Pa].
* `RH`: Relative humidity [0--1].
* `LE`: Latent heat flux [W m^-2].
* `bowen`: Bowen ratio.
* `cond_atm`: Aerodynamic conductance [m s^-1].
"""
function cond_sfc(temp, pressure, RH, LE, bowen, cond_atm)
    Delta = e_sat_deriv(temp)
    gamma = psychromet(temp, pressure, RH)
    rho_a = air_density(temp, pressure, RH)
    c_p = heat_cap_mass_moistair(temp, pressure, RH)
    D = vapor_deficit(temp, RH)
    return cond_atm * gamma / (cond_atm * rho_a * c_p * D / LE +
                               Delta * bowen - gamma)
end

end  # module


text:
"Converts a bitmatrix to plain PBM format."
function pbm(bitmatrix)
    join(
        (
            "P1",
            join(reverse(size(bitmatrix)), ' '),
            (
                join(map(b -> b ? '1' : '0', row), ' ')
                for row in Iterators.partition(transpose(bitmatrix), 34)
            )...
        ),
        '\n'
    ) * '\n'
end

pbm(filename::AbstractString, bitmatrix) = open(f -> write(f, pbm(bitmatrix)), filename, "w")

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
using PhasePlots
using Test


text:
using SpecialFunctions

include("spindown_besselzeros.jl")

function u_spindown(r, t, Re)
    # Assuming unit radius, speed and viscosity
    Omega, a = 1.0, 1.0
    nu = 1.0/Re
    S = 0.0
    ds = 0.0
    relerr = 0.0
    for n=1:length(lambda)
        ds = besselj(1, lambda[n]*r/a) / (lambda[n]*besselj(0, lambda[n]))*exp(-lambda[n]^2*nu*t/a)
        S += ds
        relerr = abs(ds/S) # If convergence is nice enough
        if relerr < eps()
            break
        end
    end
    if relerr > eps()
        warn("Spindown solution not converged, relerr=", relerr)
    end
    u = -2*Omega*a*S
end

function Lu_spindown(r, t, Re)
    # Laplacian of the above
    Omega, a = 1.0, 1.0
    nu = 1.0/Re
    S = 0.0
    ds = 0.0
    relerr = 0.0
    for n=1:length(lambda)
        b = lambda[n]/a
        ds = -b^2*besselj(1, b*r) / (lambda[n]*besselj(0, lambda[n]))*exp(-lambda[n]^2*nu*t/a)
        S += ds
        relerr = abs(ds/S) # If convergence is nice enough
        if relerr < eps()
            break
        end
    end
    if relerr > eps()
        warn("Spindown solution not converged, relerr=", relerr)
    end
    u = -2*Omega*a*S
end

function Du_spindown(r, t, Re)
    # r derivative
    Omega, a = 1.0, 1.0
    nu = 1.0/Re
    S = 0.0
    ds = 0.0
    relerr = 0.0
    for n=1:length(lambda)
        b = lambda[n]/a
        ds = b/2*(besselj(0, b*r)-besselj(2, b*r)) / (lambda[n]*besselj(0, lambda[n]))*exp(-lambda[n]^2*nu*t/a)
        S += ds
        relerr = abs(ds/S) # If convergence is nice enough
        if relerr < eps()
            break
        end
    end
    if relerr > eps()
        warn("Spindown solution not converged, relerr=", relerr)
    end
    u = -2*Omega*a*S
end

function plot_spindown(t, Re)
    r = linspace(0,1,200)
    u = map(x -> u_spindown(x, t, Re), r)
    plot(r, u, label="Reference")
end


text:
module Poisson1DTests

using Gridap
using Test

u(x) = x[1]^2 + 1
∇u(x) = ∇(u)(x)
f(x) = -Δ(u)(x)

domain = (0,3)
cells = (10,)
model = CartesianDiscreteModel(domain,cells)

order = 2
V = FESpace(model, ReferenceFE(lagrangian,Float64,order),conformity=:H1,dirichlet_tags=2)
U = TrialFESpace(V,u)

Ω = Triangulation(model)
Γ = BoundaryTriangulation(model,tags=1)
n_Γ = get_normal_vector(Γ)

degree = 2*order
dΩ = Measure(Ω,degree)
dΓ = Measure(Γ,degree)

a(u,v) = ∫( ∇(v)⋅∇(u) )*dΩ
l(v) = ∫( v*f )*dΩ + ∫( v*(n_Γ⋅∇u) )*dΓ

op = AffineFEOperator(a,l,U,V)
uh = solve(op)

e = u - uh

l2(u) = sqrt(sum( ∫( u⊙u )*dΩ ))
h1(u) = sqrt(sum( ∫( u⊙u + ∇(u)⊙∇(u) )*dΩ ))

el2 = l2(e)
eh1 = h1(e)
ul2 = l2(uh)
uh1 = h1(uh)

@test el2/ul2 < 1.e-8
@test eh1/uh1 < 1.e-7

end # module


text:
using Documenter, Epicrop, DataFrames, Dates, CSV, Plots

makedocs(
    sitename = "Epicrop.jl",
    format = Documenter.HTML(),
    modules = [Epicrop],
    doctest = true,
    pages = Any[
        "Home" => "index.md",
        "Manual" => Any[
            "hlipmodel" => "man/hlipmodel.md",
            "helper functions" => "man/helperfns.md"
        ]
    ]
)

deploydocs(
    repo = "github.com/adamhsparks/Epicrop.jl.git",
    devbranch = "main"
)


text:
using BayesianIntegral
using Distributions: Normal, MvNormal, pdf
using LinearAlgebra: diagm, Symmetric
using Statistics
using HCubature: hcubature
using Sobol
using Random
standard_normal = Normal(0.0,1.0)

# Integral of gaussian distribution should be 1.
samples = 29
dims = 1
p(x) = 1.0
X = Array{Float64,2}(undef,samples,1)
X[:,1] =  collect(-2.0:(4.0/(samples-1)):2.0) #
y = p.(X)[:,1]
cov_func_parameters = gaussian_kernel_hyperparameters(1.0, repeat([2.0] , outer = dims))
prob_means = [0.0]
covar = Array{Float64,2}(undef,1,1)
covar[1,1] = 1.0
covar = Symmetric(covar)
noise = 0.2
integ = bayesian_integral_gaussian_exponential(X, y, prob_means, covar, cov_func_parameters, noise)
(integ.expectation - 1) < 0.01

# Integral in two dimensions
samples = 25
dims = 2
p(x) = 1.0
s = SobolSeq(dims)
X = convert( Array{Float64}, hcat([next!(s, repeat([0.5] , outer = dims)     ) for i = 1:samples]...)' )
y = repeat([1.0] , outer = samples)
cov_func_parameters = gaussian_kernel_hyperparameters(1.0, repeat([10.0] , outer = dims))
prob_means = repeat([0.0] , outer = dims)
covar = Symmetric(diagm(0 => ones(dims)))
noise = 0.2
integ = bayesian_integral_gaussian_exponential(X, y, prob_means, covar, cov_func_parameters, noise)
(integ.expectation - 1) < 0.02

# In ten dimensinos
samples = 50
dims = 10
p(x) = 1.0
s = SobolSeq(dims)
X = convert(Array{Float64}, hcat([next!(s, repeat([0.5] , outer = dims)     ) for i = 1:samples]...)')
y = repeat([1.0] , outer = samples)
cov_func_parameters = gaussian_kernel_hyperparameters(1.0, repeat([20.0] , outer = dims))
prob_means = repeat([0.0] , outer = dims)
covar = Symmetric(diagm(0 => ones(dims)))
noise = 0.2
integ = bayesian_integral_gaussian_exponential(X, y, prob_means, covar, cov_func_parameters, noise)
(integ.expectation - 1) < 0.02

##### More complex cases #####
# Now looking at a more complex one.
function f(x::Array{Float64,1})
    ndims = length(x)
    total = 12.0
    for i in 1:ndims
        total = total - x[i]^2 - 2*x[i]
    end
    return total
end
function fp(x::Array{Float64,1})
    ndim = length(x)
    if ndim == 1
        return f(x) * pdf.(Ref(standard_normal), x)
    else
        dist = MvNormal(zeros(ndim),diagm(0 => ones(ndim)))
        return f(x) * pdf(dist, x)
    end
end
function fp(x)
    global counter = counter + 1
    flipped = convert(Array{Float64,1},x)
    return fp(flipped)
end

dims = 1
bayesianAttempts = 10

function compare_for_f(dims::Int, bayesianAttempts::Int = 0, paths::Int = 0;  seed::Int = 1988)
    Random.seed!(1234)
    # Traditional
    global counter = 0
    lims = 3.0
    maxIter = 5000
    numerical_val, numerical_err = hcubature(x -> fp(x), repeat([-lims], outer = dims), repeat([lims], outer = dims); maxevals = maxIter)
    cont = counter

    # Kriging
    if bayesianAttempts < 1
        samples = cont
    else
        samples = bayesianAttempts
    end
    s = SobolSeq(dims)
    X = (hcat([next!(s, repeat([0.5] , outer = dims)     ) for i = 1:samples]...)' .- 0.5) .* (lims*2.0)
    y = Array{Float64,1}(undef, samples)
    for r in 1:samples
        y[r] = f(X[r,:])
    end
    noise = 0.05
    cov_func_parameters = gaussian_kernel_hyperparameters(1.0, repeat([1.0] , outer = dims))
    prob_means = repeat([0.0], outer = dims)
    covar = Symmetric(diagm(0 => ones(dims)))
    bayesian_val, bayesian_err = bayesian_integral_gaussian_exponential(X, y, prob_means, covar, cov_func_parameters, noise)

    # Kriging with new weights
    steps = 2000
    batch_size = convert(Int, floor( samples / 4 ) )
    step_multiple = 0.02
    seed = 1988
    K_mat = K_matrix(X, gaussian_kernel, cov_func_parameters, noise)
    like = log_likelihood(y, K_mat)
    n_cov_func_parameters = calibrate_by_ML_with_SGD(X, y, cov_func_parameters, steps, batch_size, step_multiple, noise, seed)
    K_mat = K_matrix(X, gaussian_kernel, n_cov_func_parameters, noise)
    nlike = log_likelihood(y, K_mat)
    nbayesian_val, nbayesian_err = bayesian_integral_gaussian_exponential(X, y , prob_means, covar, n_cov_func_parameters, noise)

    # Kriging with new weights and RProp
    params = RProp_params(1.01,0.99,0.2,5.0,0.5)
    MaxIter = 2000
    nn_cov_func_parameters = calibrate_by_ML_with_Rprop(X, y, cov_func_parameters, MaxIter, noise, params)
    K_mat = K_matrix(X, gaussian_kernel, n_cov_func_parameters, noise)
    rproplikelihood = log_likelihood(y, K_mat)
    rpropnbayesian_val, rpropnbayesian_err = bayesian_integral_gaussian_exponential(X, y , prob_means, covar, nn_cov_func_parameters, noise)

    # MC Integration
    if paths < 1
        paths = cont
    end
    dist = MvNormal(zeros(dims),diagm(0 => ones(dims)))
    Xs = transpose(convert(Matrix{Float64}, rand(dist, paths)))
    ys = Array{Float64,1}(undef, paths)
    for r in 1:paths
        ys[r] = f(Xs[r,:])
    end
    MC_integral = mean(ys)
    MC_err = std(ys) / sqrt(paths)

    print("For ", dims, " dimensions \n")
    print("    Bayesian    Integral is               ", bayesian_val       , " with error ", bayesian_err       ," and ", samples , " evaluations \n")
    print("    Calibrated Bayesian Integral is       ", nbayesian_val      , " with error ", nbayesian_err      ," and ", samples , " evaluations \n")
    print("    Rprop Calibrated Bayesian Integral is ", rpropnbayesian_val , " with error ", rpropnbayesian_err ," and ", samples , " evaluations \n")
    print("    Traditional Integral is               ", numerical_val[1]   , " with error ", numerical_err[1]   ," and ", cont    , " evaluations \n")
    print("    MC          Integral is               ", MC_integral        , " with error ", MC_err             ," and ", paths   , " evaluations \n")
    print("        Likelihood was ", like   , " and calibrated is ", nlike, " and Rprop calibrated is ", rproplikelihood  ,"\n")
    print("        weights became ", n_cov_func_parameters.w_0   , " and  ", n_cov_func_parameters.w_i, " with SGD \n")
    print("        weights became ", nn_cov_func_parameters.w_0    , " and  ", nn_cov_func_parameters.w_i, " with RProp \n")
    return (bayesian_val = bayesian_val, bayesian_err = bayesian_err, samples = samples, traditional_value = numerical_val[1], traditional_error = numerical_err, traditional_evals = cont, MC_integral = MC_integral, MC_err = MC_err, paths = paths)
end

f1_results = compare_for_f(1, 100)
abs(f1_results.bayesian_val - f1_results.MC_integral) < 0.01
f2_results = compare_for_f(2, 300)
abs(f2_results.bayesian_val - f2_results.MC_integral) < 0.11
#f10_results = compare_for_f(10, 200, 20000)
#f20_results = compare_for_f(20, 2000, 2000)


text:
using Documenter, DistributedData

makedocs(
    modules = [DistributedData],
    clean = false,
    format = Documenter.HTML(
        prettyurls = !("local" in ARGS),
        canonical = "https://lcsb-biocore.github.io/DistributedData.jl/stable/",
        assets = ["assets/logo.ico"],
    ),
    sitename = "DistributedData.jl",
    authors = "The developers of DistributedData.jl",
    linkcheck = !("skiplinks" in ARGS),
    pages = [
        "Documentation" => "index.md",
        "Tutorial" => "tutorial.md",
        "Function reference" => "functions.md",
    ],
)

deploydocs(
    repo = "github.com/LCSB-BioCore/DistributedData.jl.git",
    target = "build",
    branch = "gh-pages",
    devbranch = "develop",
    push_preview = true,
)


text:
abstract type Component{T} end
abstract type AbstractBody{T} <: Component{T} end
abstract type AbstractConstraint{T,N} <: Component{T} end

# TODO do id differently?
CURRENTID = -1
getGlobalID() = (global CURRENTID -= 1; return CURRENTID + 1)
resetGlobalID() = (global CURRENTID = -1; return)

Base.show(io::IO, component::Component) = summary(io, component)

@inline Base.foreach(f, itr::Vector{<:Component}, arg...) = (for x in itr; f(x, arg...); end; return)


text:
using Test
using NBInclude

testdir = joinpath(dirname(@__DIR__), "test")
notebookdir = joinpath(dirname(@__DIR__), "notebooks")

function test_solvers(solvers=[MCTSPWSolver, CEMSolver, PPOSolver, TRPOSolver, RandomSearchSolver]; skip_trpo=false)
    for solver in solvers
        @test begin
            if skip_trpo && solver == TRPOSolver
                @info "Skipping TRPOSolver."
            else
                @info solver
                run_ast(solver=solver)
            end
            true
        end
    end
end

@test begin
    @info "Walk1D Jupyter notebook"
    @nbinclude(joinpath(notebookdir, "Walk1D.ipynb"))
    true
end

cd(testdir)

@test begin
    @info "Walk1D"
    include("Walk1D.jl")
    true
end

@test begin
    @info "Extra functions test"
    (planner, _, _) = run_ast(solver=MCTSPWSolver)
    d3tree = visualize(planner)
    action_trace = search!(planner; verbose=true)
    actions = online_path(planner.mdp, planner)
    x_trace = playback(planner, actions, sim->sim.x; return_trace=true)
    final_state = playback(planner, actions, sim->sim.x)
    AST.state_trace(final_state)
    AST.q_trace(final_state)
    reset_metrics!(planner.mdp)

    # Policy saving
    solver = PPOSolver(num_episodes=100, episode_length=30, save=true, verbose=true)
    planner = solve(solver, planner.mdp)
    search!(planner)

    # Policy loading
    solver.resume = true
    search!(planner)

    # No top_k action trace
    solver.resume = false
    planner.mdp.params.top_k = 0
    ast_action = search!(planner)
    true
end

# Walk1D with different solvers
test_solvers()

@test begin
    @info "Walk2D"
    include("Walk2D.jl")
    true
end
test_solvers(skip_trpo=true) # TRPO slows down the entire test suite, so skip it for Walk2D (coverd by Walk1D)

@test begin
    @info "EpisodicWalk1D"
    include("EpisodicWalk1D.jl")
    true
end

@test begin
    @info "CategoricalWalk1D"
    include("CategoricalWalk1D.jl")

    # Policy saving
    (planner, _, _) = run_ast(solver=PPOSolver)
    solver = PPOSolver(num_episodes=100, episode_length=30, save=true, verbose=true, policy_type=:discrete)
    planner = solve(solver, planner.mdp)
    search!(planner)

    # Policy loading
    solver.resume = true
    search!(planner)

    true
end


text:
using Plots
using LinearAlgebra

"
This code shows how errors for different finite difference methods converge as h->0.
"

# define forward, centered differences, D3 = some other difference approx
D1(u,x,h) = (u(x+h) - u(x))/h
D2(u,x,h) = (u(x+h) - u(x-h))/(2*h)
D3(u,x,h) = (u(x-2*h) - 6*u(x-h) + 3*u(x) + 2*u(x+h))/(6*h)

include("fdcoeffV.jl")
x = -2:2
a = fdcoeffV(1,x[3],x) # 1 = order of derivative, x̄, stencil points x
D4(u,x,h,c) = dot(c,(u(x-2*h),u(x-h),u(x),u(x+h),u(x+2*h)))/h

# function, exact derivative
u(x) = x^2 + sin(1+pi*x)
dudx(x) = 2*x + pi*cos(1+pi*x)
x = 0.0 # point at which we approximate the derivative

# measure how error changes as h decreases
hvec = 2.0 .^(-(0:8))
err = (zeros(size(hvec)),zeros(size(hvec)),zeros(size(hvec)),zeros(size(hvec)))
for (i,h) in enumerate(hvec)
    err[1][i] = abs(dudx(x) - D1(u,x,h))
    err[2][i] = abs(dudx(x) - D2(u,x,h))
    err[3][i] = abs(dudx(x) - D3(u,x,h))
    err[4][i] = abs(dudx(x) - D4(u,x,h,a))
end

plot(xaxis=:log,yaxis=:log,legend=:bottomright)
labels = ("Forward diff","Centered diff","Some other 4-point diff","5-point centered diff")
for i = 1:4
    plot!(hvec,err[i],marker=:dot,label=labels[i])
    plot!(hvec,hvec.^i,linestyle=:dash,label="Order $(i)")
end
display(plot!())


text:
isdefined(Base, :__precompile__) && __precompile__()

module Einsum

using Compat # for Array{}(undef,...)

export @einsum, @einsimd, @vielsum, @vielsimd

macro einsum(ex)
    _einsum(ex) # true, false, false
end

macro einsimd(ex)
    _einsum(ex, true, true) # false
end

macro vielsum(ex)
    _einsum(ex, true, false, true)
end

macro vielsimd(ex)
    _einsum(ex, true, true, true)
end

macro einsum_checkinbounds(ex)
    _einsum(ex, false) # false, false
end

function _einsum(expr::Expr, inbounds = true, simd = false, threads = false)
    # Get left hand side (lhs) and right hand side (rhs) of equation
    lhs = expr.args[1]
    rhs = expr.args[2]

    # Get info on the left-hand side
    lhs_arrays, lhs_indices, lhs_axis_exprs = extractindices(lhs)
    length(lhs_arrays) != 1 && throw(ArgumentError(string(
        "Left-hand side of equation contains multiple arguments. Only a single ",
        "referencing expression (e.g. @einsum A[i] = ...) should be used.")))

    # Get info on the right-hand side
    rhs_arrays, rhs_indices, rhs_axis_exprs = extractindices(rhs)

    check_index_occurrence(lhs_indices, rhs_indices)

    # remove duplicate indices on left-hand and right-hand side
    # and ensure that the array sizes match along these dimensions
    ###########################################################
    dimension_checks = Expr[]

    # remove duplicate indices on the right hand side
    for i in reverse(eachindex(rhs_indices))
        duplicated = false
        di = rhs_axis_exprs[i]

        for j = 1:(i - 1)
            if rhs_indices[j] == rhs_indices[i]
                # found a duplicate
                duplicated = true
                dj = rhs_axis_exprs[j]

                # add dimension check ensuring consistency
                push!(dimension_checks, :(@assert $dj == $di))
            end
        end

        for j = eachindex(lhs_indices)
            if lhs_indices[j] == rhs_indices[i]
                dj = lhs_axis_exprs[j]
                if Meta.isexpr(expr, :(:=))
                    # expr.head is :=
                    # infer the size of the lhs array
                    lhs_axis_exprs[j] = di
                else
                    # expr.head is =, +=, *=, etc.
                    lhs_axis_exprs[j] = :(min($dj, $di))
                end
                duplicated = true
            end
        end

        if duplicated
            deleteat!(rhs_indices, i)
            deleteat!(rhs_axis_exprs, i)
        end
    end

    # remove duplicate indices on the left hand side
    for i in reverse(eachindex(lhs_indices))
        duplicated = false
        di = lhs_axis_exprs[i]

        # don't need to check rhs, already done above

        for j = 1:(i - 1)
            if lhs_indices[j] == lhs_indices[i]
                # found a duplicate
                duplicated = true
                dj = lhs_axis_exprs[j]

                # add dimension check
                push!(dimension_checks, :(@assert $dj == $di))
            end
        end

        if duplicated
            deleteat!(lhs_indices, i)
            deleteat!(lhs_axis_exprs, i)
        end
    end

    # Create output array if specified by user
    @gensym T

    if Meta.isexpr(expr, :(:=))
        # infer type of allocated array
        #    e.g. rhs_arrays = [:A, :B]
        #    then the following line produces :(promote_type(eltype(A), eltype(B)))
        rhs_type = :(promote_type($([:(eltype($arr)) for arr in rhs_arrays]...)))

        type_definition = :(local $T = $rhs_type)

        output_definition = if length(lhs_axis_exprs) > 0
            :($(lhs_arrays[1]) = Array{$T}(undef, $(lhs_axis_exprs...)))
        else
            :($(lhs_arrays[1]) = zero($T))
        end

        assignment_op = :(=)
    else
        type_definition = :(local $T = eltype($(lhs_arrays[1])))
        output_definition = :(nothing)
        assignment_op = expr.head
    end

    if threads && !Meta.isexpr(expr, :(=)) && !Meta.isexpr(expr, :(:=))
        throw(ArgumentError(string(
            "Threaded @vielsum can only assign with = or := right now. ",
            "To use ", expr.head, " try @einsum instead.")))
        # could allow :(+=) by simply then removing $lhs = zero($T) line
    end

    if threads && length(lhs_indices) == 0
        # this won't actually cause problems, but won't use threads
        throw(ArgumentError(string(
            "Threaded @vielsum can't assign to a scalar LHS. ",
            "Try @einsum instead.")))
    end

    # copy the index expression to modify it; loop_expr is the Expr we'll build loops around
    loop_expr = unquote_offsets!(copy(expr))

    # Nest loops to iterate over the destination variables
    if length(rhs_indices) > 0
        # There are indices on rhs that do not appear in lhs.
        # We sum over these variables.

        if !threads # then use temporaries to write into, as before

            # Innermost expression has form s += rhs
            @gensym s
            loop_expr.args[1] = s
            loop_expr.head = :(+=)

            # Nest loops to iterate over the summed out variables
            loop_expr = nest_loops(loop_expr, rhs_indices, rhs_axis_exprs, simd, false)

            # Prepend with s = 0, and append with assignment
            # to the left hand side of the equation.
            lhs_assignment = Expr(assignment_op, lhs, s)

            loop_expr = quote
                local $s = zero($T)
                $loop_expr
                $lhs_assignment
            end

        else # we are threading, and thus should write directly to lhs array

            loop_expr.args[1] = lhs
            loop_expr.head = :(+=)

            loop_expr = nest_loops(loop_expr, rhs_indices, rhs_axis_exprs, simd, false)

            loop_expr = quote
                $lhs = zero($T)
                $loop_expr
            end
        end

        # Now loop over indices appearing on lhs, if any
        loop_expr = nest_loops(loop_expr, lhs_indices, lhs_axis_exprs, false, threads)
    else
        # We do not sum over any indices, only loop over lhs
        loop_expr.head = assignment_op
        loop_expr = nest_loops(loop_expr, lhs_indices, lhs_axis_exprs, simd, threads)
    end

    if inbounds
        loop_expr = :(@inbounds $loop_expr)
    end

    full_expression = quote
        $type_definition
        $output_definition
        $(dimension_checks...)

        # remove let when we drop 0.6 support -- see #31
        let $([lhs_indices; rhs_indices]...)
            $loop_expr
        end

        $(lhs_arrays[1])
    end

    return esc(full_expression)
end

function check_index_occurrence(lhs_indices, rhs_indices)
    if !issubset(lhs_indices, rhs_indices)
        missing_indices = setdiff(lhs_indices, rhs_indices)

        if length(missing_indices) == 1
            missing_string = "\"$(missing_indices[1])\""
            throw(ArgumentError(string(
                "Index ", missing_string, " is occuring only on left side")))
        else
            missing_string = join(("\"$ix\"" for ix in missing_indices),
                                  ", ", " and ")
            throw(ArgumentError(string(
                "Indices ", missing_string, " are occuring only on left side")))
        end
    end
end

"""
    nest_loops(expr, indices, axis_exprs, simd, threads) -> Expr

Construct a nested loop around `expr`, using `indices` in ranges `axis_exprs`.

# Example
```julia-repl
julia> nest_loops(:(A[i] = B[i]), [:i], [:(size(A, 1))], true, false)
quote
    local i
    @simd for i = 1:size(A, 1)
        A[i] = B[i]
    end
end
```
"""
function nest_loops(expr::Expr,
                    index_names::Vector{Symbol}, axis_expressions::Vector{Expr},
                    simd::Bool, threads::Bool)
    isempty(index_names) && return expr

    # Add @simd to the innermost loop, if required
    # and @threads to the outermost loop
    expr = nest_loop(expr, index_names[1], axis_expressions[1],
                     simd, threads && length(index_names) == 1)

    # Add remaining for loops
    for j = 2:length(index_names)
        expr = nest_loop(expr, index_names[j], axis_expressions[j],
                         false, threads && length(index_names) == j)
    end

    return expr
end

function nest_loop(expr::Expr, index_name::Symbol, axis_expression::Expr,
                   simd::Bool, threads::Bool)
    loop = :(for $index_name = 1:$axis_expression
                 $expr
             end)

    if threads
        return :(Threads.@threads $loop)
    elseif simd
        return :(@simd $loop)
    else
        return loop
    end
end

"""
    extractindices(expr) -> (array_names, index_names, axis_expressions)

Compute all `index_names` and respective axis calculations of an expression
involving the arrays with `array_names`. Everything is ordered by first
occurence in `expr`.

# Examples
```julia-repl
julia> extractindices(:(f(A[i,j,i]) + C[j]))
(Symbol[:A, :C], Symbol[:i, :j, :i, :j], Expr[:(size(A, 1)), :(size(A, 2)), :(size(A, 3)), :(size(C, 1))])
```
"""
extractindices(expr) = extractindices!(expr, Symbol[], Symbol[], Expr[])

function extractindices!(expr::Symbol,
                         array_names::Vector{Symbol},
                         index_names::Vector{Symbol},
                         axis_expressions::Vector{Expr})
    push!(array_names, expr)
    return array_names, index_names, axis_expressions
end

function extractindices!(expr::Number,
                         array_names::Vector{Symbol},
                         index_names::Vector{Symbol},
                         axis_expressions::Vector{Expr})
    return array_names, index_names, axis_expressions
end

function extractindices!(expr::Expr,
                         array_names::Vector{Symbol},
                         index_names::Vector{Symbol},
                         axis_expressions::Vector{Expr})
    if Meta.isexpr(expr, :ref) # e.g. A[i,j,k]
        array_name = expr.args[1]
        push!(array_names, array_name)

        # expr.args[2:end] are indices (e.g. [i,j,k])
        for (dimension, index_expr) in enumerate(expr.args[2:end])
            pushindex!(index_expr, array_name, dimension,
                       array_names, index_names, axis_expressions)
        end
    elseif Meta.isexpr(expr, :call)
        # e.g. 2*A[i,j], transpose(A[i,j]), or A[i] + B[j], so
        # expr.args[2:end] recursively contain the individual tensor
        # expressions (e.g. [A[i], B[j]])
        for arg in expr.args[2:end]
            extractindices!(arg, array_names, index_names, axis_expressions)
        end
    else
        throw(ArgumentError("Invalid expression head: `:$(expr.head)`"))
    end

    return return array_names, index_names, axis_expressions
end

function pushindex!(expr::Symbol, array_name::Symbol, dimension::Int,
                    array_names, index_names, axis_expressions)
    push!(index_names, expr)
    push!(axis_expressions, :(size($array_name, $dimension)))
end

function pushindex!(expr::Number, array_name::Symbol, dimension::Int,
                    array_names, index_names, axis_expressions)
    return nothing
end

function pushindex!(expr::Expr, array_name::Symbol, dimension::Int,
                    array_names, index_names, axis_expressions)
    # e.g. A[i+:offset] or A[i+5]
    #    expr is an Expr in this case
    #    We restrict it to be a Symbol (e.g. :i) followed by either
    #        a number or quoted expression.
    #    As before, push :i to index list
    #    Need to add/subtract off the offset to dimension list

    if Meta.isexpr(expr, :call) && length(expr.args) == 3
        op = expr.args[1]

        index_name = expr.args[2]
        @assert typeof(index_name) == Symbol

        offset_expr = expr.args[3]

        if offset_expr isa Integer
            offset = expr.args[3]::Integer
        elseif offset_expr isa Expr && Meta.isexpr(offset_expr, :quote)
            offset = offset_expr.args[1]
        elseif offset_expr isa QuoteNode
            offset = offset_expr.value
        else
            throw(ArgumentError("Improper expression inside reference on rhs"))
        end

        # push :i to indices we're iterating over
        push!(index_names, index_name)

        # need to invert + or - to determine iteration range
        if op == :+
            push!(axis_expressions, :(size($array_name, $dimension) - $offset))
        elseif op == :-
            push!(axis_expressions, :(size($array_name, $dimension) + $offset))
        else
            throw(ArgumentError(string("Operations inside ref on rhs are ",
                                       "limited to `+` and `-`")))
        end
    elseif Meta.isexpr(expr, :quote)
        return nothing
    else
        throw(ArgumentError("Invalid index expression: `$(expr)`"))
    end
end

function unquote_offsets!(expr::Expr, inside_ref = false)
    inside_ref |= Meta.isexpr(expr, :ref)

    for i in eachindex(expr.args)
        if expr.args[i] isa Expr
            if Meta.isexpr(expr.args[i], :quote) && inside_ref # never seems to get here
                expr.args[i] = expr.args[i].args[1]
            else
                unquote_offsets!(expr.args[i], inside_ref)
            end
        end
    end

    return expr
end

end # module


text:
pth = "C:\\Users\\p\\Desktop\\OpenStreetMapXSim.jl\\src\\osm"
path = "C:\\Users\\p\\Desktop\\OpenStreetMapXSim.jl\\src\\sim"
datapath = "C:\\Users\\p\\Desktop\\datasets\\"

datafile = "SAMPLE_WinnipegCMA_TRAFCAN2017Q1.csv"
sim_results = "counts.csv"
mapfile = "Winnipeg CMA.osm"

include(joinpath(path,"OpenStreetMapXSim.jl"))

using Main.OpenStreetMapXSim
using DataFrames
using Nanocsv
using GLM
using Plots; pyplot()
using StatPlots

function compare_data(datapath::String,mapfile::String,resultfile::String, testfile::String, counttype::Union{Nothing,String}, save::Bool = true)
    map_data = OpenStreetMapX.get_map_data(datapath,mapfile)
    frame = DataFrames.DataFrame(Node_ID = Int[],latitude = Float64[], longitude = Float64[], empirical = Int[],simulation = Int[])
    traffic = Nanocsv.read_csv(joinpath(datapath,testfile))
    if !isa(counttype,Nothing)
        traffic = traffic[traffic[:CNTTYPE1] .== counttype,:]
    end
    nodes = Dict(key => value for (key,value) in map_data.nodes if haskey(map_data.network.v, key))
    traffic_data = Dict()
    for i in 1:nrow(traffic)
        key = OpenStreetMapX.nearest_node(map_data.nodes,OpenStreetMapX.ENU(OpenStreetMapX.LLA(traffic[:LATITUDE][i],traffic[:LONGITUDE][i]),map_data.bounds),collect(keys(nodes)))
        if haskey(traffic_data, key)
            traffic_data[key][:TRAFFIC1][1] += traffic[i,:TRAFFIC1]
        else
            traffic_data[key] = traffic[i,[:LATITUDE,:LONGITUDE,:TRAFFIC1]]
        end
    end
    lines = readlines(open(joinpath(datapath,resultfile)))
    for line in lines
        x = parse.(Int,split(line))
        if haskey(traffic_data, x[1])
            push!(frame,(x[1],traffic_data[x[1]][:LATITUDE][1], traffic_data[x[1]][:LONGITUDE][1],traffic_data[x[1]][:TRAFFIC1][1], x[2]))
        end
    end
    model = GLM.lm(@formula(simulation ~ empirical - 1),frame)
    frame[:scaled_empirical] = GLM.predict(model)
    if !isa(counttype,Nothing)
        save && write_csv(joinpath(datapath,"traffic_data_comparision_$counttype.csv"),frame)
    else
        save && Nanocsv.write_csv(joinpath(datapath,"traffic_data_comparision.csv"),frame)
    end
    return frame, model
end

function get_comparision_dataset(datapath::String,mapfile::String,resultfile::String, testfile::String, save::Bool = true)
    counttypes= ["AADT", "ADT", "AWDT"]
    dataframe = DataFrames.DataFrame(Node_ID = Int[],latitude = Float64[], longitude = Float64[], empirical = Int[],simulation = Int[], scaled_empirical = Float64[], count_type = String[])
    for count in counttypes
        dframe, regression = compare_data(datapath,mapfile,resultfile,testfile, count,save)
        dframe[:count_type] = count
        append!(dataframe,dframe)
    end
    dataframe[:approximation_error] = abs.(dataframe[:scaled_empirical] .- dataframe[:simulation])./dataframe[:scaled_empirical]
    return dataframe
end

function get_grouped_hist(dataframe::DataFrames.DataFrame,nbars::Int, title::String)
    columns =  [:scaled_empirical,:simulation]
    min_val = minimum(vcat(dataframe[:scaled_empirical],dataframe[:simulation]))
    max_val = maximum(vcat(dataframe[:scaled_empirical],dataframe[:simulation]))
    tresholds = Int.(round.(range(min_val, stop = max_val, length = nbars+1))) |> collect
    categories = [UnitRange(tresholds[i],tresholds[i+1]) for i in 1:nbars]
    count_table = zeros(nbars,2)
    for i = 1:length(columns)
        for j = 1:nbars
            count_table[j,i] = sum(round(val) in categories[j] for val in dataframe[columns[i]])
        end
    end
    ctg = repeat(["Scaled Empirical","Simulation"], inner = nbars)
    nam = repeat([(minimum(c),maximum(c)) for c in categories], outer = length(columns))
    groupedbar(nam, count_table, group = ctg, xlabel = "Number of Cars at the Intersection", xrotation = 90.0, legendfont = Plots.font(18), ylabel = "Number of Intersections",
        title = title, bar_width = 0.67,
        lw = 0, framestyle = :box)

end

dataframe = get_comparision_dataset(datapath,mapfile,sim_results,datafile,false);

get_grouped_hist(dataframe,20, "Overall Traffic Count")
get_grouped_hist(dataframe[dataframe[:count_type] .== "AADT",:],20, "Traffic Count - AADT Only")
get_grouped_hist(dataframe[dataframe[:count_type] .== "ADT",:],20, "Traffic Count - ADT Only")
get_grouped_hist(dataframe[dataframe[:count_type] .== "AWDT",:],20, "Traffic Count - AWDT Only")

p = histogram(dataframe[dataframe[:approximation_error] .< 2.45,:][:approximation_error], nbins = 25, title = "Distribution of Approximation Errors After Removing Outlying Points\n - Overall Traffic", legend = :none)
p1 = histogram(dataframe[(dataframe[:approximation_error] .< 2.45) .& (dataframe[:count_type] .== "AWDT"),:][:approximation_error], nbins = 25, title = "Distribution of Approximation Errors After Removing Outlying Points\n - AWDT Points Only", legend = :none)
p2 = histogram(dataframe[(dataframe[:approximation_error] .< 2.45) .& (dataframe[:count_type] .== "AADT"),:][:approximation_error], nbins = 25, title = "Distribution of Approximation Errors After Removing Outlying Points\n - AADT Points Only", legend = :none)
p3 = histogram(dataframe[(dataframe[:approximation_error] .< 2.45) .& (dataframe[:count_type] .== "ADT"),:][:approximation_error], nbins = 25, title = "Distribution of Approximation Errors After Removing Outlying Points\n - ADT Points Only", legend = :none)
plot(p,p1,p2,p3)