text:
# ---
# title: 429. N-ary Tree Level Order Traversal
# id: problem429
# author: zhwang
# date: 2022-03-12
# difficulty: Medium
# categories: Tree, Breadth-first Search
# link: <https://leetcode.com/problems/n-ary-tree-level-order-traversal/description/>
# hidden: true
# ---
#
# Given an n-ary tree, return the _level order_ traversal of its nodes' values.
#
# _Nary-Tree input serialization  is represented in their level order traversal,
# each group of children is separated by the null value (See examples)._
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2018/10/12/narytreeexample.png)
#
#
#
#     Input: root = [1,null,3,2,4,null,5,6]
#     Output: [[1],[3,2,4],[5,6]]
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2019/11/08/sample_4_964.png)
#
#
#
#     Input: root = [1,null,2,3,4,5,null,null,6,7,null,8,null,9,10,null,null,11,null,12,null,13,null,null,14]
#     Output: [[1],[2,3,4,5],[6,7,8,9,10],[11,12,13],[14]]
#
#
#
#
# **Constraints:**
#
#   * The height of the n-ary tree is less than or equal to `1000`
#   * The total number of nodes is between `[0, 104]`
#
#
## @lc code=start
using LeetCode

## @lc code=end

## @lc test=start
@testset "429.n-ary-tree-level-order-traversal.jl" begin
    @test
end
## @lc test=end

## @lc info=start
# link: [solution to 429]()
## @lc info=end


text:
###
### Gibbs samplers / compositional samplers.
###

"""
    isgibbscomponent(alg)

Determine whether algorithm `alg` is allowed as a Gibbs component.
"""
isgibbscomponent(alg) = false

isgibbscomponent(::ESS) = true
isgibbscomponent(::GibbsConditional) = true
isgibbscomponent(::Hamiltonian) = true
isgibbscomponent(::MH) = true
isgibbscomponent(::PG) = true

"""
    Gibbs(algs...)

Compositional MCMC interface. Gibbs sampling combines one or more
sampling algorithms, each of which samples from a different set of
variables in a model.

Example:
```julia
@model gibbs_example(x) = begin
    v1 ~ Normal(0,1)
    v2 ~ Categorical(5)
end
```

# Use PG for a 'v2' variable, and use HMC for the 'v1' variable.
# Note that v2 is discrete, so the PG sampler is more appropriate
# than is HMC.
alg = Gibbs(HMC(0.2, 3, :v1), PG(20, :v2))
```

Tips:
- `HMC` and `NUTS` are fast samplers, and can throw off particle-based
methods like Particle Gibbs. You can increase the effectiveness of particle sampling by including
more particles in the particle sampler.
"""
struct Gibbs{space, A<:Tuple} <: InferenceAlgorithm
    algs::A   # component sampling algorithms

    function Gibbs{space, A}(algs::A) where {space, A<:Tuple}
        all(isgibbscomponent, algs) || error("all algorithms have to support Gibbs sampling")
        return new{space, A}(algs)
    end
end

function Gibbs(algs...)
    # obtain space of sampling algorithms
    space = Tuple(union(getspace.(algs)...))

    Gibbs{space, typeof(algs)}(algs)
end

"""
    GibbsState{V<:VarInfo, S<:Tuple{Vararg{Sampler}}}

Stores a `VarInfo` for use in sampling, and a `Tuple` of `Samplers` that
the `Gibbs` sampler iterates through for each `step!`.
"""
struct GibbsState{V<:VarInfo,S<:Tuple{Vararg{Sampler}},T}
    vi::V
    samplers::S
    states::T
end

struct GibbsTransition{T,F}
    "The parameters for any given sample."
    θ::T
    "The joint log probability for the sample's parameters."
    lp::F
end

function GibbsTransition(vi::AbstractVarInfo)
    theta = tonamedtuple(vi)
    lp = getlogp(vi)
    return GibbsTransition(theta, lp)
end

metadata(t::GibbsTransition) = (lp = t.lp,)

DynamicPPL.getlogp(t::GibbsTransition) = t.lp

# extract varinfo object from state
"""
    gibbs_varinfo(model, sampler, state)

Return the variables corresponding to the current `state` of the Gibbs component `sampler`.
"""
gibbs_varinfo(model, sampler, state) = varinfo(state)
varinfo(state) = state.vi
varinfo(state::AbstractVarInfo) = state

"""
    gibbs_state(model, sampler, state, varinfo)

Return an updated state, taking into account the variables sampled by other Gibbs components.

# Arguments
- `model`: model targeted by the Gibbs sampler.
- `sampler`: the sampler for this Gibbs component.
- `state`: the state of `sampler` computed in the previous iteration.
- `varinfo`: the variables, including the ones sampled by other Gibbs components.
"""
gibbs_state(model, sampler, state::AbstractVarInfo, varinfo::AbstractVarInfo) = varinfo

# Update state in Gibbs sampling
function gibbs_state(
    model::Model,
    spl::Sampler{<:Hamiltonian},
    state::HMCState,
    varinfo::AbstractVarInfo,
)
    # Update hamiltonian
    θ_old = varinfo[spl]
    hamiltonian = get_hamiltonian(model, spl, varinfo, state, length(θ_old))

    # TODO: Avoid mutation
    resize!(state.z.θ, length(θ_old))
    state.z.θ .= θ_old
    z = state.z

    return HMCState(varinfo, state.i, state.traj, hamiltonian, z, state.adaptor)
end

"""
    gibbs_rerun(prev_alg, alg)

Check if the model should be rerun to recompute the log density before sampling with the
Gibbs component `alg` and after sampling from Gibbs component `prev_alg`.

By default, the function returns `true`.
"""
gibbs_rerun(prev_alg, alg) = true

# `vi.logp` already contains the log joint probability if the previous sampler
# used a `GibbsConditional` or one of the standard `Hamiltonian` algorithms
gibbs_rerun(::GibbsConditional, ::MH) = false
gibbs_rerun(::Hamiltonian, ::MH) = false

# `vi.logp` already contains the log joint probability if the previous sampler
# used a `GibbsConditional` or a `MH` algorithm
gibbs_rerun(::MH, ::Hamiltonian) = false
gibbs_rerun(::GibbsConditional, ::Hamiltonian) = false

# do not have to recompute `vi.logp` since it is not used in `step`
gibbs_rerun(prev_alg, ::GibbsConditional) = false

# Do not recompute `vi.logp` since it is reset anyway in `step`
gibbs_rerun(prev_alg, ::PG) = false

# Initialize the Gibbs sampler.
function DynamicPPL.initialstep(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:Gibbs},
    vi::AbstractVarInfo;
    kwargs...
)
    # Create tuple of samplers
    algs = spl.alg.algs
    i = 0
    samplers = map(algs) do alg
        i += 1
        if i == 1
            prev_alg = algs[end]
        else
            prev_alg = algs[i-1]
        end
        rerun = gibbs_rerun(prev_alg, alg)
        selector = DynamicPPL.Selector(Symbol(typeof(alg)), rerun)
        Sampler(alg, model, selector)
    end

    # Add Gibbs to gids for all variables.
    for sym in keys(vi.metadata)
        vns = getfield(vi.metadata, sym).vns

        for vn in vns
            # update the gid for the Gibbs sampler
            DynamicPPL.updategid!(vi, vn, spl)

            # try to store each subsampler's gid in the VarInfo
            for local_spl in samplers
                DynamicPPL.updategid!(vi, vn, local_spl)
            end
        end
    end

    # Compute initial states of the local samplers.
    states = map(samplers) do local_spl
        # Recompute `vi.logp` if needed.
        if local_spl.selector.rerun
            model(rng, vi, local_spl)
        end

        # Compute initial state.
        _, state = DynamicPPL.initialstep(rng, model, local_spl, vi; kwargs...)

        # Update `VarInfo` object.
        vi = gibbs_varinfo(model, local_spl, state)

        return state
    end

    # Compute initial transition and state.
    transition = GibbsTransition(vi)
    state = GibbsState(vi, samplers, states)

    return transition, state
end

# Subsequent steps
function AbstractMCMC.step(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:Gibbs},
    state::GibbsState;
    kwargs...
)
    # Iterate through each of the samplers.
    vi = state.vi
    samplers = state.samplers
    states = map(samplers, state.states) do _sampler, _state
        # Recompute `vi.logp` if needed.
        if _sampler.selector.rerun
            model(rng, vi, _sampler)
        end

        # Update state of current sampler with updated `VarInfo` object.
        current_state = gibbs_state(model, _sampler, _state, vi)

        # Step through the local sampler.
        _, newstate = AbstractMCMC.step(rng, model, _sampler, current_state; kwargs...)

        # Update `VarInfo` object.
        vi = gibbs_varinfo(model, _sampler, newstate)

        return newstate
    end

    return GibbsTransition(vi), GibbsState(vi, samplers, states)
end

ext: jl


text:
export
    LowLevelSearchHeuristic,
    NullHeuristic,
    PerfectHeuristic,
    DefaultPerfectHeuristic,
    MultiStagePerfectHeuristic,
    ConflictTableHeuristic,
    HardConflictHeuristic,
    SoftConflictHeuristic,
    CompositeHeuristic,
        construct_composite_heuristic

"""
    LowLevelSearchHeuristic{C} <: AbstractCostModel{C}

Abstract type of a heuristic that returns cost of type `C`.
"""
abstract type LowLevelSearchHeuristic{C} <: AbstractCostModel{C} end
# function get_heuristic_cost(env::AbstractLowLevelEnv,m::AbstractCostModel,args...)
#     get_heuristic_cost(m,args...)
# end
################################################################################
############################## CompositeHeuristic ##############################
################################################################################
struct CompositeHeuristic{M<:Tuple,T<:Tuple} <: LowLevelSearchHeuristic{T}
    cost_models::M
end
function construct_composite_heuristic(args...)
    models = Tuple(args)
    for m in models
        @assert typeof(m) <: LowLevelSearchHeuristic
    end
    cost_types = map(m->cost_type(m),models)
    CompositeHeuristic{typeof(models),Tuple{cost_types...}}(models)
end
# cost_type(h::H) where {T,M,H<:CompositeHeuristic{T,M}} = T
function get_heuristic_cost(model::H,args...) where {T,M,H<:CompositeHeuristic{M,T}}
    T(map(h->get_heuristic_cost(h,args...), model.cost_models))
end
# function get_heuristic_cost(env,model::H,args...) where {T,M,H<:CompositeHeuristic{M,T}}
#     T(map(m->get_heuristic_cost(env,m,args...), model.cost_models))
# end

################################################################################
################################ NullHeuristic #################################
################################################################################
"""
    `NullHeuristic`
"""
struct NullHeuristic <: LowLevelSearchHeuristic{Float64} end
get_heuristic_cost(h::NullHeuristic,args...) = 0.0
get_heuristic_model(env) = NullHeuristic()
# get_heuristic_cost(h::NullHeuristic,env::AbstractLowLevelEnv,args...) = get_heuristic_cost(h,args...)

################################################################################
############################### PerfectHeuristic ###############################
################################################################################
"""
    `PerfectHeuristic`

    The Perfect Heuristic stores the exact distance between any vertex in a
    graph to all goal vertices specified during construction. The distances are
    stored in `dists`, a dictionary which maps a goal vertex `v` to a vector of
    distances from any other vertex in `1:nv(G)` to `v`.

    An example of how to access the distance from vertex `v1` to goal `v2`:
    `get_heuristic_cost(h,v1,v2) = h.dists[v2][v1]`
"""
@with_kw struct PerfectHeuristic{F} <: LowLevelSearchHeuristic{Float64}
    # dists::Dict{Int,Vector{Float64}} = Dict{Int,Vector{Float64}}()
    dists::F = (v1,v2)->0.0
end
function PerfectHeuristic(dists::Dict{Int,Vector{R}}) where {R<:Real}
    PerfectHeuristic((v1,v2)->get(get(dists,v1,Int[]),v2,0)*1.0)
end
get_heuristic_cost(h::PerfectHeuristic,goal_vtx::Int,vtx::Int) = h.dists(goal_vtx,vtx)
function PerfectHeuristic(graph,starts::Vector{Int},goals::Vector{Int})
    dists = Dict(v => gdistances(graph,v) for v in goals if 0 < v <= nv(graph))
    PerfectHeuristic(dists)
end
function PerfectHeuristic(dist_matrix::AbstractMatrix)
    dists = Dict(v => dist_matrix[:,v] for v in 1:size(dist_matrix,2))
    PerfectHeuristic(dists)
end

struct DefaultPerfectHeuristic <: LowLevelSearchHeuristic{Float64}
    h::PerfectHeuristic
end
get_heuristic_cost(h::DefaultPerfectHeuristic,goal_vtx::Int,vtx::Int) = haskey(h.h.dists,goal_vtx) ? get_heuristic_cost(h.h,goal_vtx,vtx) : 0.0

export EnvDistanceHeuristic
"""
    EnvDistanceHeuristic

A convenience struct that allows for the distance matrix to be stored in env
instead of the heuristic struct.
"""
struct EnvDistanceHeuristic <: LowLevelSearchHeuristic{Float64} end

