text:
export  Problem,
            result, step_cost, goal_test, actions, heuristic, state_value, successor_states,
        Node,
            solution, failure, isless,
        SearchAlgorithm,
            BreadthFirstSearch,
                execute,
            UniformCostSearch,
            DepthLimitedSearch,
            IterativeDeepeningSearch,
            RecursiveBestFirstSearch,
            GraphSearch,
                GraphSearchDepth,
                GraphSearchBreadth,
                GraphSearchUniformCost,
                GraphSearchBestFirst,
                GraphSearchAStar,
            HillClimbingSearch,
            SimulatedAnnealingSearch,
            GeneticAlgorithmSearch,
                mutate, reproduce, random_state,
        AIMASearchFIFO,
        AIMASearchLIFO,
        AIMASearchPQ,
        AIMASearchSet

using Compat

using DataStructures

import Base: ==, isless, in, eltype, length, empty!

@compat abstract type Problem end

result(problem::Problem, state::State, action::Action) = error(E_ABSTRACT)
step_cost(problem::Problem, state::State, action::Action) = error(E_ABSTRACT)
goal_test(problem::Problem, state::State) = error(E_ABSTRACT)
actions(problem::Problem,state::State) = error(E_ABSTRACT)
heuristic(problem::Problem, state::State) = error(E_ABSTRACT)
state_value(problem::Problem, state::State) = error(E_ABSTRACT)
successor_states(problem::Problem, state::State) = error(E_ABSTRACT)

search(problem::Problem) = execute(problem.search_algorithm, problem)

mutable struct NodeT{S<:State, T}
    state::S
    parent::Union{NodeT{S, T}, Void}
    action::Action
    path_cost::Number
    f::Number
end

const Node{S <: State} = NodeT{S, :none}

NodeT(state::S, s::Symbol=:none,
      action::Action=Action_NoOp, path_cost::Number=0, f::Number=0) where {S<: State} =
    NodeT{S, s}(state, nothing, action, path_cost, f)

NodeT(state::S, parent::NodeT{S,T},
    action::Action=Action_NoOp,
    path_cost::Number=zero(parent.path_cost),
    f::Number=zero(parent.f)) where {S <:State, T} =
    NodeT{S, T}(state, parent, action, path_cost, f)

isless(n1::NodeT, n2::NodeT) = isless(n1.f, n2.f)

make_node{S<:State}(state::S, s::Symbol) = NodeT(state, s)
make_node{S<:State}(state::S) = make_node(state, :none)

function child_node{S <:State}(problem, parent::NodeT{S, :none}, action)
    state = result(problem, parent.state, action)
    path_cost = parent.path_cost + step_cost(problem, parent.state, action)
    return NodeT(state, parent, action, path_cost, path_cost)
end

function child_node{S <:State}(problem, parent::NodeT{S, :greedy}, action)
    state = result(problem, parent.state, action)
    path_cost = parent.path_cost + step_cost(problem, parent.state, action)
    return NodeT(state, parent, action, path_cost, heuristic(problem, state))
end

function child_node{S <:State}(problem, parent::NodeT{S, :astar}, action)
    state = result(problem, parent.state, action)
    path_cost = parent.path_cost + step_cost(problem, parent.state, action)
    return NodeT(state, parent, action, path_cost, path_cost + heuristic(problem, state))
end

"""
The default implementation of *solution* is returning the sequence of nodes.
This may not be ideal for many implementations. They may provide a more
elaborate reporting as well.
"""
function solution{S<:State, T}(node::NodeT{S, T})
    nodes = Vector{NodeT{S, T}}()
    while true
        unshift!(nodes, node)
        node.parent === nothing && break
        node = node.parent
    end
    return nodes
end

"""
The default implementation of *failure* is throwing an error
"""
failure(node::NodeT) = error(E_NODE_NOT_REACHABLE)

@compat abstract type SearchAlgorithm end

const AIMASearchSet{S<:State} = Set{S}

mutable struct AIMASearchQueue{S<:State, T}
     s::Vector{Node{S}}
     AIMASearchQueue{S, T}() where {S<:State, T} = new(Vector{Node{S}}())
end
eltype(q::AIMASearchQueue) = eltype(q.s)
empty(q::AIMASearchQueue) = empty!(q.s)

const AIMASearchFIFO{S<:State} = AIMASearchQueue{S, :fifo}
const AIMASearchLIFO{S<:State} = AIMASearchQueue{S, :lifo}

isempty(q::AIMASearchQueue) = isempty(q.s)
insert{S<:State, T}(q::AIMASearchQueue{S, T}, d::Node{S}) = push!(q.s, d)
length(q::AIMASearchQueue) = length(q.s)

pop(q::AIMASearchFIFO) = shift!(q.s)
pop(q::AIMASearchLIFO) = pop!(q.s)

function insert_or_replace{S <: State, T}(f::Function,
    q::AIMASearchQueue{S, T}, d::Node{S})
    state = start(q.s)
    while !done(q.s, state)
        ps = state
        (i, state) = next(q.s, state)
        i.state == d.state && return f(i) ? splice!(q.s, ps, [d]) : nothing
    end
    return insert(q.s, d)
end

in{S <: State, T}(state::S, q::AIMASearchQueue{S, T}) = any(x -> x.state == state, q.s)

"""
*BreadthFirstSearch* An uninformed graph search technique to reach goal.

Cost of search is completely ignored.

*frontier* is a FIFO
*explored* is a Set

pg. 82 Fig. 3.11 AIMA 3ed
"""
mutable struct BreadthFirstSearch{FIFO, SET, S<:State} <: SearchAlgorithm
    frontier::FIFO
    explored::SET
    BreadthFirstSearch{FIFO, SET, S}() where{FIFO, SET, S<:State} =
        new(has_trait_queue(FIFO{S}()), has_trait_set(SET{S}()))
