text:
"""
This function creates the model for a full themal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(canonical::Canonical, sys::PSY.System,
                           model::DeviceModel{T, D},
                           ::Type{S};
                           kwargs...) where {T<:PSY.ThermalGen,
                                             D<:AbstractThermalFormulation,
                                             S<:PM.AbstractPowerModel}

    devices = PSY.get_components(T, sys)

    if validate_available_devices(devices, T)
        return
    end

    #Variables
    activepower_variables!(canonical, devices)

    reactivepower_variables!(canonical, devices)

    commitment_variables!(canonical, devices)

    #Initial Conditions

    initial_conditions!(canonical, devices, D)

    #Constraints
    activepower_constraints!(canonical, devices, D, S)

    reactivepower_constraints!(canonical, devices, D, S)

    commitment_constraints!(canonical, devices, D, S)

    ramp_constraints!(canonical, devices, D, S)

    time_constraints!(canonical, devices, D, S)

    feedforward!(canonical, T, model.feedforward)

    #Cost Function
    cost_function(canonical, devices, D, S)

    return

end

"""
This function creates the model for a full themal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(canonical::Canonical, sys::PSY.System,
                           model::DeviceModel{T, D},
                           ::Type{S};
                           kwargs...) where {T<:PSY.ThermalGen,
                                             D<:AbstractThermalFormulation,
                                             S<:PM.AbstractActivePowerModel}

    devices = PSY.get_components(T, sys)

    if validate_available_devices(devices, T)
        return
    end

    #Variables
    activepower_variables!(canonical, devices)

    commitment_variables!(canonical, devices)

    #Initial Conditions

    initial_conditions!(canonical, devices, D)

    #Constraints
    activepower_constraints!(canonical, devices, D, S)

    commitment_constraints!(canonical, devices, D, S)

    ramp_constraints!(canonical, devices, D, S)

    time_constraints!(canonical, devices, D, S)

    feedforward!(canonical, T, model.feedforward)

    #Cost Function
    cost_function(canonical, devices, D, S)

    return

end

"""
This function creates the model for a full themal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(canonical::Canonical, sys::PSY.System,
                           model::DeviceModel{T, ThermalBasicUnitCommitment},
                           ::Type{S};
                           kwargs...) where {T<:PSY.ThermalGen,
                                             S<:PM.AbstractPowerModel}

    devices = PSY.get_components(T, sys)

    if validate_available_devices(devices, T)
        return
    end

    #Variables
    activepower_variables!(canonical, devices)

    reactivepower_variables!(canonical, devices)

    commitment_variables!(canonical, devices)

    #Initial Conditions

    initial_conditions!(canonical, devices, model.formulation)

    #Constraints
    activepower_constraints!(canonical, devices, model.formulation, S)

    reactivepower_constraints!(canonical, devices, model.formulation, S)

    commitment_constraints!(canonical, devices, model.formulation, S)

    feedforward!(canonical, T, model.feedforward)

    #Cost Function
    cost_function(canonical, devices, model.formulation, S)

    return

end

"""
This function creates the model for a full themal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(canonical::Canonical, sys::PSY.System,
                           model::DeviceModel{T, ThermalBasicUnitCommitment},
                           ::Type{S};
                           kwargs...) where {T<:PSY.ThermalGen,
                                             S<:PM.AbstractActivePowerModel}

    devices = PSY.get_components(T, sys)

    if validate_available_devices(devices, T)
        return
    end

    #Variables
    activepower_variables!(canonical, devices)

    commitment_variables!(canonical, devices)

    #Initial Conditions

    initial_conditions!(canonical, devices, model.formulation)

    #Constraints
    activepower_constraints!(canonical, devices, model.formulation, S)

    commitment_constraints!(canonical, devices, model.formulation, S)

    feedforward!(canonical, T, model.feedforward)

    #Cost Function
    cost_function(canonical, devices, model.formulation, S)

    return

end

"""
This function creates the model for a full themal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(canonical::Canonical, sys::PSY.System,
                           model::DeviceModel{T, ThermalRampLimited},
                           ::Type{S};
                           kwargs...) where {T<:PSY.ThermalGen,
                                             S<:PM.AbstractPowerModel}

    devices = PSY.get_components(T, sys)

    if validate_available_devices(devices, T)
        return
    end

    #Variables
    activepower_variables!(canonical, devices)

    reactivepower_variables!(canonical, devices)

    #Initial Conditions

    initial_conditions!(canonical, devices, model.formulation)

    #Constraints
    if !(isa(model.feedforward, SemiContinuousFF))
        activepower_constraints!(canonical, devices, ThermalRampLimited, S)
    end

    reactivepower_constraints!(canonical, devices, model.formulation, S)

    ramp_constraints!(canonical, devices, model.formulation, S)

    feedforward!(canonical, T, model.feedforward)

    #Cost Function
    cost_function(canonical, devices, model.formulation, S)

    return

end

"""
This function creates the model for a full themal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(canonical::Canonical, sys::PSY.System,
                           model::DeviceModel{T, ThermalRampLimited},
                           ::Type{S};
                           kwargs...) where {T<:PSY.ThermalGen,
                                             S<:PM.AbstractActivePowerModel}

    devices = PSY.get_components(T, sys)

    if validate_available_devices(devices, T)
        return
    end

    #Variables
    activepower_variables!(canonical, devices)

    #Initial Conditions

    initial_conditions!(canonical, devices, model.formulation)

    #Constraints
    if !(isa(model.feedforward, SemiContinuousFF))
        activepower_constraints!(canonical, devices, ThermalRampLimited, S)
    end

    ramp_constraints!(canonical, devices, model.formulation, S)

    feedforward!(canonical, T, model.feedforward)

    #Cost Function
    cost_function(canonical, devices, model.formulation, S)

    return

end

function construct_device!(canonical::Canonical, sys::PSY.System,
                           model::DeviceModel{T, D},
                           ::Type{S};
                           kwargs...) where {T<:PSY.ThermalGen,
                                             D<:AbstractThermalDispatchFormulation,
                                             S<:PM.AbstractPowerModel}

    devices = PSY.get_components(T, sys)

    if validate_available_devices(devices, T)
        return
    end

    #Variables
    activepower_variables!(canonical, devices)

    reactivepower_variables!(canonical, devices)

    #Initial Conditions

    #Constraints
    if !(isa(model.feedforward, SemiContinuousFF))
        activepower_constraints!(canonical, devices, D, S)
    end

    reactivepower_constraints!(canonical, devices, D, S)

    feedforward!(canonical, T, model.feedforward)

    #Cost Function
    cost_function(canonical, devices, D, S)

    return

end

function construct_device!(canonical::Canonical, sys::PSY.System,
                           model::DeviceModel{T, D},
                           ::Type{S};
                           kwargs...) where {T<:PSY.ThermalGen,
                                             D<:AbstractThermalDispatchFormulation,
                                             S<:PM.AbstractActivePowerModel}

    devices = PSY.get_components(T, sys)

    if validate_available_devices(devices, T)
        return
    end

    #Variables
    activepower_variables!(canonical, devices)

    #Initial Conditions

    #Constraints
    # Slighly hacky for now
    if !(isa(model.feedforward, SemiContinuousFF))
        activepower_constraints!(canonical, devices, D, S)
    end

    feedforward!(canonical, T, model.feedforward)

    #Cost Function
    cost_function(canonical, devices, D, S)

    return

end


text:
# <hide>
function printsum(a)
    println(summary(a), ": ", repr(a))
end
# </hide>

# dicts can be initialised directly:
a1 = Dict(1=>"one", 2=>"two")
printsum(a1)
#> Dict{Int64,String} with 2 entries: Dict(2=>"two",1=>"one")

# then added to:
a1[3]="three"
printsum(a1)
#> Dict{Int64,String} with 3 entries: Dict(2=>"two",3=>"three",1=>"one")
# (note dicts cannot be assumed to keep their original order)

# dicts may also be created with the type explicitly set
a2 = Dict{Int64, AbstractString}()
a2[0]="zero"
printsum(a2)
#> Dict{Int64,AbstractString} with 1 entry: Dict{Int64,AbstractString}(0=>"zero")

# dicts, like arrays, may also be created from [comprehensions](https://docs.julialang.org/en/v1/manual/arrays/#Comprehensions-1)
using Printf
a3 = Dict([i => @sprintf("%d", i) for i = 1:10])
printsum(a3)
#> Dict{Int64,String} with 10 entries: Dict(7=>"7",4=>"4",9=>"9",10=>"10",
#>  2=>"2",3=>"3",5=>"5",8=>"8",6=>"6",1=>"1")

# as you would expect, Julia comes with all the normal helper functions
# for dicts, e.g., [haskey](https://docs.julialang.org/en/v1/base/collections/#Base.haskey)
println(haskey(a1,1)) #> true

# which is equivalent to
println(1 in keys(a1)) #> true
# where [keys](https://docs.julialang.org/en/v1/base/collections/#Base.keys) creates an iterator over the keys of the dictionary

# similar to keys, [values](https://docs.julialang.org/en/v1/base/collections/#Base.values) get iterators over the dict's values:
printsum(values(a1))
#> Base.ValueIterator for a Dict{Int64,String} with 3 entries: ["two", "three", "one"]

# use [collect](https://docs.julialang.org/en/v1/base/collections/#Base.collect-Tuple{Any}) to get an array:
printsum(collect(values(a1)))
#> 3-element Array{String,1}: ["two", "three", "one"]

ext: jl

max_stars_repo_name: talgatomarov/JuliaByExample

max_issues_repo_name: talgatomarov/JuliaByExample

max_forks_repo_name: talgatomarov/JuliaByExample


text:
struct Fuselage{T <: Real} <: AbstractAircraft
    length   :: T
    weights  :: Vector{T}
    radii    :: Vector{T}
    position :: SVector{3,T}
end

function Fuselage(L :: T, weights :: AbstractVector{T}, radii :: AbstractVector{T}, position = zeros(3)) where T <: Real
    @assert all(x -> 0 <= x <= 1., weights) "The distribution of weights must be lie in [0,1] ⊂ ℝ."
    perms = sortperm(weights); sort!(weights)
    Fuselage{T}(L, weights, radii[perms], position)
end

projected_area(fuse :: Fuselage) = forward_sum(fuse.radii) / 2 .* fuse.weights
Base.length(fuse :: Fuselage) = fuse.length

function coordinates(fuse :: Fuselage, n)
    ws_rads = cosine_spacing(fuse, n)

    [ ws_rads[:,1] .* fuse.length ws_rads[:,2] ]
end

function cosine_spacing(fuse :: Fuselage, n)
    xs = fuse.weights
    ys = fuse.radii
    cosine_interp([ xs ys ], n)
end


text:
module FirePOMDP

using POMDPs
using POMDPModelTools
using POMDPModels
using POMDPSimulators
using POMDPPolicies
using MCTS
using D3Trees
using Random
using StaticArrays
using Parameters
using DiscreteValueIteration
using LinearAlgebra
using Combinatorics
using Distances
using ParticleFilters
using BeliefUpdaters
using StatsBase
using BasicPOMCP
using Distributions
import Distributions: Uniform
using POMCPOW
using JLD

export
    FireWorld,
    FireState,
    FireObs

# make initial state
const RAND_INIT = 264
rng = MersenneTwister(RAND_INIT);

# initial fuel level
const DEFAULT_FUEL = 5

# initial burn map
# criteria: only a % burning - set maximum number
const BURN_PERC = 0.1
const BURN_THRESHOLD = 0.6
function make_burn_size(total_size::Int64, burn_perc::Float64)
    burn_size = floor(Int, total_size * burn_perc)
    if burn_size in [0, 1]
        burn_size = 2
    end
    not_burn_size = total_size - burn_size
    return [burn_size, not_burn_size]
end

function burn_map(burn_threshold::Float64, burns_size::Array{Int64,1}, rng::AbstractRNG)
    burn_prob_array = rand(rng, Uniform(burn_threshold,1), burns_size[1])
    # println("burn_prob_array ", burn_prob_array)
    not_burn_prob_array = rand(rng, Uniform(0, burn_threshold), burns_size[2])
    burn_map = vcat(burn_prob_array, not_burn_prob_array)
    return shuffle!(rng, burn_map)
end

