text:
import Base: ==, length, in, eltype;

export Problem, InstrumentedProblem,
        actions, get_result, goal_test, path_cost, value,
        format_instrumented_results,
        Node, expand, child_node, solution, path, ==,
        search,
        GAState, mate, mutate,
        tree_search, graph_search,
        breadth_first_tree_search, depth_first_tree_search, depth_first_graph_search,
        breadth_first_search, best_first_graph_search, uniform_cost_search,
        recursive_dls, depth_limited_search, iterative_deepening_search,
        greedy_best_first_graph_search,
        Graph, make_undirected, connect_nodes, get_linked_nodes, get_nodes,
        UndirectedGraph, RandomGraph,
        GraphProblem,
        astar_search, recursive_best_first_search,
        hill_climbing, exp_schedule, simulated_annealing,
        or_search, and_search, and_or_graph_search,
        OnlineDFSAgentProgram, update_state, execute,
        OnlineSearchProblem, LRTAStarAgentProgram,
        learning_realtime_astar_cost,
        genetic_search, genetic_algorithm,
        NQueensProblem, conflict, conflicted,
        random_boggle, print_boggle, boggle_neighbors, int_sqrt,
        WordList, lookup, length, in,
        BoggleFinder, set_board, find, words, score,
        boggle_hill_climbing, mutate_boggle,
        execute_searcher, compare_searchers, beautify_node;

#=

    Problem is a abstract problem that contains a initial state and goal state.

=#
mutable struct Problem <: AbstractProblem
    initial::String
    goal::Union{Nothing, String}

    function Problem(initial_state::String; goal_state::Union{Nothing, String}=nothing)
        return new(initial_state, goal_state);
    end
end

"""
    actions(ap::T, state::String) where {T <: AbstractProblem}

Return an array of possible actions that can be executed in the given state 'state'.
"""
function actions(ap::T, state::String) where {T <: AbstractProblem}
    println("actions() is not implemented yet for ", typeof(ap), "!");
    nothing;
end

"""
    get_result(ap::T, state::String, action::String) where {T <: AbstractProblem}

Return the resulting state from executing the given action 'action' in the given state 'state'.
"""
function get_result(ap::T, state::String, action::String) where {T <: AbstractProblem}
    println("get_result() is not implemented yet for ", typeof(ap), "!");
    nothing;
end

"""
    goal_test(ap::T, state::String) where {T <: AbstractProblem}

Return a boolean value representing whether the given state 'state' is a goal state in the given
problem 'ap'.
"""
function goal_test(ap::T, state::String) where {T <: AbstractProblem}
    return ap.goal == state;
end

"""
    path_cost(ap::T, cost::Float64, state1::String, action::String, state2::String) where {T <: AbstractProblem}
    path_cost(ap::T, cost::Float64, state1::AbstractVector, action::Int64, state2::AbstractVector) where {T <: AbstractProblem}

Return the cost of a solution path arriving at 'state2' from 'state1' with the given action 'action' and
cost 'cost' to arrive at 'state1'. The default path_cost() method costs 1 for every step in a path.
"""
function path_cost(ap::T, cost::Float64, state1::String, action::String, state2::String) where {T <: AbstractProblem}
    return cost + 1;
end

function path_cost(ap::T, cost::Float64, state1::AbstractVector, action::Int64, state2::AbstractVector) where {T <: AbstractProblem}
    return cost + 1;
end

"""
    value(ap::T, state::String) where {T <: AbstractProblem}

Return a value for the given state 'state' in the given problem 'ap'.

This value is used in optimization problems such as hill climbing or simulated annealing.
"""
function value(ap::T, state::String) where {T <: AbstractProblem}
    println("value() is not implemented yet for ", typeof(ap), "!");
    nothing;
end

#=

    InstrumentedProblem is a AbstractProblem implementation that wraps another AbstractProblem

    implementation and tracks the number of function calls made. This problem is used in

    compare_searchers() and execute_searcher().

=#
mutable struct InstrumentedProblem <: AbstractProblem
    problem::AbstractProblem
    actions::Int64
    results::Int64
    goal_tests::Int64
    found   # can be any DataType, but check for Nothing DataType later

    function InstrumentedProblem(ap::T) where {T <: AbstractProblem}
        return new(ap, Int64(0), Int64(0), Int64(0), nothing);
    end
end

function actions(ap::InstrumentedProblem, state::AbstractVector)
    ap.actions = ap.actions + 1;
    return actions(ap.problem, state);
end

function actions(ap::InstrumentedProblem, state::String)
    ap.actions = ap.actions + 1;
    return actions(ap.problem, state);
end

function get_result(ap::InstrumentedProblem, state::String, action::String)
    ap.results = ap.results + 1;
    return get_result(ap.problem, state, action);
end

function get_result(ap::InstrumentedProblem, state::AbstractVector, action::Int64)
    ap.results = ap.results + 1;
    return get_result(ap.problem, state, action);
end

function goal_test(ap::InstrumentedProblem, state::String)
    ap.goal_tests = ap.goal_tests + 1;
    local result::Bool = goal_test(ap.problem, state);
    if (result)
        ap.found = state;
    end
    return result;
end

function goal_test(ap::InstrumentedProblem, state::AbstractVector)
    ap.goal_tests = ap.goal_tests + 1;
    local result::Bool = goal_test(ap.problem, state);
    if (result)
        ap.found = state;
    end
    return result;
end

function path_cost(ap::InstrumentedProblem, cost::Float64, state1::String, action::String, state2::String)
    return path_cost(ap.problem, cost, state1, action, state2);
end

function path_cost(ap::InstrumentedProblem, cost::Float64, state1::AbstractVector, action::Int64, state2::AbstractVector)
    return path_cost(ap.problem, cost, state1, action, state2);
end

function value(ap::InstrumentedProblem, state::String)
    return value(ap.problem, state);
end

function value(ap::InstrumentedProblem, state::AbstractVector)
    return value(ap.problem, state);
end

function format_instrumented_results(ap::InstrumentedProblem)
    return @sprintf("<%4d/%4d/%4d/%s>", ap.actions, ap.goal_tests, ap.results, string(ap.found));
end

# A node should not exist without a state.
mutable struct Node{T}
    state::T
    path_cost::Float64
    depth::UInt32
    action::Union{Nothing, String, Int64, Tuple}
    parent::Union{Nothing, Node}
    f::Float64

    function Node{T}(state::T; parent::Union{Nothing, Node}=nothing, action::Union{Nothing, String, Int64, Tuple}=nothing, path_cost::Float64=0.0, f::Union{Nothing, Float64}=nothing) where T
        nn = new(state, path_cost, UInt32(0), action, parent);
        if (typeof(parent) <: Node)
            nn.depth = UInt32(parent.depth + 1);
        end
        if (typeof(f) <: Float64)
            nn.f = f;
        end
        return nn;
    end
end

"""
    expand(n::Node, ap::T) where {T <: AbstractProblem}

Return an array of nodes reachable by 1 step from the given node 'n' in the problem 'ap'.
"""
function expand(n::Node, ap::T) where {T <: AbstractProblem}
    return collect(child_node(n, ap, act) for act in actions(ap, n.state));
end

"""
    child_node(n::Node, ap::T, action::String) where {T <: AbstractProblem}

Return a child node for the given node 'n' in problem 'ap' after executing the action 'action' (Fig. 3.10).
"""
function child_node(n::Node, ap::T, action::String) where {T <: AbstractProblem}
    local next_node = get_result(ap, n.state, action);
    return Node{typeof(next_node)}(next_node, parent=n, action=action, path_cost=path_cost(ap, n.path_cost, n.state, action, next_node));
end

function child_node(n::Node, ap::T, action::Int64) where {T <: AbstractProblem}
    local next_node = get_result(ap, n.state, action);
    return Node{typeof(next_node)}(next_node, parent=n, action=action, path_cost=path_cost(ap, n.path_cost, n.state, action, next_node));
end

function child_node(n::Node, ap::T, action::Tuple) where {T <: AbstractProblem}
    local next_node = get_result(ap, n.state, action);
    return Node{typeof(next_node)}(next_node, parent=n, action=action, path_cost=path_cost(ap, n.path_cost, n.state, action, next_node));
end

"""
    solution(n::Node)

Return an array of actions to get from the root node of node 'n' to the given node 'n'.
"""
function solution(n::Node)
    local path_sequence = path(n);
    return [node.action for node in path_sequence[2:length(path_sequence)]];
end

"""
    path(n::Node)

Return the path between the root node of node 'n' to the given node 'n' as an array of nodes.
"""
function path(n::Node)
    local node = n;
    local path_back = [];
    while true
        push!(path_back, node);
        if (!(node.parent === nothing))
            node = node.parent;
        else
            # The root node does not have a parent node.
            break;
        end
    end
    path_back = reverse(path_back);
    return path_back;
end

function ==(n1::Node, n2::Node)
    return (n1.state == n2.state);
end

#=

    SimpleProblemSolvingAgentProgram is a abstract problem solving agent (Fig. 3.1).

=#
mutable struct SimpleProblemSolvingAgentProgram <: AgentProgram
    state::Union{Nothing, String}
    goal::Union{Nothing, String}
    seq::Array{String, 1}
    problem::Union{Nothing, Problem}

    function SimpleProblemSolvingAgentProgram(;initial_state::Union{Nothing, String}=nothing)
        return new(initial_state, nothing, Array{String, 1}(), nothing);
    end
end

function execute(spsap::SimpleProblemSolvingAgentProgram, percept::Tuple{Any, Any})
    spsap.state = update_state(spsap, spsap.state, percept);
    if (length(spsap.seq) == 0)
        spsap.goal = formulate_problem(spsap, spsap.state);
        spsap.problem = forumate_problem(spsap, spsap.state, spsap.goal);
        spsap.seq = search(spsap, spsap.problem);
        if (length(spsap.seq) == 0)
            return Nothing;
        end
    end
    local action = popfirst!(spsap.seq);
    return action;
end

function update_state(spsap::SimpleProblemSolvingAgentProgram, state::String, percept::Tuple{Any, Any})
    println("update_state() is not implemented yet for ", typeof(spsap), "!");
    nothing;
end

function formulate_goal(spsap::SimpleProblemSolvingAgentProgram, state::String)
    println("formulate_goal() is not implemented yet for ", typeof(spsap), "!");
    nothing;
end

function formulate_problem(spsap::SimpleProblemSolvingAgentProgram, state::String, goal::String)
    println("formulate_problem() is not implemented yet for ", typeof(spsap), "!");
    nothing;
end

function search(spsap::SimpleProblemSolvingAgentProgram, problem::T) where {T <: AbstractProblem}
    println("search() is not implemented yet for ", typeof(spsap), "!");
    nothing;
end

struct GAState
    genes::Array{Any, 1}

    function GAState(genes::Array{Any, 1})
        return new(Array{Any,1}(deepcopy(genes)));
    end
end

function mate(ga_state::T, other::T) where {T <: GAState}
    local c = rand(RandomDeviceInstance, range(1, stop=length(ga_state.genes)));
    local new_ga_state = deepcopy(ga_state[1:c]);
    for element in other.genes[(c + 1):length(other.genes)]
        push!(new_ga_state, element);
    end
    return new_ga_state;
end

function mutate(ga_state::T) where {T <: GAState}
    println("mutate() is not implemented yet for ", typeof(ga_state), "!");
    nothing;
end

"""
    tree_search{T1 <: AbstractProblem, T2 <: Queue}(problem::T1, frontier::T2)

Search the given problem by using the general tree search algorithm (Fig. 3.7) and return the node solution.
"""
function tree_search(problem::T1, frontier::T2) where {T1 <: AbstractProblem, T2 <: Queue}
    push!(frontier, Node{typeof(problem.initial)}(problem.initial));
    while (length(frontier) != 0)
        local node = pop!(frontier);
        if (goal_test(problem, node.state))
            return node;
        end
        extend!(frontier, expand(node, problem));
    end
    return nothing;
end

function tree_search(problem::InstrumentedProblem, frontier::T) where {T <: Queue}
    push!(frontier, Node{typeof(problem.problem.initial)}(problem.problem.initial));
    while (length(frontier) != 0)
        local node = pop!(frontier);
        if (goal_test(problem, node.state))
            return node;
        end
        extend!(frontier, expand(node, problem));
    end
    return nothing;
end

"""
    graph_search{T1 <: AbstractProblem, T2 <: Queue}(problem::T1, frontier::T2)

Search the given problem by using the general graph search algorithm (Fig. 3.7) and return the node solution.

The uniform cost algorithm (Fig. 3.14) should be used when the frontier is a priority queue.
"""
function graph_search(problem::T1, frontier::T2) where {T1 <: AbstractProblem, T2 <: Queue}
    local explored::Set;
    if (typeof(problem.initial) <: Tuple)
        explored = Set{NTuple}();
    else
        explored = Set{typeof(problem.initial)}();
    end
    push!(frontier, Node{typeof(problem.initial)}(problem.initial));
    while (length(frontier) != 0)
        local node = pop!(frontier);
        if (goal_test(problem, node.state))
            return node;
        end
        push!(explored, node.state);
        extend!(frontier, collect(child_node for child_node in expand(node, problem)
                                if (!(child_node.state in explored) && !(child_node in frontier))));
    end
    return nothing;
end