end

BreadthFirstSearch{S<:State}(::S) = BreadthFirstSearch{AIMASearchFIFO, AIMASearchSet, S}()

function execute(search::BreadthFirstSearch, problem::Problem)
    empty(search.frontier)
    empty(search.explored)
    node = NodeT(problem.initial_state, :none)
    goal_test(problem, node.state) && return solution(node)
    insert(search.frontier, node)

    while(true)
        isempty(search.frontier) && return failure(node)
        node = pop(search.frontier)
        append(search.explored, node.state)

        for action in actions(problem, node.state)
            child = child_node(problem, node, action)
            !(child.state in search.explored) && !(child.state in search.frontier) &&
            goal_test(problem, child.state) && return solution(child)
            insert(search.frontier, child)
        end
    end
end

const AIMASearchPQ_Base{S<:State, NT} = PriorityQueue{S, NodeT{S, NT}}

pop(pq::AIMASearchPQ_Base) = ((key,val) = peek(pq); dequeue!(pq); val)
peekv(pq::AIMASearchPQ_Base) = ((key,val) = peek(pq); val)
empty(pq::AIMASearchPQ_Base) = while(!isempty(pq));dequeue!(pq);end
insert{S <: State, NT}(pq::AIMASearchPQ_Base{S, NT}, node::NodeT{S, NT}) =
    enqueue!(pq, node.state, node)
replace{S <: State, NT}(pq::AIMASearchPQ_Base{S, NT}, node::NodeT{S, NT}) =
    (pq[node.state] = node)
in{S <: State, NT}(state::S, pq::AIMASearchPQ_Base{S, NT}) = haskey(pq, state)

function insert_or_replace{S <: State, NT}(f::Function,
    pq::AIMASearchPQ_Base{S, NT}, node::NodeT{S, NT})
    !haskey(pq, node.state) && return enqueue!(pq, node.state, node)
    f(pq[node.state]) && return (pq[node.state] = node)
end

const AIMASearchPQ{S <: State}   = AIMASearchPQ_Base{S, :none}
const AIMASearchPQ_G{S <: State} = AIMASearchPQ_Base{S, :greedy}
const AIMASearchPQ_A{S <: State} = AIMASearchPQ_Base{S, :astar}

"""
*UniformCostSearch* An uninformed graph search technique to reach goal.

Cost of search is used in taking decisions.

*frontier* is a priority queue
*explored* is a Set

pg. 82 Fig. 3.11 AIMA 3ed
"""
mutable struct UniformCostSearch{PQ, SET, S<:State} <: SearchAlgorithm
    frontier::PQ
    explored::SET
    UniformCostSearch{PQ, SET, S}() where {PQ, SET, S<:State} =
        has_trait_queue(PQ{S}, Node{S}) &&
        method_exists(replace, (PQ{S}, Node{S})) &&
        method_exists(in, (S, PQ{S})) &&
        method_exists(get, (PQ{S}, S, Node{S})) &&
        new(PQ{S}(), has_trait_set(SET{S}()))
end

UniformCostSearch{S<:State}(::S) = UniformCostSearch{AIMASearchPQ, AIMASearchSet, S}()

function execute(search::UniformCostSearch, problem) #returns a solution, or failure
    empty(search.frontier)
    empty(search.explored)
    node = make_node(problem.initial_state)
    insert(search.frontier, node)

    while(true)
        isempty(search.frontier) && return failure(node)
        node = pop(search.frontier)
        goal_test(problem, node.state) && return solution(node)
        append(search.explored, node.state)

        for action in actions(problem, node.state)
            child = child_node(problem, node, action)
            if !(child.state in search.explored) && !(child.state in search.frontier)
                insert(search.frontier, child)
            elseif (child.state in search.frontier) &&
                    (child.path_cost < get(search.frontier, child.state,
                                           make_node(problem.initial_state)).path_cost)
                replace(search.frontier, child)
            end
        end
    end
end

"""
*DepthLimitedSearch* An uninformed graph search technique to reach goal.

This is a depth first search limited by the depth of the search.

*limit* depth to which search should be conducted.

pg. 88 Fig. 3.17 AIMA 3ed
"""
struct DepthLimitedSearch <: SearchAlgorithm
    limit::Int
end

execute(search::DepthLimitedSearch, problem) =
    recursive_DLS(make_node(problem.initial_state), problem, search.limit)

function recursive_DLS(node, problem, limit)
    goal_test(problem, node.state) && return solution(node)
    limit == 0 && return :cutoff
    cutoff_occured = false
    for action in actions(problem, node.state)
        child = child_node(problem, node, action)
        result = recursive_DLS(child, problem, limit-1)
        result != :failure && result != :cutoff && return result
        cutoff_occured = (result == :cutoff)
    end
    return cutoff_occured ? :cutoff : :failure
end

"""
*IterativeDeepeningSearch* An uninformed graph search technique to reach goal.

Iteratively *DepthLimitedSearch* is used by incrementing the depth limit till a point comes
when there are no `:cutoff` return value.

pg. 89 Fig. 3.18 AIMA 3ed
"""
struct IterativeDeepeningSearch <: SearchAlgorithm
end

function execute(search::IterativeDeepeningSearch, problem)
    depth = 0
    while true
        result = execute(DepthLimitedSearch(depth), problem)
        result != :cutoff && return result
        depth += 1
    end
end

mutable struct GraphSearch{SQ, SET, S<:State, NT} <: SearchAlgorithm
    frontier::SQ
    explored::SET
    GraphSearch{SQ, SET, S, NT}() where {SQ, SET, S<:State, NT} =
        has_trait_queue(SQ{S}, NodeT{S, NT}) &&
        method_exists(insert_or_replace, (Function, SQ{S}, NodeT{S, NT})) &&
        new(SQ{S}(), has_trait_set(SET{S}()))