"""
    MultiStagePerfectHeuristic

Stores multiple lookup tables corresponding to different stages of a Path-
Finding search. Each stage has a different goal. The heuristic value at a
particular stage must reflect not just the distance to the next goal but the
length of the path through all remaining goals.
"""
@with_kw struct MultiStagePerfectHeuristic <: LowLevelSearchHeuristic{Float64}
    dists::Dict{Int,Vector{Vector{Float64}}} = Dict{Int,Vector{Vector{Float64}}}()
end
get_heuristic_cost(h::MultiStagePerfectHeuristic,agent_idx::Int,stage::Int,vtx::Int) = h.dists[agent_idx][stage][vtx]
function construct_multi_stage_distance_array(G,goals)
    if length(goals) > 0
        dists = map(v->Float64.(gdistances(G,v)),goals)
        d = 0
        for i in reverse(1:length(goals)-1)
            d = d + dists[i][goals[i+1]]
            dists[i][:] = dists[i] .+ d
        end
        return dists
    end
    return Vector{Vector{Float64}}()
end
function MultiStagePerfectHeuristic(graph,goals::Vector{Vector{Int}})
    dists = Dict(idx => construct_multi_stage_distance_array(graph,g) for (idx,g) in enumerate(goals))
    MultiStagePerfectHeuristic(dists)
end

export MultiStageEnvDistanceHeuristic
"""
    MultiStageEnvDistanceHeuristic <: LowLevelSearchHeuristic{Float64}

`m.dists[i][stage]` stores the distance from that stage's goal to the final
stage's goal for agent `i`.
The heuristic cost is computed as follows:

`h = get_heuristic_cost(m.h,env,s) + cost_from_stage()`
"""
@with_kw struct MultiStageEnvDistanceHeuristic <: LowLevelSearchHeuristic{Float64}
    h::EnvDistanceHeuristic          = EnvDistanceHeuristic()
    dists::Dict{Int,Vector{Float64}} = Dict{Int,Vector{Float64}}()
end
function cost_from_stage(h::MultiStageEnvDistanceHeuristic,agent_idx::Int,stage::Int)
    h.dists[agent_idx][stage]
end
function construct_multi_stage_distance_list(env,goals)
    dists = Vector{Float64}()
    if length(goals) > 0
        d = 0
        g = goals[end]
        for v in reverse(goals)
            d = get_distance(env,v,g) + d
            push!(dists,d)
            g = v
        end
    end
    return reverse(dists)
end
function construct_multi_stage_env_distance_heuristic(env,goal_dict)
    MultiStageEnvDistanceHeuristic(
        EnvDistanceHeuristic(),
        Dict(idx => construct_multi_stage_distance_list(env,goals) for (idx,goals) in goal_dict)
    )
end

################################################################################
############################# ConflictTableHeuristic ###########################
################################################################################
"""
    ConflictTableHeuristic{T<:Union{HardConflictTable,SoftConflictTable}} <: LowLevelSearchHeuristic{Float64}

Heuristic model based on conflicts between paths.
"""
struct ConflictTableHeuristic{T<:Union{HardConflictTable,SoftConflictTable}} <: LowLevelSearchHeuristic{Float64}
    table::T
end
get_heuristic_cost(h::H,args...) where {H<:ConflictTableHeuristic} = get_conflict_value(h.table, args...)

HardConflictHeuristic(args...) = ConflictTableHeuristic(HardConflictTable(args...))
SoftConflictHeuristic(args...) = ConflictTableHeuristic(SoftConflictTable(args...))

get_time_horizon(h::H) where {H<:ConflictTableHeuristic} = get_time_horizon(h.table)
get_planned_vtx(h::H,args...) where {T<:HardConflictTable,H<:ConflictTableHeuristic}  = get_planned_vtx(h.table,args...)
reset_path!(h::H,args...) where {T<:HardConflictTable,H<:ConflictTableHeuristic{T}}   = reset_path!(h.table,args...)
set_path!(h::H,args...) where {T<:HardConflictTable,H<:ConflictTableHeuristic{T}}     = set_path!(h.table,args...)
set_path!(h::H,args...) where {H<:LowLevelSearchHeuristic} = nothing
function set_path!(h::H,args...) where {H<:CompositeHeuristic}
    for m in h.cost_models
        set_path!(m,args...)
    end
end


text:
"""
    Table([kwargs...])

`Table`s can be initialized with field-column pairs by passing the latter as
keyword arguments.

Examples:

* `Table(a = [1, 2], b = [3, 4])`
* `Table(ID = collect(1:10))`

Notes: Columns are not coerced upon `Table` initialization.
"""
function (::Type{Table})(; kwargs...)
    res = Table(Dict{Symbol,Int}(), [])
    for (k, v) in kwargs
        res[k] = v
    end
    return res
end


text:
#' # Solving Day 23 with Julia

#' [Day 23](https://adventofcode.com/2018/day/23) of the [Advent of Code](https://adventofcode.com/) this year
#' was a tough problem I had trouble solving. The basic idea is that a
#' thousand tiny nanobots, each with a 3D location in space and a range representing
#' the maximum [manhattan distance](https://en.wikipedia.org/wiki/Taxicab_geometry)
#' they could interact with, were going to teleport me as long as I was in range of the
#' nanobot that was doing the teleporting. Not knowing for sure which bot that would be,
#' participants had to find the point in 3D space that had the most overlapping nanobots' ranges.
#' My first brute-force approach was to get the cube containing all the nanobots and check
#' each point within. If two or more points have the number of nanobots' overlapping ranges, choose
#' the point closest to 0,0,0 and the answer is the manhattan distance between this point and 0,0,0.
#'
#' After realizing my input had nanobots' locations listed as pos=<-66538252,24214519,54774103>, r=94247941
#' I knew the bruteforce approach wouldn't work. We're talking about a cube thats about
BigInt(100000000)^3
#' ((100 million cubed = 1 septillion) units I'll need to check! Way too many. Some people came up with clever ways to do this. I decided to try
#' [BlackBoxOptim.jl](https://github.com/robertfeldt/BlackBoxOptim.jl).

#' ## Setup
#' After making a structure for my nanobots
struct Nanobot
    x::Int64
    y::Int64
    z::Int64
    r::Int64
end

#' and reading the input:
function getlines(path)::Array{String}
    retval = []
    open(path) do file
        for line in eachline(file)
            push!(retval, line)
        end
    end
    return retval
end
#' I parsed the lines to get positions for my nanobots:
function parse_input()::Array{Nanobot,1}
    lines = getlines("23.txt")
    re = r"^pos=<(-?\d+),(-?\d+),(-?\d+)>, r=(\d+)$"
    nanobots = []
    for line in lines
        matches = match(re, line)
        try
            parts = [parse(Int64, matches[i]) for i in 1:4]
            nanobot = Nanobot(parts...)
            push!(nanobots, nanobot)
        catch
            println(line)
        end
    end
    return nanobots
end
nanobots = parse_input()

#' Helper to calculate Manhattan distance:
function mdist(x,y,z,x2,y2,z2)::Int64
    return abs(x-x2)+abs(y-y2)+abs(z-z2)
end

#' Scoring function to get score of a point in 3D space (this actually returns a function. I didn't want
#' to be doing the nanobot parse_input() function every time):
function score()
    nanobots = parse_input()
    return (w) -> begin
        x = round(w[1])
        y = round(w[2])
        z = round(w[3])
        score = 0
        for n in nanobots
            if mdist(x,y,z,n.x,n.y,n.z) <= n.r
                score += 1
            end
        end
        return float(1000-score)+(float(mdist(x,y,z,0,0,0))/848312125.0)
    end
end
score_func = score()
#' *Note: I will be passing in a 3D point as an array of floating point numbers since
#' that is what BlackBoxOptim will do. In the function I split it into x, y, and z components
#' before rounding them to their integer values and using those. The final score I return is mostly
#' just the score of that point (the number of nanobots it is in range of) with a small part (less than 1 total)
#' to break ties so that two points with the same score will tell me the points that are closer to 0,0,0.*
#' ## The Magic
#' Next up we will be
using BlackBoxOptim
#'
#' BlackBoxOptim needs a range that it will check in. From some other code (not copied here) I found the minimum
#' and maximum values for any of the x, y, and z points for any of the nanobots, including the range they could communicate with.
#' My minimum will be
range_min = -300000000.0
#' and my maximum will be
range_max = 400000000.0
#' (again, both floating point numbers since that is what BlackBoxOptim needs)
#'
#' At this point all the hard work is done. Now you just need to make the optimizer and run it:
res = bboptimize(score_func;
    Method=:adaptive_de_rand_1_bin,
    SearchRange = (range_min,range_max),
    NumDimensions = 3,
    MaxTime = 60.0,
    TraceInterval = 10.0,
    PopulationSize=1000,
    TraceMode=:compact)
#'
#' The first argument is the function to optimize. Given a vector (in this case 3 dimensions) with each
#' value in the search range find the vector that returns a minimum score from the function.
#' The method I have here was just suggested to me by the BlackBoxOptim docs. I have no idea what it's actually doing.
#' I set the max time to 10 minutes and population size to 1000 because I was about to run off to lunch. Even though I
#' got back before the ten minutes was up I was able to see it had converged on a solution and stopped it early,
#' using the best intermediate results it had at that point.
#'
point = map(x->convert(Int64,round(x)), best_candidate(res))
println(point)
println(score_func([point[1],point[2],point[3]]))
#' Pop that through the mdist function to get the distance from 0,0,0 and I had my answer.
mdist(0,0,0,point...)


text:
function deflate_leading_zeros(ps::Vector{T}) where {T}
    ## trim any 0s from the end of ps
    N = findlast(!iszero, ps)
    K = findfirst(!iszero, ps)

    N == 0 && return(zeros(T,0), length(ps))
    ps = ps[K:N]
    ps, K-1
end

## take poly [p0, p1, ..., pn] and return
## [q_m-1, q_m-2, ..., q0], k
## where we trim of k roots of 0, and then make p monic, then reverese
## monomial x^5
function reverse_poly(ps::Vector{T}) where {T}
    # assume we have called deflate_leading_zeros
    qs = reverse(ps./ps[end])[2:end]
    qs
end

#
function quadratic_equation(a::T, b::T, c::T) where {T <: Real}
    qdrtc(a, -(0.5)*b, c)
end

## make more robust
function quadratic_equation(a::Complex{T}, b::Complex{T}, c::Complex{T}) where {T}
    d = sqrt(b^2 - 4*a*c)
    e1 = (-b + d)/(2a); e2 = (-b-d)/(2a)
    return (real(e1), imag(e1), real(e2), imag(e2))

end

## Kahan quadratic equation with fma
##  https://people.eecs.berkeley.edu/~wkahan/Qdrtcs.pdf

## solve ax^2 - 2bx + c
function qdrtc(a::T, b::T, c::T) where {T <: Real}
    # z1, z2 roots of ax^2 - 2bx + c
    d = discr(a,b,c)  # (b^2 - a*c), as 2 removes 4

    if d <= 0
        r = b/a  # real
        s = sqrt(-d)/a #imag
        return (r,s,r,-s)
    else
        r = sqrt(d) * (sign(b) + iszero(b)) + b
        return (r/a, zero(T), c/r, zero(T))
    end
end

## more work could be done here.
function discr(a::T,b::T,c::T) where {T}
    pie = 3.0 # depends on 53 or 64 bit...
    d = b*b - a*c
    e = b*b + a*c

    pie*abs(d) > e && return d

    p = b*b
    dp = muladd(b,b,-p)
    q = a*c
    dq = muladd(a,c,-q)

    (p-q) + (dp - dq)
end

##
## solve degree 2 or less case
function solve_simple_cases(ps::Vector{T}) where {T}
    S = T <: Complex ? T : Complex{T}

    N = length(ps)

    if N <= 1
        return S[]
    elseif N == 2
        return S[-ps[1]/ps[2]]
    elseif N == 3
        c,b,a = ps
        r1,i1, r2,i2 = quadratic_equation(a,b,c)
        S[complex(r1, i1), complex(r2, i2)]
    end
end


text:
# RUN: julia -e "import Brutus; Brutus.lit(:emit_lowered)" --startup-file=no %s 2>&1 | FileCheck %s

sle_int(x, y) = Base.sle_int(x, y)
emit(sle_int, Int64, Int64)

# CHECK: Core.MethodMatch(Tuple{typeof(Main.Main.sle_int), Int64, Int64}, svec(), sle_int(x, y) in Main.Main at /{{.*}}/test/Codegen/lower/sle.jl:3, true)after translating to MLIR in JLIR dialect:module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Main.sle_int), Int64, Int64}"(%arg0: !jlir<"typeof(Main.sle_int)">, %arg1: !jlir.Int64, %arg2: !jlir.Int64) -> !jlir.Bool attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     "jlir.goto"()[^bb1] : () -> ()
# CHECK-NEXT:   ^bb1:  // pred: ^bb0
# CHECK-NEXT:     %0 = "jlir.constant"() {value = #jlir<"#<intrinsic #29 sle_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %1 = "jlir.call"(%0, %arg1, %arg2) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     "jlir.return"(%1) : (!jlir.Bool) -> ()
# CHECK-NEXT:   }
# CHECK-NEXT: }