function graph_search(problem::InstrumentedProblem, frontier::T) where {T <: Queue}
    local explored::Set;
    if (typeof(problem.problem.initial) <: Tuple)
        explored = Set{NTuple}();
    else
        explored = Set{typeof(problem.problem.initial)}();
    end
    push!(frontier, Node{typeof(problem.problem.initial)}(problem.problem.initial));
    while (length(frontier) != 0)
        local node = pop!(frontier);
        if (goal_test(problem, node.state))
            return node;
        end
        push!(explored, node.state);
        extend!(frontier, collect(child_node for child_node in expand(node, problem)
                                if (!(child_node.state in explored) && !(child_node in frontier))));
    end
    return nothing;
end

"""
    breadth_first_tree_search(problem::T) where {T <: AbstractProblem}

Search the shallowest nodes in the search tree first.
"""
function breadth_first_tree_search(problem::T) where {T <: AbstractProblem}
    return tree_search(problem, FIFOQueue());
end

"""
    depth_first_tree_search(problem::T) where {T <: AbstractProblem}

Search the deepest nodes in the search tree first.
"""
function depth_first_tree_search(problem::T) where {T <: AbstractProblem}
    return tree_search(problem, Stack());
end

"""
    depth_first_graph_search(problem::T) where {T <: AbstractProblem}

Search the deepest nodes in the search tree first.
"""
function depth_first_graph_search(problem::T) where {T <: AbstractProblem}
    return graph_search(problem, Stack());
end

"""
    breadth_first_search(problem::T) where {T <: AbstractProblem}
    breadth_first_search(problem::InstrumentedProblem)

Return a solution by using the breadth-first search algorithm (Fig. 3.11)
on the given problem 'problem'. Otherwise, return 'nothing' on failure.
"""
function breadth_first_search(problem::T) where {T <: AbstractProblem}
    local node = Node{typeof(problem.initial)}(problem.initial);
    if (goal_test(problem, node.state))
        return node;
    end
    local frontier = FIFOQueue();
    push!(frontier, node);
    local explored = Set{String}();
    while (length(frontier) != 0)
        node = pop!(frontier);
        push!(explored, node.state);
        for child_node in expand(node, problem)
            if (!(child_node.state in explored) && !(child_node in frontier))
                if (goal_test(problem, child_node.state))
                    return child_node;
                end
                push!(frontier, child_node);
            end
        end
    end
    return nothing;
end

function breadth_first_search(problem::InstrumentedProblem)
    local node = Node{typeof(problem.problem.initial)}(problem.problem.initial);
    if (goal_test(problem, node.state))
        return node;
    end
    local frontier = FIFOQueue();
    push!(frontier, node);
    local explored = Set{String}();
    while (length(frontier) != 0)
        node = pop!(frontier);
        push!(explored, node.state);
        for child_node in expand(node, problem)
            if (!(child_node.state in explored) && !(child_node in frontier))
                if (goal_test(problem, child_node.state))
                    return child_node;
                end
                push!(frontier, child_node);
            end
        end
    end
    return nothing;
end

"""
    best_first_graph_search(problem::T, f::Function) where {T <: AbstractProblem}

Search the nodes in the given problem 'problem' by visiting the nodes with the lowest
scores returned by f(). If f() is a heuristics estimate function to the goal state, then
this function becomes greedy best first search. If f() is a function that gets the node's
depth, then this function becomes breadth-first search.

Returns a solution if found, otherwise returns 'nothing' on failure.

This function uses f as a Function, because using f as an MemoizedFunction exhibits unusual
behavior when relying on MemoizedFunction by producing unexpected results.
"""
function best_first_graph_search(problem::T, f::Function) where {T <: AbstractProblem}
    local node = Node{typeof(problem.initial)}(problem.initial);
    if (goal_test(problem, node.state))
        return node;
    end
    local frontier = PQueue();
    push!(frontier, node, f);
    local explored = Set{typeof(problem.initial)}();
    while (length(frontier) != 0)
        node = pop!(frontier);
        if (goal_test(problem, node.state))
            return node;
        end
        push!(explored, node.state);
        for child_node in expand(node, problem)
            if (!(child_node.state in explored) &&
                !(child_node in collect(getindex(x, 2) for x in frontier.array)))
                push!(frontier, child_node, f);
            elseif (child_node in [getindex(x, 2) for x in frontier.array])
                # Recall that Nodes can share the same state and different values for other fields.
                local existing_node = pop!(collect(getindex(x, 2)
                                                    for x in frontier.array
                                                    if (getindex(x, 2) == child_node)));
                if (f(child_node) < f(existing_node))
                    delete!(frontier, existing_node);
                    push!(frontier, child_node, f);
                end
            end
        end
    end
    return nothing;
end

"""
    uniform_cost_search(problem::T) where {T <: AbstractProblem}

Search the given problem by using the uniform cost algorithm (Fig. 3.14) and return the node solution.

solution() can be used on the node solution to reconstruct the path taken to the solution.
"""
function uniform_cost_search(problem::T) where {T <: AbstractProblem}
    return best_first_graph_search(problem, (function(n::Node)return n.path_cost;end));
end

function recursive_dls(node::Node, problem::T, limit::Int64) where {T <: AbstractProblem}
    if (goal_test(problem, node.state))
        return node;
    elseif (node.depth == limit)
        return "cutoff";
    else
        local cutoff_occurred = false;
        for child_node in expand(node, problem)
            local result = recursive_dls(child_node, problem, limit);
            if (result == "cutoff")
                cutoff_occurred = true;
            elseif (!(typeof(result) <: Nothing))
                return result;
            end
        end
        return if_(cutoff_occurred, "cutoff", nothing);
    end
end;

"""
    depth_limited_search(problem::T; limit::Int64) where {T <: AbstractProblem}

Search the given problem by using the depth limited tree search algorithm (Fig. 3.17)
and return the node solution if a solution was found. Otherwise, this function returns 'nothing'.

solution() can be used on the node solution to reconstruct the path taken to the solution.
"""
function depth_limited_search(problem::T; limit::Int64=50) where {T <: AbstractProblem}
    return recursive_dls(Node{typeof(problem.initial)}(problem.initial), problem, limit);
end

function depth_limited_search(problem::InstrumentedProblem; limit::Int64=50)
    return recursive_dls(Node{typeof(problem.problem.initial)}(problem.problem.initial), problem, limit);
end

"""
    iterative_deepening_search(problem::T) where {T <: AbstractProblem}

Search the given problem by using the iterative deepening search algorithm (Fig. 3.18)
and return the node solution if a solution was found. Otherwise, this function returns 'nothing'.

solution() can be used on the node solution to reconstruct the path taken to the solution.
"""
function iterative_deepening_search(problem::T) where {T <: AbstractProblem}
    for depth in 1:typemax(Int64)
        local result = depth_limited_search(problem, limit=depth)
        if (result != "cutoff")
            return result;
        end
    end
    return nothing;
end

const greedy_best_first_graph_search = best_first_graph_search;

#=

    Graph is a graph that consists of nodes (vertices) and edges (links).

    The Graph constructor uses the keyword 'directed' to specify if the graph

    is directed or undirected.

    For an example Graph instance:

        Graph(dict=Dict([("A", Dict([("B", 1), ("C", 2)]))]))

    The example Graph is a directed graph with 3 vertices ("A", "B", and "C")

    and link "A"=>"B" (length 1) and "A"=>"C" (length 2).

=#
struct Graph{N}
    dict::Dict{N, Any}
    locations::Dict{N, Tuple{Any, Any}}
    directed::Bool

    function Graph{N}(;dict::Union{Nothing, Dict{N, }}=nothing, locations::Union{Nothing, Dict{N, Tuple{Any, Any}}}=nothing, directed::Bool=true) where N
        local ng::Graph;
        if ((typeof(dict) <: Nothing) && (typeof(locations) <: Nothing))
            ng = new(Dict{Any, Any}(), Dict{Any, Tuple{Any, Any}}(), Bool(directed));
        elseif (typeof(locations) <: Nothing)
            ng = new(Dict{eltype(dict.keys), Any}(dict), Dict{Any, Tuple{Any, Any}}(), Bool(directed));
        else
            ng = new(Dict{eltype(dict.keys), Any}(dict), Dict{eltype(locations.keys), Tuple{Any, Any}}(locations), Bool(directed));
        end
        if (!ng.directed)
            make_undirected(ng);
        end
        return ng;
    end

    function Graph{N}(graph::Graph{N}) where N
        return new(Dict{Any, Any}(graph.dict), Dict{String, Tuple{Any, Any}}(graph.locations), Bool(graph.directed));
    end
end

eltype(::Type{<:Graph{T}}) where {T} = T

function make_undirected(graph::Graph)
    for location_A in keys(graph.dict)
        for (location_B, d) in graph.dict[location_A]
            connect_nodes(graph, location_B, location_A, distance=d);
        end
    end
end

"""
    connect_nodes(graph::Graph{N}, A::N, B::N; distance::Int64=Int64(1)) where N

Add a link between Node 'A' to Node 'B'. If the graph is undirected, then add
the inverse link from Node 'B' to Node 'A'.
"""
function connect_nodes(graph::Graph{N}, A::N, B::N; distance::Int64=Int64(1)) where N
    get!(graph.dict, A, Dict{String, Int64}())[B]=distance;
    if (!graph.directed)
        get!(graph.dict, B, Dict{String, Int64}())[A]=distance;
    end
    nothing;
end

"""
    get_linked_nodes(graph::Graph{N}, a::N; b::Union{Nothing, N}=nothing) where N

Return a dictionary of nodes and their distances if the 'b' keyword is not given.
Otherwise, return the distance between 'a' and 'b'.
"""
function get_linked_nodes(graph::Graph{N}, a::N; b::Union{Nothing, N}=nothing) where N
    local linked = get!(graph.dict, a, Dict{Any, Any}());
    if (typeof(b) <: Nothing)
        return linked;
    else
        return get(linked, b, nothing);
    end
end

function get_nodes(graph::Graph)
    return collect(keys(graph.dict));
end

"""
    UndirectedGraph(dict::Dict{T, }, locations::Dict{T, Tuple{Any, Any}}) where T
    UndirectedGraph()

Return an undirected graph from the given dictionary of links 'dict' and dictionary
of locations 'locations' if given.
"""
function UndirectedGraph(dict::Dict{T, }, locations::Dict{T, Tuple{Any, Any}}) where T
    return Graph{eltype(dict.keys)}(dict=dict, locations=locations, directed=false);
end

function UndirectedGraph()
    return Graph{Any}(directed=false);
end

"""
    RandomGraph()

Return a random graph with the specified nodes and number of links.
"""
function RandomGraph(;nodes::UnitRange=1:10,
                    min_links::Int64=2,
                    width::Int64=400,
                    height::Int64=300,
                    curvature::Function=(function()
                                            return (0.4*rand(RandomDeviceInstance)) + 1.1;
                                        end))
    local g = UndirectedGraph();
    for node in nodes
        g.locations[node] = Tuple((rand(RandomDeviceInstance, 1:width), rand(RandomDeviceInstance, 1:height)));
    end
    for i in 1:min_link
        for node in nodes
            if (get_linked_nodes(g, node) < min_links)
                local here = g.locations[node];
                local neighbor = argmin(nodes, (function(n, ; graph::Graph=g, current_node::Node=node, current_location::Tuple=here)
                                                    if (n == current_node || get_linked_nodes(graph, current_node, n) != nothing)
                                                        return Inf;
                                                    end
                                                    return distance(g.locations[n], current_location);
                                                end));
                local d = distance(g.locations[neighbor], here) * curvature();
                connect(g, node, neighbor, Int64(floor(d)));
            end
        end
    end
    return g;
end

#=

    GraphProblem is the problem of searching a graph from one node to another node.

=#
struct GraphProblem <: AbstractProblem
    initial::String
    goal::String
    graph::Graph
    h::MemoizedFunction

    function GraphProblem(initial_state::String, goal_state::String, graph::Graph)
        return new(initial_state, goal_state, Graph{eltype(graph)}(graph), MemoizedFunction(initial_to_goal_distance));
    end
end

function actions(gp::GraphProblem, loc::String)
    return collect(keys(get_linked_nodes(gp.graph,loc)));
end

function get_result(gp::GraphProblem, state::String, action::String)
    return action;
end

function path_cost(gp::GraphProblem, current_cost::Float64, location_A::String, action::String, location_B::String)
    local AB_distance::Float64;
    if (haskey(gp.graph.dict, location_A) && haskey(gp.graph.dict[location_A], location_B))
        AB_distance= Float64(get_linked_nodes(gp.graph,location_A, b=location_B));
    else
        AB_distance = Float64(Inf);
    end
    return current_cost + AB_distance;
end

"""
    initial_to_goal_distance(gp::GraphProblem, n::Node)

Compute the straight line distance between the initial state and goal state.
"""
function initial_to_goal_distance(gp::GraphProblem, n::Node)
    local locations = gp.graph.locations;
    if (isempty(locations))
        return Inf;
    else
        return Float64(floor(distance(locations[n.state], locations[gp.goal])));
    end
end

function initial_to_goal_distance(gp::InstrumentedProblem, n::Node)
    local locations = gp.problem.graph.locations;
    if (isempty(locations))
        return Inf;
    else
        return Float64(floor(distance(locations[n.state], locations[gp.problem.goal])));
    end
end

"""
    astar_search(problem::GraphProblem; h::Union{Nothing, Function}=nothing)

Apply the A* search (best-first graph search with f(n)=g(n)+h(n)) to the given problem 'problem'.
If the 'h' keyword is not used, this function uses the function problem.h.

This function uses mh as a Function, because using mh as an MemoizedFunction exhibits unusual
behavior when relying on MemoizedFunction by producing unexpected results.
"""
function astar_search(problem::GraphProblem; h::Union{Nothing, Function}=nothing)
    local mh::Function;
    if (!(typeof(h) <: Nothing))
        mh = h;
    else
        mh = problem.h.f;
    end
    return best_first_graph_search(problem,
                                    (function(node::Node; h::Function=mh, prob::GraphProblem=problem)
                                        return node.path_cost + h(prob, node);
                                    end));
end

"""
    RBFS(problem::T1, node::T2, flmt::Float64, h::MemoizedFunction) where {T1 <: AbstractProblem, T2 <: Node}

Recursively calls RBFS() with a new 'flmt' value and returns its solution to recursive_best_first_search().
"""
function RBFS(problem::T1, node::T2, flmt::Float64, h::MemoizedFunction) where {T1 <: AbstractProblem, T2 <: Node}
    if (goal_test(problem, node.state))
        return node, 0.0;
    end
    local successors = expand(node, problem);
    if (length(successors) == 0);
        return node, Inf;
    end
    for successor in successors
        successor.f = max(successor.path_cost + eval_memoized_function(h, problem, successor), node.f);
    end
    while (true)
        sort!(successors, lt=(function(n1::Node, n2::Node)return isless(n1.f, n2.f);end));
        local best::Node = successors[1];
        if (best.f > flmt)
            return nothing, best.f;
        end
        local alternative::Float64;
        if (length(successors) > 1)
            alternative = successors[1].f;
        else
            alternative = Inf;
        end
        result, best.f = RBFS(problem, best, min(flmt, alternative), h);
        if (!(result === nothing))
            return result, best.f;
        end
    end
end

"""
    recursive_best_first_search(problem::T; h::Union{Nothing, MemoizedFunction}) where {T <: AbstractProblem}

Search the given problem by using the recursive best first search algorithm (Fig. 3.26)
and return the node solution.

solution() can be used on the node solution to reconstruct the path taken to the solution.
"""
function recursive_best_first_search(problem::T; h::Union{Nothing, MemoizedFunction}=nothing) where {T <: AbstractProblem}
    local mh::MemoizedFunction; #memoized h(n) function
    if (!(typeof(h) <: Nothing))
        mh = MemoizedFunction(h);
    else
        mh = problem.h;
    end

    local node = Node{typeof(problem.initial)}(problem.initial);
    node.f = eval_memoized_function(mh, problem, node);
    result, bestf = RBFS(problem, node, Inf, mh);
    return result;
end

function recursive_best_first_search(problem::InstrumentedProblem; h::Union{Nothing, MemoizedFunction}=nothing)
    local mh::MemoizedFunction; #memoized h(n) function
    if (!(typeof(h) <: Nothing))
        mh = MemoizedFunction(h);
    else
        mh = problem.problem.h;
    end

    local node = Node{typeof(problem.problem.initial)}(problem.problem.initial);
    node.f = eval_memoized_function(mh, problem, node);
    result, bestf = RBFS(problem, node, Inf, mh);
    return result;
end

"""
    hill_climbing(problem::T) where {T <: AbstractProblem}

Return a state that is a local maximum for the given problem 'problem' by using
the hill-climbing search algorithm (Fig. 4.2) on the initial state of the problem.
"""
function hill_climbing(problem::T) where {T <: AbstractProblem}
    local current_node = Node{typeof(problem.initial)}(problem.initial);
    while (true)
        local neighbors = expand(current_node, problem);
        if (length(neighbors) == 0)
            break;
        end
        local neighbor = argmax_random_tie(neighbors,
                                            (function(n::Node,; p::AbstractProblem=problem)
                                                return value(p, n.state);
                                            end));
        if (value(problem, neighbor.state) <= value(problem, current_node.state))
            break;
        end
        current_node = neighbor;
    end
    return current_node.state;
end

"""
    exp_schedule(;kvar::Int64=20, delta::Float64=0.005, lmt::Int64=100)

Return a scheduled time for simulated annealing.
"""
function exp_schedule(;kvar::Int64=20, delta::Float64=0.005, lmt::Int64=100)
    return (function(t::Real; k=kvar, d=delta, limit=lmt)
                return if_((t < limit), (k * exp(-d * t)), 0);
            end);
end

"""
    simulated_annealing(problem::T; schedule::Function=exp_schedule()) where {T <: AbstractProblem}

Return the solution node by applying the simulated annealing algorithm (Fig. 4.5) on the given
problem 'problem' and schedule function 'schedule'. If a solution node can't be found,
this function returns 'nothing' on failure.
"""
function simulated_annealing(problem::T; schedule::Function=exp_schedule()) where {T <: AbstractProblem}
    local current_node = Node{typeof(problem.initial)}(problem.initial);
    for t in 0:(typemax(Int64) - 1)
        local temperature::Float64 = schedule(t);
        if (temperature == 0)
            return current_node;
        end
        local neighbors = expand(current_node, problem);
        if (length(neighbors) == 0)
            return current_node;
        end
        local next_node = rand(RandomDeviceInstance, neighbors);
        delta_e = value(problem, next_node.state) - value(problem, current_node.state);
        if ((delta_e > 0) || (exp(delta_e/temperature) > rand(RandomDeviceInstance)))
            current_node = next_node;
        end
    end
    return nothing;
end

#=

    and_search() and or_search() are used by and_or_graph_search().

=#
function or_search(problem::T, state::AbstractVector, path::AbstractVector) where {T <: AbstractProblem}
    if (goal_test(problem, state))
        return [];
    end
    if (state in path)
        return nothing;
    end
    for action in actions(problem, state)
        local plan = and_search(get_result(problem, state, action), vcat(path, [state,]));
        if (plan != nothing)
            return [action, plan];
        end
    end
    return nothing;
end

function or_search(problem::T, state::String, path::AbstractVector) where {T <: AbstractProblem}
    if (goal_test(problem, state))
        return [];
    end
    if (state in path)
        return nothing;
    end
    for action in actions(problem, state)
        local plan = and_search(problem, get_result(problem, state, action), vcat(path, [state,]));
        if (plan != nothing)
            return [action, plan];
        end
    end
    return nothing;
end

function and_search(problem::T, states::AbstractVector, path::AbstractVector) where {T <: AbstractVector}
    local plan = Dict{Any, Any}();
    for state in states
        plan[state] = or_search(problem, state, path);
        if (plan[state] == nothing)
            return nothing;
        end
    end
    return plan;
end

"""
    and_or_graph_search(problem::T) where {T <: AbstractProblem}

Return a conditional plan by using the algorithm for searching and-or graphs (Fig. 4.11)
on the given problem 'problem'. This function returns 'nothing' on failure.
"""
function and_or_graph_search(problem::T) where {T <: AbstractProblem}
    return or_search(problem, problem.initial, []);
end

#=

    OnlineDFSAgentProgram is a online depth first search agent (Fig. 4.21)

    implementation of AgentProgram.

=#
mutable struct OnlineDFSAgentProgram <: AgentProgram
    result::Dict
    untried::Dict
    unbacktracked::Dict
    state::Union{Nothing, String}
    action::Union{Nothing, String}
    problem::AbstractProblem

    function OnlineDFSAgentProgram(problem::T) where {T <: AbstractProblem}
        return new(Dict(), Dict(), Dict(), nothing, nothing, problem);
    end
end

function update_state(odfsap::OnlineDFSAgentProgram, percept::String)
    return percept;
end

function execute(odfsap::OnlineDFSAgentProgram, percept::String)
    local s_prime::String = update_state(odfsap, percept);
    if (goal_test(odfsap.problem, s_prime))
        odfsap.action = nothing;
    else
        if (!(s_prime in keys(odfsap.untried)))
            odfsap.untried[s_prime] = actions(odfsap.problem, s_prime);
        end
        if (!(odfsap.state === nothing))
            if (haskey(odfsap.result, (odfsap.state, odfsap.action)))
                if (s_prime != odfsap.result[(odfsap.state, odfsap.action)])
                    odfsap.result[(odfsap.state, odfsap.action)] = s_prime;
                    pushfirst!(odfsap.unbacktracked[s_prime], odfsap.state);
                end
            else
                if (s_prime != [])
                    odfsap.result[(odfsap.state, odfsap.action)] = s_prime;
                    pushfirst!(odfsap.unbacktracked[s_prime], odfsap.state);
                end
            end
        end
        if (length(odfsap.untried[s_prime]) == 0)
            if (length(odfsap.unbacktracked[s_prime]) == 0)
                odfsap.action = nothing;
            else
                first_item = popfirst!(odfsap.unbacktracked[s_prime]);
                for (state, b) in keys(odfsap.result)
                    if (odfsap.result[(state, b)] == first_item)
                        odfsap.action = b;
                        break;
                    end
                end
            end
        else
            odfsap.action = popfirst!(odfsap.untried[s_prime]);
        end
    end
    odfsap.state = s_prime;
    return odfsap.action;
end

#=

    OnlineSearchProblem is a AbstractProblem implementation of a online search problem

    that can be solved by a online search agent.

=#
struct OnlineSearchProblem <: AbstractProblem
    initial::String
    goal::String
    graph::Graph
    least_costs::Dict
    h::Function

    function OnlineSearchProblem(initial::String, goal::String, graph::Graph, least_costs::Dict)
        return new(initial, goal, graph, least_costs, online_search_least_cost);
    end
end

function actions(osp::OnlineSearchProblem, state::String)
    return collect(keys(osp.graph.dict[state]));
end

function get_result(osp::OnlineSearchProblem, state::String, action::String)
    return osp.graph.dict[state][action];
end

function online_search_least_cost(osp::OnlineSearchProblem, state::String)
    return osp.least_costs[state];
end

function path_cost(osp::OnlineSearchProblem, state1::String, action::String, state2::String)
    return 1;
end

function goal_test(osp::OnlineSearchProblem, state::String)
    if (state == osp.goal)
        return true;
    else
        return false;
    end
end

#=

    LRTAStarAgentProgram is an AgentProgram implementation of LRTA*-Agent (Fig. 4.24).

    The 'result' field is not necessary as the given problem contains the results table.

=#
mutable struct LRTAStarAgentProgram <: AgentProgram
    H::Dict
    state::Union{Nothing, String}
    action::Union{Nothing, String}
    problem::AbstractProblem

    function LRTAStarAgentProgram(problem::T) where {T <: AbstractProblem}
        return new(Dict(), nothing, nothing, problem);
    end
end

function learning_realtime_astar_cost(lrtaap::LRTAStarAgentProgram, state::String, action::String, s_prime::String, H::Dict)
    if (haskey(lrtaap.H, s_prime))
        return path_cost(lrtaap.problem, state, action, s_prime) + lrtaap.H[s_prime];
    else
        return path_cost(lrtaap.problem, state, action, s_prime) + lrtaap.problem.h(lrtaap.problem, s_prime);
    end
end

"""
    execute(lrtaap::LRTAStarAgentProgram, s_prime::String)

Return an action given the percept 's_prime' and by using the LRTA*-Agent
program (Fig. 4.24). If the current state of the agent is at the goal
state, return 'nothing'.
"""
function execute(lrtaap::LRTAStarAgentProgram, s_prime::String)
    if (goal_test(lrtaap.problem, s_prime))
        lrtaap.action = nothing;
        return nothing;
    else
        if (!haskey(lrtaap.H, s_prime))
            lrtaap.H[s_prime] = lrtaap.problem.h(lrtaap.problem, s_prime);
        end
        if (!(lrtaap.state === nothing))
            lrtaap.H[lrtaap.state] = reduce(min, learning_realtime_astar_cost(lrtaap,
                                                                        lrtaap.state,
                                                                        b,
                                                                        get_result(lrtaap.problem, lrtaap.state, b),
                                                                        lrtaap.H)
                                        for b in actions(lrtaap.problem, lrtaap.state));
        end
        lrtaap.action = argmin(actions(lrtaap.problem, s_prime),
                                (function(b::String)
                                    return learning_realtime_astar_cost(lrtaap,
                                                                        s_prime,
                                                                        b,
                                                                        get_result(lrtaap.problem, s_prime, b),
                                                                        lrtaap.H);
                                end));
        lrtaap.state = s_prime;
        return lrtaap.action;
    end
end

function genetic_search(problem::T; ngen::Int64=1000, pmut::Float64=0.1, n::Int64=20) where {T <: AbstractProblem}
    local s = problem.initial;
    local states = [result(s, action) for action in actions(problem, s)];
    shuffle!(RandomDeviceInstance, states);
    if (length(states) < n)
        n = length(states);
    end
    return genetic_algorithm(states[1:n], value, ngen=ngen, pmut=pmut);
end

function genetic_algorithm(population::T, fitness::Function; ngen::Int64=1000, pmut::Float64=0.1) where {T <: AbstractVector}
    for i in 1:ngen
        local new_population = Array{Any, 1}();
        for j in 1:length(population)
            local fitnesses = map(fitness, population);
            p1, p2 = weighted_sample_with_replacement(population, fitnesses, 2);
            local child = mate(p1, p2);
            if (rand(RandomDeviceInstance) < pmut)
                mutate(child);
            end
            push!(new_population, child);
        end
        population = new_population;
    end
    return argmax(population, fitness);
end

# Simplified road map of Romania example (Fig. 3.2)
romania = UndirectedGraph(Dict(
                            Pair("A", Dict("Z"=>75, "S"=>140, "T"=>118)),
                            Pair("B", Dict("U"=>85, "P"=>101, "G"=>90, "F"=>211)),
                            Pair("C", Dict("D"=>120, "R"=>146, "P"=>138)),
                            Pair("D", Dict("M"=>75)),
                            Pair("E", Dict("H"=>86)),
                            Pair("F", Dict("S"=>99)),
                            Pair("H", Dict("U"=>98)),
                            Pair("I", Dict("V"=>92, "N"=>87)),
                            Pair("L", Dict("T"=>111, "M"=>70)),
                            Pair("O", Dict("Z"=>71, "S"=>151)),
                            Pair("P", Dict("R"=>97)),
                            Pair("R", Dict("S"=>80)),
                            Pair("U", Dict("V"=>142)),
                                ),
                            Dict{String, Tuple{Any, Any}}(
                                "A"=>( 91, 492), "B"=>(400, 327), "C"=>(253, 288), "D"=>(165, 299),
                                "E"=>(562, 293), "F"=>(305, 449), "G"=>(375, 270), "H"=>(534, 350),
                                "I"=>(473, 506), "L"=>(165, 379), "M"=>(168, 339), "N"=>(406, 537),
                                "O"=>(131, 571), "P"=>(320, 368), "R"=>(233, 410), "S"=>(207, 457),
                                "T"=>( 94, 410), "U"=>(456, 350), "V"=>(509, 444), "Z"=>(108, 531),
                                )
                            );

# One-dimensional state space example (Fig. 4.23)
one_dim_state_space = Graph{String}(dict=Dict{String, Dict{String, String}}([Pair("State_1", Dict([Pair("Right", "State_2")])),
                                        Pair("State_2", Dict([Pair("Right", "State_3"),
                                                            Pair("Left", "State_1")])),
                                        Pair("State_3", Dict([Pair("Right", "State_4"),
                                                            Pair("Left", "State_2")])),
                                        Pair("State_4", Dict([Pair("Right", "State_5"),
                                                            Pair("Left", "State_3")])),
                                        Pair("State_5", Dict([Pair("Right", "State_6"),
                                                            Pair("Left", "State_4")])),
                                        Pair("State_6", Dict([Pair("Left", "State_5")]))]));

one_dim_state_space_least_costs = Dict([Pair("State_1", 8),
                                        Pair("State_2", 9),
                                        Pair("State_3", 2),
                                        Pair("State_4", 2),
                                        Pair("State_5", 4),
                                        Pair("State_6", 3)]);

# Principal states and territories of Australia example (Fig. 6.1)
australia = UndirectedGraph(Dict(
                                Pair("T",   Dict()),
                                Pair("SA",  Dict("WA"=>1, "NT"=>1, "Q"=>1, "NSW"=>1, "V"=>1)),
                                Pair("NT",  Dict("WA"=>1, "Q"=>1)),
                                Pair("NSW", Dict("Q"=>1, "V"=>1)),
                            ),
                            Dict{String, Tuple{Any, Any}}("WA"=>(120, 24), "NT"=>(135, 20), "SA"=>(135, 30),
                                "Q"=>(145, 20), "NSW"=>(145, 32), "T"=>(145, 42), "V"=>(145, 37),
                                )
                            );

#=

    NQueensProblem is the problem of placing 'N' non-attacking queens on a 'N'x'N' chess board.

    Each state is represented as an 'N' element array where the value of 'r' at index 'c' implies

    that a queen occupies the position at row 'r' and column 'c'. The columns are values are filled

    from left to right.

=#
struct NQueensProblem <: AbstractProblem
    N::Int64
    initial::Array{Union{Nothing, Int64}, 1}

    function NQueensProblem(n::Int64)
        return new(n, fill(nothing, n));
    end
end

function actions(problem::NQueensProblem, state::AbstractVector)
    if (!(state[length(state)] === nothing))
        return Array{Any, 1}([]);
    else
        local col = utils.null_index(state);
        return collect(row for row in 1:problem.N if (!conflicted(problem, state, row, col)));
    end
end

function conflict(problem::NQueensProblem, row1::Int64, col1::Int64, row2::Int64, col2::Int64)
    return ((row1 == row2) ||
            (col1 == col2) ||
            (row1 - col1 == row2 - col2) ||
            (row1 + col1 == row2 + col2));
end

function conflict(problem::NQueensProblem, row1::Int64, col1::Int64, row2::Nothing, col2::Int64)
    error("conflict(): 'row2' is not initialized!");
end

function conflicted(problem::NQueensProblem, state::AbstractVector, row::Int64, col::Int64)
    return any(conflict(problem, row, col, state[i], i) for i in 1:(col-1));
end

function get_result(problem::NQueensProblem, state::AbstractVector, row::Int64)
    local col = utils.null_index(state);
    local new_result = deepcopy(state);
    new_result[col] = row;
    return new_result;
end

function goal_test(problem::NQueensProblem, state::AbstractVector)
    if ((state[length(state)] === nothing))
        return false;
    end
    return !any(conflicted(problem, state, state[col], col) for col in 1:length(state));
end

capital_case_alphabet = "ABCDEFGHIJKLMNOPQRSTUVWXYZ";

cubes16 = Array{String, 1}(["FORIXB", "MOQABJ", "GURILW", "SETUPL",
                            "CMPDAE", "ACITAO", "SLCRAE", "ROMASH",
                            "NODESW", "HEFIYE", "ONUDTK", "TEVIGN",
                            "ANEDVZ", "PINESH", "ABILYT", "GKYLEU"]);

function random_boggle(;n::Int64=4)
    local cubes = collect(cubes16[(i % 16) + 1] for i in 0:((n * n) - 1));
    shuffle!(RandomDeviceInstance, cubes);
    return map((function(array::String)
                    return rand(RandomDeviceInstance, collect(array));
                end), cubes);
end

boyan_best = collect("RSTCSDEIAEGNLRPEATESMSSID");

function print_boggle(board::Array{Char, 1})
    local nn::Int64 = length(board);
    local n::Int64 = int_sqrt(nn);
    local board_str::String = "";
    for i in 0:(nn - 1)
        if ((i % n == 0) && (i > 0))
            board_str = board_str * "\n";
        end
        if (board[i + 1] == 'Q')
            board_str = board_str * "Qu ";
        else
            board_str = board_str * String([board[i + 1]]) * "  ";
        end
    end
    print(board_str);
    nothing;
end

function boggle_neighbors(nn::Int64; cache::Dict=Dict{Any, Any}())
    if haskey(cache, nn)
        return cache[nn];
    end
    local n::Int64 = int_sqrt(nn)
    local neighbors::AbstractVector = Array{Any, 1}(undef, nn);
    for i in 0:(nn - 1)
        neighbors[i + 1] = Array{Int64, 1}([]);
        on_top::Bool = (i < n);
        on_bottom::Bool = (i >= (nn - n));
        on_left::Bool = (i % n == 0);
        on_right::Bool = ((i + 1) % n == 0);
        if (!on_top)
            push!(neighbors[i + 1], (i + 1) - n);
            if (!on_left)
                push!(neighbors[i + 1], (i + 1) - n - 1);
            end
            if (!on_right)
                push!(neighbors[i + 1], (i + 1) - n + 1);
            end
        end
        if (!on_bottom)
            push!(neighbors[i + 1], (i + 1) + n);
            if (!on_left)
                push!(neighbors[i + 1], (i + 1) + n - 1);
            end
            if (!on_right)
                push!(neighbors[i + 1], (i + 1) + n + 1);
            end
        end
        if (!on_left)
            push!(neighbors[i + 1], (i + 1) - 1);
        end
        if (!on_right)
            push!(neighbors[i + 1], (i + 1) + 1);
        end
    end
    cache[nn] = neighbors;
    return neighbors;
end

function int_sqrt(n::Number)
    return Int64(sqrt(n));
end

#=

    WordList contains an array of words.

=#
struct WordList
    words::Array{String, 1}
    bounds::Dict{Char, Tuple{Any, Any}}

    function WordList(filename::String; min_len::Int64=3)
        local wlba = read(filename);
        local wls = uppercase(String(wlba));
        local wlsa = sort(map(strip, split(wls, '\n')));
        local wlsa_filtered = collect(s for s in wlsa if (length(s) >= min_len));
        nwl::WordList = new(wlsa_filtered, Dict{Char, Tuple{Any, Any}}());
        for c in capital_case_alphabet
            fc::Char = c + 1;   #following character
            nwl.bounds[c] = (searchsorted(nwl.words, String([c]), 1, length(nwl.words), Base.Order.Forward).stop + 1,
                        searchsorted(nwl.words, String([fc]), 1, length(nwl.words), Base.Order.Forward).stop + 1);
        end
        return nwl;
    end
end

function lookup(wl::WordList, prefix::String; lo::Int64=1, hi::Union{Nothing, Int64}=nothing)
    local words = wl.words;
    if (typeof(hi) <: Nothing)
        hi = length(words);
    end
    local i::Int64 = searchsorted(words, prefix, lo, hi, Base.Order.Forward).start;

    #'i' is only larger than length of words when the returned index is not in WordList.
    if (i <= length(words) && startswith(words[i], prefix))
        return i, (words[i] == prefix);
    else
        return nothing, false;
    end
end

length(wl::WordList) = length(wl.words);

in(prefix::String, wl::WordList) = getindex(lookup(wl, prefix), 2);

#=

    BoggleFinder contains the words found on the Boggle board

    and the array of possible words for the Boggle board.

=#
mutable struct BoggleFinder
    wordlist::WordList
    scores::AbstractVector
    found::Dict
    board::AbstractVector
    neighbors::AbstractVector

    function BoggleFinder(;board::Union{Nothing, AbstractVector}=nothing, fn::Union{Nothing, String}=nothing)
        local wlfn::String;
        if (typeof(fn) <: Nothing)
            if (is_windows())
                wlfn = "..\\aima-data\\EN-text\\wordlist.txt";
            elseif (is_apple() || is_unix())
                wlfn = "../aima-data/EN-text/wordlist.txt";
            end
        else
            wlfn = fn;
        end
        nbf = new(WordList(wlfn),
                vcat([0, 0, 0, 0, 1, 2, 3, 5], fill(11, 100)),
                Dict{Any, Any}())
        if (!(typeof(board) <: Nothing))
            set_board(nbf, board=board);
        end
        return nbf;
    end
end

function set_board(bf::BoggleFinder; board::Union{Nothing, AbstractVector}=nothing)
    if (typeof(board) <: Nothing)
        board = random_boggle();
    end
    bf.board = board;
    bf.neighbors = boggle_neighbors(length(board));
    bf.found = Dict{Any, Any}();
    for i in 1:length(board)
        lo::Int64, hi::Int64 = bf.wordlist.bounds[board[i]];
        find(bf, lo, hi, i, [], "");
    end
    return bf;
end

function find(bf::BoggleFinder, lo::Int64, hi::Int64, i::Int64, visited::AbstractVector, prefix::String)
    if i in visited
        return nothing;
    end
    wordpos, is_word::Bool = lookup(bf.wordlist, prefix, lo=lo, hi=hi);
    if (!(typeof(wordpos) <: Nothing))
        if (is_word)
            bf.found[prefix] = true;
        end
        push!(visited, i);
        local c = bf.board[i];
        if (c == 'Q')
            c = "QU";
        else
            c = String([c]);
        end
        prefix = prefix * c;
        for j in bf.neighbors[i]
            find(bf, wordpos, hi, j, visited, prefix);
        end
        pop!(visited);
    end
    return nothing;
end

function words(bf::BoggleFinder)
    return collect(keys(bf.found));
end

function score(bf::BoggleFinder)
    return sum(collect(scores[len(w)] for w in words(bf)));
end

length(bf::BoggleFinder) = length(bf.found);

"""
    boggle_hill_climbing(;board::Union{Nothing, AbstractVector}=nothing, ntimes::Int64=100, verbose::Bool=true)

Solve the inverse Boggle by using hill climbing (initially use a random Boggle board and changing it).

Return the best Boggle board and its length.
"""
function boggle_hill_climbing(;board::Union{Nothing, AbstractVector}=nothing, ntimes::Int64=100, verbose::Bool=true)
    finder = BoggleFinder();
    if (typeof(board) <: Nothing)
        board = random_boggle();
    end
    local best_length::Int64 = length(set_board(finder, board=board));
    for t in 1:ntimes
        i, old_char = mutate_boggle(board);
        local new_length::Int64 = length(set_board(finder, board=board));
        if (new_length > best_length)
            best_length = new_length;
            if (verbose)
                println(best_length, " ", t, " ", board);
            end
        else
            board[i] = old_char;
        end
    end
    if (verbose)
        print_boggle(board);
    end
    return board, best_length;
end