# initial state
function make_initial_state(GRID_SIZE::Int64, rng::AbstractRNG)
    burns_size = make_burn_size(GRID_SIZE * GRID_SIZE, BURN_PERC)
    burn_prob_map = burn_map(BURN_THRESHOLD, burns_size, rng)
    init_burn = burn_prob_map .> BURN_THRESHOLD
    init_fuels = ones(Int, GRID_SIZE * GRID_SIZE) * DEFAULT_FUEL
    return FireState(init_burn, burn_prob_map, init_fuels)
end

# @show init_state = FireState(init_burn, burn_prob_map, init_fuels)

# make initial cost map
# high_cost_perc = 0.2 # according to the 2010 Census,
#                      # 95 percent of Californians live on just 5.3 percent of the land in the state
#                      # Less than 17 percent of the San Francisco Bay Area is developed
# mid_cost_perc = 0.3 # percentage of California forest coverage
# low_cost_perc = 1 - high_cost_perc - mid_cost_perc # 0.5

const COSTS_PERC = [0.2, 0.3, 0.5]
const COSTS_ARRAY = [-10, -5, -1]
# costs = rand(costs_array, total_size);

function make_cost_size(total_size::Int64, costs_perc::Array{Float64,1})
    high_cost_size = floor(Int, total_size * costs_perc[1])
    mid_cost_size = floor(Int, total_size * costs_perc[2])
    low_cost_size = total_size - high_cost_size - mid_cost_size
    return [high_cost_size, mid_cost_size, low_cost_size]
end

function cost_map(costs_array::Array{Int64,1}, costs_size::Array{Int64,1}, rng::AbstractRNG)
    high_cost_array = repeat([costs_array[1]], costs_size[1])
    mid_cost_array = repeat([costs_array[2]], costs_size[2])
    low_cost_array = repeat([costs_array[3]], costs_size[3])
    costs = vcat(high_cost_array, mid_cost_array, low_cost_array)
    return shuffle!(rng, costs)
end

function make_cost_map(GRID_SIZE::Int64, COSTS_PERC::Array{Float64,1}, COSTS_ARRAY::Array{Int64,1}, rng::AbstractRNG)
    costs_size = make_cost_size(GRID_SIZE * GRID_SIZE, COSTS_PERC)
    COSTS = cost_map(COSTS_ARRAY, costs_size, rng)
    return COSTS
end

# initial weather condition
const WIND = [1, 1, 5]

# POMDP State
struct FireState
    burning::BitArray{1} # an array maintaining what cells are burning
    burn_probs::Array{Float64,1} # an array of probability of cells burning
    fuels::Array{Int64,1} # an array maintaining each cell's fuel level
end

# POMDP Observation: only observes burning or not
struct FireObs
    burning::BitArray{1} # an array what cells are seen to be burning
end

# POMDP{State, Action, Observation}
@with_kw struct FireWorld <: POMDP{FireState, Array{Int64, 1}, FireObs}
    # size of grid
    grid_size::Int64 = 4
    # initial state
    state::FireState = make_initial_state(grid_size, rng)
    # fire total size, for extension of non-square shape
    map_size::Tuple{Int,Int} = (grid_size, grid_size)
    # cost map
    costs::Array{Int64,1} = make_cost_map(grid_size, COSTS_PERC, COSTS_ARRAY, rng)
    # 1-sensitivity; probability of an observation of no burning when it actually is, e.g. burning but sensor does not detect
    bprob_fn::Float64 = 0.2
    # 1-specificity; probability of an observation of burning when it isn't, false positive; due to lag cells are likely to be burning (more likely if p_small is lower)
    bprob_fp::Float64 = 0.1
    # probability of successfully putting out a fire
    tprob::Float64 = 1.0
    # discount factor
    discount::Float64 = 0.95
    # consists of strength, acceleration, and direction
    wind::Array{Int64,1} = WIND
end

# Discount factor
POMDPs.discount(pomdp::FireWorld) = pomdp.discount;

# Terminal condition
POMDPs.isterminal(pomdp::FireWorld, state::FireState) = sum(state.burning) == 0

# Equal condition
POMDPs.isequal(s1::FireState, s2::FireState) = s1.burning == s2.burning && s1.fuels == s2.fuels && s1.burn_probs == s2.burn_probs

include("updater.jl")
include("states.jl")
include("actions.jl")
include("transition.jl")
include("observations.jl")
include("reward.jl")
# include("visualization.jl")

end # module


text:
"""$(TYPEDEF) Solver options for the ALTRO solver.
$(FIELDS)
"""
@with_kw mutable struct ALTROSolverOptions{T} <: AbstractSolverOptions{T}

    verbose::Bool=false

    "Augmented Lagrangian solver options."
    opts_al::AugmentedLagrangianSolverOptions=AugmentedLagrangianSolverOptions{T}()

    # Infeasible Start
    "infeasible control constraint tolerance."
    constraint_tolerance_infeasible::T = 1.0e-5

    "regularization term for infeasible controls."
    R_inf::T = 1.0

    "project infeasible results to feasible space using TVLQR."
    dynamically_feasible_projection::Bool = true

    "resolve feasible problem after infeasible solve."
    resolve_feasible_problem::Bool = true

    "initial penalty term for infeasible controls."
    penalty_initial_infeasible::T = 1.0

    "penalty update rate for infeasible controls."
    penalty_scaling_infeasible::T = 10.0

    # Minimum Time
    "regularization term for dt."
    R_minimum_time::T = 1.0

    "maximum allowable dt."
    dt_max::T = 1.0

    "minimum allowable dt."
    dt_min::T = 1.0e-3

    "initial penalty term for minimum time bounds constraints."
    penalty_initial_minimum_time_inequality::T = 1.0

    "initial penalty term for minimum time equality constraints."
    penalty_initial_minimum_time_equality::T = 1.0

    "penalty update rate for minimum time bounds constraints."
    penalty_scaling_minimum_time_inequality::T = 1.0

    "penalty update rate for minimum time equality constraints."
    penalty_scaling_minimum_time_equality::T = 1.0

    # Projected Newton
    "finish with a projecte newton solve."
    projected_newton::Bool = false

    "options for projected newton solver."
    opts_pn::ProjectedNewtonSolverOptions{T} = ProjectedNewtonSolverOptions{T}()

    "constraint satisfaction tolerance that triggers the projected newton solver.
    If set to a non-positive number it will kick out when the maximum penalty is reached."
    projected_newton_tolerance::T = 1.0e-3

end

"""$(TYPEDEF)
Augmented Lagrangian Trajectory Optimizer (ALTRO) is a solver developed by the Robotic Exploration Lab at Stanford University.
    The solver is special-cased to solve Markov Decision Processes by leveraging the internal problem structure.

ALTRO consists of two "phases":
1) AL-iLQR: iLQR is used with an Augmented Lagrangian framework to solve the problem quickly to rough constraint satisfaction
2) Projected Newton: A collocation-flavored active-set solver projects the solution from AL-iLQR onto the feasible subspace to achieve machine-precision constraint satisfaction.
"""
struct ALTROSolver{T} <: AbstractSolver{T}
    opts::ALTROSolverOptions{T}
    stats::Dict{Symbol,Any}
    solver_al::AugmentedLagrangianSolver{T}
    solver_pn::ProjectedNewtonSolver{T}
end

ALTROSolver(prob::Problem, opts::ALTROSolverOptions) = AbstractSolver(prob, opts)

function AbstractSolver(prob::Problem{T,D},opts::ALTROSolverOptions{T}) where {T<:AbstractFloat,D<:DynamicsType}
    to = TimerOutput()
    stats = Dict{Symbol,Any}(:timer=>to)
    solver_al = AugmentedLagrangianSolver(prob, opts.opts_al)
    solver_pn = ProjectedNewtonSolver(prob, opts.opts_pn)
    solver_al.stats[:timer] = to
    solver_al.solver_uncon.stats[:timer] = to
    solver_pn.stats[:timer] = to
    ALTROSolver{T}(opts,stats,solver_al,solver_pn)
end

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
#=== platform apis ===#
@ocl_func(clGetPlatformIDs, CL_int,
              (CL_uint, Ptr{CL_platform_id}, Ptr{CL_uint}))

@ocl_func(clGetPlatformInfo,
              CL_int, (CL_platform_id, CL_platform_info, Csize_t, Ptr{Void}, Ptr{Csize_t}))

#=== device apis ===#
@ocl_func(clGetDeviceIDs, CL_int,
              (CL_platform_id, CL_device_type, CL_uint, Ptr{CL_device_id}, Ptr{CL_uint}))

@ocl_func(clGetDeviceInfo, CL_int,
              (CL_device_id, CL_device_info, Csize_t, Ptr{Void}, Ptr{Csize_t}))

#=== context apis ===#
#TODO: pass user data as Any type
@ocl_func(clCreateContext, CL_context,
              (Ptr{CL_context_properties}, CL_uint, Ptr{CL_device_id}, CL_callback, CL_user_data, Ptr{CL_int}))

@ocl_func(clCreateContextFromType, CL_context,
              (Ptr{CL_context_properties}, CL_device_type, CL_callback, CL_user_data, Ptr{CL_int}))

@ocl_func(clRetainContext, CL_int, (CL_context,))

@ocl_func(clReleaseContext, CL_int, (CL_context,))

@ocl_func(clGetContextInfo, CL_int,
              (CL_context, CL_context_info, Csize_t, Ptr{Void}, Ptr{Csize_t}))

#=== command queue apis ===#
@ocl_func(clCreateCommandQueue, CL_command_queue,
              (CL_context, CL_device_id, CL_command_queue_properties, Ptr{CL_int}))

@ocl_func(clRetainCommandQueue, CL_int, (CL_command_queue,))

@ocl_func(clReleaseCommandQueue, CL_int, (CL_command_queue,))

@ocl_func(clGetCommandQueueInfo, CL_int,
              (CL_command_queue, CL_command_queue_info, Csize_t, Ptr{Void}, Ptr{Csize_t}))

#=== memory object apis ===#
@ocl_func(clCreateBuffer, CL_mem,
              (CL_context, CL_mem_flags, Csize_t, Ptr{Void}, Ptr{CL_int}))

@ocl_func(clRetainMemObject, CL_int, (CL_mem,))

@ocl_func(clReleaseMemObject, CL_int, (CL_mem,))

@ocl_func(clGetSupportedImageFormats, CL_int,
              (CL_context, CL_mem_flags, CL_mem_object_type, CL_uint, Ptr{CL_image_format}, Ptr{CL_uint}))

@ocl_func(clGetMemObjectInfo, CL_mem,
              (CL_mem, CL_mem_info, Csize_t, Ptr{Void}, Ptr{Csize_t}))

@ocl_func(clGetImageInfo, CL_mem,
              (CL_mem, CL_image_info, Csize_t, Ptr{Void}, Ptr{Csize_t}))

#=== sampler apis ===#
@ocl_func(clCreateSampler, CL_sampler,
              (CL_context, CL_bool, CL_addressing_mode, CL_filter_mode, Ptr{CL_int}))

@ocl_func(clRetainSampler, CL_int, (CL_sampler,))

@ocl_func(clReleaseSampler, CL_int, (CL_sampler,))

@ocl_func(clGetSamplerInfo, CL_int,
              (CL_sampler, CL_sampler_info, Csize_t, Ptr{Void}, Ptr{Csize_t}))

#=== program object apis ===#
@ocl_func(clCreateProgramWithSource, CL_program,
              (CL_context, CL_uint, Ptr{Ptr{Cchar}}, Ptr{Csize_t}, Ptr{CL_int}))

@ocl_func(clCreateProgramWithBinary, CL_program,
              (CL_context, CL_uint, Ptr{CL_device_id}, Ptr{Csize_t},
               Ptr{Ptr{Cuchar}}, Ptr{CL_int}, Ptr{CL_int}))

@ocl_func(clRetainProgram, CL_int, (CL_program,))

@ocl_func(clReleaseProgram, CL_int, (CL_program,))

@ocl_func(clBuildProgram, CL_int,
              (CL_program, CL_uint, Ptr{CL_device_id}, Ptr{Cchar}, CL_callback, Ptr{Void}))

@ocl_func(clGetProgramBuildInfo, CL_int,
              (CL_program, CL_device_id, CL_program_build_info, Csize_t, Ptr{Void}, Ptr{Csize_t}))

#=== kernel object apis ===#
@ocl_func(clCreateKernel, CL_kernel,
              (CL_program, Ptr{Cchar}, Ptr{CL_int}))

@ocl_func(clCreateKernelsInProgram, CL_int,
              (CL_program, CL_uint, Ptr{CL_kernel}, Ptr{CL_uint}))