end

const GraphSearchBreadth{S <: State} = GraphSearch{AIMASearchFIFO, AIMASearchSet, S, :none}
const GraphSearchDepth{S <: State}   = GraphSearch{AIMASearchLIFO, AIMASearchSet, S, :none}
const GraphSearchUniformCost{S <: State} =
    GraphSearch{AIMASearchPQ, AIMASearchSet, S, :none}

const GraphSearchBestFirst{S <: State} =
    GraphSearch{AIMASearchPQ_G, AIMASearchSet, S, :greedy}
const GraphSearchAStar{S <: State} =
        GraphSearch{AIMASearchPQ_A, AIMASearchSet, S, :astar}

GraphSearchBreadth{S<:State}(::S)      = GraphSearchBreadth{S}()
GraphSearchDepth{S<:State}(::S)        = GraphSearchDepth{S}()
GraphSearchUniformCost{S<:State}(::S)  = GraphSearchUniformCost{S}()
GraphSearchBestFirst{S<:State}(::S)    = GraphSearchBestFirst{S}()
GraphSearchAStar{S<:State}(::S)        = GraphSearchAStar{S}()

function execute{SQ, SET, S<:State, T}(search::GraphSearch{SQ, SET, S, T}, problem)
    empty(search.frontier)
    empty(search.explored)
    node = NodeT(problem.initial_state, T)
    insert(search.frontier, node)

    while(true)
        isempty(search.frontier) && return failure(node)
        node = pop(search.frontier)
        goal_test(problem, node.state) && return solution(node)
        append(search.explored, node.state)
        for action in actions(problem, node.state)
            child = child_node(problem, node, action)
            if !(child.state in search.explored)
                insert_or_replace(x-> child < x, search.frontier, child)
            end
        end
    end
end

const AIMASearchSequence{S<:State} = Vector{Node{S}}
sort(sq::AIMASearchSequence) = sort!(sq)

struct RecursiveBestFirstSearch{SQ, S<:State} <: SearchAlgorithm
    SQ_t::Type
    RecursiveBestFirstSearch{SQ, S}() where {SQ, S<:State} =
        has_trait_sequence(SQ{S}, Node{S}) && new(Type(SQ{S}))
end

RecursiveBestFirstSearch{S<:State}(::S) = RecursiveBestFirstSearch{AIMASearchSequence, S}()

function execute(search::RecursiveBestFirstSearch, problem)
    node = make_node(problem.initial_state)
    status, val = RBFS(search, problem, node, Infinite(node.f))
    return status == :failure ? :failure : val
end

function RBFS(search, problem, node, f_limit)
    goal_test(problem, node.state) && return :success, solution(node)
    successors = search.SQ_t()

    for action in actions(problem, node.state)
        child = child_node(problem, node, action)
        insert(successors, child)
    end
    isempty(successors) && return :failure, Infinite(node.f)

    for s in successors
        s.f = max(s.path_cost + heuristic(problem, s.state), node.f)
    end

    while(true)
        sort(successors)
        best = successors[1]
        best.f > f_limit && return :failure, best.f
        alternative = length(successors) > 1 ? successors[2].f : Infinite(node.f)
        result = RBFS(search, problem, best, min(f_limit, alternative))
        result[1] == :success && return result
        best.f = result[2]
    end
end

struct StateNode{S<:State}
    state::S
    value::Number
end

StateNode(state::S) where S <: State = StateNode(state, 0)

isless(n1::StateNode, n2::StateNode)  = isless(n1.value, n2.value)
==(n1::StateNode, n2::StateNode) = n1.value == n2.value

successor_node(problem, state::State) = StateNode(state, state_value(problem, state))

const AIMASearchStateSequence{S<:State} = Vector{StateNode{S}}
sort(sq::AIMASearchStateSequence) = sort!(sq, rev=true)

struct HillClimbingSearch{SQ, S<:State} <: SearchAlgorithm
    SQ_t::Type
    plateau::Int
    HillClimbingSearch{SQ, S}(plateau::Int=typemax(Int)) where {SQ, S<:State} =
        has_trait_sequence(SQ{S}, StateNode{S}) && new(Type(SQ{S}), plateau)
end

HillClimbingSearch{S<:State}(::S, plateau::Int=typemax(Int)) =
    HillClimbingSearch{AIMASearchStateSequence, S}(plateau)

function execute(search::HillClimbingSearch, problem)
    node = StateNode(problem.initial_state, state_value(problem, problem.initial_state))
    cnt = 0
    while true
        successors = search.SQ_t()
        for state in successor_states(problem, node.state)
            successor = successor_node(problem, state)
            insert(successors, successor)
        end
        sort(successors)
        neighbor = successors[1]
        neighbor < node && return node.state
        if neighbor == node
            cnt >= search.plateau && return node.state
            cnt += 1
        elseif neighbor > node
            cnt = 0
        end
        node = neighbor
    end
end

struct SimulatedAnnealingSearch{SQ, S<:State} <: SearchAlgorithm
    SQ_t::Type
    schedule::Function
    SimulatedAnnealingSearch{SQ, S}(schedule::Function) where {SQ, S<:State} =
        has_trait_sequence(SQ{S}, Node{S}) && new(Type(SQ{S}), schedule)
end

SimulatedAnnealingSearch{S<:State}(::S,
    schedule::Function=(x-> x > 10000.0 ? 0.0 : 1.0/x)) =
    SimulatedAnnealingSearch{AIMASearchStateSequence, S}(schedule)

function execute(search::SimulatedAnnealingSearch, problem)
    node = StateNode(problem.initial_state, state_value(problem, problem.initial_state))
    schedule = search.schedule
    t = 1.0
    while true
        T = schedule(t)
        T <= 0.0000001 && return node.state
        successors = search.SQ_t()
        for state in successor_states(problem, node.state)
            successor = successor_node(problem, state)
            insert(successors, successor)
        end
        next = rand(successors)
        de = next.value - node.value
        p = de > 0 ? 1.0 : exp(de/T)
        if rand(Float64) < p
            node = next
        end
        t += 1.0
    end
end

mutable struct Population{SQN, SQF}
    nodes::SQN
    cweights::SQF
    Population{SQN, SQF}() where {SQN, SQF} = new(SQN(), SQF())
end

length(p::Population) = length(p.nodes)

function add(p::Population{SQN, SQF}, s::S, v::Number) where {SQN, SQF, S <:State }
    append(p.nodes, StateNode(s, v))
    append(p.cweights, isempty(p.cweights) ? Float64(v) : Float64(v) + p.cweights[end])
end

function random_selection(p::Population)
    wt = rand()*p.cweights[end]
    i = findfirst(x -> x > wt, p.cweights)
    return p.nodes[i].state
end

function initialize(p::Population{SQN, SQF},
    st::S, fitnessFN::Function, size::Int) where {SQN, SQF, S <:State}
    p.nodes = SQN()
    p.cweights = SQF()
    for i = 1:size
        s = random_state(st)
        add(p, s, fitnessFN(s))
    end
    return p
end

best_individual(p::Population) = findmax(p.nodes)

mutate(s::S) where {S <: State} = error(E_ABSTRACT)
reproduce(s::S) where {S <: State} = error(E_ABSTRACT)
random_state(s::S) where {S <: State} = error(E_ABSTRACT)

struct GeneticAlgorithmSearch{SQ, S<:State} <: SearchAlgorithm
    SQ_t::Type
    fitnessFN::Function
    niter::Int
    pop_size::Int
    p_mutate::Float64
    GeneticAlgorithmSearch{SQ, S}(fitnessFN::Function,
        niter::Int=1000, pop_size::Int=1000, p_mutate::Float64=0.1) where {SQ, S<:State} =
        has_trait_sequence(SQ{S}, S) &&
            new(Type(SQ), fitnessFN, niter, pop_size, p_mutate)
end

GeneticAlgorithmSearch{S<:State}(::S, fitnessFN::Function) =
    GeneticAlgorithmSearch{AIMASequence, S}(fitnessFN)

function execute(search::GeneticAlgorithmSearch, problem)
    SQN = search.SQ_t{StateNode}
    SQF = search.SQ_t{Float64}
    population = Population{SQN, SQF}()
    initialize(population, problem.initial_state, search.fitnessFN, search.pop_size)

    citer = 0
    while true
        new_population = Population{SQN, SQF}()
        for i = 1:length(population)
            x = random_selection(population)
            y = random_selection(population)
            child = reproduce(x , y)
            if (rand() < search.p_mutate)
                child = mutate(child)
            end
            add(new_population, child, search.fitnessFN(child))
        end
        population = new_population
        citer += 1
        result, idx = best_individual(population)
        if (citer > search.niter || goal_test(problem, result.state))
            return result.state
        end
    end
end
struct AND_ORGraphSearch{SQ,S<:State} <: SearchAlgorithm
    SQ_t::Type
    #AND_Search::SearchAlgorithm
    AND_ORGraphSearch{SQ,S}() where {SQ,S<:State} =

end
function OR_Search(::


text:
using LinearAlgebra
using Test

using CompScienceMeshes
using SauterSchwabQuadrature
using StaticArrays

pI = point(1,5,3)
pII = point(2,5,3)
pIII = point(7,1,0)
pIV = point(5,1,-3)

Sourcechart = simplex(pI,pIII,pII)
Testchart = simplex(pI,pIV,pII)

Accuracy = 12
ce = CommonEdge(SauterSchwabQuadrature._legendre(Accuracy,0.0,1.0))

function integrand(x,y)
			return(((x-pI)'*(y-pII))*exp(-im*1*norm(x-y))/(4pi*norm(x-y)))
end

function INTEGRAND(û,v̂)
	n1 = neighborhood(Testchart, û)
	n2 = neighborhood(Sourcechart, v̂)
	x = cartesian(n1)
	y = cartesian(n2)
	output = integrand(x,y)*jacobian(n1)*jacobian(n2)

return(output)
end

result = sauterschwab_parameterized(INTEGRAND, ce)-
		   verifintegral2(Sourcechart, Testchart, integrand, Accuracy)

@test norm(result) < 1.e-3

kernel(x,y) = 1/norm(cartesian(x)-cartesian(y))

t1 = simplex(
	@SVector[0.180878, -0.941848, -0.283207],
	@SVector[0.0, -0.92388, -0.382683],
 	@SVector[0.0, -0.980785, -0.19509])
t2 = simplex(
	@SVector[0.180878, -0.941848, -0.283207],
	@SVector[0.0, -0.92388, -0.382683],
	@SVector[0.158174, -0.881178, -0.44554])

@test indexin(t1.vertices, t2.vertices) == [1, 2, nothing]

rt = RTRefSpace{Float64}()
igd = generate_integrand_uv(kernel, rt, rt, t1, t2)

i5 = sauterschwab_parameterized(igd, CommonEdge(SauterSchwabQuadrature._legendre(5,0.0,1.0)))
i10 = sauterschwab_parameterized(igd, CommonEdge(SauterSchwabQuadrature._legendre(10,0.0,1.0)))
i15 = sauterschwab_parameterized(igd, CommonEdge(SauterSchwabQuadrature._legendre(15,0.0,1.0)))

# brute numerical approach
q1 = quadpoints(t1, 10)
q2 = quadpoints(t2, 10)

M = N = numfunctions(rt)
iref = zero(i5)
for (x,w1) in q1
    f = rt(x)
    for (y,w2) in q2
        g = rt(y)
        G = kernel(x,y)
        ds = w1*w2
        global iref += SMatrix{M,N}([dot(f[i][1], G*g[j][1])*ds for i=1:M, j=1:N])
    end
end

include(joinpath(dirname(@__FILE__,),"numquad.jl"))
ibf = numquad(kernel, rt, rt, t1, t2, zero(i5))

@test i5  ≈ iref atol=1e-3
@test i10 ≈ iref atol=1e-3
@test i10 ≈ ibf  atol=1e-3
@test i10 ≈ i15  atol=1e-5

# Test the more (or less) singular case of the second kind kernel
function kernel2nd(x,y)

	r = cartesian(x) - cartesian(y)
	R = norm(r)

	gradgreen = - r / R^3

	@SMatrix [
		0             -gradgreen[3]  gradgreen[2]
		 gradgreen[3]             0 -gradgreen[1]
		-gradgreen[2]  gradgreen[1]             0 ]
end

igd = generate_integrand_uv(kernel2nd, rt, rt, t1, t2)
i10 = sauterschwab_parameterized(igd, CommonEdge(SauterSchwabQuadrature._legendre(10,0.0,1.0)))
i15 = sauterschwab_parameterized(igd, CommonEdge(SauterSchwabQuadrature._legendre(15,0.0,1.0)))
i20 = sauterschwab_parameterized(igd, CommonEdge(SauterSchwabQuadrature._legendre(20,0.0,1.0)))

iref = numquad(kernel2nd, rt, rt, t1, t2, zero(i15))

# # Compare to BEAST:
# tqd = CompScienceMeshes.quadpoints(rt, [t1], (12,))
# bqd = CompScienceMeshes.quadpoints(rt, [t2], (13,))
#
# SE_strategy = BEAST.WiltonSEStrategy(
#   tqd[1,1],
#   BEAST.DoubleQuadStrategy(
# 	tqd[1,1],
# 	bqd[1,1]))
#
# op = BEAST.MWDoubleLayer3D(0.0)
# z2 = zeros(3,3)
# BEAST.momintegrals!(op, rt, rt, t1, t2, z2, SE_strategy)

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ 3c57ceb0-f934-11ea-1bef-c16b272d0c6b
using Pkg, DrWatson

# ╔═╡ 3c580a7e-f934-11ea-3b0c-4d5d7d5a7382
begin
	@quickactivate "StatisticalRethinkingTuring"
	using Turing
	using StatisticalRethinking
end

# ╔═╡ 83ad95c0-f933-11ea-1552-23ac7681803c
md"## Pluto-helpers-namedtuple.jl"

# ╔═╡ 3c9b1db4-f934-11ea-3d97-f51c2acc9037
md"## snippet 2.6"

# ╔═╡ 3ca3ac90-f934-11ea-1a8e-edfa2c664f1a
@model function ppl2_0(W, L)
    p ~ Uniform(0, 1)
    W ~ Binomial(W + L, p)
end

# ╔═╡ 3caf45dc-f934-11ea-2d1d-796707d89526
m2_0t = ppl2_0(6, 3);

# ╔═╡ 3cb18732-f934-11ea-3ebf-991c7b602d12
q2_0t = quap(m2_0t)

# ╔═╡ 3cbf6b30-f934-11ea-3b10-433589e36648
begin
	quap2_0t_df = DataFrame(:μ => rand(q2_0t.distr, 4000))
	Text(precis(quap2_0t_df; io=String))
end

# ╔═╡ d1670910-2c0d-11eb-12d0-616c518e22da
begin
	post2_0t = sample(m2_0t, NUTS(), 5000)
	Text(sprint(show, "text/plain", post2_0t))
end

# ╔═╡ 0604d19a-2cd4-11eb-3797-975be539c612
QM(q2_0t)

# ╔═╡ 57af25e0-2d97-11eb-0c3e-bb88e690e5ae
q2_0t

# ╔═╡ 3ce617ae-f934-11ea-0e4a-ff1ec9ee2610
md"## End pluto_helpers-namedtuple.jl"

# ╔═╡ Cell order:
# ╟─83ad95c0-f933-11ea-1552-23ac7681803c
# ╠═3c57ceb0-f934-11ea-1bef-c16b272d0c6b
# ╠═3c580a7e-f934-11ea-3b0c-4d5d7d5a7382
# ╟─3c9b1db4-f934-11ea-3d97-f51c2acc9037
# ╠═3ca3ac90-f934-11ea-1a8e-edfa2c664f1a
# ╠═3caf45dc-f934-11ea-2d1d-796707d89526
# ╠═3cb18732-f934-11ea-3ebf-991c7b602d12
# ╠═3cbf6b30-f934-11ea-3b10-433589e36648
# ╠═d1670910-2c0d-11eb-12d0-616c518e22da
# ╠═0604d19a-2cd4-11eb-3797-975be539c612
# ╠═57af25e0-2d97-11eb-0c3e-bb88e690e5ae
# ╟─3ce617ae-f934-11ea-0e4a-ff1ec9ee2610


text:
using DeltaSugar, Test, Flux

function reverse_grad_tester(d0,d1,d2)
    return all(-(d0.weight - d1.weight) .≈ (d0.weight - d2.weight)) &&
        all(-(d0.bias - d1.bias) .≈ (d0.bias - d2.bias))
end

loss( x1, x2, m ) = sum( abs2, ( m( x1 ) - x2 ) .^ 2 )
mloss( x1, x2, m ) = sum( abs2, ( m( x1 ) - x2 ) .^ 2 ) / length(x2)

@testset "TenseLayer" begin
    a = randn( 10, 6 )
    b = randn( 10, 3 )
    model = TenseLayer( 6, 3 )
    @test all(size( model( a ) ) .== ( 10, 3 ) )

    #ensure the gradient moves
    w0 = deepcopy(model.weight)
    train_light!(a, b, model, Flux.params(model), loss, Descent(0.01))
    @test sum(abs, w0 - model.weight) > 0.0
end

@testset "TenseLayer - Reverse Gradient Check" begin
    a,b = randn( 10, 6 ), randn( 10, 3 )
    original      = TenseLayer( 6, 3 )
    model_forward = deepcopy(original)
    model_reverse = RevGrad( deepcopy(original) )
    @test all(size( model_forward( a ) ) .== ( 10, 3 ) )
    @test all(size( model_reverse( a ) ) .== ( 10, 3 ) )
    #ensure the gradient moves
    ps = Flux.params( model_forward )
    train_light!(a,b, model_forward, ps, loss, Descent(1.))
    ps = Flux.params( model_reverse )
    train_light!(a,b, model_reverse, ps, loss, Descent(1.))
    @test reverse_grad_tester(original, model_forward, model_reverse)
    #now a 2 layer model check...
    a, b = randn( 10, 6 ), randn( 10 )
    layer1, layer2 = TenseLayer( 6, 3 ), TenseLayer( 3, 1 )
    model_forward = Chain( deepcopy(layer1), deepcopy(layer2) )
    model_reverse = Chain( RevGrad(layer1), RevGrad(layer2) )
    #ensure the gradient moves
    ps = Flux.params( model_forward )
    train_light!(a,b, model_forward, ps, loss, Descent(1.))
    ps = Flux.params( model_reverse )
    train_light!(a,b, model_reverse, ps, loss, Descent(1.))
    @test reverse_grad_tester(layer1, model_forward[1], model_reverse[1])
    @test reverse_grad_tester(layer2, model_forward[2], model_reverse[2])

    #test that the chain operator works
    layer1, layer2 = TenseLayer( 6, 3 ), TenseLayer( 3, 1 )
    model_forward = Chain( deepcopy(layer1), deepcopy(layer2) )
    model_reverse = RevGrad( model_forward )
    ps = Flux.params( model_forward )
    train_light!(a,b, model_forward, ps, loss, Descent(1.))
    ps = Flux.params( model_reverse )
    train_light!(a,b, model_reverse, ps, loss, Descent(1.))
    @test reverse_grad_tester(layer1, model_forward[1], model_reverse[1])
    @test reverse_grad_tester(layer2, model_forward[2], model_reverse[2])
end

@testset "DenseLayer" begin
    a = randn( 6, 10 )
    b = randn( 3, 10 )
    model = DenseLayer( 6, 3 )
    @test all(size( model( a ) ) .== ( 3, 10 ) )
    #ensure the gradient moves
    w0 = deepcopy(model.weight)
    train_light!(a, b, model, Flux.params(model), loss, Descent(0.01))
    @test sum(abs, w0 - model.weight) > 0.0
end

@testset "DenseLayer - Reverse Gradient Check" begin
    a = randn( 6, 10 )
    b = randn( 3, 10 )
    original      = DenseLayer( 6, 3 )
    model_forward = deepcopy(original)
    model_reverse = RevGrad( deepcopy(original) )
    @test all(size( model_forward( a ) ) .== ( 3, 10 ) )
    @test all(size( model_reverse( a ) ) .== ( 3, 10 ) )
    #ensure the gradient moves
    ps = Flux.params( model_forward )
    train_light!(a,b, model_forward, ps, loss, Descent(1.))
    ps = Flux.params( model_reverse )
    train_light!(a,b, model_reverse, ps, loss, Descent(1.))
    @test reverse_grad_tester(original,model_forward,model_reverse)

    #now a 2 layer model check...
    a, b = randn( 6, 10 ), randn( 1, 10 )
    layer1, layer2 = DenseLayer( 6, 3 ), DenseLayer( 3, 1 )
    model_forward = Chain( deepcopy(layer1), deepcopy(layer2) )
    model_reverse = Chain( RevGrad(layer1), RevGrad(layer2) )
    #ensure the gradient moves
    ps = Flux.params( model_forward )
    train_light!(a,b, model_forward, ps, loss, Descent(1.))
    ps = Flux.params( model_reverse )
    train_light!(a,b, model_reverse, ps, loss, Descent(1.))
    @test reverse_grad_tester(layer1, model_forward[1], model_reverse[1])
    @test reverse_grad_tester(layer2, model_forward[2], model_reverse[2])

    #test that the chain operator works
    layer1, layer2 = DenseLayer( 6, 3 ), DenseLayer( 3, 1 )
    model_forward = Chain( deepcopy(layer1), deepcopy(layer2) )
    model_reverse = RevGrad( model_forward )
    ps = Flux.params( model_forward )
    train_light!(a,b, model_forward, ps, loss, Descent(1.))
    ps = Flux.params( model_reverse )
    train_light!(a,b, model_reverse, ps, loss, Descent(1.))
    @test reverse_grad_tester(layer1, model_forward[1], model_reverse[1])
    @test reverse_grad_tester(layer2, model_forward[2], model_reverse[2])
end

@testset "Gradient Clip Test" begin
    a = randn( 6, 10 )
    b = randn( 3, 10 )
    ϵ = 1e-2
    model = GradientClamp( DenseLayer( 6, 3 ), ϵ )
    @test all(size( model( a ) ) .== ( 3, 10 ) )
    #ensure the gradient is clipped
    ps = Flux.Params(Flux.params(model))
    gs = Flux.gradient(ps) do
        loss(a, b, model)
    end
    @test all( [ reduce(max, g) <= ϵ for (_,g) in gs.grads ])
    @test all( [ reduce(min, g) >= -ϵ for (_,g) in gs.grads ])

    ϵ = 100.0
    model = GradientClamp( DenseLayer( 6, 3 ), ϵ )
    @test all(size( model( a ) ) .== ( 3, 10 ) )
    #ensure the gradient is clipped
    ps = Flux.Params(Flux.params(model))
    gs = Flux.gradient(ps) do
        loss(a, b, model)
    end
    #println.([ g for (_,g) in gs.grads ])
    @test all( [ reduce(max, g) <= ϵ for (_,g) in gs.grads ])
    @test all( [ reduce(min, g) >= -ϵ for (_,g) in gs.grads ])
end

@testset "Ridge Test" begin
    a = randn( 10, 20 )
    b = randn( 1, 20 )
    a[5,:]  = 0.5 * b[:]
    model = L2( DenseLayer( 10, 1 ), 0.02 )
    weight_i = deepcopy( model.weight )
    @test all(size( model( a ) ) .== ( 1, 20 ) )
    ps = Flux.params( model )
    for _ in 1:100
        train_light!( a, b, model, ps, mloss, Descent( 0.01 ) )
    end
    @test sum(abs2, model.weight ) < sum(abs2, weight_i )
end

@testset "FF Chain Sugar Tests" begin
    z = PlainChain( 5,4,3,2,1 ; σ = identity )
    @test all(size(z(randn(5,10))) .== (1,10))

    x = TChain( 5,4,3,2,1 ; σ = identity )
    @test all(size(x(randn(10,5))) .== (10,1))
end

@testset "LASSO Test" begin
    """
    This one won't always work, c'est la vie...
    """
    a = randn( 10, 20 )
    b = randn( 1, 20 )
    a[5,:]  = 0.5 * b[:]
    model = L1( DenseLayer( 10, 1 ), 0.00001 )
    #model.weight[1,5] .+= 0.1
    @test all(size( model( a ) ) .== ( 1, 20 ) )
    ps = Flux.params( model )
    for _ in 1:10000
        train_light!( a, b, model, ps, mloss, Descent( 0.001 ) )
    end
    #println( model.weight )
    @test argmax( model.weight[:] ) == 5
end


text:
using DifferentialEquations
using ProgressMeter
include("Yfind.jl")

using PyCall
pygui(:qt)
using PyPlot
pygui(true)

@everywhere function PSS_function_parallel(K)
    Q_start=1e-10
    Q_end=1
    P_start=0
    P_end=.3

    n_iter_P=200
    n_iter_Q=200

    delta_P=(P_end-P_start)/n_iter_P
    delta_Q=(Q_end-Q_start)/n_iter_Q

    ArrP=P_start: delta_P:P_end
    ArrQ=Q_start: delta_Q: Q_end

    Q=ArrQ[Int(ceil(K/n_iter_Q))]
    P=ArrP[mod(K,n_iter_P)+1]
    max_hit=20
    Energy=.25
    t_end=1e4
    H=(2*Energy)^2

    Y=Yfind(Q,P,H)

    if ~isempty(Y) && Y>0

        u0=zeros(5)
        u0[1]=0 #X
        u0[2]=P #P
        u0[3]=Q #Q
        u0[4]=Y #Y
        u0[5]=0 # counter variable
        #are the initial conditions right?
        prob = ODEProblem(Eq_of_M,u0,(0., t_end))
        sol=solve(prob, Vern9(),maxiters=1e20, reltol=1e-13,abstol=1e-15,callback=cb,save_start=true,save_end=true,save_everystep=false)
        # sol=solve(prob, Vern9(),maxiters=1e20, reltol=1e-10,abstol=1e-12,callback=cb,save_start=true,save_end=true,save_everystep=false)
        # sol=solve(prob, RK4(),reltol=1e-8,abstol=1e-10,callback=cb,save_start=true,save_end=true,save_everystep=false)

        uf=zeros(5)
        uf[1]=sol[1,end] #X
        uf[2]=sol[2,end] #P
        uf[3]=sol[3,end] #Q
        uf[4]=sol[4,end] #Y
        uf[5]=0
        dH=abs(H_test(uf)-H)
        #
        # prob_test= ODEProblem(Eq_of_M,u0,(0., 10.))
        # sol_test=solve(prob_test, RK4(),maxiters=1e10 )
        # S=sign(sol_test.u[2][1])
        # if dH<1e-8 && S<0
        #     return sol[:,1:end-1][3,:],sol[:,1:end-1][2,:]
        if dH<1e-8 #&& S>0
            return sol[:,2:end-1][3,:],sol[:,2:end-1][2,:]
        end
    end
end


text:
#=
example/fft.jl

1D FFT example

Because FFT returns complex numbers, we must use T=ComplexF32 here
for LinearMap to work properly.

Because the usual `fft` is not the unitary FFT, we need a factor of N
for the adjoint.
=#

using FFTW
using Test: @test

N = 8
A = LinearMapAA(fft, y -> N*ifft(y), (N, N), (name="fft",), T=ComplexF32)
@show A[:,2]

@test Matrix(A') == Matrix(A)' # test the adjoint

# timing tests to confirm that LinearMapAA does not have overhead

using BenchmarkTools

N = 2^11
L = LinearMap{ComplexF32}(fft, y -> N*ifft(y), N, N)
A = LinearMapAA(fft, y -> N*ifft(y), (N, N), (name="fft",), T=ComplexF32)

x = rand(N)
y = rand(Float32, N) + 1im * rand(Float32, N)
y2 = rand(Float32, N) + 1im * rand(Float32, N)

if true
	mul!(y, L, x)
	mul!(y2, A, x)
	@test isapprox(y, y2)
end

if true # essentially identical
	@btime x = ($L)' * $y
	@btime x = ($A)' * $y
end

if true # essentially identical
	@btime y = $L * $x
	@btime y = $A * $x
end

if true # essentially identical
	@btime mul!($y, $L, $x)
	@btime mul!($y, $A, $x)
end

nothing


text:
module MemorableUniqueIdentifier

export create, animal, validate, difficulty, search, mine_until, MiningResult

import JSON
import SHA
using Pkg.Artifacts

# Load the animals
corpus = JSON.parsefile(joinpath(artifact"animals", "animals.json"), dicttype=Dict, inttype=Int8)

"""
    bhash(key)

Hash the passed key using SHA256 and hex encoded
the result and return the first 32 bytes.

"""
function bhash(key)
    SubString(bytes2hex(SHA.sha256(key)), 1, 32)
end

"""
    animal("40db1854b0d8664b1f1baf3f6ddb6377")

Return the animal name from the specified key.

# Examples
```julia-repl
julia> animal("40db1854b0d8664b1f1baf3f6ddb6377")
"Eme Cod"
```
"""
function animal(key)
    search(bhash(key))
end

"""
    validate("40db1854b0d8664b1f1baf3f6ddb6377")

Check to see if the key is a memorable unique identifier

# Examples
```julia-repl
julia> validate("40db1854b0d8664b1f1baf3f6ddb6377")
true
```
"""
function validate(key)
    animal(key) !== nothing
end

"""
    difficulty("40db1854b0d8664b1f1baf3f6ddb6377")

Return the difficulty or length of the passed key, zero
if the key is not valid.

# Examples
```julia-repl
julia> difficulty("40db1854b0d8664b1f1baf3f6ddb6377")
6
```
"""
function difficulty(key)
    a = animal(key)
    if a === nothing
        0
    else
        length(replace(a, " " => ""))
    end
end

"""
    struct MiningResult
        key::String
        hash::String
        length::UInt
        pretty::String
    end

A struct containing the result of mining for a key.
"""
struct MiningResult
    "The key found from mining"
    key::String
    "The hash value of the key"
    hash::String
    "The length of the key that was mined"
    length::UInt
    "The pretty name of the key that was mined"
    pretty::String

    function MiningResult(key, hash, length, pretty)
        new(key, hash, length, pretty)
    end
end

"""
    mine_until(6, 10)

Mine new memorable unique identifiers until the quota is
statisfied, return an array of [`MiningResult`](@ref) structs.

# Examples
```julia-repl
julia> mine_until(6, 3)
3-element Array{MiningResult,1}:
 MiningResult("fd4482e12aa6f335f6c7cd7bf6e25fbd", "6a6c0d0240f64a044a57e8dcbd05b30f", 0x0000000000000006, "Hah Cod")
 MiningResult("dd540c82c49a983254c341b8e75f6290", "a1bf14daf63cc7e6cd4b0d17925b6320", 0x0000000000000006, "Alb Fly")
 MiningResult("4f97d7d6b65c8d4debe0baf3ec7695ad", "10bf08107ba1414cc4882f4e1b61277e", 0x0000000000000006, "Lob Fox")
```
"""
function mine_until(difficulty, quota)
    results = MiningResult[]
    while length(results) < quota
        key = string(rand(UInt128), base=16)
        hashed = bhash(key)

        short = SubString(String(hashed), 1, difficulty)

        if haskey(corpus, short)
            push!(results, MiningResult(key, hashed, difficulty, pretty(hashed, corpus[short]...)))
        end
    end
    results
end

"""
    create(6)

Create a new memorable unique identifier with the specified
difficulty.  Returns a [`MiningResult`](@ref) struct.

# Examples
```julia-repl
julia> create(6)
MiningResult("6278483d0d444cf841e500ace6d5d9cb", "0ddd09f61175c875bc53953759683e9c", 0x0000000000000006, "Odd Dog")
```
"""
function create(difficulty)
    mine_until(difficulty, 1)[1]
end

"""
    search(code)

Return the spirit animal from the given public identity.

# Examples
```julia-repl
julia> search("e3ec0d444d8d4b7ed9762380de7618de")
"Eme Cod"
```
"""
function search(code)
    for k = 15:-1:6
        short = SubString(code, 1, k)
        if haskey(corpus, short)
            return pretty(short, corpus[short]...)
        end
    end
    return nothing
end

"""
    pretty(code, length1, length2)

Return the animal name of code using the lengths from the corpus

"""
function pretty(code, k1, k2)
    w1 = to_readable_hex(SubString(code, 1, k1))
    w2 = to_readable_hex(SubString(code, k1 + 1, k2 + k1))
    titlecase(join([w1, w2], " "))
end

# A lookup table from hex values to characters
to_hex_dict = Dict(
"0" => "o",
"1" => "l",
"2" => "z",
"3" => "m",
"4" => "y",
"5" => "s",
"6" => "h",
"7" => "t",
"8" => "x",
"9" => "g")

from_hex_dict = Dict(value => key for (key, value) in to_hex_dict)

to_readable_hex(word) = replace(word, r"[0-9]" => c -> to_hex_dict[c])
from_readable_hex(word) = replace(word, r"[olzmyshtxg]", c -> from_hex_dict[c])

end # module

max_stars_repo_name: rustyconover/muid-julia

max_issues_repo_name: rustyconover/muid-julia

max_forks_repo_name: rustyconover/MemorableUniqueIdentifier


text:
using Pkg

Pkg.update()
Pkg.add(Pkg.PackageSpec(url="https://github.com/kandread/Sads.jl"))
Pkg.add("PackageCompiler")
Pkg.add("Distributions")
Pkg.add("NCDatasets")