function mutate_boggle(board::AbstractArray)
    local i::Int64 = rand(RandomDeviceInstance, 1:length(board));
    local old_char::Char = board[i];
    board[i] = rand(RandomDeviceInstance, collect(rand(RandomDeviceInstance, cubes16)));
    return i, old_char;
end

function execute_searcher(searcher::Function, problem::T) where {T <: AbstractProblem}
    local p = InstrumentedProblem(problem);
    searcher(p);
    return p;
end

function compare_searchers(problems::Array{T, 1},
                            header::Array{String, 1};
                            searchers::Array{Function, 1}=[breadth_first_tree_search,
                                                        breadth_first_search,
                                                        depth_first_graph_search,
                                                        iterative_deepening_search,
                                                        depth_limited_search,
                                                        recursive_best_first_search]) where {T <: AbstractProblem}
    local table = vcat(permutedims(hcat(header), [2, 1]),
                        hcat(map(string, searchers),
                            permutedims(reduce(hcat,
                                collect(
                                    collect(format_instrumented_results(execute_searcher(s, p)) for p in problems)
                                for s in searchers)),
                            [2,1])));
    return table;
end

function beautify_node(n::Node)
    return @sprintf("%s%s%s", "<Node ", string(n.state), ">");
end

max_stars_repo_name: mikhail-j/aima-julia-aimacode-fork

max_issues_repo_name: mikhail-j/aima-julia-aimacode-fork

max_forks_repo_name: mikhail-j/aima-julia-aimacode-fork


text:
@testset "Peters' Finite State" begin

    submodel = Submodel(Peters{4}())

    # test provided jacobians
    submodel_jacobian_tests(submodel)

    # test consistency of input/output functions
    submodel_io_tests(submodel)
end

ext: jl


text:
"""
check/get β from contour or explicitly given value
"""
function get_beta(grid::TimeGrid, β)
  im_b = get_branch(grid.contour, imaginary_branch)
  @assert (im_b === nothing && β !== nothing) || (im_b !== nothing && β === nothing)
  β === nothing ? length(im_b) : β
end

function gf_1level(t1::BranchPoint, t2::BranchPoint; ϵ, β)
    -1.0im * (θ(t1, t2) - fermi(ϵ, β)) * exp(-1.0im * (t1.val - t2.val) * ϵ)
end

function gf_1level(grid::TimeGrid; ϵ, β=nothing)
  β = get_beta(grid, β)
  TimeGF(grid) do t1, t2
    gf_1level(t1.val, t2.val, ϵ=ϵ, β=β)
  end
end

function dos2gf(dos, t1::BranchPoint, t2::BranchPoint; β, integrator = dos_integrator)
    theta = θ(t1, t2)
    Δt = t1.val - t2.val
    f = ω -> (ω > 0.0 ? exp(-1.0im * ω * (Δt - 1.0im * (1.0 - theta) * β)) / (exp(-β * ω) + 1) :
                        exp(-1.0im * ω * (Δt + 1.0im * theta * β)) / (exp(β * ω) + 1))
    return -1.0im * (2 * theta - 1) * integrator(f, dos)
end

function dos2gf(dos, grid::TimeGrid; β=nothing, integrator=dos_integrator)
  β = get_beta(grid, β)
  TimeGF(grid, time_invariant=true) do t1, t2
    dos2gf(dos, t1.val, t2.val, β=β, integrator=integrator)
  end
end


text:
export inverselines, singlet_block

"""
    inverselines(nbit::Int; n_reg::Int=nbit) -> ChainBlock

inverse first `n_reg` lines

TODO:
deprecate this function, it is not used.
"""
function inverselines(nbit::Int; n_reg::Int=nbit)
    c = chain(nbit)
    for i = 1:(n_reg ÷ 2)
        push!(c, swap(i,(n_reg-i+1)))
    end
    c
end

function singlet_block(::Type{T}, nbit::Int, i::Int, j::Int) where T
    unit = chain(nbit)
    push!(unit, put(nbit, i=>chain(XGate{T}(), HGate{T}())))
    push!(unit, control(nbit, -i, j=>XGate{T}()))
end

singlet_block(nbit::Int, i::Int, j::Int) = singlet_block(ComplexF64, nbit, i, j)
singlet_block() = singlet_block(2,1,2)

Yao.mat(ρ::DensityMatrix{1}) = dropdims(state(ρ), dims=3)


text:
using Documenter
using SkyCoords

DocMeta.setdocmeta!(SkyCoords, :DocTestSetup, :(using SkyCoords); recursive = true)

makedocs(
    modules = [SkyCoords],
    sitename = "SkyCoords.jl",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        canonical = "https://juliaastro.github.io/SkyCoords.jl",
    ),
    pages = ["Home" => "index.md", "API/Reference" => "api.md"],
    strict = true,
)

deploydocs(repo = "github.com/JuliaAstro/SkyCoords.jl", push_preview=true)


text:
#Function that builds global variable to hold authentication keys
#TODO: Figure out how to do OAuth authentication directly, rather than user putting in credentials directly
function twitterauth(consumer_key::String, consumer_secret::String, oauth_token::String, oauth_secret::String)
    #Create a global variable to hold return from this function
    global TWITTERCRED

    return TWITTERCRED = TWCredential(consumer_key, consumer_secret, oauth_token, oauth_secret)

end

#General functions for OAuth authenticated GET/POST request
get_oauth(endpoint::String, options::Dict) = oauth_request_resource(endpoint, "GET", options, TWITTERCRED.consumer_key, TWITTERCRED.consumer_secret, TWITTERCRED.oauth_token, TWITTERCRED.oauth_secret)

post_oauth(endpoint::String, options::Dict) = oauth_request_resource(endpoint, "POST", options, TWITTERCRED.consumer_key, TWITTERCRED.consumer_secret, TWITTERCRED.oauth_token, TWITTERCRED.oauth_secret)


text:
using Test

using LinearAlgebra, Statistics
using OrdinaryDiffEq
using TransportBasedInference
using TransportBasedInference: evaluate, optimize
using ForwardDiff
using FastGaussQuadrature
using SpecialFunctions
using QuadGK
using Optim
using NLsolve
using MLDataUtils
using Quadrature
using FiniteDiff
using Distributions
using Random

# Test for Tools
# include("tools/tools.jl")
# include("tools/normal.jl")
# include("tools/transform.jl")
#
# # Tools for data assimilation
# include("DA/inflation.jl")
#
#
# # Tools for state-space models
# include("statespace/system.jl")
#
#
# # Functions to manage margins
# include("margin/reducedmargin.jl")
# include("margin/totalorder.jl")
#
# # Test for Hermite polynomials and functions
# include("hermitefunction/phypolyhermite.jl")
# include("hermitefunction/propolyhermite.jl")
#
# include("hermitefunction/phyhermite.jl")
# include("hermitefunction/prohermite.jl")

# Test tools for Basis, MultiBasis, ExpandedFunction
# include("hermitemap/rectifier.jl")
# include("hermitemap/basis.jl")
# include("hermitemap/expandedfunction.jl")
# include("hermitemap/reduced.jl")
# include("hermitemap/expandedfunction2.jl")
# include("hermitemap/multifunction.jl")
#
# # Test tools for integrated function
# include("hermitemap/integratedfunction.jl")
# include("hermitemap/storage.jl")

# Test tools for HermiteMap component
# include("hermitemap/hermitemapcomponent.jl")
# include("hermitemap/linhermitemapcomponent.jl")
#
# # Test greedy procedure
# include("hermitemap/greedyfit.jl")
#
# # Test optimization of HermiteMap component
# include("hermitemap/qr.jl")
# include("hermitemap/qraccelerated.jl")
# include("hermitemap/optimize.jl")
#
# include("hermitemap/inverse.jl")
# # include("hermitemap/hybridinverse.jl")
# include("hermitemap/hermitemap.jl")

# Test for Radial Maps
# include("radialmap/function.jl")
# include("radialmap/separablecomponent.jl")
# include("radialmap/mapcomponent.jl")
# include("radialmap/sparsemapcomponent.jl")
# include("radialmap/map.jl")
# include("radialmap/weights.jl")
# include("radialmap/quantile.jl")
# include("radialmap/cost.jl")
# include("radialmap/optimize.jl")
# include("radialmap/optimize_multithreading.jl")
# include("radialmap/greedyfit.jl")
# include("radialmap/iterativeoptimize.jl")
# include("radialmap/inverse.jl")


text:
using JuMP, BilevelJuMP
using Random

function bench_toll(nodes, optimizer, mode, seed = 1234)

    rng = Random.MersenneTwister(seed)

    Nodes = nodes

    # g, d = LightGraphs.euclidean_graph(Nodes, 2, seed = seed)
    # costs = Dict(k => 10*d[k]^2 for k in keys(d))

    costs = rand(rng, Nodes, Nodes) # complete graph

    cap = min(3, Nodes)

    Fares = 3

    b = zeros(Nodes)
    b[1] = -cap
    b[end] = +cap

    ctrs = []
    ctrs_e = []
    vars = []

    # MOI.empty!(optimizer)
    model = BilevelModel(optimizer, mode = mode)
    try
        JuMP.set_time_limit_sec(model, MAX_TIME)
    catch e
        @show e
        @show "failed to set limit time"
    end

    v = @variable(Upper(model), 0 <= F[n=1:Nodes,w=1:Nodes,s=1:Fares;n>w] <= cap)
    for vv in v; push!(vars, vv); end

    c = @constraint(Upper(model),[n=1:Nodes,w=1:Nodes;n>w],
        sum(F[n,w,s] for s in 1:Fares) <= cap)
    for cc in c; push!(ctrs, cc); end

    #=
        lower model
    =#

    v = @variable(Lower(model), 0 <= f[n=1:Nodes,w=1:Nodes,s=1:Fares;n>w] <= cap)
    for vv in v; push!(vars, vv); end

    v = @variable(Lower(model), 0 <= a[n=1:Nodes,w=1:Nodes,s=1:Fares;n>w] <= 100)
    for vv in v; push!(vars, vv); end

    c = @constraint(Lower(model),[n = 1:Nodes],
         sum(f[n,w,s] for w in 1:Nodes, s in 1:Fares if n>w)
        -sum(f[w,n,s] for w in 1:Nodes, s in 1:Fares if w>n) == b[n])
    for cc in c; push!(ctrs_e, cc); end

    c = @constraint(Lower(model),[n=1:Nodes,w=1:Nodes,s=1:Fares;n>w],
         s * costs[n,w] * f[n,w,s] <= a[n,w,s])
    for cc in c; push!(ctrs, cc); end

    c = @constraint(Lower(model),[n=1:Nodes,w=1:Nodes,s=1:Fares;n>w],
        f[n,w,s] <= F[n,w,s])
    for cc in c; push!(ctrs, cc); end

    @objective(Lower(model), Min,
        sum(costs[n,w] * f[n,w,s] + a[n,w,s]
            for n in 1:Nodes, w in 1:Nodes, s in 1:Fares if n>w)
        )

    #=
        upper model
    =#

    @objective(Upper(model), Max,
        sum(a[n,w,s]
            for n in 1:Nodes, w in 1:Nodes, s in 1:Fares if n>w)
    )

    #=
        Optimize
    =#

    optimize!(model)

    #=
    for v in vars
        val = value(v) # >= 0
    end

    for c in ctrs
        val = value(c) - normalized_rhs(c) # >= 0
    end
    =#

    @show primal_st = primal_status(model)
    @show term_st = termination_status(model) #in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED, MOI.ALMOST_LOCALLY_SOLVED]

    solve_t = JuMP.solve_time(model)
    build_t = BilevelJuMP.build_time(model)

    obj_l = try
        objective_value(Lower(model))
    catch
        NaN
    end
    obj_u = try
        objective_value(Upper(model))
    catch
         NaN
    end
    gap = try
        bound = objective_bound(Upper(model))
        abs(obj_u - bound)/max(abs(bound), 1e-8)
    catch
        NaN
    end

    return primal_st, term_st, solve_t, build_t, obj_l, obj_u, gap
end


text:
################################################
# Description:
# Black-Scholes option pricing
#
# Author: Asaf Manela
# Date: 2015-08-02
################################################
module BlackScholes

using Distributions
N(x) = cdf(Normal(),x)
Nx(x) = pdf(Normal(),x)
d1(S,K,σ,r,T,δ) = (log(S/K) + (r-δ+σ^2/2)*T)/(σ*sqrt(T))
d2(S,K,σ,r,T,δ) = d1(S,K,σ,r,T,δ) - σ*sqrt(T)

BSForward(S,K,r,T,δ) = S*exp(-δ*T) - K*exp(-r*T)
BSCall(S,K,σ,r,T,δ) = S*exp(-δ*T)*N(d1(S,K,σ,r,T,δ)) - K*exp(-r*T)*N(d2(S,K,σ,r,T,δ))
BSPut(S,K,σ,r,T,δ) = K*exp(-r*T)*N(-d2(S,K,σ,r,T,δ)) - S*exp(-δ*T)*N(-d1(S,K,σ,r,T,δ))
BSPremium(S,K,σ,r,T,δ,iscall) = iscall ? BSCall(S,K,σ,r,T,δ) : BSPut(S,K,σ,r,T,δ)

BSCallVega(S,K,σ,r,T,δ) = sqrt(T)*S*exp(-δ*T)*Nx(d1(S,K,σ,r,T,δ))