@ocl_func(clRetainKernel, CL_int, (CL_kernel,))

@ocl_func(clReleaseKernel, CL_int, (CL_kernel,))

@ocl_func(clSetKernelArg, CL_int,
              (CL_kernel, CL_uint, Csize_t, Ptr{Void}))

@ocl_func(clGetKernelInfo, CL_int,
              (CL_kernel, CL_kernel_info, Csize_t, Ptr{Void}, Ptr{Csize_t}))

@ocl_func(clGetKernelWorkGroupInfo, CL_int,
              (CL_kernel, CL_device_id, CL_kernel_work_group_info, Csize_t, Ptr{Void}, Ptr{Csize_t}))

#=== event object apis ===#
@ocl_func(clWaitForEvents, CL_int,
              (CL_uint, Ptr{CL_event_info}))

@ocl_func(clGetEventInfo, CL_int,
              (CL_event, CL_event_info, Csize_t, Ptr{Void}, Ptr{Csize_t}))

@ocl_func(clRetainEvent, CL_int, (CL_event,))

@ocl_func(clReleaseEvent, CL_int, (CL_event,))

#=== profiling apis ===#
@ocl_func(clGetEventProfilingInfo, CL_int,
              (CL_event, CL_profiling_info, Csize_t, Ptr{Void}, Ptr{Csize_t}))

#=== flush and finish apis ===#
@ocl_func(clFlush, CL_int, (CL_command_queue,))

@ocl_func(clFinish, CL_int, (CL_command_queue,))

#=== enqueued commands apis ===#
@ocl_func(clEnqueueReadBuffer, CL_int,
              (CL_command_queue, CL_mem, CL_bool, Csize_t, Csize_t, Ptr{Void},
               CL_uint, Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clEnqueueWriteBuffer, CL_int,
              (CL_command_queue, CL_mem, CL_bool,
               Csize_t, Csize_t, Ptr{Void}, CL_uint,
               Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clEnqueueCopyBuffer, CL_int,
              (CL_command_queue, CL_mem, CL_mem,
               Csize_t, Csize_t, Csize_t, CL_uint, Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clEnqueueReadImage, CL_int,
              (CL_command_queue, CL_mem, CL_bool,
               Ptr{Csize_t}, Ptr{Csize_t}, Csize_t, Csize_t,
               Ptr{Void}, CL_uint, Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clEnqueueWriteImage, CL_int,
              (CL_command_queue, CL_mem, CL_bool, Ptr{Csize_t}, Ptr{Csize_t},
               Csize_t, Csize_t, Ptr{Void}, CL_uint, Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clEnqueueCopyImage, CL_int,
              (CL_command_queue, CL_mem, CL_mem, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Csize_t},
               CL_uint, Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clEnqueueCopyImageToBuffer, CL_int,
               (CL_command_queue, CL_mem, CL_mem, Ptr{Csize_t}, Ptr{Csize_t},
                Csize_t, CL_uint, Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clEnqueueCopyBufferToImage, CL_int,
               (CL_command_queue, CL_mem, CL_mem, Csize_t, Ptr{Csize_t}, Ptr{Csize_t},
                CL_uint, Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clEnqueueMapBuffer, Ptr{Void},
              (CL_command_queue, CL_mem, CL_bool, CL_map_flags, Csize_t, Csize_t,
               CL_uint, Ptr{CL_event}, Ptr{CL_event}, Ptr{CL_int}))

@ocl_func(clEnqueueMapImage, Ptr{Void},
              (CL_command_queue, CL_mem, CL_bool, CL_map_flags,
               Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Csize_t},
               CL_uint, Ptr{CL_event}, Ptr{CL_event}, Ptr{CL_int}))

@ocl_func(clEnqueueUnmapMemObject, CL_int,
              (CL_command_queue, CL_mem, Ptr{Void}, CL_uint, Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clEnqueueNDRangeKernel, CL_int,
              (CL_command_queue, CL_kernel, CL_uint,
               Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Csize_t},
               CL_uint, Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clEnqueueTask, CL_int,
              (CL_command_queue, CL_kernel, CL_uint, Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clEnqueueNativeKernel, CL_int,
              (CL_command_queue, Ptr{Void}, Csize_t, CL_uint,
               Ptr{CL_mem}, Ptr{Ptr{Void}}, CL_uint, Ptr{CL_event}, Ptr{CL_event}))

#== opengl interop functions ==#

@ocl_func(clEnqueueAcquireGLObjects, CL_int,
              (CL_command_queue, CL_uint, Ptr{CL_mem}, CL_uint, Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clEnqueueReleaseGLObjects, CL_int,
              (CL_command_queue, CL_uint, Ptr{CL_mem}, CL_uint, Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clCreateFromGLBuffer, CL_mem,
              (CL_context, CL_mem_flags, GL_uint, Ptr{CL_int}))

@ocl_func(clCreateFromGLRenderbuffer, CL_mem,
              (CL_context, CL_mem_flags, GL_uint, Ptr{CL_int}))

@ocl_func(clCreateFromGLTexture2D, CL_mem,
              (CL_context, CL_mem_flags, GL_enum, GL_int, GL_uint, Ptr{CL_int}))

@ocl_func(clCreateFromGLTexture3D, CL_mem,
              (CL_context, CL_mem_flags, GL_enum, GL_int, GL_uint, Ptr{CL_int}))

@ocl_func(clGetGLObjectInfo, CL_int,
              (CL_mem, Ptr{CL_GL_object_type}, Ptr{GL_uint}))

@ocl_func(clGetGLTextureInfo, CL_int,
              (CL_mem, CL_GL_texture_info, Csize_t, Ptr{Void}, Ptr{Csize_t}))


text:
using SymbolEquations
using Test

@testset "SymbolEquations.jl" begin
    # Write your tests here.
end


text:
using Terming

function println_animation(str::String; delay::Real=0.05)
    # save current position of cursor
    Terming.csave()

    for c in str
        Terming.print(c)
        sleep(delay)
        delay *= 0.98
    end

    # restore position of cursor and move down
    Terming.crestore(); Terming.cmove_down()
end

function main()
    # set term size and clear
    Terming.displaysize(20, 80); Terming.clear()
    # move cursor to (row=2, col=2)
    Terming.cmove(2, 2)

    # +----------------+
    # | without buffer |
    # +----------------+
    discription = "The following string is blocked bue to the time-consuming calculations:"
    println_animation(discription)

    # save current position of cursor
    Terming.csave()

    str = "This string will be finished printing once the calculations are..."
    Terming.print(str)
    sleep(1) # fake time consuming calculation
    Terming.println(" done!!")

    # restore position of cursor and move down 2 row
    Terming.crestore(); Terming.cmove_down(2)

    # +-------------+
    # | with buffer |
    # +-------------+
    discription = "The following string is not blocked by the time-consuming calculations:"
    println_animation(discription)

    # save current position of cursor
    Terming.csave()

    Terming.buffered() do buffer
        str = "This string will be finished printing once the calculations are"
        Terming.print(buffer, str)
        sleep(1) # fake time consuming calculation
        Terming.println(buffer, " done!!")
    end

    # restore position of cursor and move down 2 row
    Terming.crestore(); Terming.cmove_down(2)

    return
end

main()


text:
module Bifurcations

# Re-export methods from DifferentialEquations extended here:
export init, solve, solve!, step!
using DiffEqBase: init, solve, solve!, step!
# see: continuations/solver.jl

export solving!

using FillArrays: Eye

include("utils/utils.jl")

include("continuations/continuations.jl")
using .Continuations: AbstractContinuationProblem, AbstractProblemCache, solving!
import .Continuations: get_prob_cache, get_u0, residual!, residual_jacobian!,
    residual, isindomain

include("base/base.jl")
# TimeKind trait:
using .BifurcationsBase: timekind, Continuous, Discrete
import .BifurcationsBase: TimeKind
# StateKind trait:
using .BifurcationsBase: statekind, MutableState, ImmutableState
import .BifurcationsBase: StateKind

using .BifurcationsBase: BifurcationProblem
export BifurcationProblem

include("codim1/codim1.jl")

include("fixedpoint.jl")
include("api.jl")
include("diffeq.jl")

include("codim1lc/codim1lc.jl")
using .Codim1LimitCycle: LimitCycleProblem

include("codim2/codim2.jl")
include("codim2lc/codim2lc.jl")

include("examples/examples.jl")

# using Requires
# @require RecipesBase include("plotting.jl")
include("plotting.jl")
export plot_state_space, plot_state_space!

end # module


text:
cd(@__DIR__)
using Pkg
Pkg.activate(@__DIR__)
CI = get(ENV, "CI", nothing) == "true" || get(ENV, "GITHUB_TOKEN", nothing) !== nothing
CI && Pkg.instantiate()

using InteractiveDynamics
using DynamicalSystems, DynamicalBilliards, Agents
using Documenter
using DocumenterTools: Themes

# %%
# download the themes
for file in ("juliadynamics-lightdefs.scss", "juliadynamics-darkdefs.scss", "juliadynamics-style.scss")
    download("https://raw.githubusercontent.com/JuliaDynamics/doctheme/master/$file", joinpath(@__DIR__, file))
end
# create the themes
for w in ("light", "dark")
    header = read(joinpath(@__DIR__, "juliadynamics-style.scss"), String)
    theme = read(joinpath(@__DIR__, "juliadynamics-$(w)defs.scss"), String)
    write(joinpath(@__DIR__, "juliadynamics-$(w).scss"), header*"\n"*theme)
end
# compile the themes
Themes.compile(joinpath(@__DIR__, "juliadynamics-light.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-light.css"))
Themes.compile(joinpath(@__DIR__, "juliadynamics-dark.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-dark.css"))

makedocs(
modules=[InteractiveDynamics, DynamicalSystems, DynamicalBilliards, Agents],
doctest=false,
sitename= "InteractiveDynamics",
root = @__DIR__,
format = Documenter.HTML(
    prettyurls = CI,
    assets = [
        asset("https://fonts.googleapis.com/css?family=Montserrat|Source+Code+Pro&display=swap", class=:css),
        ],
    ),
pages = [
    "Introduction" => "index.md",
    "Trajectory Evolution" => "trajectory.md",
    "Orbit Diagram" => "od.md",
    "Cobweb Diagram" => "cobweb.md",
    "Poincaré Surface of Section" => "psos.md",
    "Brainscanning PSOS" => "brainscan.md",
    "Billiards" => "billiards.md",
    "Agent Based Models" => "agents.md",
    # "Trajectory Highlighter" => "highlight.md",
],
)

if CI
    deploydocs(
        repo = "github.com/JuliaDynamics/InteractiveDynamics.jl.git",
        target = "build",
        push_preview = true
    )
end


text:
struct Disks
	R::Vector{Float64}
	PA::Vector{Float64}
	I::Vector{Float64}
	V::Vector{Float64}
end

# struct Geometry
# 	R::Vector{Float64}
# 	Rs::Vector{Float64}
# 	PA::Vector{Float64}
# 	I::Vector{Float64}
# 	M::Vector{Float64}
# 	X::Vector{Float64}
# 	Y::Vector{Float64}
# 	V::Vector{Float64}
# 	VC::Vector{Float64}
# 	dV::Vector{Float64}
# 	Z::Vector{Float64}
# 	Φ::Vector{Float64}
# 	dT::Vector{Float64}
# 	T::Vector{Float64}
# end

"""
Rotation Matrix
"""
function rot2d(angle::Float64)::Array{Float64,2}
    s, c = sin(angle), cos(angle)
    return [c -s ; s c]
end
#
"""
Minimum Distance between an ellipse and a point
algorithm from https://blog.chatfield.io/simple-method-for-distance-to-ellipse/

Input
======
- ``a`` major axis of the ellipse
- ``inclination`` inclination of the ellipse
- ``pa`` position angle of the ellipse
- ``p`` xy vector of the point

Parameters
=======
- ``debug`` plot

Output
====
distance, xy vector of the closest point on the ellipse
"""
function distf(a::Float64, inclination::Float64, pa::Float64, p::Vector{Float64}; debug=false,tx=0.70710678118654752,ty=0.70710678118654752)
    angle=pa+pi/2.0
    p_rot = rot2d(-angle)*p
    p_abs=abs.(p_rot)

    b= abs(cos(inclination)*a)
    ss_sub = a^2.0 - b^2.0
    #local x=0.0 # =Float64()
    #y =Float64()
    for i in 1:3
        @fastmath x = tx*a
        @fastmath y = ty*b
        @fastmath eM = ss_sub*tx^3.0/a
        @fastmath em = -ss_sub*ty^3.0/b
        #eM = ss_sub*cos(tx)^3.0/semi_major
        #em = -ss_sub*cos(ty)^3.0/semi_minor
        rx=x-eM
        ry=y-em

        qx=p_abs[1]-eM
        qy=p_abs[2]-em

        r=hypot(rx,ry)
        q=hypot(qx,qy)

        # tx= clamp((qx*r/q+eM)/semi_major,0.0,1.0)
  #       ty= clamp((qy*r/q+em)/semi_minor,0.0,1.0)
        @fastmath tx= (qx*r/q+eM)/a
        @fastmath ty= (qy*r/q+em)/b
        tn=hypot(tx,ty)
        tx /= tn
        ty /= tn
        #@show x,y,tx,ty
    end
    x=a*tx
    y=b*ty

    x_edge=copysign(x, p_rot[1])
    y_edge=copysign(y, p_rot[2])
    p_edge = [x_edge,y_edge]
    dst = norm(p_edge .- p_rot, 2)
    p_edge_rot = rot2d(angle)*p_edge
    if debug
        fig=Figure()
        ax = Axis(fig[1, 1])#, xlabel = "x label", ylabel = "y label",title = "Title")
        ax.aspect = DataAspect()
        phis=0.0:0.01:2*pi
        R=a ./ sqrt.(1.0 .+ tan(inclination)^2.0 .* cos.(phis .-pa).^2.0)
        xd=R.*cos.(phis)
        yd=R .*sin.(phis)
        lines!(ax,xd,yd,ratio=1,label=false)
        GLMakie.scatter!(ax,[p[1]],[p[2]],label=false)
        #scatter!(pl,[p_rot[1]],[p_rot[2]],label="p_rot")
        GLMakie.scatter!(ax,[p_edge_rot[1]],[p_edge_rot[2]],label=false)
        lines!(ax,[p[1],p_edge_rot[1]],[p[2],p_edge_rot[2]])
        @show norm(p-p_edge_rot,2)
        @show dst, p_edge_rot
        display(fig)
    end
    return dst, p_edge_rot
end

"""
Disk Creation routine

Input
======
- ``Rd`` vector of radii
- ``P`` Parameters struct which contains function ppi(R; I params),ppa(R; PA params),Vcirc(R;V params)

Output
====
Disks struct (R,PA(R),I(R),V(R))
"""
function make_disks(Rd::Vector{Float64},P)::Disks

	PA = [getfield(P,s) for s in P.PA]
	I = [getfield(P,s) for s in P.I]
	V = [getfield(P,s) for s in P.V]

	PAd=deg2rad.(P.ppi(P.R_i,PA).(Rd))
	Id=deg2rad.(P.ppa(P.R_pa,I).(Rd))
	Vd=P.Vcirc.(Rd,V...)

	return Disks(Rd,PAd,Id,Vd)
end

function distrv(x::Float64,y::Float64,v::Float64,disks::Disks,dR::Float64;dvi_min=300.0)::Float64
	dVi=dvi_min
	for (rd,i_d,phi_d,V_d) in zip(disks.R,disks.I,disks.PA,disks.V)
		d,xy_ellipse=distf(rd, i_d, phi_d, [x,y])
		if d<=dR
			phi_ell=atan(y,x)
			Vd_sky= V_d* sin(i_d)*sin(phi_ell-phi_d)/sqrt(1.0 +tan(i_d)^2.0*cos(phi_ell-phi_d)^2.0)
			dvi=abs(Vd_sky-v)
			if dvi<dvi_min
				dvi_min=dvi
				dVi=Vd_sky-v
			end
		end
	end
	return dVi
end
#X::Vector{Float64}, Y::Vector{Float64}, V::Vector{Float64}
function cloud_find(disks::Disks,clouds::Clouds,P;dV0=300.0)::Vector{Float64}
	dV=fill(dV0,length(clouds.Xc))
    Threads.@threads for i in 1:length(clouds.Xc)
		x=clouds.Xc[i]
		y=clouds.Yc[i]
		v=clouds.Vc[i]
		@inbounds dV[i]=distrv(x,y,v,disks,P.dR;dvi_min=300.0)
    end
	return dV
end

function update_parameters!(P,pars::Vector{Float64})
	for (p,s) in zip(pars,P.Fit)
		setfield!(P,s,p)
	end
end

function likelihood(clouds::Clouds,P; dV0=300.0,dv=30.0,dv_hi=30.0):Float64
	dV=cloud_find(make_disks(P.Rd,P),clouds,P)
	return sum(dV .^2.0) /dv^2.0
end

function logprior(P,priors)::Float64
	ps=0.0
	for p in fieldnames(typeof(priors))
		pr=getfield(priors,p)
		ps += (getfield(P,p) - pr[1] )^2.0/ pr[2]^2.0
	end
	return ps
end

function cloud_properties(x::Float64,y::Float64,v::Float64,disks::Disks,dR::Float64;τ=nothing,dvi_min=300.0)::NTuple{8, Float64}
	dvc=dvi_min

	zd=NaN
	ic=NaN
	vc=NaN
    rc=NaN
	phic=NaN
	dT=NaN
	Tnet=NaN
    k_best=0
    found=false
	for (k,(rd,i_d,phi_d,V_d)) in enumerate(zip(disks.R,disks.I,disks.PA,disks.V))
		d,xy_ellipse=distf(rd, i_d, phi_d, [x,y])
		if d<=dR
			phi_ell=atan(y,x)
			Vd_sky= V_d* sin(i_d)*sin(phi_ell-phi_d)/sqrt(1.0 +tan(i_d)^2.0*cos(phi_ell-phi_d)^2.0)

			dvi=abs(v-Vd_sky)
			if dvi<dvi_min
				dvi_min=dvi
				dvc=v-Vd_sky
				vc=V_d
				zd = -tan(i_d)*(xy_ellipse[1]*cos(phi_d)+xy_ellipse[2]*sin(phi_d))
				ic=i_d
				phic=phi_d
				rc=rd
                k_best=k
				found =true
			end
		end
	end
	if found
        if isnothing(τ)
            return dvc,vc,zd,ic,phic,NaN,NaN,k_best
        else
            Wd_ = [sin(ic)*cos(phic),sin(ic)*sin(phic),cos(ic)]
            Nd_= Wd_ ./norm(Wd_)
            Id_= [-sin(phic),cos(phic),0.0]
            Id_=Id_ ./norm(Id_)
            Md_=[cos(ic)*cos(phic),cos(ic)*sin(phic),-sin(ic)]#cross(Nd_,Id_)
            Md_=Md_ ./norm(Md_)

            dT=dot(τ([x,y,zd]),Nd_)
            τ_net(ϕ)=dot(τ(rc .*cos(ϕ) .* Id_ .+ rc .* sin(ϕ) .*Md_),Nd_)
            Tnet=rc*3.0857e16  /vc * quadgk(τ_net,0.0,2*pi;atol=0.1)[1]
            return dvc,vc,zd,ic,phic,dT,Tnet,k_best
        end
	else
		return NaN,NaN,NaN,NaN,NaN,NaN,NaN,0
	end
end

function cloud_geometry!(clouds::Clouds,disks::Disks,P;dV0=300.0,τ=nothing)#::Geometry

    X=clouds.Xp
    Y=clouds.Yp
    V=clouds.Vp
    M=clouds.Ip

	VC=fill(NaN,length(X))
    Z=fill(NaN,length(X))
	K=fill(NaN,length(X))
	dT=fill(NaN,length(X))
	Tnet=fill(NaN,length(X))
    PA=fill(NaN,length(X))
    I=fill(NaN,length(X))
    K=zeros(Int64,length(X))
    FF=Bool.(zeros(length(X)))
	dV=fill(dV0,length(X))
    Threads.@threads for j in 1:length(X)
		x=X[j]
		y=Y[j]
		v=V[j]
        dvi_min=dV0
		dV[j],VC[j],Z[j],I[j],PA[j],dT[j],Tnet[j],K[j]=cloud_properties(x,y,v,disks,P.dR;τ=τ,dvi_min=dV0)
    end

    clouds.r=sqrt.(X .^2.0 .+Y .^2.0 )
    clouds.R=sqrt.(X .^2.0 .+Y .^2.0 .+Z .^2.0)
    clouds.Φ=atan.(Y,X)
    clouds.Z=Z
    clouds.inc=I
    clouds.pa=PA
    clouds.dT=dT
    clouds.Tnet=Tnet
    clouds.VC=VC
    clouds.dV=dV
    clouds.k=K
end

function manipulate_cube(ix::Int64,x::Float64,iy::Int64,y::Float64,cube::Array{Float64, 3},model_cube::Array{Float64, 3},data::Dict{String, Any},disks::Disks,search_v::Float64,Sd::Function,Sdpars::Vector{Float64},Rd::Vector{Float64},dR::Float64,lowest_sig::Float64)
    for (rd,i_d,phi_d,V_d) in zip(disks.R,disks.I,disks.PA,disks.V)
        d,xy_ellipse=distf(rd, i_d, phi_d, [x,y])
        if d<=dR
            phi_ell=atan(y,x)#xy_ellipse[2],xy_ellipse[1])
            Vd_sky=V_d*sin(i_d)*sin(phi_ell-phi_d)/sqrt(1.0 +tan(i_d)^2.0*cos(phi_ell-phi_d)^2.0)
            loc=abs.(Vd_sky.-data["V"]) .<search_v
            if any(loc)
                lsic=cube[ix,iy,:]
                local_sig=lsic[loc]
                Vlocal=data["V"][loc]
                Fdata=maximum(local_sig)
                Vmax=Vlocal[argmax(local_sig)]
                if Fdata > lowest_sig*data["rms"]
                    msig=Fdata .*exp.(-(data["V"] .-Vmax) .^2.0 ./(2.0*Sd(rd,Sdpars...)^2.0))
                    model_cube[ix,iy,:]+=msig
                    cube[ix,iy,:] -=msig
                end
            end
        end
    end
    return cube,model_cube
end

function create_model_cube(data::Dict{String, Any}, IX::Vector{Int64}, XX::Vector{Float64},IY::Vector{Int64},YY::Vector{Float64},Sd::Function,Sdpars::Vector{Float64},Rd::Vector{Float64},dR::Float64,disks::Disks;search_v=1,beam=nothing,lowest_sig=3.0)::Array{Float64, 3}
	cube=copy(data["data"])
	model_cube=fill(0.0,size(cube))
	Threads.@threads for n in 1:length(IX)
		ix=IX[n]
		x=XX[n]
		iy=IY[n]
		y=YY[n]
        cube,model_cube=manipulate_cube(ix,x,iy,y,cube,model_cube,data,disks,search_v,Sd,Sdpars,Rd,dR,lowest_sig)
	end
	if ~isnothing(beam)
		for iv in eachindex(data["V"])
			model_cube[:,:,iv]=imfilter(model_cube[:,:,iv],beam)
		end
	end
	return model_cube
end

# function create_model_cube(data::Dict{String, Any}, IX::Vector{Int64}, XX::Vector{Float64},IY::Vector{Int64},YY::Vector{Float64},Sd::Function,Sdpars::Vector{Float64},Rd::Vector{Float64},dR::Float64,disks::Disks;search_v=1,beam=nothing,lowest_sig=3.0)::Array{Float64, 3}
# 	cube=copy(data["data"])
# 	model_cube=fill(0.0,size(cube))
# 	@progress for n in 1:length(IX)
# 		ix=IX[n]
# 		x=XX[n]
# 		iy=IY[n]
# 		y=YY[n]
# 		for (rd,i_d,phi_d,V_d) in zip(disks.R,disks.I,disks.PA,disks.V)
# 			d,xy_ellipse=distf(rd, i_d, phi_d, [x,y])
# 			if d<=dR
# 				phi_ell=atan(y,x)#xy_ellipse[2],xy_ellipse[1])
# 				Vd_sky=V_d*sin(i_d)*sin(phi_ell-phi_d)/sqrt(1.0 +tan(i_d)^2.0*cos(phi_ell-phi_d)^2.0)
# 				loc=abs.(Vd_sky.-data["V"]) .<search_v
# 				if any(loc)
# 					lsic=cube[ix,iy,:]
# 					local_sig=lsic[loc]
# 					Vlocal=data["V"][loc]
# 					Fdata=maximum(local_sig)
# 					Vmax=Vlocal[argmax(local_sig)]
# 					if Fdata > lowest_sig*data["rms"]
# 						msig=Fdata .*exp.(-(data["V"] .-Vmax) .^2.0 ./(2.0*Sd(rd,Sdpars...)^2.0))
# 						model_cube[ix,iy,:]+=msig
# 						cube[ix,iy,:] -=msig
# 					end
# 				end
# 			end
# 		end
# 	end
# 	if ~isnothing(beam)
# 		for iv in eachindex(data["V"])
# 			model_cube[:,:,iv]=imfilter(model_cube[:,:,iv],beam)
# 		end
# 	end
# 	return model_cube
# end


text:
using ProtoGrad
using Test

@testset "Fine-tuning" begin
    @testset "$(T)" for T in [Float32, Float64]
        input_size, hidden_size, output_size = 3, 2, 1
        batch_size = 4

        m1 = ProtoGrad.Linear(T, input_size=>hidden_size)
        m3 = ProtoGrad.Linear(T, hidden_size=>output_size)

        function get_complete_model(b)
            ProtoGrad.Compose(
                m1,
                ProtoGrad.ReLU(),
                b,
                x -> ProtoGrad.relu.(x),
                m3
            )
        end

        m2 = ProtoGrad.Linear(T, hidden_size=>hidden_size)
        m2_orig = copy(m2)

        W_true = randn(T, output_size, input_size)
        b_true = randn(T, output_size)

        x = randn(T, input_size, batch_size)
        y = W_true * x .+ b_true + randn(T, output_size, batch_size)

        data_iter = Iterators.repeated((x, y))
        f = ProtoGrad.SupervisedObjective(ProtoGrad.mse, data_iter) ∘ get_complete_model

        optimizer = ProtoGrad.GradientDescent(stepsize=0.1)

        m2_final = nothing

        for output in Iterators.take(optimizer(m2, f), 10)
            @test length(ProtoGrad.overlap(output.solution, m2)) == 0
            @test length(ProtoGrad.overlap(output.solution, m2_orig)) == 0
            m2_final = output.solution
        end

        @test length(ProtoGrad.overlap(m2_final, m2_orig)) == 0
        @test_skip m2_final.W != m2_orig.W
        @test_skip m2_final.b != m2_orig.b
    end
end


text:
@inline function setnextprevself!(x, tag::Val)
    setnext!(x, x, tag)
    setprev!(x, x, tag)
    nothing
end

@inline function circlist_link_before!(x::T, node::T, tag::Val) where {T}
    prev = getprev(x, tag)::T
    setnext!(prev, node, tag)
    setprev!(node, prev, tag)
    setnext!(node, x, tag)
    setprev!(x, node, tag)
    nothing
end

@inline function circlist_unlink!(x::T, tag::Val) where {T}
    prev = getprev(x, tag)::T
    next = getnext(x, tag)::T
    setnext!(prev, next, tag)
    setprev!(next, prev, tag)
    next
end

@inline function circslist_prev(node::T, start::T, tag::Val) where {T}
    prev = start
    x = getnext(start, tag)::T
    while x !== node
        prev = x
        x = getnext(x, tag)::T
    end
    prev
end


text:
#!/usr/bin/env julia

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "config"
    help = "a YAML configuration file"
    default = "config.yaml"
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

# Make some diagnostic plots to show what the model looks like in analytic terms

using DiskJockey.hmodel
using DiskJockey.constants
import PyPlot.plt
using LaTeXStrings
using Dierckx

# velocity structure
function plot_vel(pars::Parameters, grid)

    vel = Array(Float64, (nz,nr))
    for j=1:nz
        for i=1:nr
            # convert from cm/s to km/s
            vel[j, i] = DiskJockey.hmodel.velocity(rs[i], zs[j], pars) .* 1e-5
        end
    end

    fig = plt[:figure]()
    ax = fig[:add_subplot](111)

    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.77)

    img = ax[:contourf](xx./AU, yy./AU, vel) #, levels=levels)

    cax = fig[:add_axes]([0.82, 0.22, 0.03, 0.65])
    cb = fig[:colorbar](img, cax=cax)

    ax[:set_ylabel](L"$z$ [AU]")
    ax[:set_xlabel](L"$r$ [AU]")
    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.85)

    plt[:savefig]("velocity.png")
end

# temperature structure
function plot_temp(pars::Parameters, grid)

    Ts = Array(Float64, (nz,nr));
    for j=1:nz
        for i=1:nr
            Ts[j, i] = DiskJockey.hmodel.temperature(rs[i], zs[j], pars)
        end
    end

    fig = plt[:figure](figsize=(6,4))
    ax = fig[:add_subplot](111)

    levels = logspace(0, log10(200), 50);
    img = ax[:contourf](xx./AU, yy./AU, Ts, cmap=plt[:get_cmap]("Blues_r"), levels=levels)
    # ax[:plot](rs./AU, z_qs./AU, "r--")

    cax = fig[:add_axes]([0.82, 0.22, 0.03, 0.65])
    cb = fig[:colorbar](img, cax=cax)
    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.77)

    ax[:set_ylabel](L"$z$ [AU]")
    ax[:set_xlabel](L"$r$ [AU]")

    fig[:savefig]("temperature.png")

end

# scale height
function plot_height(pars::Parameters, grid)

    heights = DiskJockey.hmodel.Hp(grid.rs, pars) ./ AU

    fig = plt[:figure]()
    ax = fig[:add_subplot](111)
    ax[:plot](rr, heights)
    ax[:set_ylabel](L"$H_p$ [AU]")
    ax[:set_xlabel](L"$r$ [AU]")
    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.85)

    plt[:savefig]("scale_height.png")
end

# Atmospheric height
function plot_atm_height(pars::Parameters, grid)
    heights = DiskJockey.hmodel.z_q(grid.rs, pars) ./ AU

    fig = plt[:figure]()
    ax = fig[:add_subplot](111)
    ax[:plot](rr, heights)
    ax[:set_ylabel](L"$z_q$ [AU]")
    ax[:set_xlabel](L"$r$ [AU]")
    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.85)

    plt[:savefig]("atm_height.png")
end

function plot_surface_density(pars::Parameters, grid)

    Sigmas = DiskJockey.hmodel.Sigma(grid.rs, pars)

    fig = plt[:figure]()
    ax = fig[:add_subplot](111)
    ax[:loglog](rr, Sigmas)

    # Now, go overlay small grey lines vertically
    for cell_edge in grid.Rs/AU
        ax[:axvline](cell_edge, color="0.5", lw=0.4)
    end

    ax[:set_ylabel](L"$\Sigma\, [\mathrm{g/cm}^2]$")
    ax[:set_xlabel](L"$r$ [AU]")
    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.85)

    plt[:savefig]("surface_density.png")
end

function plot_density_slice(pars::Parameters)

    # Plot this at a few radii.
    rs = Float64[10., 20., 50., 100.]

    fig = plt[:figure]()
    ax = fig[:add_subplot](111)

    for r in rs
        zs, slice = DiskJockey.hmodel.density_slice(r * AU, pars)
        # ax[:semilogy](zs ./AU, slice)
        ax[:semilogy](zs ./AU, slice/(mu_gas * m_H))
    end
    ax[:set_xlabel](L"$z$ [AU]")
    ax[:set_ylabel](L"$\rho$ [$g/{cm{^3}}$]")
    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.85)

    plt[:savefig]("density_slices.png")

end

function plot_density_gradient(pars::Parameters)

    nz = 500
    zs = linspace(0, 5 * AU, nz)
    R = 0.5 * AU
    y = Array(Float64, nz)
    for i=1:nz
        y[i] = DiskJockey.hmodel.dlnrho(R, zs[i], pars)
    end

    fig = plt.figure()
    ax = fig[:add_subplot](111)

    ax[:plot](zs/AU, y * AU)

    ax[:set_ylabel](L"$\partial \ln \rho_\mathrm{gas} / \partial z$ [1/AU]")
    ax[:set_xlabel](L"$z$ [AU]")
    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.85)

    plt.savefig("density_gradient.png")

end

function plot_density_1D_unnormed(pars::Parameters)
    nz = 50
    zs = linspace(0, 14 * AU, nz)
    R = 10 * AU
    y = Array(Float64, nz)
    for i=1:nz
        y[i] = DiskJockey.hmodel.un_lnrho(R, zs[i], pars)
    end

    fig = plt.figure()
    ax = fig[:add_subplot](111)

    ax[:plot](zs/AU, y)

    ax[:set_ylabel](L"$\ln \rho_\mathrm{gas}$ unnormalized")
    ax[:set_xlabel](L"$z$ [AU]")
    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.85)

    plt.savefig("un_ln_gas_slice.png")

end

function plot_dlnrho(pars::Parameters)
    nz = 50
    zs = linspace(0, 200 * AU, nz)
    R = 10 * AU
    dlnrhos = Array(Float64, nz)
    for i=1:nz
        dlnrhos[i] = DiskJockey.hmodel.dlnrho(R, zs[i], pars)
    end

    fig = plt.figure()
    ax = fig[:add_subplot](111)

    ax[:plot](zs/AU, dlnrhos)

    ax[:set_ylabel](L"$d \ln \rho$")
    ax[:set_xlabel](L"$z$ [AU]")
    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.85)

    fig[:savefig]("dlnrhos.png")
end

function plot_density_1D(pars::Parameters)
    nz = 50
    zs = linspace(0, 30 * AU, nz)
    R = 0.5 * AU
    y = Array(Float64, nz)
    for i=1:nz
        y[i] = DiskJockey.hmodel.un_lnrho(R, zs[i], pars)
    end

    # Now calculate the correction factor here
    lncor = exp(DiskJockey.hmodel.lncorrection_factor(R, pars))

    lngas = lncor +  y - log(mu_gas * m_H)

    fig = plt.figure()
    ax = fig[:add_subplot](111)

    ax[:plot](zs/AU, lngas)

    ax[:set_ylabel](L"$\ln n_\mathrm{gas}$")
    ax[:set_xlabel](L"$z$ [AU]")
    ax[:set_xlim](0,30)
    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.85)

    plt.savefig("ln_gas_slice.png")
end

# Go through some trial points and plot the interpolator to see how it's doing.
function plot_interpolator(pars::Parameters)

    nr = grid.nr
    rs = grid.rs

    lncors = Array(Float64, nr)

    for i=1:nr
        lncors[i] = DiskJockey.hmodel.lncorrection_factor(rs[i], pars)
    end

    # Now, let's use the interpolator to make finer samples and see if they line up.
    spl = DiskJockey.hmodel.make_correction_interpolator(pars, grid)

    # Now, come up with a finer spaceing of grid points
    n_fine = 100
    rs_fine = logspace(log10(rs[1]), log10(rs[end]), n_fine)
    lncors_fine = Array(Float64, n_fine)
    for i=1:n_fine
        lncors_fine[i] = evaluate(spl, rs_fine[i])
    end

    fig = plt.figure()
    ax = fig[:add_subplot](111)

    ax[:plot](rs/AU, lncors)
    ax[:plot](rs_fine/AU, lncors_fine, "o")

    ax[:set_xlabel](L"$r$ [AU]")
    ax[:set_ylabel]("ln correction factor")

    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.85)

    # ax[:set_xlim](0, 20)

    plt.savefig("cor_factor.png")

end

# density structure
function plot_dens(pars::Parameters, grid)

    nr = grid.nr
    nz = 64
    rs = grid.rs
    zq = DiskJockey.hmodel.z_q(rs[end], pars)

    zs = cat(1, [0], logspace(log10(0.1 * AU), log10(1.5 * zq), nz-1))

    xx = Array(Float64, (nz, nr));
    yy = Array(Float64, (nz, nr));
    for i=1:nz
        xx[i, :] = rs
    end

    for j=1:nr
        yy[:, j] = zs
    end

    println("Making correction interpolator")
    tic()
    spl = DiskJockey.hmodel.make_correction_interpolator(pars, grid)
    toc()
    println("Finished correction interpolator")

    rhos = Array(Float64, (nz,nr))
    for j=1:nz
        for i=1:nr
            temp = DiskJockey.hmodel.temperature(rs[i], zs[j], pars)
            XF = DiskJockey.hmodel.X_freeze(temp, pars)
            rhos[j, i] = XF * DiskJockey.hmodel.rho_gas(rs[i], zs[j], pars, spl)
        end
    end

    nlog = log10(rhos/(mu_gas * m_H))
    #Add a tiny bit to nn to prevent log10(0.0)

    levels = Float64[4.0, 5, 6, 7, 8, 9]
    # levels = linspace(4, 12., 8)

    fig = plt.figure()
    ax = fig[:add_subplot](111)
    ax[:set_ylabel](L"$z$ [AU]")
    ax[:set_xlabel](L"$r$ [AU]")
    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.77)

    img = ax[:contourf](xx./AU, yy./AU, nlog, levels=levels)

    cax = fig[:add_axes]([0.82, 0.22, 0.03, 0.65])
    cb = fig[:colorbar](img, cax=cax)

    println("Minimum ", minimum(nlog))
    println("Maximum ", maximum(nlog))

    #ticks = np.linspace(0, np.max(cov), num=6)
    #cb.set_ticks(ticks)
    #cb.set_ticks(MaxNLocator(nbins=5))

    #Plot the contoured density in cylindrical coordinates, then plot the spherical grid on top of it?
    # Do this by plotting a bunch of lines
    # First, the radial lines

    # Basically, each of these originates from x = 0, z = 0, and has a slope of theta
    # for theta in (pi/2 - grid.Thetas)
    #     slope = tan(theta)
    #     ax[:plot](rr, slope .* rr, "k", lw=0.1)
    # end
    #
    # # Arcs
    # xs = Array(Float64, (grid.ntheta + 1))
    # ys = Array(Float64, (grid.ntheta + 1))
    #
    # for r in rr
    #     for (i, theta) in enumerate(pi/2 - grid.Thetas)
    #         xs[i] = cos(theta) * r
    #         ys[i] = sin(theta) * r
    #     end
    #     ax[:plot](xs, ys, "k", lw=0.1)
    # end

    # ax[:set_ylim](0, maximum(zz))

    plt.savefig("density.png")

end

pp = config["parameters"]
params = ["M_star", "r_c", "r_in", "r_out", "T_10m", "q_m", "T_10a", "q_a", "T_freeze", "X_freeze", "sigma_s", "gamma", "h", "delta", "logM_gas", "delta_gas", "r_cav", "ksi", "dpc", "incl", "PA", "vel", "mu_RA", "mu_DEC"]
nparam = length(params)
starting_param = Array(Float64, nparam)

for i=1:nparam
    starting_param[i] = pp[params[i]][1]
end

# Convert logM_gas to M_gas
starting_param[15] = 10^starting_param[15]

pars = Parameters(starting_param...)

grd = config["grid"]
grid = Grid(grd["nr"], grd["ntheta"], pars.r_in, pars.r_out, true)
const global rr = grid.rs ./ AU # convert to AU
const global nr = grid.nr

global nr = grid.nr
global nz = 64
global rs = grid.rs
global zq = DiskJockey.hmodel.z_q(rs[end], pars)

global zs = cat(1, [0], logspace(log10(0.1 * AU), log10(zq), nz-1))

global xx = Array(Float64, (nz, nr));
global yy = Array(Float64, (nz, nr));

for i=1:nz
    xx[i, :] = rs
end

for j=1:nr
    yy[:, j] = zs
end

# These work
# plot_vel(pars, grid)
# plot_temp(pars, grid)
# plot_height(pars, grid)
# plot_atm_height(pars, grid)
plot_surface_density(pars, grid)
plot_density_slice(pars)

# These have yet to be tested, or are for the old "lncorrection_factor" setup
# # plot_density_gradient(pars)
# plot_density_1D_unnormed(pars)
# plot_density_1D(pars)
# plot_interpolator(pars)
# plot_dens(pars, grid)

# plot_dlnrho(pars)

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
using Comonicon

"""
my command line interface.

# Arguments

- `arg`: an argument

# Options

- `-o, --option <name>`: an option that has short option.

# Flags

- `-f, --flag`: a flag that has short flag.
"""
@main function mycmd(arg; option = "Sam", flag::Bool = false)
    @show arg
    @show option
    @show flag
end


text:
__precompile__()

module ManpowerPlanningIO

    using DataFrames
    using Dates
    using ManpowerPlanning
    using Plots
    using SQLite
    # using WebIO
    using XLSX

    MP = ManpowerPlanning
    version = v"0.7.0"

    export versionMPIO
    versionMPIO() = @info string( "Running version ", version,
        " of ManpowerPlanningIO module (ManpowerPlanning v", MP.version,
        ") in Julia v", VERSION )

    MPsim = ManpowerSimulation
    CR = XLSX.CellRange
    WS = XLSX.Worksheet
    XF = XLSX.XLSXFile
    privPath = "private"
    plotly()

    include( "configSim.jl" )
    include( joinpath( privPath, "base.jl" ) )
    include( joinpath( privPath, "attrition.jl" ) )
    include( joinpath( privPath, "attribute.jl" ) )
    include( joinpath( privPath, "basenode.jl" ) )
    include( joinpath( privPath, "compoundnode.jl" ) )
    include( joinpath( privPath, "recruitment.jl" ) )
    include( joinpath( privPath, "transition.jl" ) )
    include( joinpath( privPath, "retirement.jl" ) )
    include( joinpath( privPath, "initpop.jl" ) )
    include( "reporting/excelreports.jl" )
    include( "plotting/simulationplots.jl" )
    include( "excelrequest/exceloutputrequests.jl" )

end # module ManpowerPlanningIO


text:
"""
    LobattoLegendreBasis([RealT=Float64,] polydeg::Integer)

Create a nodal Lobatto-Legendre basis for polynomials of degree `polydeg`.
"""
struct LobattoLegendreBasis{RealT<:Real, NNODES,
                            InverseVandermondeLegendre<:AbstractMatrix{RealT},
                            BoundaryMatrix<:AbstractMatrix{RealT},
                            DerivativeMatrix<:AbstractMatrix{RealT}} <: AbstractBasisSBP{RealT}
  nodes          ::SVector{NNODES, RealT}
  weights        ::SVector{NNODES, RealT}
  inverse_weights::SVector{NNODES, RealT}

  inverse_vandermonde_legendre::InverseVandermondeLegendre
  boundary_interpolation      ::BoundaryMatrix # lhat

  derivative_matrix         ::DerivativeMatrix # dsplit
  derivative_split          ::DerivativeMatrix # dsplit
  derivative_split_transpose::DerivativeMatrix # dsplit_transposed
  derivative_dhat           ::DerivativeMatrix # dhat, neg. adjoint wrt the SBP dot product
end