# CHECK: module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Main.sle_int), Int64, Int64}"(%arg0: !jlir<"typeof(Main.sle_int)">, %arg1: i64, %arg2: i64) -> i1 attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     %0 = cmpi sle, %arg1, %arg2 : i64
# CHECK-NEXT:     return %0 : i1
# CHECK-NEXT:   }
# CHECK-NEXT: }

# CHECK:   llvm.func @"Tuple{typeof(Main.sle_int), Int64, Int64}"(%arg0: !llvm.ptr<struct<"struct_jl_value_type", opaque>>, %arg1: i64, %arg2: i64) -> i1 attributes {llvm.emit_c_interface, sym_visibility = "nested"} {
# CHECK-NEXT:     %0 = llvm.icmp "sle" %arg1, %arg2 : i64
# CHECK-NEXT:     llvm.return %0 : i1
# CHECK-NEXT:   }
# CHECK-NEXT:   llvm.func @"_mlir_ciface_Tuple{typeof(Main.sle_int), Int64, Int64}"(%arg0: !llvm.ptr<struct<"struct_jl_value_type", opaque>>, %arg1: i64, %arg2: i64) -> i1 attributes {llvm.emit_c_interface, sym_visibility = "nested"} {
# CHECK-NEXT:     %0 = llvm.call @"Tuple{typeof(Main.sle_int), Int64, Int64}"(%arg0, %arg1, %arg2) : (!llvm.ptr<struct<"struct_jl_value_type", opaque>>, i64, i64) -> i1
# CHECK-NEXT:     llvm.return %0 : i1
# CHECK-NEXT:   }
# CHECK-NEXT: }

max_stars_repo_name: stephenneuendorffer/brutus

max_issues_repo_name: stephenneuendorffer/brutus

max_forks_repo_name: stephenneuendorffer/brutus


text:
# Autogenerated wrapper script for ONNXRuntime_jll for powerpc64le-linux-gnu-cxx11
export libonnxruntime

JLLWrappers.@generate_wrapper_header("ONNXRuntime")
JLLWrappers.@declare_library_product(libonnxruntime, "libonnxruntime.so.1.10.0")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libonnxruntime,
        "lib/libonnxruntime.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\kernelEvolutions.jl")
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\structs\\BasicStructs.jl")
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\utils\\CUDAGpuUtils.jl")
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\utils\\IterationUtils.jl")
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\utils\\ReductionUtils.jl")
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\utils\\MemoryUtils.jl")
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\distanceMetrics\\MeansMahalinobis.jl")

using ..BasicPreds, ..CUDAGpuUtils , ..MeansMahalinobis, ..IterationUtils,..ReductionUtils , ..MemoryUtils

function getExampleKernelArgs()
    nx=512 ; ny=512 ; nz=317
    #first we initialize the metrics on CPU so we will modify them easier
    goldBoolCPU= zeros(Float32,nx,ny,nz); #mimicks gold standard mask
    segmBoolCPU= ones(Float32,nx,ny,nz); #mimicks other  mask
    cartTrueGold =  CartesianIndices(zeros(3,3,5) ).+CartesianIndex(5,5,5);
    cartTrueSegm =  CartesianIndices(zeros(150,150,150) ).+CartesianIndex(100,100,100);

    goldBoolCPU[cartTrueGold].=Float32(1.0);
    segmBoolCPU[cartTrueSegm].=Float32(1.0);

    numberToLooFor= Float32(1.0);
    goldBoolGPU = CuArray(goldBoolCPU);
    segmBoolGPU = CuArray(segmBoolCPU);
    #we will fill it after we work with launch configuration
    loopXdim = UInt32(1);loopYdim = UInt32(1) ;loopZdim = UInt32(1) ;
    sizz = size(goldBoolCPU);maxX = UInt32(sizz[1]);maxY = UInt32(sizz[2]);maxZ = UInt32(sizz[3])
    threads = (32,20)
    blocks = 50
    loopXdim = UInt32(fld(maxX, threads[1]))
    loopYdim = UInt32(fld(maxY, threads[2]))
    loopZdim = UInt32(fld(maxZ,blocks ))

    totalCountGold= CuArray([0]);
    totalCountSegm= CuArray([0]);

    return (goldBoolGPU,segmBoolGPU,numberToLooFor
    ,loopYdim,loopXdim,loopZdim
    ,(maxX, maxY,maxZ)
    ,totalCountGold
    ,totalCountSegm)
end

# @testset "iter3d" begin

    args = getExampleKernelArgs()
    goldBoolGPU,segmBoolGPU,numberToLooFor,loopYdim,loopXdim,loopZdim,maxes,totalCountGold,totalCountSegm= args

    function testLoopKernelA(goldArr,segmArr
        ,numberToLooFor
        ,loopYdim::UInt32
        ,loopXdim::UInt32
        ,loopZdim::UInt32
        ,arrDims::Tuple{UInt32,UInt32,UInt32}
        ,totalCountGold,totalCountSegm)

        countGold::UInt32 = UInt32(0)
        countSegm::UInt32 = UInt32(0)
        shmemSum= @cuStaticSharedMem(UInt32, (32,2))
        clearSharedMemWarpLong(shmemSum, UInt8(2))

        @iter3d(arrDims,  loopXdim,loopYdim,loopZdim,if(  @inbounds(goldArr[x,y,z])  ==numberToLooFor)
            #updating variables needed to calculate means
            countGold+=UInt32(1)
        end)#if bool in arr

        sync_threads()

        @iter3dWithVal(segmArr,arrDims, loopXdim,loopYdim,loopZdim, if(  value  ==numberToLooFor)
            #updating variables needed to calculate means
            countSegm+=UInt32(1)
        end)#if bool in arr

        offsetIter = UInt8(1)
        @redWitAct(offsetIter,shmemSum,  countGold,+, countSegm,+  )
        @addAtomic(shmemSum,totalCountGold, totalCountSegm)

        return
    end

    @cuda threads=(32,20) blocks=50 testLoopKernelA(args...)

    @test totalCountGold[1]==3*3*5
    nx=512 ; ny=512 ; nz=317
    @test totalCountSegm[1]==(nx*ny*nz)

    using Revise, Parameters, Logging, Test
    using CUDA
    includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
    using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils, ..BasicStructs
    using ..ResultListUtils, ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils, ..ScanForDuplicates
    using ..MainOverlap, ..RandIndex , ..ProbabilisticMetrics , ..VolumeMetric ,..InformationTheorhetic
    using ..CUDAAtomicUtils, ..TpfpfnKernel, ..InterClassCorrKernel

    dimx = 529
    dimy = 556
    dimz = 339
    # dimx = 333
    # dimy = 33
    # dimz = 33
    arrGold = CUDA.ones(UInt8,(dimx,dimy,dimz))
    arrAlgo =  CUDA.ones(UInt8,(dimx,dimy,dimz))
    sizz= size(arrGold)
    ##### load tests ...
    conf = ConfigurtationStruct(trues(12)...)
    numberToLooFor = UInt8(1)

    args,threads,blocks,metricsTuplGlobal= TpfpfnKernel.prepareForconfusionTableMetricsNoSliceWise(numberToLooFor  ,conf)

    # iterLoop= args[9]
    # pixPerSlice = Int64(args[10])
    singleSum= CUDA.zeros(UInt64,1)
    function testIter(goldGPU#goldBoolGPU
        , segmGPU#segmBoolGPU
        ,sliceMetricsTupl
        ,tp,tn,fp,fn#tp,tn,fp,fn
        ,arrDims,totalNumbOfVoxels,iterLoop,pixPerSlice
        ,numberToLooFor#numberToLooFor
        #,metricsTuplGlobal
        ,conf,singleSum)
        @iterateLinearlyMultipleBlocks(iterLoop,pixPerSlice,totalNumbOfVoxels,begin
           CUDA.@atomic singleSum[1]+=1
            #CUDA.@cuprint "i $(i) \n"
        end)

    return
    end

    @cuda threads=threads blocks=blocks testIter(arrGold, arrAlgo,args..., singleSum) #args[8][3]  is number of slices ...

    @test Int64(singleSum[1]) ==dimx*dimy*dimz

    Int64(singleSum[1])-dimx*dimy*dimz

# end#testset


text:
# packing and unpacking of Graph related types

using IncrementalInference
using DistributedFactorGraphs
using Manifolds

using Test

##

@testset "Serialization of SamplableBelief types" begin
##

td = Uniform()

ptd = convert(String, td) # TODO, PackedSamplableBelief
utd = convert(SamplableBelief, td)

@test td.a - utd.a |> abs < 1e-10
@test td.b - utd.b |> abs < 1e-10

##
end
##

fg = initfg()
N=100

# doors = zeros(1,4)
T = ContinuousScalar
doors = [[-100.0;],[0.0;],[100.0;],[300.0;]]
pd = manikde!(T, doors, [3.0;])
pd = resample(pd,N);
bws = getBW(pd)[:,1]
doors2 = getPoints(pd);
v1 = addVariable!(fg,:x1, ContinuousScalar,N=N)
f1  = addFactor!(fg,[:x1], Prior(manikde!(T, doors2, bws)))

v2 = addVariable!(fg,:x2, ContinuousScalar, N=N)
lc = LinearRelative(  Normal(50.0, 2.0) )
f2 = addFactor!(fg, [:x1; :x2], lc)

##

@testset "Testing conversion to packed function node data structure and back" begin
##

topack = getSolverData(f1)
dd = convert(PackedFunctionNodeData{PackedPrior},topack)
upd = reconstFactorData(fg, [:x1;], FunctionNodeData{CommonConvWrapper{Prior}}, dd)

@test compare(topack, upd)

topack = getSolverData(f2)
dd =  convert(IncrementalInference.PackedFunctionNodeData{PackedLinearRelative},topack)
upd = reconstFactorData(fg, [:x1;:x2], IncrementalInference.FunctionNodeData{CommonConvWrapper{LinearRelative}}, dd)

@test compare(topack, upd)

##
end

@testset "Testing conversion to packed variable node data structure and back" begin
##

dat = getSolverData(getVariable(fg,:x1))

# dat.BayesNetVertID

pd = packVariableNodeData(fg, dat)
unpckd = unpackVariableNodeData(fg, pd)

@test compareFields(dat, unpckd, skip=[:variableType])
@test compareFields(getVariableType(dat), getVariableType(unpckd))
@test isa(getVariableType(dat), ContinuousScalar)
@test isa(getVariableType(unpckd), ContinuousScalar)

##
end

@testset "test serialization of ManifoldKernelDensity" begin
##

# create a basic manifoldkerneldensity
mkd = manikde!(TranslationGroup(2), [randn(2) for _ in 1:100])

# convert up and down
st = convert(String, mkd) # TODO, PackedSamplableBelief
upk = convert(SamplableBelief, st)

# and check the basics
@test isapprox( getPoints(mkd)[1], getPoints(upk)[1])
@test isapprox( getPoints(mkd)[end], getPoints(upk)[end])

@test mkd.manifold == upk.manifold
@test mkd._partial == upk._partial
@test mkd.infoPerCoord == upk.infoPerCoord

##
end

#


text:
using PointSampler
using Test

#%%
P1 = collect(MonteCarloThDesign(; dims = 1, n_points = 10))
@test [0.0] ∈ P1
@test length(P1) == 10
#%%

P1_2 = collect( MonteCarloThDesign(; dims = 1, lb = [1.0], ub = [5.0], n_points = 10 ) )
@test [1.0] ∈ P1_2
@test min( P1_2... )[end] >= 1.0
@test max( P1_2... )[end] <= 5.0

#%%
seeds = [ rand(2) for i = 1:5 ]
P2_1 = collect(MonteCarloThDesign(dims = 2, n_points = 3, seeds = seeds ))
@test length( P2_1 ) == 3
@test all( seeds[1:3] .≈ P2_1 )

#%%
P2_2 = collect(MonteCarloThDesign(dims = 2, n_points = 3 ,
    lb = 2.0 .* ones(2),
    ub = 3.0 .* ones(2),
    seeds = seeds )
)   # clean_seeds = true (default)
@test !( seeds ⊆ P2_2 )

#%%
push!(seeds, [2.5; 2.5])
P2_3 = collect(MonteCarloThDesign(dims = 2, n_points = 3 ,
    lb = 2.0 .* ones(2),
    ub = 3.0 .* ones(2),
    seeds = seeds )
)
@test [2.5; 2.5] ∈ P2_3

#%%
P2_4 = collect(MonteCarloThDesign(dims = 2, n_points = 3 ,
    lb = 2.0 .* ones(2),
    ub = 3.0 .* ones(2),
    seeds = seeds,
    clean_seeds = false,
))
@test all( seeds[1:3] .≈ P2_4 )
#%%

P2_5 = collect(MonteCarloThDesign(dims = 2, n_points = 10, seeds = seeds))
@test length(P2_5) == 10
@test all( seeds[1:5] .≈ P2_5[1:5] )
@test !( seeds[6] ≈ P2_5[6])    # sixt seed should have been filtered out


text:
module Ejemplo

# Write your package code here.
export f_xy,f_loga
include("funciones.jl")

end