MAX_ITERATIONS = 100
function BSImpliedVol(S,K,r,T,δ,premium,iscall)
  if iscall
    C = premium
  else # convert using put-call parity
    C = premium + BSForward(S,K,r,T,δ)
  end
  σ = 0.0
  σnext = sqrt(2*pi/T)*C/S
  i = 0
  while abs(σnext - σ) > eps(σnext) && i < MAX_ITERATIONS
    σ = σnext
    σnext = σ - (BSCall(S,K,σ,r,T,δ) - C)/BSCallVega(S,K,σ,r,T,δ)
    i = i+1
  end
  return σnext
end

export BSForward, BSPut, BSCall, BSImpliedVol, BSCallVega, BSPremium, N, Nx, d1, d2
end #module


text:
"""
Connection Genes
"""

mutable struct Connection
    input::Int
    output::Int
    weight::Float64
    disabled::Bool

    function Connection(input::Int, output::Int, weight::Float64=rand(), disabled::Bool=false)
        new(input, output, weight, disabled)
    end
end

Base.copy(c::Connection) = Connection(c.input, c.output, c.weight, c.disabled)


text:
using RationalExtensions
using Base.Test

# write your own tests here
@test Rad(0,0,1) == 0
@test (1 + 2*Sqrt(3))//(2+3*Sqrt(3)) == 16//23 - Sqrt(3)//23
@test Rad(-big(3)//4,1//2,5) == -3//4 + Sqrt(5)//2
@test Rad(big(1),2,3) == Rad(1,2,3)
@test (r = Sqrt(200); (r.a,r.b,r.n) == (0,10,2))
@test 3//4 * Sqrt(3)//4 + Sqrt(3) == 19Sqrt(3)//16
@test true*(1-2Sqrt(5)) == 1 - 2Sqrt(5)
@test Sqrt(3)/Sqrt(5) ≈ 0.7745966692414833
@test isa(Sqrt(3)/Sqrt(5),Float64)
@test 2/Sqrt(5) ≈ 0.8944271909999159
@test -(1+2Sqrt(3)) == -1 - 2Sqrt(3)
@test big(3//4) + 2//3*Sqrt(1) == 17//12
@test convert(Rad,3//4) == 3//4
@test isa(convert(Rad,3//4),Rad{Rational{Int64},Int64})
@test repr(1-Sqrt(5)//2) == "1 - √(5)//2"


text:
# Contraction

function Contract(left::ArrayEvaluable, right::ArrayEvaluable, linds::Vector{Int}, rinds::Vector{Int}, tinds::Vector{Int})
    Contract(ArrayEvaluable[left, right], [linds, rinds], tinds)
end

function Contract(left::Contract, right::ArrayEvaluable, linds::Vector{Int}, rinds::Vector{Int}, tinds::Vector{Int})
    rename = Dict{Int,Int}(zip(linds, left.target))
    new = newindex(left) - 1

    for ind in Iterators.flatten((rinds, tinds))
        get!(rename, ind) do
            new += 1
        end
    end

    newindices = copy(left.indices)
    push!(newindices, [rename[i] for i in rinds])
    newterms = copy(left.terms)
    push!(newterms, right)

    Contract(newterms, newindices, [rename[i] for i in tinds])
end

Contract(left::ArrayEvaluable, right::Contract, linds::Vector{Int}, rinds::Vector{Int}, tinds::Vector{Int}) =
    Contract(right, left, rinds, linds, tinds)

function Contract(left::Zeros, right::ArrayEvaluable, linds::Vector{Int}, rinds::Vector{Int}, tinds::Vector{Int})
    lmap = Dict(k => v for (v,k) in enumerate(linds))
    rmap = Dict(k => v for (v,k) in enumerate(rinds))
    new_shape = Tuple(axid in keys(lmap) ? size(left, lmap[axid]) : size(right, rmap[axid]) for axid in tinds)
    new_type = promote_type(eltype(left), eltype(right))
    Zeros(new_type, new_shape...)
end

Contract(left::ArrayEvaluable, right::Zeros, linds::Vector{Int}, rinds::Vector{Int}, tinds::Vector{Int}) =
    Contract(right, left, rinds, linds, tinds)

# Ensure that Inflate commutes past Contract
function Contract(left::Inflate, right::ArrayEvaluable, linds::Vector{Int}, rinds::Vector{Int}, tinds::Vector{Int})
    # Find the index expressions associated with the axes that are contracted over,
    # apply a corresponding getindex followed by a contraction
    contr_indices = Dict(axid => get(left.map, i, :) for (i, axid) in enumerate(linds))
    new_right = right[(get(contr_indices, axid, :) for axid in rinds)...]
    new_contraction = Contract(left.data, new_right, linds, rinds, tinds)

    # Mappings from name to axis index for both the left and right
    lmap = Dict(k => v for (v,k) in enumerate(linds))
    rmap = Dict(k => v for (v,k) in enumerate(rinds))

    # Inflate the output
    new_inds = IndexMap((k => left.map[lmap[k]] for k in tinds if k in keys(lmap))...)
    new_shape = Tuple(axid in keys(lmap) ? size(left, lmap[axid]) : size(right, rmap[axid]) for axid in tinds)
    Inflate(new_contraction, new_shape, new_inds)
end

Contract(left::ArrayEvaluable, right::Inflate, linds::Vector{Int}, rinds::Vector{Int}, tinds::Vector{Int}) =
    Contract(right, left, rinds, linds, tinds)

# Convenience contraction constructors

# Contract laxis and raxis, leaving other axes as before (left before right)
function Contract(left::ArrayEvaluable, right::ArrayEvaluable, laxis::Int, raxis::Int)
    (nl, nr) = ndims(left), ndims(right)
    tinds = collect(Int, 1:nl+nr-2)
    linds = collect(Int, 1:nl-1)
    rinds = collect(Int, nl:nl+nr-2)
    insert!(linds, laxis, nl+nr-1)
    insert!(rinds, raxis, nl+nr-1)
    Contract(left, right, linds, rinds, tinds)
end

# Matrix multiplication onto a specific axis, from right or left
function Contract(left::ArrayEvaluable, right::ArrayEvaluable, onto::Int; fromright::Bool=false)
    @assert ndims(left) == 2
    tinds = collect(Int, 1:ndims(right))
    rinds = collect(Int, 1:ndims(right))

    z = ndims(right) + 1
    rinds[onto] = z
    linds = fromright ? [z, onto] : [onto, z]

    Contract(left, right, linds, rinds, tinds)
end

# Contract last axis on the left with the first axis on the right
Contract(left::ArrayEvaluable, right::ArrayEvaluable) = Contract(left, right, ndims(left), 1)

Base.:*(left::ArrayEvaluable, right::ArrayEvaluable) = Contract(left, right)
Base.:*(left::ArrayEvaluable, right) = Contract(left, asarray(right))
Base.:*(left, right::ArrayEvaluable) = Contract(asarray(left), right)

# GetIndex

Base.getindex(self::ArrayEvaluable, indices...) = GetIndex(self, IndexMap(indices...))

function Base.getindex(self::Constant, indices...)
    indices = Index[isa(ix, Colon) ? ix : asarray(ix) for ix in indices]
    if all(isa(ix, Colon) || isa(ix, Constant) for ix in indices)
        Constant(self.value[(isa(i, Colon) ? i : i.value for i in indices)...])
    else
        GetIndex(self, IndexMap(indices...))
    end
end

# Ensure that Inflate commutes past GetIndex
function Base.getindex(self::Inflate, indices...)
    map = IndexMap(indices...)
    @assert isempty(keys(self.map) ∩ keys(map))
    @assert dimcheck(map, 0, 1)

    new_inflatemap = Pair[]
    new_shape = Int[]
    skipped_axes = 0
    for i in 1:ndims(self)
        i ∉ keys(self.map) && i ∉ keys(map) && (push!(new_shape, size(self, i)); continue)
        if i in keys(self.map)
            push!(new_inflatemap, (i - skipped_axes) => self.map[i])
            push!(new_shape, size(self, i))
        elseif ndims(map, i) == 1
            push!(new_shape, size(map, i, 1))
        elseif ndims(map, i) == 0
            skipped_axes += 1
        end
    end

    new_data = isempty(map) ? self.data : map[self.data]
    Inflate(new_data, Tuple(new_shape), IndexMap(new_inflatemap...))
end

# Inflate

Inflate(source::ArrayEvaluable, shape::Shape, indices::Vector{Index}) = Inflate(source, shape, IndexMap(indices...))
Inflate(source::ArrayEvaluable, shape::Shape, indices...) = Inflate(source, shape, IndexMap(indices...))

# InsertAxis

InsertAxis(self::ArrayEvaluable, axes::Tuple{Vararg{Int}}) = InsertAxis(self, collect(Int, axes))
InsertAxis(self::ArrayEvaluable, axes::Int...) = InsertAxis(self, collect(Int, axes))

function InsertAxis(self::Constant, axes::Vector{Int})
    isempty(axes) && return self
    newshape = collect(Int, size(self))
    insertmany!(newshape, axes, 1)
    Constant(reshape(self.value, newshape...))
end

# Ensure that Inflate commutes past InsertAxis
function InsertAxis(self::Inflate, axes::Vector{Int})
    isempty(axes) && return self
    newdata = InsertAxis(self.data, axes)
    newindices = IndexMap((k + sum(axes .<= k) => v for (k, v) in self.map)...)
    newshape = collect(Int, size(self))
    insertmany!(newshape, axes, 1)
    Inflate(newdata, Tuple(newshape), newindices)
end

expandleft(self::ArrayEvaluable, tdims::Int) = InsertAxis(self, fill(1, tdims - ndims(self)))
expandright(self::ArrayEvaluable, tdims::Int) = InsertAxis(self, fill(ndims(self) + 1, tdims - ndims(self)))

# Inv

Base.inv(self::ArrayEvaluable{T,2}) where T = Inv(self)

# Derivatives

function grad(self::ArrayEvaluable, geom::ArrayEvaluable)
    @assert ndims(geom) == 1
    dims = size(geom, 1)
    grad(self, dims) * inv(grad(geom, dims))
end

function grad(self::ApplyTransform, d::Int)
    self.dims == d || error("Inconsistent number of dimension")
    ApplyTransformGrad(self.trans, self.arg, self.dims)
end

function grad(self::Contract, d::Int)
    new = newindex(self)
    newtarget = [self.target; new]

    terms = ArrayEvaluable[]
    for i in 1:length(self.terms)
        newterms, newinds = copy(self.terms), copy(self.indices)
        newterms[i] = grad(newterms[i], d)
        newinds[i] = [newinds[i]; new]
        push!(terms, Contract(newterms, newinds, newtarget))
    end

    .+(terms...)
end

function grad(self::Monomials, d::Int)
    self.degree == 0 && return Zeros(eltype(self), size(self)..., d)
    newmono = Monomials(self.points, self.degree-1, self.padding+1)
    gradpts = InsertAxis(grad(self.points, d), 1)
    scale = Constant(Matrix(Diagonal(append!(fill(0.0, (self.padding+1,)), 1:self.degree))))
    scale * (newmono .* gradpts)
end

function grad(self::Product, d::Int)
    # Since the gradient dimension comes last, we need to explicitly broadcast first
    maxndims = maximum(ndims(term) for term in self.terms)
    vterms = [expandright(term, maxndims) for term in self.terms]
    gterms = [grad(term, d) for term in vterms]

    ret = gterms[1]
    for (i, (vt, gt)) in enumerate(zip(vterms[2:end], gterms[2:end]))
        ret = .+(ret .* vt, .*(vterms[1:i]..., gt))
    end
    ret
end

function grad(self::Add, d::Int)
    # Since the gradient dimension comes last, we need to explicitly broadcast first
    maxndims = maximum(ndims(term) for term in self.terms)
    .+((grad(expandright(term, maxndims), d) for term in self.terms)...)
end

grad(self::Constant, d::Int) = Zeros(eltype(self), size(self)..., d)
grad(self::GetIndex, d::Int) = getindex(grad(self.value, d), self.indices..., :)
grad(self::Inflate, d::Int) = Inflate(grad(self.data, d), (size(self)..., d), self.map)
grad(self::InsertAxis, d::Int) = InsertAxis(grad(self.source, d), self.axes)
grad(self::Inv, d::Int) = -Contract(self, self * grad(self.source, d), 2; fromright=true)
grad(self::Neg, d::Int) = -grad(self.source, d)
grad(self::Normalize, d::Int) = grad(self.source, d)
grad(self::Reshape, d::Int) = reshape(grad(self.source, d), self.newshape..., d)
grad(self::Zeros, d::Int) = Zeros(eltype(self), size(self)..., d)

# Reshape

function Base.reshape(self::ArrayEvaluable, newshape::Int...)
    newshape == size(self) && return self
    Reshape(self, newshape)
end

# Multiplication

Product(left::ArrayEvaluable, right::ArrayEvaluable) = Product((left, right))
Product(left::Constant, right::Constant) = Constant(left.value .* right.value)

function Product(left::ArrayEvaluable, right::Zeros)
    newtype = promote_type(eltype(left), eltype(right))
    newshape = broadcast_shape(size(left), size(right))
    newdims = length(newshape)
    Zeros(newtype, newshape...)
end
Product(left::Zeros, right::ArrayEvaluable) = Product(right, left)

# Ensure that Inflate commutes past Product
function Product(left::Inflate, right::Inflate)
    isempty(keys(left.map) ∩ keys(right.map)) || error("Joining inflations with overlapping axes not supported")
    ndims(left) == ndims(right) || error("Joining inflations with different number of dimensions not supported")

    newdata = left.data .* right.data
    newmap = IndexMap(left.map..., right.map...)
    newshape = broadcast_shape(size(left), size(right))
    Inflate(newdata, newshape, newmap)
end

Broadcast.broadcasted(::typeof(*), left::ArrayEvaluable) = left
Broadcast.broadcasted(::typeof(*), left::ArrayEvaluable, right::ArrayEvaluable) = Product(left, right)
Broadcast.broadcasted(::typeof(*), left::ArrayEvaluable, right) = Product(left, asarray(right))
Broadcast.broadcasted(::typeof(*), left, right::ArrayEvaluable) = Product(asarray(left), right)

# Summation and negation

Add(left::ArrayEvaluable, right::ArrayEvaluable) = Add((left, right))
Add(left::Constant, right::Constant) = Constant(left.value .+ right.value)

# TODO: Apply broadcasting during optimization phase?
function Add(left::ArrayEvaluable, right::Zeros)
    any(s != 1 for s in size(right)[ndims(left)+1:end]) && return Add((left, right))
    expandright(left, ndims(right))
end
Add(left::Zeros, right::ArrayEvaluable) = Add(right, left)

Broadcast.broadcasted(::typeof(+), left::ArrayEvaluable) = left
Broadcast.broadcasted(::typeof(+), left::ArrayEvaluable, right::ArrayEvaluable) = Add(left, right)
Broadcast.broadcasted(::typeof(+), left::ArrayEvaluable, right) = Add(left, asarray(right))
Broadcast.broadcasted(::typeof(+), left, right::ArrayEvaluable) = Add(asarray(left), right)

Neg(self::Neg) = self.source

Base.:-(self::ArrayEvaluable) = Neg(self)
Broadcast.broadcasted(::typeof(-), self::ArrayEvaluable) = Neg(self)

# Dimension reduction

Base.dropdims(self::ArrayEvaluable, dims) = DropDims(self, Tuple(dims))

# Ensure that Inflate commutes past DropDims
function Base.dropdims(self::Inflate, dims)
    dims = collect(dims)
    newdata = dropdims(self.data, dims)
    newmap = IndexMap((k + sum(dims .< k) => v for (k, v) in self.map if k ∉ dims)...)
    newshape = Tuple(k for (i, k) in enumerate(size(self)) if i ∉ dims)
    Inflate(newdata, newshape, newmap)
end

Base.sum(self::ArrayEvaluable, dims) = Sum(self, Tuple(dims))

# Ensure that Inflate commutes past Sum
function Base.sum(self::Inflate, dims)
    dims = collect(dims)
    newdata = sum(self.data, dims)
    newmap = IndexMap((k => v for (k, v) in self.map if k ∉ dims)...)
    newshape = Tuple(i in dims ? 1 : k for (i, k) in enumerate(size(self)))
    Inflate(newdata, newshape, newmap)
end

# Miscellaneous

const element = Argument{Element}(:element, false, true)
const elemindex = ArrayArgument{Int,0}(:(fill(element.index, ())), false, true, ())
const dimtrans = Argument{TransformChain}(:(element.transform[1]), false, true)
const fulltrans = Argument{TransformChain}(:(element.transform), false, true)
const _point = Argument{Vector{Float64}}(:point, false, false)
point(n::Int) = ApplyTransform(dimtrans, _point, n)

outer(left::ArrayEvaluable, right::ArrayEvaluable) = InsertAxis(left, 1) .* InsertAxis(right, 2)
outer(self::ArrayEvaluable) = outer(self, self)
rootcoords(ndims::Int) = ApplyTransform(fulltrans, _point, ndims)

max_stars_repo_name: TheBB/Jutils

max_issues_repo_name: TheBB/Jutils

max_forks_repo_name: TheBB/Jutils


text:
function ChainRulesCore.rrule(::typeof(getindex), x::ITensor, I...)
  y = getindex(x, I...)
  function getindex_pullback(ȳ)
    # TODO: add definition `ITensor(::Tuple{}) = ITensor()`
    # to ITensors.jl so no splatting is needed here.
    x̄ = ITensor(inds(x)...)
    x̄[I...] = unthunk(ȳ)
    Ī = broadcast_notangent(I)
    return (NoTangent(), x̄, Ī...)
  end
  return y, getindex_pullback
end

# Specialized version in order to avoid call to `setindex!`
# within the pullback, should be better for taking higher order
# derivatives in Zygote.
function ChainRulesCore.rrule(::typeof(getindex), x::ITensor)
  y = x[]
  function getindex_pullback(ȳ)
    x̄ = ITensor(unthunk(ȳ))
    return (NoTangent(), x̄)
  end
  return y, getindex_pullback
end

function setinds_pullback(ȳ, x, a...)
  x̄ = ITensors.setinds(ȳ, inds(x))
  ā = broadcast_notangent(a)
  return (NoTangent(), x̄, ā...)
end

function inv_op(f::Function, args...; kwargs...)
  return error(
    "Trying to differentiate `$f` but the inverse of the operation (`inv_op`) `$f` with arguments $args and keyword arguments $kwargs is not defined.",
  )
end

function inv_op(::typeof(prime), x, n::Integer=1; kwargs...)
  return prime(x, -n; kwargs...)
end

function inv_op(::typeof(replaceprime), x, n1n2::Pair; kwargs...)
  return replaceprime(x, reverse(n1n2); kwargs...)
end

function inv_op(::typeof(swapprime), x, n1n2::Pair; kwargs...)
  return swapprime(x, reverse(n1n2); kwargs...)
end

function inv_op(::typeof(addtags), x, args...; kwargs...)
  return removetags(x, args...; kwargs...)
end

function inv_op(::typeof(removetags), x, args...; kwargs...)
  return addtags(x, args...; kwargs...)
end

function inv_op(::typeof(replacetags), x, n1n2::Pair; kwargs...)
  return replacetags(x, reverse(n1n2); kwargs...)
end

function inv_op(::typeof(swaptags), x, n1n2::Pair; kwargs...)
  return swaptags(x, reverse(n1n2); kwargs...)
end

function inv_op(::typeof(replaceind), x, n1n2::Pair; kwargs...)
  return replaceind(x, reverse(n1n2); kwargs...)
end

function inv_op(::typeof(replaceinds), x, n1n2::Pair; kwargs...)
  return replaceinds(x, reverse(n1n2); kwargs...)
end

function inv_op(::typeof(swapind), x, args...; kwargs...)
  return swapind(x, reverse(args)...; kwargs...)
end

function inv_op(::typeof(swapinds), x, args...; kwargs...)
  return swapinds(x, reverse(args)...; kwargs...)
end

_check_inds(x::ITensor, y::ITensor) = hassameinds(x, y)
_check_inds(x::MPS, y::MPS) = hassameinds(siteinds, x, y)
_check_inds(x::MPO, y::MPO) = hassameinds(siteinds, x, y)

for fname in (
  :prime,
  :setprime,
  :noprime,
  :replaceprime,
  :swapprime,
  :addtags,
  :removetags,
  :replacetags,
  :settags,
  :swaptags,
  :replaceind,
  :replaceinds,
  :swapind,
  :swapinds,
)
  @eval begin
    function ChainRulesCore.rrule(
      f::typeof($fname), x::Union{ITensor,MPS,MPO}, a...; kwargs...
    )
      y = f(x, a...; kwargs...)
      function f_pullback(ȳ)
        x̄ = inv_op(f, unthunk(ȳ), a...; kwargs...)
        if !_check_inds(x, x̄)
          error(
            "Trying to differentiate function `$f` with arguments $a and keyword arguments $kwargs. The forward pass indices $(inds(x)) do not match the reverse pass indices $(inds(x̄)). Likely this is because the priming/tagging operation you tried to perform is not invertible. Please write your code in a way where the index manipulation operation you are performing is invertible. For example, `prime(A::ITensor)` is invertible, with an inverse `prime(A, -1)`. However, `noprime(A)` is in general not invertible since the information about the prime levels of the original tensor are lost. Instead, you might try `prime(A, -1)` or `replaceprime(A, 1 => 0)` which are invertible.",
          )
        end
        ā = broadcast_notangent(a)
        return (NoTangent(), x̄, ā...)
      end
      return y, f_pullback
    end
  end
end

@non_differentiable permute(::Indices, ::Indices)


text:
using EscapeString
using Base.Test

s  = "
Hi there 'Sam'!
"
@test escapestring(s, :html_text) == "<div><b>Hi there 'Sam'!</b></div>"

# EOF


text:
include("cleanup.jl")

if isdir(joinpath("..","build"))
  @info "build dir exists"
else
  @info "build dir does not exist"
  mkdir(joinpath("..","build"))
end

if isdir(joinpath("..", "build", "guides"))
  @info "guides dir exist"
else
  @info "guides dir does not exist"
  mkdir(joinpath("..", "build", "guides"))
end

if isdir(joinpath("..", "build", "tutorials"))
  @info "tutorials dir exist"
else
  @info "tutorials dir does not exist"
  mkdir(joinpath("..", "build", "tutorials"))
end

t1 = `julia --project -e 'using Pkg; Pkg.instantiate(); include("runtest_tutorial.jl")'` |> run

if Base.process_exited(t1)
  Base.kill(t1)
end

t2 = `julia --project -e 'using Pkg; Pkg.instantiate(); include("runtest_tutorial_breaking.jl")'` |> run

if Base.process_exited(t2)
  Base.kill(t2)
end

t3 = `julia --project -e 'using Pkg; Pkg.instantiate(); include("runtest_guide.jl")'` |> run

if Base.process_exited(t3)
  Base.kill(t3)
end

t4 = `julia --project -e 'using Pkg; Pkg.instantiate(); include("runtest_guide_breaking.jl")'` |> run

if Base.process_exited(t4)
  Base.kill(t4)
end


text:
using MathPhysicalConstants, Measurements, SpecialFunctions, QuadGK, Calculus, Unitful
import MathPhysicalConstants.SI: α, O_M, c, e, ε0, h, µ0, ℯ, PlanckConstantHbar, ƛe
using Test

import Base: isapprox
import Measurements: value, uncertainty

if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@testset "Comparisons and Tests" begin
    # isapprox is overloaded only for tests.  Make sure it works as expected.
    @test measurement(big"0.75", 0.01) == 3//4 ± 1//100
    @test measurement(big(π)) ≠ π
    @test ℯ ≠ measurement(Float32(ℯ))
    @test 3//4 == measurement(Float32(0.75), Float32(0)) ≠ 4//3
    @test MathPhysicalConstants.SI.PlanckConstantHbar*MathPhysicalConstants.SI.c == 3.1615262057012976e-26
    @test MathPhysicalConstants.SI.PlanckConstantHbar/(MathPhysicalConstants.SI.MassElectron*MathPhysicalConstants.SI.c) == 3.861592046068738e-13
    @test MathPhysicalConstants.SI.GravitationalConstant^2 == 4.45433438464e-21
    @test MathPhysicalConstants.SI.MassElectron*MathPhysicalConstants.SI.MassProton*MathPhysicalConstants.SI.MassHydrogen == 2.550178071499845e-84
    @test (2*MathPhysicalConstants.SI.PlanckConstantHbar^2)/(MathPhysicalConstants.SI.SanchezGravitationalConstant*MathPhysicalConstants.SI.MassElectron*MathPhysicalConstants.SI.MassProton*MathPhysicalConstants.SI.MassNeutron) == 1.3056254088279212e26
    @test MathPhysicalConstants.SI.KotovCosmicPeriodicity == 9600.061
    @test MathPhysicalConstants.SI.SanchezGravitationalConstant == 6.675453818e-11
    @test MathPhysicalConstants.SI.O_M == 80801742479451287588645990496171075700575436800000000
    @test isinteger(-3 ± 0)
    @test !isinteger(3 ± 1)
    @test !isinteger(pi ± 0)
    @test iszero(0 ± 0)
    @test iszero(zero(Measurement))
    @test !iszero(0 ± 1)
    @test !iszero(1 ± 0)
    @test isone(1 ± 0)
    @test isone(one(Measurement))
    @test !isone(1 ± 1)
    @test !isone(0 ± 0)
end

@testset "Base" begin
#    @test ustrip(big(h)) == big"6.626070040e-34"
    @test setprecision(BigFloat, 768) do; precision(ustrip(big(c))) end == 768
    @test measurement(h) === measurement(h)
    @test iszero(measurement(α) - measurement(α))
#    @test isone(measurement(BigFloat, ƛe) / measurement(BigFloat, ƛe))
#    @test iszero(measurement(BigFloat, PlanckConstantHbar) - (measurement(BigFloat, h) / 2big(pi)))
#    @test isone(measurement(BigFloat, PlanckConstantHbar) / (measurement(BigFloat, h) / 2big(pi)))
end

@testset "Utils" begin
    @test c === ustrip(float(c))
    @test_throws ErrorException h.foo
end

#@testset "Promotion" begin
#    x = @inferred(inv(μ0 * c ^ 2))
#    T = promote_type(typeof(ε0), typeof(x))
#    u = u"A^2 * kg^-1 * m^-3 * s^4"
#    @test promote_type(typeof(α), BigInt) === BigFloat
#    @test promote_type(typeof(α), typeof(1u"m/cm")) === Float64
#    @test T === Quantity{Float64, dimension(x), typeof(u)}
#    @test convert(T, ε0) === (ε0 / unit(ε0)) * u
#    @test convert(Float32, α) === float(Float32, α)
#    @test uconvert(u"cm/s", c) === uconvert(u"cm/s", float(c))
#end
@testset "Maths" begin
#    @test α ≈ @inferred(ℯ^2/(4 * pi * ε0 * PlanckConstantHbar * c))
    @test @inferred(α + 2) ≈ 2 + float(α)
    @test @inferred(5 + α) ≈ float(α) + 5
    @test @inferred(α + 2.718) ≈ 2.718 + float(α)
    @test @inferred(-3.14 + α) ≈ float(α) - 3.14
    @test ε0 ≈ @inferred(1 / (μ0 * c ^ 2))
    @test @inferred(big(0) + α) == big(α)
    @test @inferred(α * 1.0) == float(α)
end

@testset "Show" begin
    @test repr(c) ==
        "2.99792458e8"

end
# using Base.Test
#include("../src/RunTests.jl")


text:
"""
    contrast(img::AbstractMatrix) -> AbstractFloat

计算图像`img`的对比度
"""
contrast(img::AbstractMatrix) = std(img) / mean(img)

"""
    ssi(; noised::AbstractMatrix{T}, filtered::AbstractMatrix{T})) -> AbstractFloat  where T

计算图像`img`的SSI值

# Arguments
- `noised`: 噪声图像
- `filtered`: 降噪后的图像
"""
ssi(; noised::AbstractMatrix{T}, filtered::AbstractMatrix{T}) where T = contrast(filtered) / contrast(noised)

"""
    smpi(; noised::AbstractMatrix{T}, filtered::AbstractMatrix{T})) -> AbstractFloat  where T

计算图像`img`的SMPI值

# Arguments
- `noised`: 噪声图像
- `filtered`: 降噪后的图像
"""
smpi(; noised::AbstractMatrix{T}, filtered::AbstractMatrix{T}) where T = (1 + abs(mean(filtered) - mean(noised))) * (std(filtered) / std(noised))

"""
    enl(img::AbstractMatrix) -> AbstractFloat

计算图像`img`的`ENL`值
"""
enl(img::AbstractMatrix) = (mean(img) / std(img)) ^2


text:
const DD = DimensionalData
const TYPE_dimname = Union{Symbol,AbstractString}
const dimnames_default = ["x", "y", "z", "t"]

# https://github.com/rafaqz/DimensionalData.jl/issues/370

# function make_dims(vals::Vector, names::Vector{<:TYPE_dimname})
#   (; zip(Symbol.(names), vals)...)
# end

# function make_dims(val, name::TYPE_dimname)
#   make_tuple([val], [name])
# end

# vals = [1:i for i in Size] |> collect
# make_tuple(vals, dimnames[1:n])

function make_dims(array::AbstractArray, dimnames::Vector{<:TYPE_dimname} = dimnames_default)
  n = size(array) |> length
  length(dimnames) < n && error("The length of `names` is short than `dims`")
  Size = size(array)
  dimnames = Symbol.(dimnames)

  [Dim{dimnames[i]}(1:Size[i]) for i in 1:length(Size)]
end

function DimensionalData.DimArray(
  array::AbstractArray, dimnames::Vector{<:TYPE_dimname} = dimnames_default; kw...)

  dims = make_dims(array, dimnames)
  DimArray(array, Tuple(dims), kw...) #
end

const TYPE_DIM = Union{AbstractString,Symbol,Integer,Type{<:Dim}}

which_dim(d, dim::AbstractString) = findall(string.(name(d.dims)) .== dim)[1]
which_dim(d, dim::Symbol) = findall(name(d.dims) .== dim)[1]
which_dim(d, dim::Type{<:Dim}) = findall(name(d.dims) .== name(dim))[1]
which_dim(d, dim::Integer) = dim

# not passed test
dimnum2(d, dims::Vector{<:TYPE_DIM}) = begin
  nums = [which_dim(d, dim) for dim in dims]
  tuple(nums...)
end

dimnum2(d, dims::TYPE_DIM) = begin
  dimnum2(d, [dims])
end

export DimArray, which_dim, which_dims, dimnum2


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

################################## INSTALLATION ##################################
# (1) Enter Package Manager via     ]
# (2) Install FMI via               add "https://github.com/ThummeTo/FMI.jl"
# (3) Install FMIFlux via           add "https://github.com/ThummeTo/FMIFlux.jl"
################################ END INSTALLATION ################################

using FMI
using FMIFlux
using Flux
using DifferentialEquations: Tsit5
import Plots

modelFMUPath = joinpath(dirname(@__FILE__), "../model/SpringPendulum1D.fmu")
realFMUPath = joinpath(dirname(@__FILE__), "../model/SpringFrictionPendulum1D.fmu")

t_start = 0.0
t_step = 0.01
t_stop = 5.0

myFMU = fmiLoad(realFMUPath)
fmiInstantiate!(myFMU; loggingOn=false)
fmiSetupExperiment(myFMU, t_start, t_stop)

fmiEnterInitializationMode(myFMU)
fmiExitInitializationMode(myFMU)

x0 = fmi2GetContinuousStates(myFMU)

realSimData = fmi2SimulateCS(myFMU, t_step, t_start, t_stop, ["mass.s", "mass.v", "mass.f", "mass.a"], false)
fmiUnload(myFMU)

fmiPlot(realSimData)

displacement = 0.1
myFMU = fmiLoad(modelFMUPath)

fmiInstantiate!(myFMU; loggingOn=false)
fmuSimData = fmiSimulate(myFMU, t_step, t_start, t_stop, ["mass.s", "mass.v", "mass.a"])
fmiReset(myFMU)
fmiSetupExperiment(myFMU, 0.0)

fmi2SetReal(myFMU, "fixed.s0", displacement)
fmi2SetReal(myFMU, "mass_s0", 0.5 + displacement) # das sollte wurst sein oder?

fmiEnterInitializationMode(myFMU)
fmiExitInitializationMode(myFMU)
#x0 = fmi2GetContinuousStates(myFMU)
tData = t_start:t_step:t_stop
posData = fmi2SimulationResultGetValues(realSimData, "mass.s")
velData = fmi2SimulationResultGetValues(realSimData, "mass.v")

# loss function for training
function losssum()
    solution = problem(t_start, x0)

    posNet = collect(data[2] for data in solution.u)
    velNet = collect(data[3] for data in solution.u)

    Flux.Losses.mse(posNet, posData) + Flux.Losses.mse(velNet, velData)
end

# callback function for training
global iterCB = 0
function callb()
    global iterCB += 1

    # freeze first layer parameters (2,4,6) for velocity -> (static) direct feed trough for velocity
    # parameters for position (1,3,5) are learned
    p_net[1][2] = 0.0
    p_net[1][4] = 1.0
    p_net[1][6] = 0.0

    if iterCB % 10 == 1
        avg_ls = losssum()
        display("Loss: $(round(avg_ls, digits=5))   Avg displacement in data: $(round(sqrt(avg_ls / 2.0), digits=5))   Weight/Scale: $(p_net[1][1])   Bias/Offset: $(p_net[1][5])")
    end

end

# NeuralFMU setup
numStates = fmiGetNumberOfStates(myFMU)

net = Chain(inputs -> NeuralFMUInputLayer(myFMU, inputs),
            Dense(numStates, numStates, identity; initW = (out, in) -> [[1.0, 0.0] [0.0, 1.0]], initb = out -> zeros(out)),
            inputs -> fmi2DoStepME(myFMU, 0.0, inputs),
            Dense(numStates, 16, tanh),
            Dense(16, 16, tanh),
            Dense(16, numStates),
            inputs -> NeuralFMUOutputLayer(inputs))

problem = NeuralFMU(net, (t_start, t_stop), Tsit5(), tData)
problem.fmu = myFMU
solutionBefore = problem(t_start, x0)
fmiPlot(problem)

# train it ...
p_net = Flux.params(problem)

optim = ADAM()
for i in 1:3
    display("epoch: $i/3")
    Flux.train!(losssum, p_net, Iterators.repeated((), 1000), optim; cb=callb)
end

###### plot results s
solutionAfter = problem(t_start, x0)
fig = Plots.plot(xlabel="t [s]", ylabel="mass position [m]", linewidth=2,
    xtickfontsize=12, ytickfontsize=12,
    xguidefontsize=12, yguidefontsize=12,
    legendfontsize=12, legend=:bottomright)
Plots.plot!(fig, tData, fmi2SimulationResultGetValues(fmuSimData, "mass.s"), label="FMU", linewidth=2)
Plots.plot!(fig, tData, posData, label="reference", linewidth=2)
Plots.plot!(fig, tData, collect(data[2] for data in solutionAfter.u), label="NeuralFMU", linewidth=2)
Plots.savefig(fig, "exampleResult_s.pdf")

###### plot results v
solutionAfter = problem(t_start, x0)
fig = Plots.plot(xlabel="t [s]", ylabel="mass velocity [m/s]", linewidth=2,
    xtickfontsize=12, ytickfontsize=12,
    xguidefontsize=12, yguidefontsize=12,
    legendfontsize=12, legend=:bottomright)
Plots.plot!(fig, tData, fmi2SimulationResultGetValues(fmuSimData, "mass.v"), label="FMU", linewidth=2)
Plots.plot!(fig, tData, velData, label="reference", linewidth=2)
Plots.plot!(fig, tData, collect(data[3] for data in solutionAfter.u), label="NeuralFMU", linewidth=2)
Plots.savefig(fig, "exampleResult_v.pdf")

# write training parameters *p_net* back to *net* with data offset *c*
function transferParams!(net, p_net, c=0)
    numLayers = length(net.layers)
    for l in 1:numLayers
        ni = size(net.layers[l].weight,2)
        no = size(net.layers[l].weight,1)

        w = zeros(no, ni)
        b = zeros(no)

        for i in 1:ni
            for o in 1:no
                w[o,i] = p_net[1][c + (i-1)*no + (o-1)]
            end
        end

        c += ni*no

        for o in 1:no
            b[o] = p_net[1][c + (o-1)]
        end

        c += no

        copy!(net.layers[l].weight, w)
        copy!(net.layers[l].bias, b)
    end
end

###### friction model extraction
layers = problem.neuralODE.model.layers[4:6]
net_bottom = Chain(layers...)
transferParams!(net_bottom, p_net, 7)

s_neural = collect(data[2] for data in solutionAfter.u)
v_neural = collect(data[3] for data in solutionAfter.u)

s_fmu = fmi2SimulationResultGetValues(fmuSimData, "mass.s")
v_fmu = fmi2SimulationResultGetValues(fmuSimData, "mass.v")
a_fmu = fmi2SimulationResultGetValues(fmuSimData, "mass.a")

f_real = fmi2SimulationResultGetValues(realSimData, "mass.f")
s_real = fmi2SimulationResultGetValues(realSimData, "mass.s")
v_real = fmi2SimulationResultGetValues(realSimData, "mass.v")
a_real = fmi2SimulationResultGetValues(realSimData, "mass.a")

fs = zeros(length(v_real))
for i in 1:length(v_real)
    fs[i] = -net_bottom([v_real[i], 0.0])[2]
end

fig = Plots.plot(xlabel="v [m/s]", ylabel="friction force [N]", linewidth=2,
    xtickfontsize=12, ytickfontsize=12,
    xguidefontsize=12, yguidefontsize=12,
    legendfontsize=12, legend=:bottomright)

mat = hcat(v_real, zeros(length(v_real)))
mat[sortperm(mat[:, 1]), :]
Plots.plot!(fig, mat[:,1], mat[:,2], label="FMU", linewidth=2)

mat = hcat(v_real, f_real)
mat[sortperm(mat[:, 1]), :]
Plots.plot!(fig, mat[:,1], mat[:,2], label="reference", linewidth=2)

mat = hcat(v_real, fs)
mat[sortperm(mat[:, 1]), :]
Plots.plot!(fig, mat[:,1], mat[:,2], label="NeuralFMU", linewidth=2)

Plots.savefig(fig, "frictionModel.pdf")

#########
layers = problem.neuralODE.model.layers[2:2]
net_top = Chain(layers...)
transferParams!(net_top, p_net, 1)

disp_s = zeros(length(s_real))
for i in 1:length(s_real)
    disp_s[i] = net_top([s_real[i], 0.0])[1] - s_real[i] - displacement
end

fig = Plots.plot(xlabel="t [s]", ylabel="displacement [m]", linewidth=2,
    xtickfontsize=12, ytickfontsize=12,
    xguidefontsize=12, yguidefontsize=12,
    legendfontsize=12, legend=:topright)

Plots.plot!(fig, [t_start, t_stop], [displacement, displacement], label="FMU", linewidth=2)
Plots.plot!(fig, [t_start, t_stop], [0.0, 0.0], label="reference", linewidth=2)
Plots.plot!(fig, tData, disp_s, label="NeuralFMU", linewidth=2)

Plots.savefig(fig, "displacementModel.pdf")

fmiUnload(myFMU)