function LobattoLegendreBasis(RealT, polydeg::Integer)
  nnodes_ = polydeg + 1
  nodes, weights = gauss_lobatto_nodes_weights(nnodes_)
  inverse_weights = inv.(weights)

  _, inverse_vandermonde_legendre = vandermonde_legendre(nodes)

  boundary_interpolation = zeros(nnodes_, 2)
  boundary_interpolation[:, 1] = calc_lhat(-1.0, nodes, weights)
  boundary_interpolation[:, 2] = calc_lhat( 1.0, nodes, weights)

  derivative_matrix          = polynomial_derivative_matrix(nodes)
  derivative_split           = calc_dsplit(nodes, weights)
  derivative_split_transpose = Matrix(derivative_split')
  derivative_dhat            = calc_dhat(nodes, weights)

  # type conversions to make use of StaticArrays etc.
  nodes           = SVector{nnodes_}(convert.(RealT, nodes))
  weights         = SVector{nnodes_}(convert.(RealT, weights))
  inverse_weights = SVector{nnodes_}(convert.(RealT, inverse_weights))

  inverse_vandermonde_legendre = convert.(RealT, inverse_vandermonde_legendre)
  boundary_interpolation       = SMatrix{nnodes_, 2}(convert.(RealT, boundary_interpolation))

  derivative_matrix          = SMatrix{nnodes_, nnodes_}(convert.(RealT, derivative_matrix))
  derivative_split           = SMatrix{nnodes_, nnodes_}(convert.(RealT, derivative_split))
  derivative_split_transpose = SMatrix{nnodes_, nnodes_}(convert.(RealT, derivative_split_transpose))
  derivative_dhat            = SMatrix{nnodes_, nnodes_}(convert.(RealT, derivative_dhat))

  return LobattoLegendreBasis{RealT, nnodes_, typeof(inverse_vandermonde_legendre), typeof(boundary_interpolation), typeof(derivative_matrix)}(
    nodes, weights, inverse_weights,
    inverse_vandermonde_legendre, boundary_interpolation,
    derivative_matrix, derivative_split, derivative_split_transpose, derivative_dhat
  )
end

LobattoLegendreBasis(polydeg::Integer) = LobattoLegendreBasis(Float64, polydeg)

function Base.show(io::IO, basis::LobattoLegendreBasis)
  @nospecialize basis # reduce precompilation time

  print(io, "LobattoLegendreBasis{", real(basis), "}(polydeg=", polydeg(basis), ")")
end
function Base.show(io::IO, ::MIME"text/plain", basis::LobattoLegendreBasis)
  @nospecialize basis # reduce precompilation time

  print(io, "LobattoLegendreBasis{", real(basis), "} with polynomials of degree ", polydeg(basis))
end

@inline Base.real(basis::LobattoLegendreBasis{RealT}) where {RealT} = RealT

@inline nnodes(basis::LobattoLegendreBasis{RealT, NNODES}) where {RealT, NNODES} = NNODES

@inline polydeg(basis::LobattoLegendreBasis) = nnodes(basis) - 1

struct LobattoLegendreMortarL2{RealT<:Real, NNODES, MortarMatrix<:AbstractMatrix{RealT}} <: AbstractMortarL2{RealT}
  forward_upper::MortarMatrix
  forward_lower::MortarMatrix
  reverse_upper::MortarMatrix
  reverse_lower::MortarMatrix
end

function MortarL2(basis::LobattoLegendreBasis)
  RealT = real(basis)
  NNODES = nnodes(basis)

  forward_upper = calc_forward_upper(NNODES)
  forward_lower = calc_forward_lower(NNODES)
  reverse_upper = calc_reverse_upper(NNODES, Val(:gauss))
  reverse_lower = calc_reverse_lower(NNODES, Val(:gauss))

  # type conversions to make use of StaticArrays etc.
  forward_upper = SMatrix{NNODES, NNODES}(convert.(RealT, forward_upper))
  forward_lower = SMatrix{NNODES, NNODES}(convert.(RealT, forward_lower))
  reverse_upper = SMatrix{NNODES, NNODES}(convert.(RealT, reverse_upper))
  reverse_lower = SMatrix{NNODES, NNODES}(convert.(RealT, reverse_lower))

  LobattoLegendreMortarL2{RealT, NNODES, typeof(forward_upper)}(
    forward_upper, forward_lower,
    reverse_upper, reverse_lower)
end

function Base.show(io::IO, mortar::LobattoLegendreMortarL2)
  @nospecialize mortar # reduce precompilation time

  print(io, "LobattoLegendreMortarL2{", real(mortar), "}(polydeg=", polydeg(mortar), ")")
end
function Base.show(io::IO, ::MIME"text/plain", mortar::LobattoLegendreMortarL2)
  @nospecialize mortar # reduce precompilation time

  print(io, "LobattoLegendreMortarL2{", real(mortar), "} with polynomials of degree ", polydeg(mortar))
end

@inline Base.real(mortar::LobattoLegendreMortarL2{RealT}) where {RealT} = RealT

@inline nnodes(mortar::LobattoLegendreMortarL2{RealT, NNODES}) where {RealT, NNODES} = NNODES

@inline polydeg(mortar::LobattoLegendreMortarL2) = nnodes(mortar) - 1

# TODO: We can create EC mortars along the lines of the following implementation.
# abstract type AbstractMortarEC{RealT} <: AbstractMortar{RealT} end

# struct LobattoLegendreMortarEC{RealT<:Real, NNODES, MortarMatrix<:AbstractMatrix{RealT}, SurfaceFlux} <: AbstractMortarEC{RealT}
#   forward_upper::MortarMatrix
#   forward_lower::MortarMatrix
#   reverse_upper::MortarMatrix
#   reverse_lower::MortarMatrix
#   surface_flux::SurfaceFlux
# end

# function MortarEC(basis::LobattoLegendreBasis{RealT}, surface_flux)
#   forward_upper   = calc_forward_upper(n_nodes)
#   forward_lower   = calc_forward_lower(n_nodes)
#   l2reverse_upper = calc_reverse_upper(n_nodes, Val(:gauss_lobatto))
#   l2reverse_lower = calc_reverse_lower(n_nodes, Val(:gauss_lobatto))

#   # type conversions to make use of StaticArrays etc.
#   nnodes_ = nnodes(basis)
#   forward_upper   = SMatrix{nnodes_, nnodes_}(forward_upper)
#   forward_lower   = SMatrix{nnodes_, nnodes_}(forward_lower)
#   l2reverse_upper = SMatrix{nnodes_, nnodes_}(l2reverse_upper)
#   l2reverse_lower = SMatrix{nnodes_, nnodes_}(l2reverse_lower)

#   LobattoLegendreMortarEC{RealT, nnodes_, typeof(forward_upper), typeof(surface_flux)}(
#     forward_upper, forward_lower,
#     l2reverse_upper, l2reverse_lower,
#     surface_flux
#   )
# end

# @inline nnodes(mortar::LobattoLegendreMortarEC{RealT, NNODES}) = NNODES

struct LobattoLegendreAnalyzer{RealT<:Real, NNODES, Vandermonde<:AbstractMatrix{RealT}} <: SolutionAnalyzer{RealT}
  nodes  ::SVector{NNODES, RealT}
  weights::SVector{NNODES, RealT}
  vandermonde::Vandermonde
end

function SolutionAnalyzer(basis::LobattoLegendreBasis{RealT}; analysis_polydeg=2*polydeg(basis)) where {RealT}
  nnodes_ = analysis_polydeg + 1
  nodes, weights = gauss_lobatto_nodes_weights(nnodes_)

  vandermonde = polynomial_interpolation_matrix(basis.nodes, nodes)

  # type conversions to make use of StaticArrays etc.
  nodes   = SVector{nnodes_}(convert.(RealT, nodes))
  weights = SVector{nnodes_}(convert.(RealT, weights))

  vandermonde = convert.(RealT, vandermonde)

  return LobattoLegendreAnalyzer{RealT, nnodes_, typeof(vandermonde)}(
    nodes, weights, vandermonde)
end

function Base.show(io::IO, analyzer::LobattoLegendreAnalyzer)
  @nospecialize analyzer # reduce precompilation time

  print(io, "LobattoLegendreAnalyzer{", real(analyzer), "}(polydeg=", polydeg(analyzer), ")")
end
function Base.show(io::IO, ::MIME"text/plain", analyzer::LobattoLegendreAnalyzer)
  @nospecialize analyzer # reduce precompilation time

  print(io, "LobattoLegendreAnalyzer{", real(analyzer), "} with polynomials of degree ", polydeg(analyzer))
end

@inline Base.real(analyzer::LobattoLegendreAnalyzer{RealT}) where {RealT} = RealT

@inline nnodes(analyzer::LobattoLegendreAnalyzer{RealT, NNODES}) where {RealT, NNODES} = NNODES
@inline eachnode(analyzer::LobattoLegendreAnalyzer) = Base.OneTo(nnodes(analyzer))

@inline polydeg(analyzer::LobattoLegendreAnalyzer) = nnodes(analyzer) - 1

struct LobattoLegendreAdaptorL2{RealT<:Real, NNODES, MortarMatrix<:AbstractMatrix{RealT}} <: AdaptorL2{RealT}
  forward_upper::MortarMatrix
  forward_lower::MortarMatrix
  reverse_upper::MortarMatrix
  reverse_lower::MortarMatrix
end

function AdaptorL2(basis::LobattoLegendreBasis{RealT}) where {RealT}
  nnodes_ = nnodes(basis)
  forward_upper   = calc_forward_upper(nnodes_)
  forward_lower   = calc_forward_lower(nnodes_)
  l2reverse_upper = calc_reverse_upper(nnodes_, Val(:gauss))
  l2reverse_lower = calc_reverse_lower(nnodes_, Val(:gauss))

  # type conversions to make use of StaticArrays etc.
  forward_upper   = SMatrix{nnodes_, nnodes_}(convert.(RealT, forward_upper))
  forward_lower   = SMatrix{nnodes_, nnodes_}(convert.(RealT, forward_lower))
  l2reverse_upper = SMatrix{nnodes_, nnodes_}(convert.(RealT, l2reverse_upper))
  l2reverse_lower = SMatrix{nnodes_, nnodes_}(convert.(RealT, l2reverse_lower))

  LobattoLegendreAdaptorL2{RealT, nnodes_, typeof(forward_upper)}(
    forward_upper, forward_lower,
    l2reverse_upper, l2reverse_lower)
end

function Base.show(io::IO, adaptor::LobattoLegendreAdaptorL2)
  @nospecialize adaptor # reduce precompilation time

  print(io, "LobattoLegendreAdaptorL2{", real(adaptor), "}(polydeg=", polydeg(adaptor), ")")
end
function Base.show(io::IO, ::MIME"text/plain", adaptor::LobattoLegendreAdaptorL2)
  @nospecialize adaptor # reduce precompilation time

  print(io, "LobattoLegendreAdaptorL2{", real(adaptor), "} with polynomials of degree ", polydeg(adaptor))
end

@inline Base.real(adaptor::LobattoLegendreAdaptorL2{RealT}) where {RealT} = RealT

@inline nnodes(adaptor::LobattoLegendreAdaptorL2{RealT, NNODES}) where {RealT, NNODES} = NNODES

@inline polydeg(adaptor::LobattoLegendreAdaptorL2) = nnodes(adaptor) - 1

###############################################################################
# Polynomial derivative and interpolation functions

# TODO: Taal refactor, allow other RealT below and adapt constructors above accordingly

# Calculate the Dhat matrix
function calc_dhat(nodes, weights)
  n_nodes = length(nodes)
  dhat = Matrix(polynomial_derivative_matrix(nodes)')

  for n in 1:n_nodes, j in 1:n_nodes
    dhat[j, n] *= -weights[n] / weights[j]
  end

  return dhat
end

# Calculate the Dsplit matrix for split-form differentiation: dplit = 2D - M⁻¹B
function calc_dsplit(nodes, weights)
  # Start with 2 x the normal D matrix
  dsplit = 2 .* polynomial_derivative_matrix(nodes)

  # Modify to account for
  dsplit[  1,   1] += 1 / weights[1]
  dsplit[end, end] -= 1 / weights[end]

  return dsplit
end

# Calculate the polynomial derivative matrix D
function polynomial_derivative_matrix(nodes)
  n_nodes = length(nodes)
  d = zeros(n_nodes, n_nodes)
  wbary = barycentric_weights(nodes)

  for i in 1:n_nodes, j in 1:n_nodes
    if j != i
      d[i, j] = wbary[j] / wbary[i] * 1 / (nodes[i] - nodes[j])
      d[i, i] -= d[i, j]
    end
  end

  return d
end

# Calculate and interpolation matrix (Vandermonde matrix) between two given sets of nodes
function polynomial_interpolation_matrix(nodes_in, nodes_out)
  n_nodes_in = length(nodes_in)
  n_nodes_out = length(nodes_out)
  wbary_in = barycentric_weights(nodes_in)
  vdm = zeros(n_nodes_out, n_nodes_in)

  for k in 1:n_nodes_out
    match = false
    for j in 1:n_nodes_in
      if isapprox(nodes_out[k], nodes_in[j], rtol=eps())
        match = true
        vdm[k, j] = 1
      end
    end

    if match == false
      s = 0.0
      for j in 1:n_nodes_in
        t = wbary_in[j] / (nodes_out[k] - nodes_in[j])
        vdm[k, j] = t
        s += t
      end
      for j in 1:n_nodes_in
        vdm[k, j] = vdm[k, j] / s
      end
    end
  end

  return vdm
end

# Calculate the barycentric weights for a given node distribution.
function barycentric_weights(nodes)
  n_nodes = length(nodes)
  weights = ones(n_nodes)

  for j = 2:n_nodes, k = 1:(j-1)
    weights[k] *= nodes[k] - nodes[j]
    weights[j] *= nodes[j] - nodes[k]
  end

  for j in 1:n_nodes
    weights[j] = 1 / weights[j]
  end

  return weights
end

# Calculate Lhat.
function calc_lhat(x, nodes, weights)
  n_nodes = length(nodes)
  wbary = barycentric_weights(nodes)

  lhat = lagrange_interpolating_polynomials(x, nodes, wbary)

  for i in 1:n_nodes
    lhat[i] /= weights[i]
  end

  return lhat
end

# Calculate Lagrange polynomials for a given node distribution.
function lagrange_interpolating_polynomials(x, nodes, wbary)
  n_nodes = length(nodes)
  polynomials = zeros(n_nodes)

  for i in 1:n_nodes
    if isapprox(x, nodes[i], rtol=eps(x))
      polynomials[i] = 1
      return polynomials
    end
  end

  for i in 1:n_nodes
    polynomials[i] = wbary[i] / (x - nodes[i])
  end
  total = sum(polynomials)

  for i in 1:n_nodes
    polynomials[i] /= total
  end

  return polynomials
end

# From FLUXO (but really from blue book by Kopriva)
function gauss_lobatto_nodes_weights(n_nodes::Integer)
  # From Kopriva's book
  n_iterations = 10
  tolerance = 1e-15

  # Initialize output
  nodes = zeros(n_nodes)
  weights = zeros(n_nodes)

  # Get polynomial degree for convenience
  N = n_nodes - 1

  # Calculate values at boundary
  nodes[1] = -1.0
  nodes[end] = 1.0
  weights[1] = 2 / (N * (N + 1))
  weights[end] = weights[1]

  # Calculate interior values
  if N > 1
    cont1 = pi/N
    cont2 = 3/(8 * N * pi)

    # Use symmetry -> only left side is computed
    for i in 1:(div(N + 1, 2) - 1)
      # Calculate node
      # Initial guess for Newton method
      nodes[i+1] = -cos(cont1*(i+0.25) - cont2/(i+0.25))

      # Newton iteration to find root of Legendre polynomial (= integration node)
      for k in 0:n_iterations
        q, qder, _ = calc_q_and_l(N, nodes[i+1])
        dx = -q/qder
        nodes[i+1] += dx
        if abs(dx) < tolerance * abs(nodes[i+1])
          break
        end
      end

      # Calculate weight
      _, _, L = calc_q_and_l(N, nodes[i+1])
      weights[i+1] = weights[1] / L^2

      # Set nodes and weights according to symmetry properties
      nodes[N+1-i] = -nodes[i+1]
      weights[N+1-i] = weights[i+1]
    end
  end

  # If odd number of nodes, set center node to origin (= 0.0) and calculate weight
  if n_nodes % 2 == 1
    _, _, L = calc_q_and_l(N, 0)
    nodes[div(N, 2) + 1] = 0.0
    weights[div(N, 2) + 1] = weights[1] / L^2
  end

  return nodes, weights
end

# From FLUXO (but really from blue book by Kopriva)
function calc_q_and_l(N::Integer, x::Float64)
  L_Nm2 = 1.0
  L_Nm1 = x
  Lder_Nm2 = 0.0
  Lder_Nm1 = 1.0

  local L
  for i in 2:N
    L = ((2 * i - 1) * x * L_Nm1 - (i - 1) * L_Nm2) / i
    Lder = Lder_Nm2 + (2 * i - 1) * L_Nm1
    L_Nm2 = L_Nm1
    L_Nm1 = L
    Lder_Nm2 = Lder_Nm1
    Lder_Nm1 = Lder
  end

  q = (2 * N + 1)/(N + 1) * (x * L - L_Nm2)
  qder = (2 * N + 1) * L

  return q, qder, L
end
calc_q_and_l(N::Integer, x::Real) = calc_q_and_l(N, convert(Float64, x))

# From FLUXO (but really from blue book by Kopriva)
function gauss_nodes_weights(n_nodes::Integer)
  # From Kopriva's book
  n_iterations = 10
  tolerance = 1e-15

  # Initialize output
  nodes = ones(n_nodes) * 1000
  weights = zeros(n_nodes)

  # Get polynomial degree for convenience
  N = n_nodes - 1
  if N == 0
    nodes .= 0.0
    weights .= 2.0
    return nodes, weights
  elseif N == 1
    nodes[1] = -sqrt(1/3)
    nodes[end] = -nodes[1]
    weights .= 1.0
    return nodes, weights
  else # N > 1
    # Use symmetry property of the roots of the Legendre polynomials
    for i in 0:(div(N + 1, 2) - 1)
      # Starting guess for Newton method
      nodes[i+1] = -cos(pi / (2 * N + 2) * (2 * i + 1))

      # Newton iteration to find root of Legendre polynomial (= integration node)
      for k in 0:n_iterations
        poly, deriv = legendre_polynomial_and_derivative(N + 1, nodes[i+1])
        dx = -poly / deriv
        nodes[i+1] += dx
        if abs(dx) < tolerance * abs(nodes[i+1])
          break
        end
      end

      # Calculate weight
      poly, deriv = legendre_polynomial_and_derivative(N + 1, nodes[i+1])
      weights[i+1] = (2 * N + 3) / ((1 - nodes[i+1]^2) * deriv^2)

      # Set nodes and weights according to symmetry properties
      nodes[N+1-i] = -nodes[i+1]
      weights[N+1-i] = weights[i+1]
    end

    # If odd number of nodes, set center node to origin (= 0.0) and calculate weight
    if n_nodes % 2 == 1
      poly, deriv = legendre_polynomial_and_derivative(N + 1, 0.0)
      nodes[div(N, 2) + 1] = 0.0
      weights[div(N, 2) + 1] = (2 * N + 3) / deriv^2
    end

    return nodes, weights
  end
end

# From FLUXO (but really from blue book by Kopriva)
function legendre_polynomial_and_derivative(N::Int, x::Real)
  if N == 0
    poly = 1.0
    deriv = 0.0
  elseif N == 1
    poly = convert(Float64, x)
    deriv = 1.0
  else
    poly_Nm2 = 1.0
    poly_Nm1 = convert(Float64, x)
    deriv_Nm2 = 0.0
    deriv_Nm1 = 1.0

    poly = 0.0
    deriv = 0.0
    for i in 2:N
      poly = ((2*i-1) * x * poly_Nm1 - (i-1) * poly_Nm2) / i
      deriv=deriv_Nm2 + (2*i-1)*poly_Nm1
      poly_Nm2=poly_Nm1
      poly_Nm1=poly
      deriv_Nm2=deriv_Nm1
      deriv_Nm1=deriv
    end
  end

  # Normalize
  poly = poly * sqrt(N+0.5)
  deriv = deriv * sqrt(N+0.5)

  return poly, deriv
end

# Calculate Legendre vandermonde matrix and its inverse
function vandermonde_legendre(nodes, N)
  n_nodes = length(nodes)
  n_modes = N + 1
  vandermonde = zeros(n_nodes, n_modes)

  for i in 1:n_nodes
    for m in 1:n_modes
      vandermonde[i, m], _ = legendre_polynomial_and_derivative(m-1, nodes[i])
    end
  end
  # for very high polynomial degree, this is not well conditioned
  inverse_vandermonde = inv(vandermonde)
  return vandermonde, inverse_vandermonde
end
vandermonde_legendre(nodes) = vandermonde_legendre(nodes, length(nodes) - 1)


text:
import Q: chkparens
@testset "parser" for (n, cmd) in [
   (0, "([{}])"),
   (-1, "([{}]"),
   (5,  "((([}])))"),
   (2, "(})"),
   (0, ""),
   (0, "q)42"),
]
   @test n == chkparens(cmd)
end


text:
# NOTE: This code is meant to be used to generate sequential rules
# from the frequent sequences generated by the SPADE algorithm.

# ==(x::SequenceRule, y::SequenceRule) = x.rule == y.rule && x.conf == y.conf
#
# function unique(v::Array{SequenceRule, 1})
#     out = Array{SequenceRule, 1}(0)
#     for i = 1:length(v)
#         if !in(v[i], out)
#             push!(out, v[i])
#         end
#     end
#     return out
# end

==(x::SeqRule, y::SeqRule) = x.prefix == y.prefix && x.postfix == y.postfix && x.conf == y.conf

function unique(v::Array{SeqRule, 1})
    out = Array{SeqRule, 1}(0)
    for i = 1:length(v)
        if !in(v[i], out)
            push!(out, v[i])
        end
    end
    return out
end

# NOTE: The gen_rules1() function is quite inefficient. In
# particular, it's runtime is O(n*m*2^k) where k is the number
# of elements in the pattern, m is the number of patterns for a
# given pattern length k, and n is the number of different pattern
# lengths. Additionally, though this corresponds to the pseudocode
# Zaki provides in his orginal paper, it is not the implementation
# that appears in the arulesSequence package.

function gen_rules1(F::Array{Array{IDList, 1}, 1}, min_conf)
    supp_cnt = count_patterns(F)
    rules = String[]

    # Check the confidence for all sub-patterns
    # from all of our frequent patterns
    for k = 1:length(F)
        for i = 1:length(F[k])
            sub_patrns = gen_combin_subpatterns(F[k][i].patrn)
            for s in sub_patrns
                if s ≠ ""
                    cnt = get(supp_cnt, s, 0)
                    conf = isfinite(cnt) ? F[k][i].supp_cnt/cnt : -Inf

                    if conf ≥ min_conf
                        push!(rules, "$s => $(pattern_string(F[k][i].patrn))")
                    end
                end
            end
        end
    end
    rules
end

# rules = gen_rules1(res, 0)

function create_children(node::PrefixNode, uniq_items::Array{String,1}, supp_cnt)
    seq_ext_children = Array{PrefixNode,1}(0)
    item_ext_children = Array{PrefixNode,1}(0)

    for item in uniq_items
        seq_patrn = sequence_extension(node.pattern, item)

        # computing support
        seq_string = pattern_string(seq_patrn)
        seq_supp = get(supp_cnt, seq_string, 0)

        # only append children with non-zero support
        if seq_supp > 0
            seq_extd_child = PrefixNode(seq_patrn)
            seq_extd_child.support = seq_supp
            push!(seq_ext_children, seq_extd_child)
        end

        itm_patrn = item_extension(node.pattern, item)
        itm_string = pattern_string(itm_patrn)
        itm_supp = get(supp_cnt, itm_string, 0)

        # We only create an item-extended child if this item
        # doesn't already appear in the last array of our our
        # parent (i.e., `node`). Otherwise, we would either be
        # duplicating an entry in the last array, or be making
        # an exact duplicate of our parent (assuming we only allow
        # unique entries in the pattern's arrays).
        if itm_supp > 0 && item ∉ node.pattern[end]
            item_extd_child = PrefixNode(itm_patrn)
            item_extd_child.support = itm_supp
            push!(item_ext_children, item_extd_child)
        end
    end
    return (seq_ext_children, item_ext_children)
end

# xroot = PrefixNode([["A"]])
# xsupp_cnt = Dict("{A} -> {B}" => 1,
#                  "{A} -> {C}" => 1,
#                  "{A} -> {D}" => 1)
# @code_warntype create_children(xroot, ["A", "B", "C", "D"], xsupp_cnt)

function growtree!(root, uniq_items, supp_cnt, depth = 0, maxdepth = 1)
    root.seq_ext_children, root.item_ext_children = create_children(root, uniq_items, supp_cnt)
    allchildren = [root.seq_ext_children; root.item_ext_children]

    for child in allchildren
        # push!(tree, child)

        if depth+1 < maxdepth
            growtree!(child, uniq_items, supp_cnt, depth+1, maxdepth)
        end
    end
end

# xroot = PrefixNode([["A"]])
# xsupp_cnt = Dict("{A} -> {B}" => 1,
#                  "{A} -> {C}" => 1,
#                  "{A} -> {D}" => 1)
# @code_warntype growtree!(xroot, ["A", "B", "C", "D"], xsupp_cnt, 0, 3)

function build_tree(F::Array{Array{IDList,1},1}, maxdepth)
    supp_cnt = count_patterns(F)
    uniq_items = String[]

    # NOTE: This loop is an embarassment. This is all to get the unique items
    # from the frequent sequences. The same unique item vector is created by
    # the spade() function. It should be passed to here.
    for k = 1:length(F)
        for i = 1:length(F[k])
            for j = 1:length(F[k][i].patrn)
                for l = 1:length(F[k][i].patrn[j])
                    if F[k][i].patrn[j][l] ∉ uniq_items
                        push!(uniq_items, F[k][i].patrn[j][l])
                    end
                end
            end
        end
    end

    root = PrefixNode([[""]])
    root.seq_ext_children = PrefixNode[]

    for itm in uniq_items
        child_node1 = PrefixNode([[itm]])
        child_node1.support = supp_cnt[pattern_string(itm)]
        push!(root.seq_ext_children, child_node1)
        growtree!(child_node1, uniq_items, supp_cnt, 1, maxdepth)
    end
    return root
end

# @time build_tree(res3, 20);
# @code_warntype build_tree(res, 20);

function generate_sr_from_tree_root!(sp_root::PrefixNode, rules, min_conf)
    if isdefined(sp_root, :seq_ext_children)
        for seq_child in sp_root.seq_ext_children
            # println(seq_child)

            generate_sr_from_subtree!(sp_root, seq_child, rules, min_conf)
        end
    end

    if isdefined(sp_root, :item_ext_children)
        for itm_child in sp_root.item_ext_children
            generate_sr_from_tree_root!(itm_child, rules, min_conf)
        end
    end
end

function generate_sr_from_subtree!(pre, seq_child, rules, min_conf, one_elem_consq = true)
    conf = seq_child.support/pre.support

    if conf ≥ min_conf
        post = postfix(pre.pattern, seq_child.pattern)
        if one_elem_consq
            if length(post) == 1
                new_rule = SeqRule(pre.pattern, post, conf)
                if new_rule ∉ rules
                    push!(rules, new_rule)
                end
            end
        elseif !one_elem_consq
            new_rule = SeqRule(pre.pattern, post, conf)
            if new_rule ∉ rules
                push!(rules, new_rule)
            end
        end

        if isdefined(seq_child, :seq_ext_children)
            if !isempty(seq_child.seq_ext_children)
                for grandchild in seq_child.seq_ext_children

                    # adding this very speculative (and based on a guess)
                    generate_sr_from_subtree!(seq_child, grandchild, rules, min_conf)

                    generate_sr_from_subtree!(pre, grandchild, rules, min_conf)
                end
            end
        end

        if isdefined(seq_child, :item_ext_children)
            if !isempty(seq_child.item_ext_children)
                for grandchild in seq_child.item_ext_children
                    generate_sr_from_subtree!(pre, grandchild, rules, min_conf)
                end
            end
        end
    end
end

# t = build_tree(res, 20)
# xrules = []
# @code_warntype generate_sr_from_subtree!(t.seq_ext_children[1], t.seq_ext_children[1].seq_ext_children[1], xrules, 0.01)
# @code_warntype generate_sr_from_tree_root!(t.seq_ext_children[1], xrules, 0.01)

# NOTE: The sequential_rules() function currently builds a
# prefix tree and then traverses the tree generating rules.
# This could be made more efficient and cache-friendly by
# building and traversing the tree in the same pass.
"""
    sequential_rules(F; min_conf, maxdepth)

This function takes an array of arrays with `IDList` objects and
returns an array of sequential rules.

### Arguments
* `F`: an array of arrays of `IDList` generated by `spade()` function
* `min_conf`: minimum level of confidence rules must have
* `maxlength`: maximum number of elements a rule must have (e.g., `maxlength = 3`
    gives rules such as <{A},{B}> => <{C}> or <{D,E}> => <{F}>.
"""
function sequential_rules(F, min_conf = 0.1, maxlength = 5)
    root = build_tree(F, maxlength)
    rules = Array{SeqRule,1}(0)

    for i = 1:length(root.seq_ext_children)
        generate_sr_from_tree_root!(root.seq_ext_children[i], rules, min_conf)
    end
    rules
end

# @code_warntype sequential_rules(res, 0.01)
# @time sequential_rules(res2, 0.01, 20);


text:
## Configuration
n = 5
adj = Float64[0 0 1 1 0;
              0 0 1 0 0;
              0 1 0 0 0;
              1 0 1 0 1;
              0 1 1 0 0]
adj ./= sum(adj,2)

## Graph algorithm
pageRank = ones(Float64,n)
N = 30
for t = 1:N
  temp = zeros(pageRank)
  for j = 1:n
    temp[j] = 0.15 + 0.85 * vecdot(pageRank,adj[:,j])
  end
  pageRank = temp
end
pageRank

## Markov Chain algorithm
comAdj = [0.85*adj 0.15*ones(n); ones(n)'/n 0.0]
comPageRank = eigvecs(comAdj')[:,1]
comPageRank = comPageRank / sum(comPageRank) * 1.15n

max_stars_repo_name: WenjieZ/PageRank

max_issues_repo_name: WenjieZ/PageRank

max_forks_repo_name: WenjieZ/PageRank