text:
int_rules_4_5_9 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.5.9*trig^m*(a+b*sec^n+c*sec^(2*n))^p =#
    @apply_utils Antiderivative((~(a') + ~(b') * sec(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * sec(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~(p'), ~x) => (1 / (4 ^ ~p * (~c) ^ ~p)) * Antiderivative((~b + 2 * ~c * sec(~d + ~e * ~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~p)))
    @apply_utils Antiderivative((~(a') + ~(b') * csc(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * csc(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~(p'), ~x) => (1 / (4 ^ ~p * (~c) ^ ~p)) * Antiderivative((~b + 2 * ~c * csc(~d + ~e * ~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~p)))
    @apply_utils Antiderivative((~(a') + ~(b') * sec(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * sec(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~p, ~x) => ((~a + ~b * sec(~d + ~e * ~x) ^ ~n + ~c * sec(~d + ~e * ~x) ^ (2 * ~n)) ^ ~p / (~b + 2 * ~c * sec(~d + ~e * ~x) ^ ~n) ^ (2 * ~p)) * Antiderivative(u * (~b + 2 * ~c * sec(~d + ~e * ~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative((~(a') + ~(b') * csc(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * csc(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~p, ~x) => ((~a + ~b * csc(~d + ~e * ~x) ^ ~n + ~c * csc(~d + ~e * ~x) ^ (2 * ~n)) ^ ~p / (~b + 2 * ~c * csc(~d + ~e * ~x) ^ ~n) ^ (2 * ~p)) * Antiderivative(u * (~b + 2 * ~c * csc(~d + ~e * ~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative(1 / (~(a') + ~(b') * sec(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * sec(~(d') + ~(e') * ~x) ^ ~(n2')), ~x) => Module([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((2 * ~c) / q) * Antiderivative(1 / ((~b - q) + 2 * ~c * sec(~d + ~e * ~x) ^ ~n), ~x) - ((2 * ~c) / q) * Antiderivative(1 / (~b + q + 2 * ~c * sec(~d + ~e * ~x) ^ ~n), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0))
    @apply_utils Antiderivative(1 / (~(a') + ~(b') * csc(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * csc(~(d') + ~(e') * ~x) ^ ~(n2')), ~x) => Module([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((2 * ~c) / q) * Antiderivative(1 / ((~b - q) + 2 * ~c * csc(~d + ~e * ~x) ^ ~n), ~x) - ((2 * ~c) / q) * Antiderivative(1 / (~b + q + 2 * ~c * csc(~d + ~e * ~x) ^ ~n), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0))
    @apply_utils Antiderivative(sin(~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * sec(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * sec(~(d') + ~(e') * ~x) ^ ~n2) ^ ~(p'), ~x) => Module([f = FreeFactors(cos(~d + ~e * ~x), ~x)], (-f / ~e) * Subst(Antiderivative(((1 - f ^ 2 * (~x) ^ 2) ^ ((~m - 1) / 2) * (~b + ~a * (f * ~x) ^ ~n) ^ ~p) / (f * ~x) ^ (~n * ~p), ~x), ~x, cos(~d + ~e * ~x) / f)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ(~n2, 2 * ~n) && (IntegerQ((~m - 1) / 2) && IntegersQ(~n, ~p)))
    @apply_utils Antiderivative(cos(~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * csc(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * csc(~(d') + ~(e') * ~x) ^ ~n2) ^ ~(p'), ~x) => Module([f = FreeFactors(sin(~d + ~e * ~x), ~x)], (f / ~e) * Subst(Antiderivative(((1 - f ^ 2 * (~x) ^ 2) ^ ((~m - 1) / 2) * (~b + ~a * (f * ~x) ^ ~n) ^ ~p) / (f * ~x) ^ (~n * ~p), ~x), ~x, sin(~d + ~e * ~x) / f)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ(~n2, 2 * ~n) && (IntegerQ((~m - 1) / 2) && IntegersQ(~n, ~p)))
    @apply_utils Antiderivative(sin(~(d') + ~(e') * ~x) ^ ~m * (~(a') + ~(b') * sec(~(d') + ~(e') * ~x) ^ ~n + ~(c') * sec(~(d') + ~(e') * ~x) ^ ~n2) ^ ~(p'), ~x) => Module([f = FreeFactors(tan(~d + ~e * ~x), ~x)], (f ^ (~m + 1) / ~e) * Subst(Antiderivative(((~x) ^ ~m * ExpandToSum(~a + ~b * (1 + f ^ 2 * (~x) ^ 2) ^ (~n / 2) + ~c * (1 + f ^ 2 * (~x) ^ 2) ^ ~n, ~x) ^ ~p) / (1 + f ^ 2 * (~x) ^ 2) ^ (~m / 2 + 1), ~x), ~x, tan(~d + ~e * ~x) / f)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (IntegerQ(~m / 2) && IntegerQ(~n / 2)))
    @apply_utils Antiderivative(cos(~(d') + ~(e') * ~x) ^ ~m * (~(a') + ~(b') * csc(~(d') + ~(e') * ~x) ^ ~n + ~(c') * csc(~(d') + ~(e') * ~x) ^ ~n2) ^ ~(p'), ~x) => Module([f = FreeFactors(cot(~d + ~e * ~x), ~x)], (-(f ^ (~m + 1)) / ~e) * Subst(Antiderivative(((~x) ^ ~m * ExpandToSum(~a + ~b * (1 + f ^ 2 * (~x) ^ 2) ^ (~n / 2) + ~c * (1 + f ^ 2 * (~x) ^ 2) ^ ~n, ~x) ^ ~p) / (1 + f ^ 2 * (~x) ^ 2) ^ (~m / 2 + 1), ~x), ~x, cot(~d + ~e * ~x) / f)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (IntegerQ(~m / 2) && IntegerQ(~n / 2)))
    @apply_utils Antiderivative(sec(~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * sec(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * sec(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~p, ~x) => (1 / (4 ^ ~p * (~c) ^ ~p)) * Antiderivative(sec(~d + ~e * ~x) ^ ~m * (~b + 2 * ~c * sec(~d + ~e * ~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~p)))
    @apply_utils Antiderivative(csc(~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * csc(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * csc(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~p, ~x) => (1 / (4 ^ ~p * (~c) ^ ~p)) * Antiderivative(csc(~d + ~e * ~x) ^ ~m * (~b + 2 * ~c * csc(~d + ~e * ~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~p)))
    @apply_utils Antiderivative(sec(~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * sec(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * sec(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~p, ~x) => ((~a + ~b * sec(~d + ~e * ~x) ^ ~n + ~c * sec(~d + ~e * ~x) ^ (2 * ~n)) ^ ~p / (~b + 2 * ~c * sec(~d + ~e * ~x) ^ ~n) ^ (2 * ~p)) * Antiderivative(sec(~d + ~e * ~x) ^ ~m * (~b + 2 * ~c * sec(~d + ~e * ~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative(csc(~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * csc(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * csc(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~p, ~x) => ((~a + ~b * csc(~d + ~e * ~x) ^ ~n + ~c * csc(~d + ~e * ~x) ^ (2 * ~n)) ^ ~p / (~b + 2 * ~c * csc(~d + ~e * ~x) ^ ~n) ^ (2 * ~p)) * Antiderivative(csc(~d + ~e * ~x) ^ ~m * (~b + 2 * ~c * csc(~d + ~e * ~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative(sec(~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * sec(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * sec(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~p, ~x) => Antiderivative(ExpandTrig(sec(~d + ~e * ~x) ^ ~m * (~a + ~b * sec(~d + ~e * ~x) ^ ~n + ~c * sec(~d + ~e * ~x) ^ (2 * ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ(~n2, 2 * ~n) && IntegersQ(~m, ~n, ~p))
    @apply_utils Antiderivative(csc(~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * csc(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * csc(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~p, ~x) => Antiderivative(ExpandTrig(csc(~d + ~e * ~x) ^ ~m * (~a + ~b * csc(~d + ~e * ~x) ^ ~n + ~c * csc(~d + ~e * ~x) ^ (2 * ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ(~n2, 2 * ~n) && IntegersQ(~m, ~n, ~p))
    @apply_utils Antiderivative(tan(~(d') + ~(e') * ~x) ^ ~(m') * (~a + ~(b') * sec(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * sec(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~(p'), ~x) => Module([f = FreeFactors(cos(~d + ~e * ~x), ~x)], (-1 / (~e * f ^ ((~m + ~n * ~p) - 1))) * Subst(Antiderivative(((1 - f ^ 2 * (~x) ^ 2) ^ ((~m - 1) / 2) * (~c + ~b * (f * ~x) ^ ~n + ~c * (f * ~x) ^ (2 * ~n)) ^ ~p) / (~x) ^ (~m + 2 * ~n * ~p), ~x), ~x, cos(~d + ~e * ~x) / f)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (IntegerQ((~m - 1) / 2) && (IntegerQ(~n) && IntegerQ(~p))))
    @apply_utils Antiderivative(cot(~(d') + ~(e') * ~x) ^ ~(m') * (~a + ~(b') * csc(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * sec(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~(p'), ~x) => Module([f = FreeFactors(sin(~d + ~e * ~x), ~x)], (1 / (~e * f ^ ((~m + ~n * ~p) - 1))) * Subst(Antiderivative(((1 - f ^ 2 * (~x) ^ 2) ^ ((~m - 1) / 2) * (~c + ~b * (f * ~x) ^ ~n + ~c * (f * ~x) ^ (2 * ~n)) ^ ~p) / (~x) ^ (~m + 2 * ~n * ~p), ~x), ~x, sin(~d + ~e * ~x) / f)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (IntegerQ((~m - 1) / 2) && (IntegerQ(~n) && IntegerQ(~p))))
    @apply_utils Antiderivative(tan(~(d') + ~(e') * ~x) ^ ~(m') * (~a + ~(b') * sec(~(d') + ~(e') * ~x) ^ ~n + ~(c') * sec(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~(p'), ~x) => Module([f = FreeFactors(tan(~d + ~e * ~x), ~x)], (f ^ (~m + 1) / ~e) * Subst(Antiderivative(((~x) ^ ~m * ExpandToSum(~a + ~b * (1 + f ^ 2 * (~x) ^ 2) ^ (~n / 2) + ~c * (1 + f ^ 2 * (~x) ^ 2) ^ ~n, ~x) ^ ~p) / (1 + f ^ 2 * (~x) ^ 2), ~x), ~x, tan(~d + ~e * ~x) / f)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ(~n2, 2 * ~n) && (IntegerQ(~m / 2) && IntegerQ(~n / 2)))
    @apply_utils Antiderivative(cot(~(d') + ~(e') * ~x) ^ ~(m') * (~a + ~(b') * csc(~(d') + ~(e') * ~x) ^ ~n + ~(c') * sec(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~(p'), ~x) => Module([f = FreeFactors(cot(~d + ~e * ~x), ~x)], (-(f ^ (~m + 1)) / ~e) * Subst(Antiderivative(((~x) ^ ~m * ExpandToSum(~a + ~b * (1 + f ^ 2 * (~x) ^ 2) ^ (~n / 2) + ~c * (1 + f ^ 2 * (~x) ^ 2) ^ ~n, ~x) ^ ~p) / (1 + f ^ 2 * (~x) ^ 2), ~x), ~x, cot(~d + ~e * ~x) / f)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ(~n2, 2 * ~n) && (IntegerQ(~m / 2) && IntegerQ(~n / 2)))
    @apply_utils Antiderivative((~A + ~(B') * sec(~(d') + ~(e') * ~x)) * (~a + ~(b') * sec(~(d') + ~(e') * ~x) + ~(c') * sec(~(d') + ~(e') * ~x) ^ 2) ^ ~n, ~x) => (1 / (4 ^ ~n * (~c) ^ ~n)) * Antiderivative((~A + ~B * sec(~d + ~e * ~x)) * (~b + 2 * ~c * sec(~d + ~e * ~x)) ^ (2 * ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B], ~x) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~n))
    @apply_utils Antiderivative((~A + ~(B') * csc(~(d') + ~(e') * ~x)) * (~a + ~(b') * csc(~(d') + ~(e') * ~x) + ~(c') * csc(~(d') + ~(e') * ~x) ^ 2) ^ ~n, ~x) => (1 / (4 ^ ~n * (~c) ^ ~n)) * Antiderivative((~A + ~B * csc(~d + ~e * ~x)) * (~b + 2 * ~c * csc(~d + ~e * ~x)) ^ (2 * ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B], ~x) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~n))
    @apply_utils Antiderivative((~A + ~(B') * sec(~(d') + ~(e') * ~x)) * (~a + ~(b') * sec(~(d') + ~(e') * ~x) + ~(c') * sec(~(d') + ~(e') * ~x) ^ 2) ^ ~n, ~x) => ((~a + ~b * sec(~d + ~e * ~x) + ~c * sec(~d + ~e * ~x) ^ 2) ^ ~n / (~b + 2 * ~c * sec(~d + ~e * ~x)) ^ (2 * ~n)) * Antiderivative((~A + ~B * sec(~d + ~e * ~x)) * (~b + 2 * ~c * sec(~d + ~e * ~x)) ^ (2 * ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B], ~x) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(~n)))
    @apply_utils Antiderivative((~A + ~(B') * csc(~(d') + ~(e') * ~x)) * (~a + ~(b') * csc(~(d') + ~(e') * ~x) + ~(c') * csc(~(d') + ~(e') * ~x) ^ 2) ^ ~n, ~x) => ((~a + ~b * csc(~d + ~e * ~x) + ~c * csc(~d + ~e * ~x) ^ 2) ^ ~n / (~b + 2 * ~c * csc(~d + ~e * ~x)) ^ (2 * ~n)) * Antiderivative((~A + ~B * csc(~d + ~e * ~x)) * (~b + 2 * ~c * csc(~d + ~e * ~x)) ^ (2 * ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B], ~x) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(~n)))
    @apply_utils Antiderivative((~A + ~(B') * sec(~(d') + ~(e') * ~x)) / (~(a') + ~(b') * sec(~(d') + ~(e') * ~x) + ~(c') * sec(~(d') + ~(e') * ~x) ^ 2), ~x) => Module([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], (~B + (~b * ~B - 2 * ~A * ~c) / q) * Antiderivative(1 / (~b + q + 2 * ~c * sec(~d + ~e * ~x)), ~x) + (~B - (~b * ~B - 2 * ~A * ~c) / q) * Antiderivative(1 / ((~b - q) + 2 * ~c * sec(~d + ~e * ~x)), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B], ~x) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative((~A + ~(B') * csc(~(d') + ~(e') * ~x)) / (~(a') + ~(b') * csc(~(d') + ~(e') * ~x) + ~(c') * csc(~(d') + ~(e') * ~x) ^ 2), ~x) => Module([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], (~B + (~b * ~B - 2 * ~A * ~c) / q) * Antiderivative(1 / (~b + q + 2 * ~c * csc(~d + ~e * ~x)), ~x) + (~B - (~b * ~B - 2 * ~A * ~c) / q) * Antiderivative(1 / ((~b - q) + 2 * ~c * csc(~d + ~e * ~x)), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B], ~x) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative((~A + ~(B') * sec(~(d') + ~(e') * ~x)) * (~(a') + ~(b') * sec(~(d') + ~(e') * ~x) + ~(c') * sec(~(d') + ~(e') * ~x) ^ 2) ^ ~n, ~x) => Antiderivative(ExpandTrig((~A + ~B * sec(~d + ~e * ~x)) * (~a + ~b * sec(~d + ~e * ~x) + ~c * sec(~d + ~e * ~x) ^ 2) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~n))
    @apply_utils Antiderivative((~A + ~(B') * csc(~(d') + ~(e') * ~x)) * (~(a') + ~(b') * csc(~(d') + ~(e') * ~x) + ~(c') * csc(~(d') + ~(e') * ~x) ^ 2) ^ ~n, ~x) => Antiderivative(ExpandTrig((~A + ~B * csc(~d + ~e * ~x)) * (~a + ~b * csc(~d + ~e * ~x) + ~c * csc(~d + ~e * ~x) ^ 2) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~n))
end

max_stars_repo_path: src/rules/4 Trig functions/4.5 Secant/4.5.9 trig^m (a+b sec^n+c sec^(2 n))^p.jl


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    @assert precompile(HyperRectangle{2,Float32}, (Int, Int, Int, Int))
    @assert precompile(==, (HyperRectangle{2,Float32}, HyperRectangle{2,Float32}))
    @assert precompile(normal_mesh, (Tesselation{3,Float32,Cylinder{3,Float32},1},))
    @assert precompile(normal_mesh, (Tesselation{3,Float32,HyperSphere{3,Float32},1},))
    @assert precompile(normal_mesh, (HyperSphere{3,Float32},))

    if Base.VERSION >= v"1.6.0-DEV.1083"
        @assert precompile(triangle_mesh, (Polygon{2, Float32, Point2f0, LineString{2, Float32, Point2f0,
                                        Base.ReinterpretArray{Line{2, Float32}, 1, Tuple{Point2f0, Point2f0}, TupleView{Tuple{Point2f0, Point2f0}, 2, 1, Vector{Point2f0}}, false}},
                                        Vector{LineString{2, Float32, Point2f0, Base.ReinterpretArray{Line{2, Float32}, 1, Tuple{Point2f0, Point2f0}, TupleView{Tuple{Point2f0, Point2f0}, 2, 1, Vector{Point2f0}}, false}}}},))
    else
        @assert precompile(triangle_mesh, (Polygon{2, Float32, Point2f0, LineString{2, Float32, Point2f0,
                                        Base.ReinterpretArray{Line{2, Float32}, 1, Tuple{Point2f0, Point2f0}, TupleView{Tuple{Point2f0, Point2f0}, 2, 1, Vector{Point2f0}}}},
                                        Vector{LineString{2, Float32, Point2f0, Base.ReinterpretArray{Line{2, Float32}, 1, Tuple{Point2f0, Point2f0}, TupleView{Tuple{Point2f0, Point2f0}, 2, 1, Vector{Point2f0}}}}}},))
    end

    @assert precompile(split_intersections, (Vector{Point2f0},))
end


text:
function GeodeticPoint(arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return GeodeticPoint((jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function equals(obj::GeodeticPoint, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_altitude(obj::GeodeticPoint)
    return jcall(obj, "getAltitude", jdouble, ())
end

function get_east(obj::GeodeticPoint)
    return jcall(obj, "getEast", Vector3D, ())
end

function get_latitude(obj::GeodeticPoint)
    return jcall(obj, "getLatitude", jdouble, ())
end

function get_longitude(obj::GeodeticPoint)
    return jcall(obj, "getLongitude", jdouble, ())
end

function get_nadir(obj::GeodeticPoint)
    return jcall(obj, "getNadir", Vector3D, ())
end

function get_north(obj::GeodeticPoint)
    return jcall(obj, "getNorth", Vector3D, ())
end

function get_south(obj::GeodeticPoint)
    return jcall(obj, "getSouth", Vector3D, ())
end

function get_west(obj::GeodeticPoint)
    return jcall(obj, "getWest", Vector3D, ())
end

function get_zenith(obj::GeodeticPoint)
    return jcall(obj, "getZenith", Vector3D, ())
end

function hash_code(obj::GeodeticPoint)
    return jcall(obj, "hashCode", jint, ())
end

function to_string(obj::GeodeticPoint)
    return jcall(obj, "toString", JString, ())
end


text:
export IntRangeDict, save

import Base: push!, in, show, foreach, collect, getindex, read, write

type IntRangeSpan{K<:Integer, V}
    lv::K
    rv::K
    data::Vector{V}
end

immutable IntRangeDict{K<:Integer, V}
    data::Vector{IntRangeSpan{K, V}}
    IntRangeDict(data::Vector=[]) = new(data)
    IntRangeDict(io::IO) = begin
        n = read(io, Int)
        data = Vector{IntRangeSpan{K, V}}(n)
        for i in 1:n
            lv = read(io, K)
            rv = read(io, K)
            len = read(io, Int)
            data[i] = IntRangeSpan{K, V}(lv, rv, read(io, V, len))
        end
        IntRangeDict{K, V}(data)
    end
end

immutable IntRangeHandler{K<:Integer, V}
    lv::K
    rv::K
    dict::IntRangeDict{K, V}
end

function find_last{K, V}(dict::IntRangeDict{K, V}, int::K)::Int
    findlast(x->x.lv <= int, dict.data)
end

function find_binary{K, V}(dict::IntRangeDict{K, V}, int::K)::Int
    isempty(dict.data) && return 0

    v(i) = dict.data[i].lv

    find_within(i, j) = if i == j
        v(i) <= int ? i : i - 1
    else
        m = (i + j + 1) ÷ 2
        v(m) <= int ? find_within(m, j) :
                      find_within(i, m-1)
    end

    find_within(1, length(dict.data))
end

function getindex{K, V}(dict::IntRangeDict{K, V}, range::UnitRange{K})
    IntRangeHandler{K, V}(range.start, range.stop, dict)
end

function push!{K, V}(handle::IntRangeHandler{K, V}, v::V)
    data, lv, rv = handle.dict.data, handle.lv, handle.rv

    if lv > rv
        return handle
    end

    i = find_last(handle.dict, lv) # use last rahter than binary to make inserting incremently faster

    if i == 0 || lv > data[i].rv
        if i+1 <= length(data) && rv >= data[i+1].lv
            insert!(data, i+1, IntRangeSpan{K, V}(lv, data[i+1].lv-1, [v]))
            push_aligned_left!(handle.dict, rv, v, i+2)
        else
            insert!(data, i+1, IntRangeSpan{K, V}(lv, rv, [v]))
        end
    else
        if lv == data[i].lv
            push_aligned_left!(handle.dict, rv, v, i)
        else
            old_rv = data[i].rv
            data[i].rv = lv-1

            if rv <= old_rv
                insert!(data, i+1, IntRangeSpan{K, V}(lv, rv, [data[i].data[:]; v]))

                if rv != old_rv
                    insert!(data, i+2, IntRangeSpan{K, V}(rv+1, old_rv, data[i].data[:]))
                end
            else
                insert!(data, i+1, IntRangeSpan{K, V}(lv, old_rv, [data[i].data[:]; v]))
                push_aligned_right!(handle.dict, rv, v, i+1)
            end
        end
    end

    handle.dict
end

"push when handle.lv == handle.dict.data[i].lv"
function push_aligned_left!{K, V}(dict::IntRangeDict{K, V}, rv::K, v::V, i::Int)
    p = dict.data[i]

    if rv < p.rv
        lv = p.lv
        p.lv = rv + 1
        insert!(dict.data, i, IntRangeSpan{K, V}(lv, rv, [p.data[:]; v]))
    else
        push!(p.data, v)
        if rv > p.rv
            push_aligned_right!(dict, rv, v, i)
        end
    end
end

"push when handle.lv == handle.dict.data[i].rv+1"
function push_aligned_right!{K, V}(dict::IntRangeDict{K, V}, rv::K, v::V, i::Int)
    p = dict.data[i]
    if i+1 <= length(dict.data)
        pnext = dict.data[i+1]
        if p.rv+1 == pnext.lv
            push_aligned_left!(dict, rv, v, i+1)
        elseif rv < pnext.lv
            insert!(dict.data, i+1, IntRangeSpan{K, V}(p.rv+1, rv, [v]))
        else
            insert!(dict.data, i+1, IntRangeSpan{K, V}(p.rv+1, pnext.lv-1, [v]))
            push_aligned_left!(dict, rv, v, i+2)
        end
    else
        push!(dict.data, IntRangeSpan{K, V}(p.rv+1, rv, [v]))
    end
end

function traverse{K, V}(f::Function, dict::IntRangeDict{K, V})::Void
    for i in dict.data
        f(i.lv:i.rv, i.data)
    end
end

function show{K, V}(io::IO, dict::IntRangeDict{K, V})::Void
    println(io, "IntRangeDict{$K, $V}:")
    if isempty(dict.data)
        println(io, "  (empty)")
    else
        traverse(dict) do range, data
            println(io, "  ", range.start, '-', range.stop, ": ", join(data, ','))
        end
    end
end

function getindex{K, V}(dict::IntRangeDict{K, V}, index::K)::Vector{V}
    i = find_binary(dict, index)
    i == 0 || dict.data[i].rv < index ? [] : dict.data[i].data
end

function save{K, V}(io::IO, dict::IntRangeDict{K, V})
    write(io, length(dict.data))
    for i in dict.data
        write(io, i.lv, i.rv, length(i.data), i.data)
    end
end


text:
module StructureLearner

using LogicCircuits
using ..Utils

using ..Probabilistic
using ..IO

export
# ChowLiuTree
learn_chow_liu_tree, parent_vector, print_tree, CLT,

# CircuitBuilder
compile_prob_circuit_from_clt, learn_probabilistic_circuit, BaseCache, ⊤, LitCache,

# PSDDInitializer
learn_struct_prob_circuit,
learn_vtree_from_clt, compile_psdd_from_clt

include("ChowLiuTree.jl")
include("CircuitBuilder.jl")
include("PSDDInitializer.jl")

end


text:
# Tests for NEP types functionality

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test,Random

@bench @testset "NEPTypes" begin

    n=3;
    Random.seed!(0);
A1 = rand(n,n); A2 = rand(n,n); A3 = rand(n,n); A4 = rand(n,n);
f1 = S -> one(S); f2 = S -> -S; f3 = S -> exp(-S); f4 = S -> sqrt(S)
nep=SPMF_NEP([A1,A2,A3,A4],[f1,f2,f3,f4]);

f1 = S -> 1
@test_logs (:warn, r".*you have not provided valid matrix-functions.*") begin
    nep = SPMF_NEP([A1,A2,A3,A4], [f1,f2,f3,f4], check_consistency=true)
end

end


text:
# Special purpose formatter to unindent multi-line arglists
function format_arglist!(wheren, nchars_moved, lcurly_pos = 0, body_indent = 4)
    # We replace the call argument
    ocall = children(wheren)[1]
    (length(children(ocall)) == 2) && return
    firstarg = children(ocall)[3]
    # Heuristic: If the first argument is on a new line, don't try to change the indentation
    if ('\n' in leading_trivia(firstarg))
        return
    end
    # Heuristic: If the arugment indents (more spefically the first indented line), is further
    # left than the original lbrace, don't unindent, it's unlikely they were aligned to the
    # lparen.
    if lcurly_pos != 0
        indents = countindent_body(ocall)
        if !isempty(indents) && (indents[1] < lcurly_pos)
            return
        end
        # Heuristic: If this formatting would the argument list past the indent
        # of the body (or 4 spaces if there is no body provided), don't perform
        # the unindent
        if any(x->(x + nchars_moved) < body_indent, indents)
            return
        end
    end
    children(wheren)[1] = format_addindent_body(ocall, nchars_moved, nothing)
    nothing
end

function format_shortfunc_body!(tree, funcdef, nchars_moved)
    # Heurisitc: In short format, it is possible that people aligned after the equal
    # sign. If the whole body is after the equal sign, indent it.
    body = children(funcdef)[3]
    indents = Int[]
    countindent_body(body, indents)
    while length(indents) != 0 && indents[end] == 0
        pop!(indents)
    end
    eq_sign_offset = line_pos(funcdef, first(children(funcdef)[2].span))
    if length(indents) != 0 && all(indents .> eq_sign_offset)
        new_body = format_addindent_body(body, textwidth(" where "), parent(body))
        children(tree)[3] = new_body
    end
    nothing
end

CSTParser.defines_function(node::OverlayNode) = CSTParser.defines_function(node.expr)
CSTParser.has_sig(node::OverlayNode) = CSTParser.has_sig(node.expr)
function CSTParser.get_id(x::OverlayNode{BinarySyntaxOpCall})
    if isexpr(children(x)[2], OPERATOR ,Tokens.ISSUBTYPE) || isexpr(children(x)[2], OPERATOR, Tokens.DECLARATION)
        return CSTParser.get_id(children(x)[1])
    else
        return x
    end
end

CSTParser.get_id(x::OverlayNode{Curly}) = CSTParser.get_id(children(x)[1])

function compute_new_line_pos(tree, node, offset = 0, consider_trailing=true)
    if isempty(children(tree))
        if consider_trailing
            ws = trailing_trivia(tree)
            if '\n' in ws
                return length(last_line(ws))
            end
        end
        x = sprint(print_replacement, tree, true, consider_trailing)
        return offset + length(x)
    end
    for (i,c) in enumerate(children(tree))
        if c == node
            return offset
        end
        is_last = i == length(children(tree))
        offset = compute_new_line_pos(c, node, offset, !isa(c, TriviaReplacementNode) &&
            !(is_last && !consider_trailing))
        if isa(c, TriviaReplacementNode) && !(!consider_trailing && is_last)
            ws = trailing_trivia(c)
            if '\n' in ws
                offset = length(last_line(ws))
            end
        end
    end
    offset
end

begin
    struct OldParametricSyntax; end
    register(OldParametricSyntax, Deprecation(
        "Parameteric syntax of the form f{T}(x::T) is deprecated and needs to be written using the `where` keyword",
        "julia",
        v"0.6.0",
        v"0.7.0-DEV.1143",
        typemax(VersionNumber)
    ))
    applies_in_macrocall(dep::OldParametricSyntax, context) = true

    function get_struct_parent(expr)
        expr.parent == nothing && return nothing
        (isexpr(expr.parent, CSTParser.Struct) ||
            isexpr(expr.parent, CSTParser.Mutable)) && return expr.parent
        # Nested function are ok
        CSTParser.defines_function(expr.parent) && return nothing
        return get_struct_parent(expr.parent)
    end

    function struct_name(expr)
        i = 1
        while isexpr(children(expr)[i], KEYWORD)
            i += 1
        end
        children(expr)[i]
    end

    function extract_identifiers(exprs)
        identifiers = Symbol[]
        for expr in exprs
            isexpr(expr, PUNCTUATION) && continue
            id = CSTParser.get_id(expr)
            @assert isexpr(id, IDENTIFIER)
            push!(identifiers, Expr(id.expr))
        end
        identifiers
    end

    function deconflict_identifiers(exprs, conflicts)
        new_exprs = map(exprs) do expr
            isexpr(expr, PUNCTUATION) && return expr
            id_expr = CSTParser.get_id(expr)
            orig_id = id = Expr(id_expr.expr)
            @assert isa(id, Symbol)
            while id in conflicts
                id = Symbol(string(id,"_"))
            end
            id == orig_id && return expr
            return replace_node(expr, id_expr, ReplacementNode(String(id), leading_trivia(id_expr), trailing_trivia(id_expr)))
        end
        new_exprs
    end

    function rewrite_param_syntax(dep, expr, resolutions)
        CSTParser.defines_function(expr) || return
        sp = get_struct_parent(expr)
        # If there's already a where expr, this is new syntax
        call = isexpr(expr, FunctionDef) ? children(expr)[2] : children(expr)[1]
        had_rt_annotation = false
        call_or_rt_expr = call
        if isexpr(call, BinarySyntaxOpCall) && isexpr(children(call)[2],
                CSTParser.OPERATOR, Tokens.DECLARATION)
            call_or_rt_expr = rt_expr = call
            call = children(call)[1]
            had_rt_annotation = true
        end
        is_where_expr(call) && return
        length(children(call)) == 0 && return
        had_curly = isexpr(children(call)[1], Curly)
        tparams = []
        if had_curly
            curly = children(call)[1]
            fname = children(curly)[1]
            # Include punctuation
            tparams = children(curly)[2:end]
        else
            fname = children(call)[1]
        end
        needs_new_curly = false
        if sp !== nothing
            needs_new_curly = isexpr(struct_name(sp), Curly) && (isexpr(fname, IDENTIFIER) && Expr(fname.expr) == Expr(children(struct_name(sp))[1].expr))
            if needs_new_curly
                # Includes puctuation
                new_curlies = struct_name(sp)[2:end]
                if isempty(tparams)
                    tparams = new_curlies
                else
                    new_curlies = deconflict_identifiers(new_curlies, extract_identifiers(tparams))
                    tparams = [tparams[1]; new_curlies[2:end-1]; ReplacementNode(",",""," "); tparams[2:end]]
                end
            end
        end
        if !had_curly && !needs_new_curly
            return
        end
        # In long form syntax, with only one parameter, strip the curly braces
        if isexpr(expr, FunctionDef) && length(tparams) == 3
            tparams = tparams[2:2]
        end
        new_tree = ChildReplacementNode(nothing, children(expr)[(isexpr(expr, FunctionDef) ? 3 : 2):end], expr)
        replace_op = CSTParser.OPERATOR(0, 0:1, Tokens.ERROR, false)
        replace_lit = CSTParser.LITERAL(0, 0:1, "", Tokens.ERROR)
        new_where = TriviaReplacementNode(new_tree, ChildReplacementNode(new_tree,
            [ReplacementNode("where"," "," "), tparams...], CSTParser.BinarySyntaxOpCall(replace_lit, replace_op, replace_lit)),
            "", trailing_trivia(call_or_rt_expr))
        pushfirst!(children(new_tree), new_where)
        call_parent = new_where
        isexpr(expr, FunctionDef) && pushfirst!(children(new_tree), children(expr)[1])
        if had_rt_annotation
            annotation = children(rt_expr)[end]
            call_parent = ChildReplacementNode(new_tree, [children(rt_expr)[2]], rt_expr)
            push!(children(call_parent), TriviaReplacementNode(call_parent, annotation, leading_trivia(annotation), ""))
            pushfirst!(children(new_where), call_parent)
        end
        new_call = TriviaReplacementNode(new_where, ChildReplacementNode(new_where, children(call)[2:end], call), "", "")
        pushfirst!(children(call_parent), new_call)
        if needs_new_curly
            new_curly = TriviaReplacementNode(new_call, ChildReplacementNode(new_call, [fname, new_curlies...], had_curly ? Curly : EXPR{Curly}(Expr[], 0, 0:1)),"","")
            pushfirst!(children(new_call), new_curly)
        else
            pushfirst!(children(new_call), fname)
        end
        nchars_moved = (needs_new_curly ? sum(expr->sum(textwidth, fullspan_text(expr)), children(new_curly)[2:end]) : 0) -
                       (had_curly ? line_pos(call, first(children(call)[2].span)) - (line_pos(curly, first(children(curly)[2].span))) : 0)
        heuristic_pos = had_curly ? line_pos(curly, first(children(curly)[2].span)) : line_pos(call, first(children(call)[2].span)-1)
        format_arglist!(new_where, nchars_moved, heuristic_pos, isexpr(expr, FunctionDef) ? 4 : 0)
        # If the parameter list is multi-line, it might need to be indented as well
        if had_curly
            new_lbrace_pos = compute_new_line_pos(new_where, children(new_where)[3])
            old_lbrace_pos = line_pos(curly, first(children(curly)[2].span))
            where_indent = new_lbrace_pos - old_lbrace_pos
            for i = 4:length(children(new_where))
                children(new_where)[i] = format_addindent_body(children(new_where)[i], where_indent, nothing)
            end
        end
        if !isexpr(expr, FunctionDef)
            format_shortfunc_body!(new_tree, expr, nchars_moved)
        end
        buf = IOBuffer()
        print_replacement(buf, new_tree, false, false)
        push!(resolutions, TextReplacement(dep, expr.span, String(take!(buf))))
    end

    match(OldParametricSyntax, CSTParser.FunctionDef) do x
        dep, expr, resolutions, context = x
        rewrite_param_syntax(dep, expr, resolutions)
    end

    match(OldParametricSyntax, CSTParser.BinarySyntaxOpCall) do x
        dep, expr, resolutions, context = x
        rewrite_param_syntax(dep, expr, resolutions)
    end
end


text:
module WeaveSupport
## Modified from
## https://github.com/SciML/SciMLTutorials.jl/blob/master/src/SciMLTutorials.jl

import Base64: base64encode
using Random

using Weave

using Mustache
import Markdown
using JSON
using Reexport
@reexport using LaTeXStrings
#using SymPy
#function Base.show(io::IO, ::MIME"text/latex", x::SymPy.SymbolicObject)
#    print(io, SymPy.as_markdown(sympy.latex(x)))
#end

using Pkg

include("formatting.jl")
include("bootstrap.jl")
include("questions.jl")
include("show-methods.jl")
include("toc.jl")

#import Plots
# just show body, not standalone
#function Plots._show(io::IO, ::MIME"text/html", plt::Plots.Plot{Plots.PlotlyBackend})
#    write(io, Plots.html_body(plt))
#end

## we have jmd files that convert to html files
## using a specialized template

const repo_directory = joinpath(@__DIR__,"..", "..")
const cssfile = joinpath(@__DIR__, "..", "..", "templates", "skeleton_css.css")
const htmlfile = joinpath(@__DIR__,"..", "..", "templates", "bootstrap.tpl")
const latexfile = joinpath(@__DIR__, "..", "..", "templates", "julia_tex.tpl")

function weave_file(folder, file; build_list=(:script,:html,:pdf,:github,:notebook), force=false, kwargs...)

    jmddir = joinpath(repo_directory,"CwJ",folder)
    tmp = joinpath(jmddir, file)
    bnm = replace(basename(tmp), r".jmd$" => "")

    if !force
        testfile = joinpath(repo_directory, "html", folder, bnm*".html")
        if isfile(testfile) && (mtime(testfile) >= mtime(tmp))
            return
        end
    end

    Pkg.activate(dirname(tmp))
    Pkg.instantiate()
    args = Dict{Symbol,String}(:folder=>folder,:file=>file)

    if :script ∈ build_list
        println("Building Script")
        dir = joinpath(repo_directory,"script",folder)
        isdir(dir) || mkpath(dir)
        args[:doctype] = "script"
        tangle(tmp;out_path=dir)
    end

    if :html ∈ build_list
        println("Building HTML")
        dir = joinpath(repo_directory,"html",folder)
        isdir(dir) || mkpath(dir)

        figdir = joinpath(jmddir,"figures")
        htmlfigdir = joinpath(dir, "figures")

        if isdir(figdir)
            isdir(htmlfigdir) && rm(htmlfigdir, recursive=true)
            cp(figdir, htmlfigdir)
        end

        Weave.set_chunk_defaults!(:wrap=>false)
        args[:doctype] = "html"
        #weave(tmp,doctype = "md2html",out_path=dir,args=args; fig_ext=".svg", css=cssfile, kwargs...)
        weave(tmp,doctype = "md2html", out_path=dir,args=args; fig_ext=".svg",
              template=htmlfile,
              fig_path=tempdir(),
              kwargs...)

        # clean up
        isdir(htmlfigdir) && rm(htmlfigdir, recursive=true)

    end

    if :pdf ∈ build_list

        eval(quote using Tectonic end) # load Tectonic; wierd testing error

        println("Building PDF")
        dir = joinpath(repo_directory,"pdf",folder)
        isdir(dir) || mkpath(dir)

        fig_path = "_figures_" * bnm
        figdir = joinpath(jmddir,"figures")
        texfigdir = joinpath(dir, "figures")

        if isdir(figdir)
            isdir(texfigdir) && rm(texfigdir, recursive=true)
            cp(figdir, texfigdir)
        end

        args[:doctype] = "pdf"
        try
            weave(tmp,doctype="md2tex",out_path=dir,args=args;
                  template=latexfile,
                  fig_path=fig_path,
                  kwargs...)

            texfile = joinpath(dir, bnm * ".tex")
            Base.invokelatest(Tectonic.tectonic, bin -> run(`$bin $texfile`))

            # clean up
            for ext in (".tex",)
                f = joinpath(dir, bnm * ext)
                isfile(f) && rm(f)
            end

        catch ex
            @warn "PDF generation failed" exception=(ex, catch_backtrace())

        end

        isdir(texfigdir) && rm(texfigdir, recursive=true)
        isdir(joinpath(dir,fig_path)) && rm(joinpath(dir,fig_path), recursive=true)
        for ext in (".aux", ".log", ".out")
            f = joinpath(dir, bnm * ext)
            isfile(f) && rm(f)
        end

    end

    if :github ∈ build_list
        println("Building Github Markdown")
        dir = joinpath(repo_directory,"markdown",folder)
        isdir(dir) || mkpath(dir)
        args[:doctype] = "github"
        weave(tmp,doctype = "github",out_path=dir, args=args;
              fig_path=tempdir(),
              kwargs...)
    end

    if :notebook ∈ build_list
        println("Building Notebook")
        dir = joinpath(repo_directory,"notebook",folder)
        isdir(dir) || mkpath(dir)
        args[:doctype] = "notebook"
        Weave.convert_doc(tmp,joinpath(dir,file[1:end-4]*".ipynb"))
    end
end

"""
    weave_all(; force=false, build_list=(:script,:html,:pdf,:github,:notebook))

Run `weave` on all source files.

* `force`: by default, only run `weave` on files with `html` file older than the source file in `CwJ`
* `build_list`: list of output types to be built. The default is all types

The files will be built as subdirectories in the package directory. This is returned by `pathof(CalculusWithJulia)`.

"""
function weave_all(;force=false, build_list=(:script,:html,:pdf,:github,:notebook))
    for folder in readdir(joinpath(repo_directory,"CwJ"))
        folder == "test.jmd" && continue
        weave_folder(folder; force=force, build_list=build_list)
    end
end

function weave_folder(folder; force=false, build_list=(:script,:html,:pdf,:github,:notebook))
    !isnothing(match(r"\.ico$", folder)) && return nothing
    for file in readdir(joinpath(repo_directory,"CwJ",folder))
        !occursin(r".jmd$", basename(file)) && continue
        println("Building $(joinpath(folder,file))")
        try
            weave_file(folder,file; force=force, build_list=build_list)
        catch
        end
    end
end

macro q_str(x)
    "`$x`"
end

export mmd
export @q_str
export ImageFile, Verbatim, Invisible, Outputonly, HTMLonly, JSXGraph
export alert, warning, note
export example, popup, table
export gif_to_data
export numericq, radioq, booleanq, yesnoq, shortq, longq, multiq

end # module


text:
module ELKI
    using ..OutlierDetectionData: ARFF
    using OutlierDetection: CLASS_NORMAL, CLASS_OUTLIER, to_categorical
    using DataFrames: DataFrame, Not, select, convert, disallowmissing!, rename!
    using DataDeps

    export list, load

    include("datasets.jl")

    const _meta = [
        (
            name = "semantic",
            md5 = "e77ad3d37baeb1e4025497496e2183d3d00c6ea2d1e603b2d388894b490fa71a",
            url = "https://www.dbs.ifi.lmu.de/research/outlier-evaluation/input/semantic.tar.gz",
        ),
        (
            name = "literature",
            md5 =  "42a69b719f846bc03e4a6a0e5776969d4cd6b61a17675035c27753e7a6e0e256",
            url =  "https://www.dbs.ifi.lmu.de/research/outlier-evaluation/input/literature.tar.gz",
        )
    ]

    list() = vcat(ELKI_SEMANTIC, ELKI_LITERATURE)
    to_name(dataset::AbstractString) = "elki-$dataset"

    function load(dataset::AbstractString)
        dep = to_dep(dataset)
        folder = match(r"^(\w+?)([_|\.]|$)", dataset, 1).captures[1]
        data = ARFF.read("$dep/$folder/$dataset.arff"; quotechar='\'', escapechar='\'') |> DataFrame
        # lowercase all columns because some datasets use :Outlier instead of :outlier
        rename!(data, lowercase.(names(data)))
        X, y = select(data, Not([:outlier, :id])), ifelse.(data.outlier .== "no", CLASS_NORMAL, CLASS_OUTLIER)
        return X, to_categorical(y)
    end

    function to_dep(dataset_name::AbstractString)
        if dataset_name in ELKI_LITERATURE
            dep = @datadep_str to_name("literature")
            return joinpath(dep, "literature")
        elseif dataset_name in ELKI_SEMANTIC
            dep = @datadep_str to_name("semantic")
            return joinpath(dep, "semantic")
        else
            throw(ArgumentError("$dataset_name was not found in the available datasets"))
        end
    end

    function __init__()
        for (name, md5, url) in _meta
            register(DataDep(
                to_name(name),
                """Dataset: $name
                Collection: On the Evaluation of Unsupervised Outlier Detection
                Authors: Campos et al.
                Website: https://www.dbs.ifi.lmu.de/research/outlier-evaluation/DAMI/""",
                url,
                md5,
                post_fetch_method = unpack
            ))
        end
    end
end


text:
get_values(x::Variable) = x.values
get_values(x) = x

function get_gy(f::DiffableFunction)
    outputs = f.grad_field.outputs
    return [output.grad for output in outputs]
end

ones_like(x::R) where {R <: Real} = one(R)

ones_like(x::AbstractArray{R}) where {R <: Real} = ones(R, size(x))

"""
    cleargrad!(x::Variable)
Reset the Variable's gradient.

```julia-repl
julia> x
name: nothing
values: 1
grad: Variable(1)
creator: User-Defined(nothing)

julia> x.grad
name: nothing
values: 1
creator: User-Defined(nothing)

julia> cleargrad!(x)

julia> x.grad === nothing
true
```
"""
function cleargrad!(x::Variable)
    x.grad = nothing
end

function as_tuple(x)
    return (x,)
end

function as_tuple(x::Tuple)
    return x
end


text:
using Plots
using MeshCat
using GeometryTypes
using CoordinateTransformations
using FileIO
using MeshIO
using Random
IMAGE_DIR = joinpath(TrajectoryOptimization.root_dir(),"development","paper","images")

Random.seed!(123)

# Solver options
N = 101 # 201
integration = :rk4
opts = SolverOptions()
opts.verbose = false
opts.square_root = true
opts.cost_tolerance = 1e-5
opts.cost_tolerance_intermediate = 1e-5
opts.constraint_tolerance = 1e-4
opts.outer_loop_update_type = :feedback

dircol_options = Dict("tol"=>opts.cost_tolerance,"constr_viol_tol"=>opts.constraint_tolerance)

# Obstacle Avoidance
model,obj_uncon = TrajectoryOptimization.Dynamics.quadrotor
r_quad = 3.0
n = model.n
m = model.m
obj_con = TrajectoryOptimization.Dynamics.quadrotor_3obs[2]
spheres = TrajectoryOptimization.Dynamics.quadrotor_3obs[3]
n_spheres = length(spheres)

solver_uncon = Solver(model,obj_uncon,integration=integration,N=N,opts=opts)
solver_con = Solver(model,obj_con,integration=integration,N=N,opts=opts)

U_hover = 0.5*9.81/4.0*ones(solver_uncon.model.m, solver_uncon.N-1)
X_hover = rollout(solver_uncon,U_hover)
@time results_uncon, stats_uncon = solve(solver_uncon,U_hover)
@time results_uncon_dircol, stats_uncon_dircol = TrajectoryOptimization.solve_dircol(solver_uncon, X_hover, U_hover, options=dircol_options)

@time results_con, stats_con = solve(solver_con,U_hover)
@time results_dircol, stats_dircol = TrajectoryOptimization.solve_dircol(solver_con, X_hover, U_hover, options=dircol_options)

# Trajectory plots
t_array = range(0,stop=solver_con.obj.tf,length=solver_con.N)
plot(t_array[1:end-1],to_array(results_con.U)',title="Quadrotor Obstacle Avoidance",xlabel="time",ylabel="control",labels="")
savefig(joinpath(IMAGE_DIR,"quadrotor_obstacles_control_traj.eps"))
plot(t_array,to_array(results_con.X)[1:3,:]',title="Quadrotor Obstacle Avoidance",xlabel="time",ylabel="position",labels=["x";"y";"z"],legend=:topleft)
savefig(joinpath(IMAGE_DIR,"quadrotor_obstacles_pos_traj.eps"))

plot(t_array[1:end-1],Array(results_uncon_dircol.U)[:,1:end-1]',title="Quadrotor Obstacle Avoidance",xlabel="time",ylabel="control",labels="dircol")
plot(t_array,Array(results_uncon_dircol.X)[1:3,:]',title="Quadrotor Obstacle Avoidance",xlabel="time",ylabel="control",labels="dircol")

t_array = range(0,stop=solver_con.obj.tf,length=solver_con.N)
plot(t_array[1:end-1],Array(results_dircol.U)[:,1:end-1]',title="Quadrotor Obstacle Avoidance",xlabel="time",ylabel="control",labels="")
savefig(joinpath(IMAGE_DIR,"quadrotor_obstacles_control_traj_dircol.eps"))
plot(t_array,Array(results_dircol.X)[1:3,:]',title="Quadrotor Obstacle Avoidance",xlabel="time",ylabel="position",labels=["x";"y";"z"],legend=:topleft)
savefig(joinpath(IMAGE_DIR,"quadrotor_obstacles_pos_traj_dircol.eps"))

@assert max_violation(results_con) <= opts.constraint_tolerance

# Constraint convergence plot
plot(stats_con["c_max"],yscale=:log10,title="Quadrotor Obstacle Avoidance",xlabel="iteration",ylabel="max constraint violation",label="ALTRO")
savefig(joinpath(IMAGE_DIR,"quadrotor_obstacles_constraint_convergence.eps"))
plot(stats_dircol["c_max"],yscale=:log10,title="Quadrotor Obstacle Avoidance",xlabel="iteration",ylabel="log(max constraint violation)",label="DIRCOL")
savefig(joinpath(IMAGE_DIR,"quadrotor_obstacles_constraint_convergence_dircol.eps"))

@show stats_con["iterations"]
@show stats_con["outer loop iterations"]
@show stats_con["c_max"][end]
@show stats_con["cost"][end]

#################
# Visualization #
#################

vis = Visualizer()
open(vis)

# Import quadrotor obj file
traj_folder = joinpath(dirname(pathof(TrajectoryOptimization)),"..")
urdf_folder = joinpath(traj_folder, "dynamics","urdf")
obj = joinpath(urdf_folder, "quadrotor_base.obj")

# color options
green_ = MeshPhongMaterial(color=RGBA(0, 1, 0, 1.0))
green_transparent = MeshPhongMaterial(color=RGBA(0, 1, 0, 0.1))
red_ = MeshPhongMaterial(color=RGBA(1, 0, 0, 1.0))
blue_ = MeshPhongMaterial(color=RGBA(0, 0, 1, 1.0))
blue_transparent = MeshPhongMaterial(color=RGBA(0, 0, 1, 0.1))

orange_ = MeshPhongMaterial(color=RGBA(233/255, 164/255, 16/255, 1.0))
orange_transparent = MeshPhongMaterial(color=RGBA(233/255, 164/255, 16/255, 0.1))
black_ = MeshPhongMaterial(color=RGBA(0, 0, 0, 1.0))
black_transparent = MeshPhongMaterial(color=RGBA(0, 0, 0, 0.1))

# geometries
robot_obj = FileIO.load(obj)
sphere_small = HyperSphere(Point3f0(0), convert(Float32,0.1*r_quad)) # trajectory points
sphere_medium = HyperSphere(Point3f0(0), convert(Float32,r_quad))

obstacles = vis["obs"]
traj = vis["traj"]
traj_uncon = vis["traj_uncon"]
target = vis["target"]
robot = vis["robot"]
robot_uncon = vis["robot_uncon"]

# Set camera location
settransform!(vis["/Cameras/default"], compose(Translation(25., 15., 20.),LinearMap(RotY(-pi/12))))

# Create and place obstacles
for i = 1:n_spheres
    setobject!(vis["obs"]["s$i"],HyperSphere(Point3f0(0), convert(Float32,spheres[i][4])),red_)
    settransform!(vis["obs"]["s$i"], Translation(spheres[i][1], spheres[i][2], spheres[i][3]))
end

# Create and place trajectory
for i = 1:N
    setobject!(vis["traj_uncon"]["t$i"],sphere_small,orange_)
    settransform!(vis["traj_uncon"]["t$i"], Translation(results_uncon.X[i][1], results_uncon.X[i][2], results_uncon.X[i][3]))
end
for i = 1:N
    setobject!(vis["traj"]["t$i"],sphere_small,blue_)
    settransform!(vis["traj"]["t$i"], Translation(results_con.X[i][1], results_con.X[i][2], results_con.X[i][3]))
end

# Create and place initial position
# setobject!(vis["robot_uncon"]["ball"],sphere_medium,orange_transparent)
# setobject!(vis["robot_uncon"]["quad"],robot_obj,black_)

# setobject!(vis["robot"]["ball"],sphere_medium,blue_transparent)
setobject!(vis["robot"]["quad"],robot_obj,black_)

# Animate quadrotor
for i = 1:N
    settransform!(vis["robot_uncon"], compose(Translation(results_uncon.X[i][1:3]...),LinearMap(Quat(results_uncon.X[i][4:7]...))))
    sleep(solver_uncon.dt)
end

for i = 1:N
    settransform!(vis["robot"], compose(Translation(results_con.X[i][1:3]...),LinearMap(Quat(results_con.X[i][4:7]...))))
    sleep(solver_con.dt)
end

# Ghose quadrotor scene
traj_idx = [1;15;25;35;50;N]
n_robots = length(traj_idx)
for i = 1:n_robots
    robot = vis["robot_$i"]
    setobject!(vis["robot_$i"]["quad"],robot_obj,black_)
    settransform!(vis["robot_$i"], compose(Translation(results_con.X[traj_idx[i]][1:3]...),LinearMap(Quat(results_con.X[traj_idx[i]][4:7]...))))
end


text:
"""
    function Simpson(f, a, b)

Calculate the integral `f` function from `x` to `y` by Simpson formula

Here's Simpson expression:

``Simpson(f, a, b) = \\frac{(b-a)}{6} * (f(a) + 4*f(\\frac{(a+b)}{2}) + f(b))`

# Arguments
- `f::Express`: integrand function
- `a::Float`: left boundary of integral interval of definite integral
- `b::Float`: right boundary of integral interval of definite integral

# Examples
```julia-repl
julia> Simpson(2*x, 1, 2)
4
```
"""
function Simpson(f, a, b)
    # Simpson Formula
    # (b-a)/6 * (f(a) + 4*f((a+b)/2) + f(b))
    h = (b - a)/2
    x = a:h:b
    sum((x[3] - x[1])/6*(f(x[3]) + 4*f(x[2]) + f(x[1])))
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

## type join (closest common ancestor, or least upper bound) ##

typejoin() = (@_pure_meta; Bottom)
typejoin(t::ANY) = (@_pure_meta; t)
typejoin(t::ANY, ts...) = (@_pure_meta; typejoin(t, typejoin(ts...)))
function typejoin(a::ANY, b::ANY)
    @_pure_meta
    if isa(a,TypeConstructor); a = a.body; end
    if isa(b,TypeConstructor); b = b.body; end
    if a <: b
        return b
    elseif b <: a
        return a
    end
    if isa(a,TypeVar)
        return typejoin(a.ub, b)
    end
    if isa(b,TypeVar)
        return typejoin(a, b.ub)
    end
    if isa(a,Union) || isa(b,Union)
        u = Union{a, b}
        if !isa(u,Union)
            return u
        end
        return reduce(typejoin, Bottom, u.types)
    end
    if a <: Tuple
        if !(b <: Tuple)
            return Any
        end
        ap, bp = a.parameters, b.parameters
        lar = length(ap)::Int; lbr = length(bp)::Int
        laf, afixed = full_va_len(ap)
        lbf, bfixed = full_va_len(bp)
        if lar==0 || lbr==0
            return Tuple
        end
        if laf < lbf
            if isvarargtype(ap[lar]) && !afixed
                c = Vector{Any}(laf)
                c[laf] = Vararg{typejoin(ap[lar].parameters[1], tailjoin(bp,laf))}
                n = laf-1
            else
                c = Vector{Any}(laf+1)
                c[laf+1] = Vararg{tailjoin(bp,laf+1)}
                n = laf
            end
        elseif lbf < laf
            if isvarargtype(bp[lbr]) && !bfixed
                c = Vector{Any}(lbf)
                c[lbf] = Vararg{typejoin(bp[lbr].parameters[1], tailjoin(ap,lbf))}
                n = lbf-1
            else
                c = Vector{Any}(lbf+1)
                c[lbf+1] = Vararg{tailjoin(ap,lbf+1)}
                n = lbf
            end
        else
            c = Vector{Any}(laf)
            n = laf
        end
        for i = 1:n
            ai = ap[min(i,lar)]; bi = bp[min(i,lbr)]
            ci = typejoin(unwrapva(ai),unwrapva(bi))
            c[i] = i == length(c) && (isvarargtype(ai) || isvarargtype(bi)) ? Vararg{ci} : ci
        end
        return Tuple{c...}
    elseif b <: Tuple
        return Any
    end
    while b !== Any
        if a <: b.name.primary
            while a.name !== b.name
                a = supertype(a)
            end
            # join on parameters
            n = length(a.parameters)
            p = Vector{Any}(n)
            for i = 1:n
                ai, bi = a.parameters[i], b.parameters[i]
                if ai === bi || (isa(ai,Type) && isa(bi,Type) && typeseq(ai,bi))
                    p[i] = ai
                else
                    p[i] = a.name.primary.parameters[i]
                end
            end
            return a.name.primary{p...}
        end
        b = supertype(b)
    end
    return Any
end

# Returns length, isfixed
function full_va_len(p)
    isempty(p) && return 0, true
    if isvarargtype(p[end])
        N = p[end].parameters[2]
        if isa(N, Integer)
            return (length(p) + N - 1)::Int, true
        end
        return length(p)::Int, false
    end
    return length(p)::Int, true
end

# reduce typejoin over A[i:end]
function tailjoin(A, i)
    if i > length(A)
        return unwrapva(A[end])
    end
    t = Bottom
    for j = i:length(A)
        t = typejoin(t, unwrapva(A[j]))
    end
    return t
end

## promotion mechanism ##

promote_type()  = (@_pure_meta; Bottom)
promote_type(T) = (@_pure_meta; T)
promote_type(T, S, U, V...) = (@_pure_meta; promote_type(T, promote_type(S, U, V...)))

promote_type(::Type{Bottom}, ::Type{Bottom}) = (@_pure_meta; Bottom)
promote_type{T}(::Type{T}, ::Type{T}) = (@_pure_meta; T)
promote_type{T}(::Type{T}, ::Type{Bottom}) = (@_pure_meta; T)
promote_type{T}(::Type{Bottom}, ::Type{T}) = (@_pure_meta; T)

# Try promote_rule in both orders. Typically only one is defined,
# and there is a fallback returning Bottom below, so the common case is
#   promote_type(T, S) =>
#   promote_result(T, S, result, Bottom) =>
#   typejoin(result, Bottom) => result
function promote_type{T,S}(::Type{T}, ::Type{S})
    @_pure_meta
    promote_result(T, S, promote_rule(T,S), promote_rule(S,T))
end

promote_rule(T, S) = (@_pure_meta; Bottom)

promote_result(t,s,T,S) = (@_pure_meta; promote_type(T,S))
# If no promote_rule is defined, both directions give Bottom. In that
# case use typejoin on the original types instead.
promote_result{T,S}(::Type{T},::Type{S},::Type{Bottom},::Type{Bottom}) = (@_pure_meta; typejoin(T, S))

promote() = ()
promote(x) = (x,)
function promote{T,S}(x::T, y::S)
    (convert(promote_type(T,S),x), convert(promote_type(T,S),y))
end
promote_typeof(x) = (@_pure_meta; typeof(x))
promote_typeof(x, xs...) = (@_pure_meta; promote_type(typeof(x), promote_typeof(xs...)))
function promote(x, y, z)
    (convert(promote_typeof(x,y,z), x),
     convert(promote_typeof(x,y,z), y),
     convert(promote_typeof(x,y,z), z))
end
function promote(x, y, zs...)
    (convert(promote_typeof(x,y,zs...), x),
     convert(promote_typeof(x,y,zs...), y),
     convert(Tuple{Vararg{promote_typeof(x,y,zs...)}}, zs)...)
end
# TODO: promote{T}(x::T, ys::T...) here to catch all circularities?

## promotions in arithmetic, etc. ##

# Because of the promoting fallback definitions for Number, we need
# a special case for undefined promote_rule on numeric types.
# Otherwise, typejoin(T,S) is called (returning Number) so no conversion
# happens, and +(promote(x,y)...) is called again, causing a stack
# overflow.
function promote_result{T<:Number,S<:Number}(::Type{T},::Type{S},::Type{Bottom},::Type{Bottom})
    @_pure_meta
    promote_to_supertype(T, S, typejoin(T,S))
end

# promote numeric types T and S to typejoin(T,S) if T<:S or S<:T
# for example this makes promote_type(Integer,Real) == Real without
# promoting arbitrary pairs of numeric types to Number.
promote_to_supertype{T<:Number          }(::Type{T}, ::Type{T}, ::Type{T}) = (@_pure_meta; T)
promote_to_supertype{T<:Number,S<:Number}(::Type{T}, ::Type{S}, ::Type{T}) = (@_pure_meta; T)
promote_to_supertype{T<:Number,S<:Number}(::Type{T}, ::Type{S}, ::Type{S}) = (@_pure_meta; S)
promote_to_supertype{T<:Number,S<:Number}(::Type{T}, ::Type{S}, ::Type) =
    error("no promotion exists for ", T, " and ", S)

+(x::Number, y::Number) = +(promote(x,y)...)
*(x::Number, y::Number) = *(promote(x,y)...)
-(x::Number, y::Number) = -(promote(x,y)...)
/(x::Number, y::Number) = /(promote(x,y)...)
^(x::Number, y::Number) = ^(promote(x,y)...)

fma(x::Number, y::Number, z::Number) = fma(promote(x,y,z)...)
muladd(x::Number, y::Number, z::Number) = muladd(promote(x,y,z)...)

(&)(x::Integer, y::Integer) = (&)(promote(x,y)...)
(|)(x::Integer, y::Integer) = (|)(promote(x,y)...)
xor(x::Integer, y::Integer) = xor(promote(x,y)...)

==(x::Number, y::Number) = (==)(promote(x,y)...)
<( x::Real, y::Real)     = (< )(promote(x,y)...)
<=(x::Real, y::Real)     = (<=)(promote(x,y)...)

div(x::Real, y::Real) = div(promote(x,y)...)
fld(x::Real, y::Real) = fld(promote(x,y)...)
cld(x::Real, y::Real) = cld(promote(x,y)...)
rem(x::Real, y::Real) = rem(promote(x,y)...)
mod(x::Real, y::Real) = mod(promote(x,y)...)

mod1(x::Real, y::Real) = mod1(promote(x,y)...)
fld1(x::Real, y::Real) = fld1(promote(x,y)...)

max(x::Real, y::Real) = max(promote(x,y)...)
min(x::Real, y::Real) = min(promote(x,y)...)
minmax(x::Real, y::Real) = minmax(promote(x, y)...)

# "Promotion" that takes a function into account and tries to preserve
# non-concrete types. These are meant to be used mainly by elementwise
# operations, so it is advised against overriding them
_default_type(T::Type) = (@_pure_meta; T)

promote_op(::Any...) = (@_pure_meta; Any)
function promote_op{S}(f, ::Type{S})
    @_inline_meta
    Z = Tuple{_default_type(S)}
    T = _default_eltype(Generator{Z, typeof(f)})
    isleaftype(S) && return isleaftype(T) ? T : Any
    return typejoin(S, T)
end
function promote_op{R,S}(f, ::Type{R}, ::Type{S})
    @_inline_meta
    Z = Iterators.Zip2{Tuple{_default_type(R)}, Tuple{_default_type(S)}}
    T = _default_eltype(Generator{Z, typeof(a -> f(a...))})
    isleaftype(R) && isleaftype(S) && return isleaftype(T) ? T : Any
    return typejoin(R, S, T)
end

## catch-alls to prevent infinite recursion when definitions are missing ##

no_op_err(name, T) = error(name," not defined for ",T)
+{T<:Number}(x::T, y::T) = no_op_err("+", T)
*{T<:Number}(x::T, y::T) = no_op_err("*", T)
-{T<:Number}(x::T, y::T) = no_op_err("-", T)
/{T<:Number}(x::T, y::T) = no_op_err("/", T)
^{T<:Number}(x::T, y::T) = no_op_err("^", T)

fma{T<:Number}(x::T, y::T, z::T) = no_op_err("fma", T)
fma(x::Integer, y::Integer, z::Integer) = x*y+z
muladd{T<:Number}(x::T, y::T, z::T) = x*y+z

(&){T<:Integer}(x::T, y::T) = no_op_err("&", T)
(|){T<:Integer}(x::T, y::T) = no_op_err("|", T)
xor{T<:Integer}(x::T, y::T) = no_op_err("xor", T)

=={T<:Number}(x::T, y::T) = x === y
 <{T<:Real}(x::T, y::T) = no_op_err("<" , T)
<={T<:Real}(x::T, y::T) = no_op_err("<=", T)

rem{T<:Real}(x::T, y::T) = no_op_err("rem", T)
mod{T<:Real}(x::T, y::T) = no_op_err("mod", T)

min(x::Real) = x
max(x::Real) = x
minmax(x::Real) = (x, x)

max{T<:Real}(x::T, y::T) = ifelse(y < x, x, y)
min{T<:Real}(x::T, y::T) = ifelse(y < x, y, x)
minmax{T<:Real}(x::T, y::T) = y < x ? (y, x) : (x, y)

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia