text:
"""Entropy MCTS solver with IPFT

Derived from MCTS-DPW (JuliaPOMDP/MCTS.jl/src/dpw_types.jl)

Fields:

    depth::Int64:
        Maximum rollout horizon and tree depth.
        default: 10

    exploration_constant::Float64:
        Specified how much the solver should explore.
        In the UCB equation, Q + c*sqrt(log(t/N)), c is the exploration constant.
        default: 1.0

    n_iterations::Int64
        Number of iterations during each action() call.
        default: 100

    max_time::Float64
        Maximum amount of CPU time spent iterating through simulations.
        default: Inf

    k_action::Float64
    alpha_action::Float64
    k_observation::Float64
    alpha_observation::Float64
        These constants control the double progressive widening. A new state
        or action will be added if the number of children is less than or equal to kN^alpha.
        defaults: k:10, alpha:0.5

    lambda::Float64
        Hyper parameter balancing reward and information reward

    keep_tree::Bool
        If true, store the tree in the planner for reuse at the next timestep (and every time it is used in the future). There is a computational cost for maintaining the state dictionary necessary for this.
        default: false

    reuse_states::Bool
        If true, reuses first sample in a tree node instead of sampling new particles in every visit

    enable_action_pw::Bool
        If true, enable progressive widening on the action space; if false just use the whole action space.
        default: true

    check_repeat_state::Bool
    check_repeat_action::Bool
        When constructing the tree, check whether a state or action has been seen before (there is a computational cost to maintaining the dictionaries necessary for this)
        default: true

    tree_in_info::Bool:
        If true, return the tree in the info dict when action_info is called. False by default because it can use a lot of memory if histories are being saved.
        default: false

    rng::AbstractRNG:
        Random number generator

    estimate_value::Any (rollout policy)
        Function, object, or number used to estimate the value at the leaf nodes.
        If this is a function `f`, `f(mdp, s, depth)` will be called to estimate the value.
        If this is an object `o`, `estimate_value(o, mdp, s, depth)` will be called.
        If this is a number, the value will be set to that number.
        default: RolloutEstimator(RandomSolver(rng))

    init_Q::Any
        Function, object, or number used to set the initial Q(s,a) value at a new node.
        If this is a function `f`, `f(mdp, s, a)` will be called to set the value.
        If this is an object `o`, `init_Q(o, mdp, s, a)` will be called.
        If this is a number, Q will always be set to that number.
        default: 0.0

    init_N::Any
        Function, object, or number used to set the initial N(s,a) value at a new node.
        If this is a function `f`, `f(mdp, s, a)` will be called to set the value.
        If this is an object `o`, `init_N(o, mdp, s, a)` will be called.
        If this is a number, N will always be set to that number.
        default: 0

    next_action::Any
        Function or object used to choose the next action to be considered for progressive widening.
        The next action is determined based on the MDP, the state, `s`, and the current `IPFTStateNode`, `snode`.
        If this is a function `f`, `f(mdp, s, snode)` will be called to set the value.
        If this is an object `o`, `next_action(o, mdp, s, snode)` will be called.
        default: RandomActionGenerator(rng)

    default_action::Any
        Function, action, or Policy used to determine the action if POMCP fails with exception `ex`.
        If this is a Function `f`, `f(pomdp, belief, ex)` will be called.
        If this is a Policy `p`, `action(p, belief)` will be called.
        If it is an object `a`, `default_action(a, pomdp, belief, ex)` will be called, and if this method is not implemented, `a` will be returned directly.
        default: `ExceptionRethrow()`
"""
mutable struct IPFTSolver <: AbstractMCTSSolver
    depth::Int
    criterion::AbstractMultiObjectiveCriterion
    final_criterion::AbstractMultiObjectiveCriterion
    # exploration_constant::Float64
    n_iterations::Int
    max_time::Float64
    k_action::Float64
    alpha_action::Float64
    k_observation::Float64
    alpha_observation::Float64
    keep_tree::Bool
    reuse_states::Bool
    enable_action_pw::Bool
    check_repeat_state::Bool
    check_repeat_action::Bool
    tree_in_info::Bool
    rng::AbstractRNG
    estimate_value::Any
    init_Q::Any
    init_N::Any
    next_action::Any
    default_action::Any
end

"""
    IPFTSolver()

Use keyword arguments to specify values for the fields
"""
function IPFTSolver(;depth::Int=10,
                    criterion::AbstractMultiObjectiveCriterion=MultiObjectiveUCB([1.0,0.0],1.0),
                    final_criterion::AbstractMultiObjectiveCriterion=MultiObjective(criterion),
                    # exploration_constant::Float64=1.0,
                    n_iterations::Int=100,
                    max_time::Float64=Inf,
                    k_action::Float64=10.0,
                    alpha_action::Float64=0.5,
                    k_observation::Float64=10.0,
                    alpha_observation::Float64=0.5,
                    keep_tree::Bool=false,
                    reuse_states::Bool=false,
                    enable_action_pw::Bool=true,
                    check_repeat_state::Bool=true,
                    check_repeat_action::Bool=true,
                    tree_in_info::Bool=false,
                    rng::AbstractRNG=Random.GLOBAL_RNG,
                    estimate_value::Any = RolloutEstimator(RandomSolver(rng)),
                    # init_Q::Any = 0.0,
                    init_Q::Any = [0.0, 0.0],
                    init_N::Any = 0,
                    next_action::Any = RandomActionGenerator(rng),
                    default_action::Any = ExceptionRethrow()
                   )
    IPFTSolver(depth, criterion, final_criterion, n_iterations, max_time, k_action, alpha_action, k_observation, alpha_observation, keep_tree, reuse_states, enable_action_pw, check_repeat_state, check_repeat_action, tree_in_info, rng, estimate_value, init_Q, init_N, next_action, default_action)
end

mutable struct IPFTree{S,A}
    # for each state node
    total_n::Vector{Int}
    children::Vector{Vector{Int}}
    s_labels::Vector{S}
    s_lookup::Dict{S, Int}

    # for each state-action node
    n::Vector{Int}
    q::Vector{Vector{Float64}}
    transitions::Vector{Vector{Tuple{Int,S,Array{Float64,1}}}}
    a_labels::Vector{A}
    a_lookup::Dict{Tuple{Int,A}, Int}

    # for tracking transitions
    n_a_children::Vector{Int}
    unique_transitions::Set{Tuple{Int,Int}}

    function IPFTree{S,A}(sz::Int=1000) where {S,A}
        sz = min(sz, 100_000)
        return new(sizehint!(Int[], sz),
                   sizehint!(Vector{Int}[], sz),
                   sizehint!(S[], sz),
                   Dict{S, Int}(),

                   sizehint!(Int[], sz),
                   sizehint!(Float64[], sz),
                   sizehint!(Vector{Tuple{Int,Float64}}[], sz),
                   sizehint!(A[], sz),
                   Dict{Tuple{Int,A}, Int}(),

                   sizehint!(Int[], sz),
                   Set{Tuple{Int,Int}}()
                  )
    end
end

function insert_state_node!(tree::IPFTree{S,A}, s::S, maintain_s_lookup=true) where {S,A}
    push!(tree.total_n, 0)
    push!(tree.children, Int[])
    push!(tree.s_labels, s)
    snode = length(tree.total_n)
    if maintain_s_lookup
        tree.s_lookup[s] = snode
    end
    return snode
end

function insert_action_node!(tree::IPFTree{S,A}, snode::Int, a::A, n0::Int, q0::Vector{Float64}, maintain_a_lookup=true) where {S,A}
    push!(tree.n, n0)
    push!(tree.q, q0)
    push!(tree.a_labels, a)
    push!(tree.transitions, Vector{Tuple{Int,Float64}}[])
    sanode = length(tree.n)
    push!(tree.children[snode], sanode)
    push!(tree.n_a_children, 0)
    if maintain_a_lookup
        tree.a_lookup[(snode, a)] = sanode
    end
    return sanode
end

Base.isempty(tree::IPFTree) = isempty(tree.n) && isempty(tree.q)

struct IPFTStateNode{S,A} <: AbstractStateNode
    tree::IPFTree{S,A}
    index::Int
end

total_n(n::IPFTStateNode) = n.tree.total_n[n.index]
children(n::IPFTStateNode) = n.tree.children[n.index]
n_children(n::IPFTStateNode) = length(children(n))
isroot(n::IPFTStateNode) = n.index == 1

mutable struct IPFTPlanner{P<:Union{MDP,POMDP}, S, A, SE, NA, RNG} <: InformationMCTSPlanner{P}
    solver::IPFTSolver
    mdp::P
    tree::Union{Nothing, IPFTree{S,A}}
    solved_estimate::SE
    next_action::NA
    rng::RNG
end

function IPFTPlanner(solver::IPFTSolver, mdp::P) where P<:Union{MDP,POMDP}
    se = convert_estimator(solver.estimate_value, solver, mdp)
    return IPFTPlanner{P,
                      statetype(P),
                      actiontype(P),
                      typeof(se),
                      typeof(solver.next_action),
                      typeof(solver.rng)}(solver,
                                            mdp,
                                            nothing,
                                            se,
                                            solver.next_action,
                                            solver.rng
    )
end

Random.seed!(p::IPFTPlanner, seed) = Random.seed!(p.rng, seed)


text:
# # optimization spaces
# abstract type Space end

# # Euclidean
# struct Euclidean <: Space
#     n::Int
# end

function conjugate(q)
	s = q[1]
	v = q[2:4]

	return [s; -v]
end

function L_multiply(q)
	s = q[1]
	v = q[2:4]

	SMatrix{4,4}([s -transpose(v);
	              v s * I + skew(v)])
end

function R_multiply(q)
	s = q[1]
	v = q[2:4]

	SMatrix{4,4}([s -transpose(v);
	              v s * I - skew(v)])
end

function multiply(q1, q2)
	L_multiply(q1) * q2
end

# eq. 14 http://roboticexplorationlab.org/papers/planning_with_attitude.pdf
function attitude_jacobian(q)
	s = q[1]
	v = q[2:4]

	[-transpose(v);
	 s * I + skew(v)]
end

function G_func(q)
	quat = q[4:7]
	[1.0 0.0 0.0 0.0 0.0 0.0;
     0.0 1.0 0.0 0.0 0.0 0.0;
	 0.0 0.0 1.0 0.0 0.0 0.0;
     zeros(4, 3) attitude_jacobian(quat)]
end

function static_quaternion(q::UnitQuaternion)
	SVector{4}([q.w, q.x, q.y, q.z])
end

function quaternion_offset(q_parent, q_child)
	L_multiply(conjugate(q_parent)) * q_child
end

function get_quaternion(q)
	return view(q, 4:7)
end

# eq. 16 http://roboticexplorationlab.org/papers/maximal_coordinate_dynamics.pdf
function ω_finite_difference(q1, q2, h)
	2.0 * multiply(conjugate(q1), (q2 - q1) ./ h)[2:4]
end

# https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation#Quaternion-derived_rotation_matrix
function quaternion_rotation_matrix(q)
	r, i, j, k  = q

	r11 = 1.0 - 2.0 * (j^2.0 + k^2.0)
	r12 = 2.0 * (i * j - k * r)
	r13 = 2.0 * (i * k + j * r)

	r21 = 2.0 * (i * j + k * r)
	r22 = 1.0 - 2.0 * (i^2.0 + k^2.0)
	r23 = 2.0 * (j * k - i * r)

	r31 = 2.0 * (i * k - j * r)
	r32 = 2.0 * (j * k + i * r)
	r33 = 1.0 - 2.0 * (i^2.0 + j^2.0)

	SMatrix{3,3}([r11 r12 r13;
	              r21 r22 r23;
				  r31 r32 r33])
end

# Cayley map
function cayley_map(ϕ)
	1.0 / sqrt(1.0 + norm(ϕ)^2.0) * [1.0; ϕ]
end

# slerp
function slerp(q0::UnitQuaternion, q1::UnitQuaternion, t::T) where T
	Rotations.params(q0 * exp(t * log(inv(q0) * q1)))
end

# Square root of quaternion: https://www.johndcook.com/blog/2021/01/06/quaternion-square-roots/
function sqrt_quat(q; ϵ=1e-16)
	r = norm(q)
	# Angle
	theta = acos(q[1]/r)
	# Axis
	u = q[2:4]
	u ./= norm(u) + ϵ
	# Half axis-angle rotation
	x = [cos(theta/2);  sin(theta/2)*u]
	# Sqrt on the norm of the quaternion
	x .*= r^0.5 # useless for unit quaternion
	return x
end

function midpoint(q0, q1)
	# Small delta rotation
	qψ = R_multiply(q0)' * q1
	# Divide delta in two
	sqψ = sqrt_quat(qψ)
	# Apply small rotation
	qmid = L_multiply(sqψ) * q0
	return qmid
end

struct RnQuaternion <: Space
	n::Int # size of the optimization space (typically = num_var - 1)
	mapping # Jacobian mapping i.e. a function of z
	r_idx # euclidean part in z
	Δr_idx # euclidean part in Δ
	quat_idx # quaternion part in z
	Δquat_idx # quaternion part in Δ
end

function rn_quaternion_space(dim, mapping, r_idx, Δr_idx, quat_idx, Δquat_idx)
	RnQuaternion(dim, mapping, r_idx, Δr_idx, quat_idx, Δquat_idx)
end

function candidate_point!(z̄::Vector{T}, s::RnQuaternion, z::Vector{T}, Δ::Vector{T}, α::T) where T
	z̄[s.r_idx] .= z[s.r_idx] - α .* Δ[s.Δr_idx]
	for i = 1:length(s.quat_idx)
		z̄[s.quat_idx[i]] .= L_multiply(z[s.quat_idx[i]]) * cayley_map(-1.0 * α .* Δ[s.Δquat_idx[i]])
	end
	return nothing
end

function mapping!(δz, s::RnQuaternion, δzs, z)
	δz .= s.mapping(z) * δzs
end

ext: jl


text:
"""
```
load_data(m::AbstractModel; try_disk::Bool = true, verbose::Symbol = :low)
```

Create a DataFrame with all data series for this model, fully transformed.

First, check the disk to see if a valid dataset is already stored in `inpath(m, "data")`. A
dataset is valid if every series in `m.data_transforms` is present and the entire sample is
contained (from `date_presample_start` to `date_mainsample_end`. If no valid dataset is
already stored, the dataset will be recreated. This check can be eliminated by passing
`try_disk=false`.

If the dataset is to be recreated, in a preliminary stage, intermediate data series, as
specified in `m.data_series`, are loaded in levels using `load_data_levels`. See
`?load_data_levels` for more details.

Then, the series in levels are transformed as specified in `m.data_transforms`. See
`?transform_data` for more details.

The resulting DataFrame is saved to disk as `data_<yymmdd>.csv` and returned to the caller.
"""
function load_data(m::AbstractModel; try_disk::Bool = true, verbose::Symbol=:low)
    recreate_data = false

    # Check if already downloaded
    if try_disk && has_saved_data(m)
        if VERBOSITY[verbose] >= VERBOSITY[:low]
            print("Reading dataset from disk...")
        end
        df = read_data(m)
        if isvalid_data(m, df)
            if VERBOSITY[verbose] >= VERBOSITY[:low]
                println("dataset from disk valid")
            end
        else
            if VERBOSITY[verbose] >= VERBOSITY[:low]
                println("dataset from disk not valid")
            end
            recreate_data = true
        end
    else
        recreate_data = true
    end

    # Download routines
    if recreate_data
        if VERBOSITY[verbose] >= VERBOSITY[:low]
            println("Creating dataset...")
        end
        df = load_data_levels(m; verbose=verbose)
        df = transform_data(m, df; verbose=verbose)

        # Ensure that only appropriate rows make it into the returned DataFrame.
        start_date = get_setting(m, :date_presample_start)
        end_date   = get_setting(m, :date_mainsample_end)
        df = df[start_date .<= df[:, :date] .<= end_date, :]

        save_data(m, df)
        if VERBOSITY[verbose] >= VERBOSITY[:low]
            println("dataset creation successful")
        end
    end

    return df
end

"""
```
load_data_levels(m::AbstractModel; verbose::Symbol=:low)
```

Load data in levels by appealing to the data sources specified for the model. Data from FRED
is loaded first, by default; then, merge other custom data sources.

Check on disk in `inpath(m, "data")` datasets, of the correct vintage, corresponding to the
ones in `keys(m.data_series)`. Load the appropriate data series (specified in
`m.data_series[source]`) for each data source.

To accomodate growth rates and other similar transformations, more rows of data may be
downloaded than otherwise specified by the date model settings. (By the end of the process,
these rows will have been dropped.)

Data from FRED (i.e. the `:fred` data source) are treated separately. These are downloaded
using `load_fred_data`. See `?load_fred_data` for more details.

Data from non-FRED data sources are read from disk, verified, and merged.
"""
function load_data_levels(m::AbstractModel; verbose::Symbol=:low)
    # Start two quarters further back than `start_date` as we need these additional
    # quarters to compute differences.
    start_date = get_setting(m, :date_presample_start) - Dates.Month(6)
    end_date = get_setting(m, :date_mainsample_end)

    # Load FRED data
    df = load_fred_data(m; start_date=firstdayofquarter(start_date), end_date=end_date)

    # Set ois series to load
    if n_anticipated_shocks(m) > 0
        m.data_series[:ois] = [symbol("ant$i") for i in 1:n_anticipated_shocks(m)]
    end

    # For each additional source, search for the file with the proper name. Open
    # it, read it in, and merge it with fred_series

    vint = data_vintage(m)

    for source in keys(m.data_series)

        # Skip FRED sources, which have already been handled
        if source == :fred
            continue
        end

        # Check that this source is actually used
        mnemonics = m.data_series[source]
        if isempty(mnemonics)
            warn("No series were specified from $(string(source))")
            continue
        end

        # Read and merge data from this source
        file = inpath(m, "data", "$(string(source))_$vint.csv")
        if isfile(file)
            if VERBOSITY[verbose] >= VERBOSITY[:low]
                println("Reading data from $file...")
            end

            # Read in dataset and check that the file contains data for the proper dates
            addl_data = readtable(file)

            # Convert dates from strings to quarter-end dates for date arithmetic
            format_dates!(:date, addl_data)

            # Warn on sources with incomplete data; missing data will be replaced with NaN
            # during merge.
            if !in(lastdayofquarter(start_date), addl_data[:date]) ||
               !in(lastdayofquarter(end_date), addl_data[:date])
                warn("$file does not contain the entire date range specified; NaNs used.")
            end

            # Make sure each mnemonic that was specified is present
            for series in mnemonics
                if !in(series, names(addl_data))
                    error("$(string(series)) is missing from $file.")
                end
            end

            # Extract just the columns and rows of the dataset we want, and merge them with
            # data
            cols = [:date; mnemonics]
            rows = start_date .<= addl_data[:date] .<= end_date

            addl_data = addl_data[rows, cols]
            df = join(df, addl_data, on=:date, kind=:outer)
        else
            # If series not found, use all NaNs
            addl_data = DataFrame(fill(NaN, (size(df,1), length(mnemonics))))
            names!(addl_data, mnemonics)
            df = hcat(df, addl_data)
            warn("$file was not found; NaNs used." )
        end
    end

    # turn NAs into NaNs
    na2nan!(df)

    sort!(df, cols = :date)

    return df
end

"""
```
save_data(m::AbstractModel, df::DataFrame)
```

Save `df` to disk as CSV. File is located in `inpath(m, "data")`.
"""
function save_data(m::AbstractModel, df::DataFrame)
    vint = data_vintage(m)
    filename = inpath(m, "data", "data_$vint.csv")
    writetable(filename, df)
end

"""
```
has_saved_data(m::AbstractModel)
```

Determine if there is a saved dataset on disk for the required vintage.
"""
function has_saved_data(m::AbstractModel)
    vint = data_vintage(m)
    filename = inpath(m, "data", "data_$vint.csv")
    isfile(filename)
end

"""
```
read_data(m::AbstractModel)
```

Read CSV from disk as DataFrame. File is located in `inpath(m, "data")`.
"""
function read_data(m::AbstractModel)
    vint     = data_vintage(m)
    filename = inpath(m, "data", "data_$vint.csv")
    df       = readtable(filename)

    # Convert date column from string to Date
    df[:date] = map(Date, df[:date])

    return df
end

"""
```
isvalid_data(m::AbstractModel, df::DataFrame)
```

Return if dataset is valid for this model, ensuring that all observables are contained and
that all quarters between the beginning of the presample and the end of the mainsample are
contained.
"""
function isvalid_data(m::AbstractModel, df::DataFrame)
    valid = true

    # Ensure that every series in m_series is present in df_series
    m_series = collect(keys(m.data_transforms))
    df_series = names(df)
    coldiff = setdiff(m_series, df_series)
    valid = valid && isempty(coldiff)
    if !isempty(coldiff)
        println("Columns of 'df' do not match expected.")
        println(coldiff)
    end

    # Ensure the dates between date_presample_start and date_mainsample_end are contained.
    actual_dates = df[:date]
    expected_dates = get_quarter_ends(get_setting(m, :date_presample_start),
                                      get_setting(m, :date_mainsample_end))
    datesdiff = setdiff(expected_dates, actual_dates)
    valid = valid && isempty(datesdiff)
    if !isempty(datesdiff)
        println("Dates of 'df' do not match expected.")
        println(datesdiff)
    end

    return valid
end

"""
```
df_to_matrix(m::AbstractModel, df::DataFrame)
```

Return `df`, converted to matrix of floats, and discard date column. Also ensure data are
sorted by date and that rows outside of sample are discarded. The output of this function is
suitable for direct use in `estimate`, `posterior`, etc.
"""
function df_to_matrix(m::AbstractModel, df::DataFrame)
    # Sort rows by date and discard rows outside of sample
    df1 = sort(df; cols=[:date])
    t_start = find(df1[:date] .== get_setting(m, :date_presample_start))[1]
    t_end = find(df1[:date] .== get_setting(m, :date_mainsample_end))[1]

    # Discard columns not used.
    #TODO sort columns as well according to observables
    datecol = df1.colindex[:date]
    colinds = setdiff(1:size(df1,2), datecol)

    df1 = df1[t_start:t_end, colinds]

    return convert(Matrix{Float64}, df1)
end


text:
using ProgressMeter
using Test

if get(ENV, "CI", "false") == "true"
    using InteractiveUtils
    display(versioninfo())   # among other things, this shows the number of threads
end

include("core.jl")
include("test.jl")
include("test_showvalues.jl")
include("test_map.jl")
include("test_float.jl")
include("test_threads.jl")

println("")
println("All tests complete")


text:
#= test/pheig.jl
=#

@testset "pheig" begin
  n = 128
  M = matrixlib(:fourier, rand(n), rand(n))
  M += M'
  opts = LRAOptions(maxdet_tol=0., sketch_randn_niter=1)

  @time for (t, s) in ((:none,                 :none ),
                      (:RandomGaussian,       :randn),
                      (:RandomSubset,         :sub  ),
                      (:SRFT,                 :srft ),
                      (:SparseRandomGaussian, :sprn ))
    opts.sketch = s
    for T in (Float32, Float64, ComplexF32, ComplexF64)
      let A
          println("  $t/$T")

          rtol = 5*eps(real(T))
          approx_rtol = 100*rtol
          opts.rtol = rtol

          A = convert(Array{T}, T <: Real ? real(M) : M)
          F = pheigfact(A, opts)
          @test norm(A - Matrix(F)) < approx_rtol*norm(A)

          s = pheigvals(A, opts, rank=F[:k], rtol=0.)
          @test norm(s - F[:values]) < approx_rtol*norm(s)

          x = rand(T, n)
          y = A  *x; @test norm(y - F  *x) < approx_rtol*norm(y)
          y = A' *x; @test norm(y - F' *x) < approx_rtol*norm(y)
          y = transpose(A)*x; @test norm(y - transpose(F)*x) < approx_rtol*norm(y)

          X = rand(T, n, n)
          C = A  *X  ; @test norm(C - F  *X  ) < approx_rtol*norm(C)
          C = A  *X' ; @test norm(C - F  *X' ) < approx_rtol*norm(C)
          C = A  *transpose(X); @test norm(C - F  *transpose(X)) < approx_rtol*norm(C)
          C = A' *X  ; @test norm(C - F' *X  ) < approx_rtol*norm(C)
          C = A' *X' ; @test norm(C - F' *X' ) < approx_rtol*norm(C)
          C = transpose(A)*X  ; @test norm(C - transpose(F)*X  ) < approx_rtol*norm(C)
          C = transpose(A)*transpose(X); @test norm(C - transpose(F)*transpose(X)) < approx_rtol*norm(C)
          C = X  *A  ; @test norm(C - X  *F  ) < approx_rtol*norm(C)
          C = X  *A' ; @test norm(C - X  *F' ) < approx_rtol*norm(C)
          C = X  *transpose(A); @test norm(C - X  *transpose(F)) < approx_rtol*norm(C)
          C = X' *A  ; @test norm(C - X' *F  ) < approx_rtol*norm(C)
          C = X' *A' ; @test norm(C - X' *F' ) < approx_rtol*norm(C)
          C = transpose(X)*A  ; @test norm(C - transpose(X)*F  ) < approx_rtol*norm(C)
          C = transpose(X)*transpose(A); @test norm(C - transpose(X)*transpose(F)) < approx_rtol*norm(C)

          x = A*rand(T, n)
          y = F\x; @test norm(x - A*y) < approx_rtol*norm(x)

          X = A*rand(T, n, n)
          C = F\X; @test norm(X - A*C) < approx_rtol*norm(X)
      end
    end
  end
end


text:
using QuantumPropagators
using Documenter

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

makedocs(;
    modules=[QuantumPropagators],
    authors="Michael Goerz <mail@michaelgoerz.net> and contributors",
    repo="https://github.com/JuliaQuantumControl/QuantumPropagators.jl/blob/{commit}{path}#{line}",
    sitename="QuantumPropagators.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://juliaquantumcontrol.github.io/QuantumPropagators.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Overview" => "overview.md",
        "Background" => "background.md",
        hide("Examples" => "examples/index.md", [
            # "Example 1" => joinpath("examples", "1.md"),
            # "Example 2" => joinpath("examples", "2.md"),
           ]
        ),
        "Howtos" => "howto.md",
        "Benchmarks" => "benchmarks.md",
        "API" => "api.md",
        "History" => "history.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaQuantumControl/QuantumPropagators.jl",
)


text:
using VisualRegressionTests
using Pkg, Test, Plots

# popup on local machine only
istravis = "TRAVIS" ∈ keys(ENV)
if !istravis
  Pkg.add("Gtk")
  using Gtk
end

# load an image for testing
import FileIO, TestImages
img = TestImages.testimage("cameraman")

# save a temporary copy
reffn = tempname() * ".png"
FileIO.save(reffn, img)

# this is the test function, which saves an image to the given location
func = fn -> FileIO.save(fn, img)

@testset "VisualRegressionTests.jl" begin
  @testset "Basic" begin
    result = test_images(VisualTest(func, reffn))
    @test isa(result, VisualTestResult)
    @test result.status == EXACT_MATCH
    @test result.err == nothing
    @test result.diff == 0
  end

  @testset "Macros" begin
    @visualtest func "VisualTest.png" !istravis

    @plottest plot([1.,2.,3.]) "PlotTest.png" !istravis

    @plottest begin
      plot([1.,2.,3.])
      plot!([3.,2.,1.])
    end "MorePlotTest.png" !istravis

    plotit() = plot([1.,2.,3.])
    @plottest plotit "FuncPlotTest.png" !istravis
  end
end


text:
module disk

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
module Chirps

using SearchLight, Faker, ChirpsValidator, Validation

export Chirp

mutable struct Chirp <: AbstractModel
  ### internals
  _table_name::String
  _id::String

  ### fields
  id::Nullable{SearchLight.DbId}
  content::String
  created_at::DateTime

  ### validator
  validator::ModelValidator

  ### relations
  # belongs_to::Vector{SearchLight.SQLRelation}
  # has_one::Vector{SearchLight.SQLRelation}
  # has_many::Vector{SearchLight.SQLRelation}

  ### callbacks
  # before_save::Function
  # after_save::Function
  # on_dehydration::Function
  # on_hydration::Function
  # after_hydration::Function

  ### scopes
  # scopes::Dict{Symbol,Vector{SearchLight.SQLWhereEntity}}

  ### constructor
  Chirp(;
    id = Nullable{SearchLight.DbId}(),
    content = "",
    created_at = Dates.now(),

    validator = ModelValidator([
      ValidationRule(:content, ChirpsValidator.not_empty),
      ValidationRule(:content, ChirpsValidator.minimum_length)
    ])

    # belongs_to = [],
    # has_one = [],
    # has_many = [],

    # before_save = (m::Chirp) -> warn("Not implemented"),
    # after_save = (m::Chirp) -> warn("Not implemented"),
    # on_dehydration = (m::Chirp, field::Symbol, value::Any) -> warn("Not implemented"),
    # on_hydration = (m::Chirp, field::Symbol, value::Any) -> warn("Not implemented")
    # after_hydration = (m::Chirp, field::Symbol, value::Any) -> warn("Not implemented")

    # scopes = Dict{Symbol,Vector{SearchLight.SQLWhereEntity}}()

  ) = new("chirps", "id",
          id, content, created_at,
          validator
          # belongs_to, has_one, has_many,
          # before_save, after_save, on_dehydration, on_hydration, after_hydration
          # scopes
          )
end

function random()
  Chirp(content = Faker.sentence())
end

end

max_stars_repo_name: essenciary/chirper

max_issues_repo_name: essenciary/chirper

max_forks_repo_name: essenciary/chirper


text:
"""is_nonbacktracking: a predicate that tests if a walk is nonbacktracking.
    That is that for no index i walk[i+2]==walk[i]
"""
function is_nonbacktracking(walk)
    n = length(walk)
    if n < 3
       return true
    end
    for i in 1:n-2
        if walk[i+2] == walk[i]
            return false
        end
    end
    return true
end

"""test_nbw: check that a walk is nonbacktracking and print it if is isn't.
    Used only for testing and debugging.
"""
function test_nbw(g, start, len)
    w = non_backtracking_randomwalk(g, start, len)
    if is_nonbacktracking(w)
        return true
    else
        print("walk was:\n  $w")
    end
    return false
end
g = PathDiGraph(10)
@test randomwalk(g, 1, 5) == [1:5;]
@test randomwalk(g, 2, 100) == [2:10;]
@test_throws BoundsError randomwalk(g, 20, 20)

g = PathGraph(10)
@test saw(g, 1, 20) == [1:10;]
@test_throws BoundsError saw(g, 20, 20)

g = PathGraph(10)
@test non_backtracking_randomwalk(g, 1, 20) == [1:10;]
@test_throws BoundsError non_backtracking_randomwalk(g, 20, 20)

g = DiGraph(PathGraph(10))
@test non_backtracking_randomwalk(g, 1, 20) == [1:10;]
@test_throws BoundsError non_backtracking_randomwalk(g, 20, 20)

g = PathDiGraph(10)
@test non_backtracking_randomwalk(g, 10, 20) == [10]

g = PathDiGraph(10)
@test non_backtracking_randomwalk(g, 1, 20) == [1:10;]

g = CycleGraph(10)
visited = non_backtracking_randomwalk(g, 1, 20)
@test visited == [1:10; 1:10;] || visited == [1; 10:-1:1; 10:-1:2;]

g = CycleDiGraph(10)
@test non_backtracking_randomwalk(g, 1, 20) == [1:10; 1:10;]

n = 10
g = CycleGraph(n)
for k = 3:n-1
    add_edge!(g, 1, k)
end

for len = 1:3*n
    @test test_nbw(g,1,len)
    @test test_nbw(g,2,len)
end
#test to make sure it works with self loops.
add_edge!(g, 1, 1)
for len = 1:3*n
    @test test_nbw(g,1,len) == true
    @test test_nbw(g,2,len) == true
end


text:
# g_tests_hybrid_opt.jl
#
# A Julia script searching for maximum cuts of 0-1 weighted benchmark
# graphs from G-set using the triangular model based solver
#
# Implements a basic iterated search with local postprocessing
#

using LightGraphs
using MatrixDepot
using GraphRecipes

include("DiceBasic.jl")

listTests = [1, 2 , 3, 4, 5, 22, 23, 24, 25, 26, 43, 44,
            45, 46, 47, 48, 49, 50, 51, 52, 53, 54]

# Common model definitions
Ks = 0.00001
tmax = 350      # the number of steps within one realization
domain = 0.4    # the magnitude of random perturbations

# How many returns
Nreruns = 200

# Array with the best cuts
bestcuts = zeros(length(listTests))
besthybridcuts = zeros(length(listTests))

for i in 1:length(listTests)
    nameG = "Gset/G"*string(listTests[i])
    println("Prosessing $nameG")
    Gcur = Graph(mdopen(nameG).A)
    Nvert = nv(Gcur)

    model = DiceBasic.Model(Gcur, DiceBasic.triangular, 20/Nvert, Ks)

    bestcut = -1
    besthybridcut = -1
    bestconf = zeros(Nvert)

    dom = domain
    for j in 1:Nreruns
        print("\rRealization $j: ")

        # Form the initial state by perturbing the seed state
        # Here, the perturbation has the form of rescaling and additive noise
        Vstart = dom.*DiceBasic.get_initial(Nvert, (-1.0, 1.0)) .+ bestconf.*0.8
        dom /=1.01
        Vend = DiceBasic.roundup(DiceBasic.propagate(model, tmax, Vstart))

        (cut, conf) = DiceBasic.get_best_configuration(Gcur, Vend)
        print("Pre-local cut = $cut,  ")
        if cut > bestcut
            bestcut = cut       # for keeping records of the "bare" performance
        end

        # the local search reorienting nodes against the majority
        conf = DiceBasic.local_search(model.graph, conf)
        print("Post-local cut = $(DiceBasic.cut(model.graph, conf)) ")

        # the local search reorienting edges against the majority
        conf = DiceBasic.local_twosearch(model.graph, conf)
        cut = DiceBasic.cut(model.graph, conf)
        print("Post-two-local cut = $cut ")

        if cut > besthybridcut
            besthybridcut = cut
            bestconf = conf     # the obtained partition makes the seed state
        end
    end
    # From https://discourse.julialang.org/t/how-clear-the-printed-content-in-terminal-and-print-to-the-same-line/19549/4
    # with the reference to https://en.wikipedia.org/wiki/ANSI_escape_code#CSI_sequences
    print("\e[2K") # clear line
    print("\e[1G") # move the cursor to column 1
    println("The best cut is $bestcut/$besthybridcut")
    bestcuts[i] = bestcut
    besthybridcuts[i] = besthybridcut
end

max_stars_repo_name: merement/Dice

max_issues_repo_name: merement/Dice

max_forks_repo_name: merement/Dice


text:
module GLIntrinsics

using StaticArrays, Sugar, GeometryTypes

import ..Transpiler: ints, floats, numbers, Numbers, Floats, int, Ints, uchar
import ..Transpiler: fixed_array_length, is_ntuple, is_fixedsize_array, GLMethod
import ..Transpiler: AbstractGLIO, ret, vecs, Vecs, vector_lengths, functions

import Sugar: typename, vecname
using SpecialFunctions: erf

immutable GLArray{T, N} <: AbstractArray{T, N} end
immutable GLTexture{T, N} <: AbstractArray{T, N} end

const GLDeviceArray = Union{GLArray, GLTexture}
const Types = Union{vecs..., numbers..., GLArray, GLTexture}
const Functions = Union{map(typeof, functions)...}

glsl_sizeof(T) = sizeof(T) * 8
# for now we disallow Float64 and map it to Float32 -> super hack alert!!!!
glsl_sizeof(::Type{Float64}) = 32
glsl_length{T <: Number}(::Type{T}) = 1
glsl_length(T) = length(T)
prescripts = Dict(
    Float32 => "",
    Float64 => "", # ignore float64 for now
    Int => "i",
    Int32 => "i",
    UInt => "u",
    Bool => "b"
)

function typename{T, N}(io::AbstractGLIO, x::Type{GLTexture{T, N}})
    string(prescripts[eltype(T)], "sampler", N, "D")
end
function typename{T, N}(io::AbstractGLIO, x::Type{GLArray{T, N}})
    if !(N in (1, 2, 3))
        # TODO, fake ND arrays with 1D array
        error("GPUArray can't have more than 3 dimensions for now")
    end
    sz = glsl_sizeof(T)
    len = glsl_length(T)
    if true
        qualifiers = [image_format(x)]
        string("layout (", join(qualifiers, ", "), ") image$(N)D")
    else
        "image$(N)D$(len)x$(sz)_bindless"
    end
end
function typename{T <: SMatrix}(io::AbstractGLIO, ::Type{T})
    M, N = size(T)
    string(prescripts[eltype(T)], "mat", M == N ? M : string(M, "x", N))
end

function vecname{T}(io::AbstractGLIO, t::Type{T})
    N = fixed_array_length(T)
    ET = eltype(T)
    return string(prescripts[ET], "vec", N)
end

imageStore{T}(x::GLArray{T, 1}, i::int, val::NTuple{4, T}) = nothing
imageStore{T}(x::GLArray{T, 2}, i::NTuple{2, int}, val::NTuple{4, T}) = nothing

imageLoad{T}(x::GLArray{T, 1}, i::int) = ret(NTuple{4, T})
imageLoad{T}(x::GLArray{T, 2}, i::NTuple{2, int}) = ret(NTuple{4, T})

texelFetch{T}(x::GLTexture{T, 1}, i::int, lod::int) = ret(NTuple{4, T})
texelFetch{T}(x::GLTexture{T, 2}, i::NTuple{2, int}, lod::int) = ret(NTuple{4, T})

texture{T}(x::GLTexture{T, 1}, i::Float32) = ret(NTuple{4, T})
texture{T}(x::GLTexture{T, 2}, i::Vec2f0) = ret(NTuple{4, T})
texture{T}(x::GLTexture{T, 3}, i::Vec3f0) = ret(NTuple{4, T})

imageSize{T, N}(x::GLArray{T, N}) = ret(NTuple{N, int})
textureSize{T, N}(::GLTexture{T, N}) = ret(NTuple{N, int})

"""
Gradient in x direction
This is sadly a bit hard to implement for a pure CPU versions, since it's pretty much backed into the GPU hardware.
How it seems to work is, that it takes the values from neighboring registers, which work in parallel on the pixels
of the triangle, so they actually do hold the neighboring values needed to calculate the gradient.
"""
dFdx{T}(value::T) = T(0.001) # just default to a small gradient if it's called on the CPU
dFdy{T}(value::T) = T(0.001) # just default to a small gradient if it's called on the CPU

EmitVertex() = nothing
EndPrimitive() = nothing

const gl_GlobalInvocationID = Vec3f0(0, 0, 0)
const gl_FragCoord = Vec4f0(0, 0, 0, 0)

function glintrinsic{F <: Function, T <: Tuple}(f::F, types::Type{T})
    glintrinsic(f, Sugar.to_tuple(types))
end
function glintrinsic{F <: Function}(f::F, types::Tuple)
    # we rewrite Ntuples as glsl arrays, so getindex becomes inbuild
    if f == broadcast
        BF = types[1]
        if BF <: cli.Functions && all(T-> T <: Types || is_fixedsize_array(T), types[2:end])
            return true
        end
    end
    if f == getindex && length(types) == 2 && first(types) <: NTuple && last(types) <: Integer
        return true
    end
    # matmul
    if f == (*)  && length(types) == 2 && all(T-> T <: StaticArray, types)
        return true
    end
    m = methods(f)
    isempty(m) && return false
    sym = first(m).name
    (F <: cli.Functions && all(T-> T <: Types || is_fixedsize_array(T), types)) || (
        # if any intrinsic funtion stub matches
        isdefined(GLIntrinsics, sym) &&
        Base.binding_module(GLIntrinsics, sym) == GLIntrinsics &&
        length(methods(f, types)) == 1
    )
end

end # end GLIntrinsics

using .GLIntrinsics
using GeometryTypes
const gli = GLIntrinsics
import .gli: glintrinsic, GLArray, GLDeviceArray, GLTexture

import Sugar.isintrinsic

function gli.glintrinsic{T}(x::Type{T})
    T <: gli.Types ||
    is_fixedsize_array(T) ||
    T <: Tuple{Numbers} ||
    (T <: NTuple{1} && is_fixedsize_array(eltype(T))) ||
    (T <: SMatrix && all(x-> x <= 4, size(T))) ||
    T <: uchar # uchar in ints makes 0.6 segfault -.-
end
function isintrinsic(x::GLMethod)
    if isfunction(x)
        isintrinsic(Sugar.getfunction(x)) ||
        gli.glintrinsic(x.signature...)
    else
        gli.glintrinsic(x.signature)
    end
end

# copied from rewriting. TODO share implementation!

# Make constructors inbuild for now. TODO, only make default constructors inbuild
glintrinsic{T}(f::Type{T}, types::ANY) = true

# homogenous tuples, translated to glsl array
function glintrinsic{T, I}(
        f::typeof(getindex), types::Type{Tuple{T, I}}
    )
    return is_fixedsize_array(T) && I <: Union{StaticArray, Integer}
end

function glintrinsic{T <: GLDeviceArray, Val, I <: Integer}(
        f::typeof(setindex!), types::Type{Tuple{T, Val, I}}
    )
    return true
end

glintrinsic(f::typeof(tuple), types::Tuple) = true

import Base: getindex, setindex!, size

GlobalInvocationID() = gli.gl_GlobalInvocationID
FragCoord() = gli.gl_FragCoord

size{T, N}(x::gli.GLArray{T, N}) = gli.imageSize(x)
size{T, N}(x::gli.GLTexture{T, N}) = gli.textureSize(x)

getindex{T}(x::gli.GLTexture{T, 1}, i::Integer) = gli.texelFetch(x, i, 0)
getindex{T}(x::gli.GLTexture{T, 1}, i::Integer, j::Integer) = gli.texelFetch(x, (i, j), 0)

getindex{T}(x::gli.GLTexture{T, 1}, i::AbstractFloat) = gli.texture(x, i)
getindex{T}(x::gli.GLTexture{T, 2}, i::AbstractFloat, j::AbstractFloat) = gli.texture(x, Vec2f0(i, j))
getindex{T}(x::gli.GLTexture{T, 2}, idx::Vec2f0) = gli.texture(x, idx)
getindex{T}(x::gli.GLTexture{T, 3}, idx::Vec3f0) = gli.texture(x, idx)

getindex{T}(x::gli.GLArray{T, 1}, i::Integer) = gli.imageLoad(x, i)
function Base.getindex{T}(x::gli.GLArray{T, 2}, i::Integer, j::Integer)
    getindex(x, (i, j))
end
function Base.getindex{T <: Number}(x::gli.GLArray{T, 2}, idx::NTuple{2, int})
    gli.imageLoad(x, idx)[1]
end
function Base.setindex!{T}(x::gli.GLArray{T, 1}, val::T, i::Integer)
    gli.imageStore(x, i, (val, val, val, val))
end
function Base.setindex!{T}(x::gli.GLArray{T, 2}, val::T, i::Integer, j::Integer)
    setindex!(x, (val, val, val, val), (i, j))
end
function Base.setindex!{T}(x::gli.GLArray{T, 2}, val::T, idx::NTuple{2, int})
    gli.imageStore(x, idx, (val, val, val, val))
end
function Base.setindex!{T}(x::gli.GLArray{T, 2}, val::NTuple{4, T}, idx::NTuple{2, int})
    gli.imageStore(x, idx, val)
end
function Base.setindex!{T}(x::gli.GLArray{T, 1}, val::NTuple{4, T}, i::Integer)
    gli.imageStore(x, i, val)
end

function gl_erf{T <: AbstractFloat}(x::T)
    # constants
    a1 =  T(0.254829592)
    a2 = T(-0.284496736)
    a3 =  T(1.421413741)
    a4 = T(-1.453152027)
    a5 =  T(1.061405429)
    p  =  T(0.3275911)

    # Save the sign of x
    sign = 1
    if (x < 0)
        sign = -1
    end
    xabs = abs(x)
    # A&S formula 7.1.26
    t = T(1.0) / (T(1.0) + p * xabs)
    y = T(1.0) - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-xabs*xabs)

    return sign * y
end
gl_erfc{T <: AbstractFloat}(x::T) = T(1.0) - gl_erf(x)

# FMA is only supported in opengl >= 4.0 . To keep it simple, we just use this fallback for now
# until we propagate opengl versions and emit code accordingly
gl_fma{T <: AbstractFloat}(a::T, b::T, c::T) = a * b + c


text:
function main(args)
    return args
end
# Load inline arguments for worker unit
num_workers, this = main(ARGS)
num_workers = parse(Int64, num_workers)
this        = parse(Int64, this)

#
# Load packages and functions
include("../../src/init.jl")
include("../../src/fct.jl")
include("../../src/worker_fct.jl")
# using Distributed
#
# Output paths and filter settings
plots_base_path = "../../../waveforms/baselines/plots/raw/"
base_path_raw   = "../../../waveforms/baselines/raw_wf/"
base_path       = "../../../waveforms/baselines/wf/"
# addprocs(3)
# nprocs()

cal = JSON.parsefile("../../dicts/A_cal.json");
channels = 0:1:36
channels = get_share_for_worker(channels, num_workers, this)
for ch in channels
    files = glob(base_path_raw * "*/" * lpad(ch,2,"0") * "-" * channel_to_name[ch] * "/*.h5");
#     IJulia.clear_output(true)
    Base.run(`clear`)
    @info("Channel Ch" * lpad(ch, 2, "0") * " in progress")
    @showprogress for file in files
        filename = base_path * split(split(file, "/raw_wf/")[2], "-bl-raw.h5")[1] * ".h5"
        run_str = split(basename(file), "-")[2]
        run = parse(Int64, split(run_str, "run")[end])
        if haskey(cal, channel_to_name[ch])
            if haskey(cal[channel_to_name[ch]], run_str) && !isfile(filename)
                data = HDF5.h5open(file, "r") do h5f
                    LegendHDF5IO.readdata(h5f, "data")
                end
                cal_pulses = []
                map(x->push!(cal_pulses, RDWaveform(x.time, (float.(x.value) .- mean(x.value)) ./ cal[channel_to_name[ch]][run_str][1] )), data.waveform);

                !isdir(dirname(filename)) ? mkpath(dirname(filename)) : ""
                HDF5.h5open(filename, "w") do h5f
                    LegendHDF5IO.writedata(h5f, "data", Table(
                            energy = data.energy,
                            multiplicity = data.multiplicity,
                            timestamp = data.timestamp,
                            run = data.run,
                            datasetID = data.datasetID,
                            channel = data.channel,
                            waveform = StructArray{RDWaveform}(Array{RDWaveform,1}(cal_pulses))
                        )
                    )
                end
            end
        end
    end
end


text:
module Config

import JSON3
using Memento

struct AppConfig
    loglevel::String
end

function getconfig()::AppConfig
    loglevel = lowercase(get(ENV, "LOG_LEVEL", "warn"))
    config = AppConfig(loglevel)
    configure_logger(config)

    return config
end

function configure_logger(config::AppConfig)
    handler = DefaultHandler(
        stderr,
        DictFormatter(Dict(:time => :date, :msg => :msg, :level => :level), JSON3.write))

    logger = Memento.getlogger()
    logger.handlers["console"] = handler
    setlevel!(logger, config.loglevel; recursive = true)
    Memento.substitute!()
end

end

max_stars_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia

max_issues_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia

max_forks_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia


text:
#!/usr/bin/env julia

using Seq

fibs = cat([0, 1], @lazyseq map(+, rest(fibs), fibs))

for k in take(10, fibs)
    println(k)
end


text:
### A Nbdev script file ###
### Autogenerated file. Don't modify. ###

module Export
#export
using ReusePatterns

#export
import Pluto: Notebook, Cell, load_notebook

#export
function processMdChunk(mdChunk::String)
	chopMdMarker(mdChunk)
end

#export
begin
Base.@kwdef mutable struct Nucleus
	code::String=""
	parsedcode::Union{Nothing,Expr}=nothing
	end

code(cell::Cell) = cell.code
parsedcode(cell::Cell) = cell.parsedcode

Nucleus(code, parsedcode) = Nucleus(code = code, parsedcode = parsedcode)
end

#export
begin
mutable struct Nb
	nuclei::Array{Nucleus,1}
	name::String
end

nuclei(nuclei::Array) = nuclei
name(name::String) = name
end

#export
begin
const _header = "### A Nbdev script file ###"
const _sub_header = "### Autogenerated file. Don't modify. ###"
end

#export
function collect_nuclei(notebook::Notebook, marker::String)
	notebook_cells=notebook.cells
	nuclei=[]
	for i in 1:length(notebook_cells)
		raw_code=code(notebook_cells[i])
		parsed_code=parsedcode(notebook_cells[i])
		nucleus=Nucleus(raw_code, parsed_code)

		if startswith(raw_code, marker)
	        push!(nuclei,nucleus)
        end

	end
	nuclei
end

#export
begin

load_nb(filename::String, marker::String) = _load_notebook(filename, marker)

function _load_notebook(filename::String, marker::String)
		notebook=load_notebook(filename)
		collected_nuclei=nuclei(collect_nuclei(notebook, marker))
		Nb(collected_nuclei, filename)
end

end

#export
strip=(x,y) -> replace(x, y=>"")

#export
begin

function save_nb(io, notebook)
    println(io, _header)
    println(io, _sub_header)
    println(io, "")
    println(io, "module $(uppercasefirst(strip(strip(notebook.name, r"[0-9_]"), r".jl")))")
	notebook
	print(io, "end")
end

function save_nb(notebook::Nb, path::String)
	open(path, "w") do io
        save_nb(io, notebook)
    end
end
end

#export
begin
function readfilenames()
	files=readdir()
	files
end

function readfilenames(dir::String)
	files=cd(readdir, dir)
	files
end
end

#export
function export_Code(file::String, code_cells)
	io=IO
    open(file, "w") do io
		for code in code_cells
            print(io, code*"\n\n")
        end
	end
end

#export
begin
function export_file(file, path::String, marker::String)
	    code_cells=collect_codecells(file, marker)
		file_path=joinpath(path,file)
        export_Code(file_path, code_cells)
end
end

#export
function export_content(files::AbstractVector, path::String, marker::String)
	map(file->export_file(file, path, marker), files)
end

#export
function getfile_extension(filename)
    return filename[findlast(isequal('.'),filename):end]
end

#export
#TODO: Too many loops. Need to clear this either
#by broadcasting or something more efficient
#TODO: User should eb able to get messages about which
#files are updated
#TODO: currently output inside loop is not shown in pluto notebooks.
#so, the show_output displays all stuffs in the console.
#maybe a progress bar would be better
function notebook2script()
	files=[file for file in readfilenames() if getfile_extension(file)== ".jl"]
	export_content(files, "../src", "#export")
end

end


text:
## The following exercises use the ThinkJulia module, described in Case Study: Interface Design:
## Exercise 5-6
## Read the following function and see if you can figure out what it does (see the examples in Case Study: Interface Design). Then run it and see if you got it right.
using ThinkJulia

function draw(t, length, n)
    if n == 0
        return
    end
    angle = 50
    forward(t, length*n)
    turn(t, -angle)
    draw(t, length, n-1)
    turn(t, 2*angle)
    draw(t, length, n-1)
    turn(t, -angle)
    forward(t, -length*n)
end

println("Ans: ")

# Ans: It should recursive pattern
turtle = Turtle()

@svg begin
    draw(turtle, 10, 5)
end

println("End.")


text:
function ultra2cheb(c_ultra::AbstractVector{T},λ::T,plan::ChebyshevUltrasphericalPlan{FORWARD,T}) where T<:AbstractFloat
    M,N,nM₀,αN,K = getconstants(plan)
    i₁,i₂,j₁,j₂ = getindices(plan)
    p₁,p₂,rp,c₁,c₂,um,vm,θ,tempsin,tempsin2,tempsinλ,tempsinλm,tempmindices,cnλ,cnmλ = getplan(plan)

    v_cheb = zero(c_ultra)

    k=0
    for j=1:K
        k = j
        copyto!(cnmλ,cnλ)
        copyto!(tempsinλm,tempsinλ)

        j₁[k] < nM₀ && (k-=1; break)

        @inbounds for i=i₁[k]:i₁[k+1]-1 v_cheb[i] += backward_recurrence(c_ultra,j₂[k],θ[i],tempsin2[N+2-i],tempsin2[i],rp) end
        @inbounds for i=i₂[k+1]+1:i₂[k] v_cheb[i] += backward_recurrence(c_ultra,j₂[k],θ[i],tempsin2[N+2-i],tempsin2[i],rp) end

        for m=0:M-1
            # Initialize diagonal N-scaling multiplied by local coefficients and zero out excess
            init_c₁c₂!(c₁,c₂.parent,cnmλ,c_ultra,j₁[k],j₂[k])

            # Apply planned DCT-I and DST-I in-place
            applyTN!(c₁,p₁);applyUN!(c₂,p₂)

            # Compute u_m(θ) and v_m(θ)
            compute_umvm!(um,vm,λ,tempsinλm,m,θ,i₁[k+1]:i₂[k+1])

            # Multiply point-wise by u_m(θ) and v_m(θ) for valid indices
            @inbounds for i=i₁[k+1]:i₂[k+1] v_cheb[i] += um[i]*c₁[i] + vm[i]*c₂.parent[i] end

            # Update C_{n,m}^λ by recurrence in m
            cst = (λ+m)/2/(m+1)*(m+1-λ)
            @inbounds for i=1:N+1 cnmλ[i] *= cst/(i+λ+m) end

            # Update sin^{λ+m}(θ)
            @inbounds for i=1:N+1 tempsinλm[i] *= tempsin[i] end
        end
    end

    # Finish off recurrence
    @inbounds for i=i₁[k+1]:i₂[k+1] v_cheb[i] += backward_recurrence(c_ultra,j₂[k+1],θ[i],tempsin2[N+2-i],tempsin2[i],rp) end

    # perform IDCT-I
    applyTNinv!(v_cheb,p₁)
end


text:
using DistributionsAD, AdvancedVI

########################################################################
function bpca(X, Σ; maxiter=1, S=1)
########################################################################

    Q = 1

    JITTER = 1e-8

    D, N = size(X)

    @printf("BPCA: %d number of data items of dimension %d, projected to %d components\n", N, D, Q)

    #---------------------------------------------
    function logPrior(W, μ)
    #---------------------------------------------

        -1e-10*(sum(W.^2) + sum(μ.^2))

    end

    #---------------------------------------------
    function marginalLogLikelihood(W, μ)
    #---------------------------------------------

        local ℓ = zero(eltype(W))

        @inbounds for n in 1:N

            local C = W*W' + Matrix(Diagonal(Σ[:,n])) + JITTER*I

            ℓ +=  -0.5*logdet(C)-0.5*sum((X[:,n].-μ)'*(C\(X[:,n].-μ)))

        end

        return ℓ

    end

    #---------------------------------------------
    function unpack(param)
    #---------------------------------------------

        local MARK = 0

        local W = reshape(param[MARK+1:MARK+D*Q], D, Q)

        MARK += D*Q

        local μ = reshape(param[MARK+1:MARK+D], D)

        MARK += D

        @assert(MARK == length(param))

        return W, μ

    end

    #---------------------------------------------
    function objective(param)
    #---------------------------------------------

        @assert(length(param) == D*Q + D)

        local W, μ = unpack(param)

        return marginalLogLikelihood(W, μ) + logPrior(W, μ)

    end

    #---------------------------------------------
    # Initialise with ppca estimates
    #---------------------------------------------

    W, μ, = ppca(X)

    #---------------------------------------------
    # Run optimiser
    #---------------------------------------------

    # options for initialisation
    opt0 = Optim.Options(show_trace = false, show_every = 5, iterations = 1)

    # options for optimisations
    opt  = Optim.Options(show_trace = false, show_every = 5, iterations = maxiter, g_tol=1e-6)

    mu, C, elbof = VIoverfitting(objective, [vec(W); vec(μ)], S = S, optimiser=LBFGS(), optimoptions = opt, initoptions = opt0)

    # @show elbof(mu, C, [randn(2D) for s=1:20_000])

    makeC(x) = (local C=reshape(x,2D,2D); cholesky(C*C'))
    getq(θ)  = TuringDenseMvNormal(θ[1:2D], makeC(θ[2D+1:end]))
    advi     = ADVI(30, 10000)
    qadvi    = vi(objective, advi, getq, randn(2D + 2D*2D))

    @show AdvancedVI.elbo(advi, qadvi, objective, 20_000)

    @show elbof(mu, C, [randn(2D) for s=1:20_000])

    return MvNormal(mu, C*C'),  elbof#, historyparameter

end


text:
## methods for core

core_object_methods = (:as_poly, :atoms,
                       :compare, :compare_pretty,
                       :doit, :dummy_eq, # :count,
                       :has,
                       :sort_key,
                       :args_cnc,
                       :as_coeff_Add, :as_coeff_Mul, :as_coeff_add,
                       :as_coeff_exponent, :as_coeff_factors, :as_coeff_mul,
                       :as_coeff_terms, :as_coefficient, :as_coefficient_dict,
                       :as_content_primitive, :as_expr, :as_independent,
                       :as_leading_term, :as_numer_denom, :as_ordered_factors,
                       :as_ordered_terms, :as_powers_dicts, :as_real_imag,
                       :as_terms,
                       :coeff, :compute_leading_term, :could_extract_minus_sign,
                       :count_ops, :equals,  :extract_additively,
                       :extract_branch_factor, :extract_multiplicatively,
                       :getO, :getn,
                       ## :is_integer called as x[:is_integer]
                       :is_constant, :is_polynomial, :is_rational_function,
                       :integer_nthroot,
                       :leadterm,
                       :primitive,
                       :removeO,
                       :series
                       )

core_object_properties = (:assumptions0,
                          :is_even, :is_odd,
                          :is_number, :is_integer, :is_real, :is_complex, :is_rational,
                          :is_commutative
                          )

"""
Return a vector of free symbols in an expression
"""
free_symbols(ex::Sym) =  convert(Vector{Sym}, collect(ex.x[:free_symbols]))
function free_symbols{T<:SymbolicObject}(exs::Vector{T})
    as = map(free_symbols, exs)
    out = as[1]
    if length(as) > 1
        for j in 2:length(as)
            for u in as[j]
                u in out || push!(out, u)
            end
        end
    end
    out
end
function free_symbols(exs::Tuple)
    as = map(free_symbols, exs)
    out = as[1]
    if length(as) > 1
        for j in 2:length(as)
            for u in as[j]
                u in out || push!(out, u)
            end
        end
    end
    out
end
export free_symbols

## From relational
core_sympy_methods = (:sympify,
                      :Dummy,
                      :Mod, :Rel,
                      :Eq, :Ne, :Lt, :Le, :Gt, :Ge,
                      :Equality, :Unequality,
                      :GreaterThan, :LessThan, :StrictGreaterThan, :StrictLessThan,
                      :PoleError,
                      :count_ops,
                      :gcdex, :half_gcdex,
                      :igcd, :ilcm
                      )

"""
Extract left and right hand side of a relation, parts of a relation.

(These are properties in SymPy, functions in SymPy.jl)

Examples:
```
x = Sym("x")
Eq(x, sin(x)) |> rhs  ## sin(x)
```
"""
rhs(ex::Sym, args...; kwargs...) = ex.x[:rhs]
lhs(ex::Sym, args...; kwargs...) = ex.x[:lhs]

"""

Returns a tuple of arguments

cf. [args](http://docs.sympy.org/latest/modules/core.html#sympy.core.basic.Basic.args)

(args is a property in SymPy, a function call in SymPy.jl.)

Examples
```
Eq(x, x^2) |> args ## (x, x^2)
sin(x) |> args ## (x,)
```
"""
args(ex::Sym) = ex.x[:args]

## need to import these
core_sympy_methods_base = (:factorial,
                           :gcd, :lcm,
                           :isqrt
                           )
for meth in core_sympy_methods_base
    meth_name = string(meth)
    @eval begin
                @doc """
`$($meth_name)`: a SymPy function.
The SymPy documentation can be found through: http://docs.sympy.org/latest/search.html?q=$($meth_name)
""" ->
        ($meth)(ex::Sym, args...; kwargs...) =
            sympy_meth($meth_name, ex, args...; kwargs...)
    end
end


text:
struct SLFN{T <: Number}
    n_features::Int
    n_outputs::Int
    hidden_layer::ImmutableHiddenLayer{T}
    output_weights::Matrix{T}
end

function predict(model::T₁,
                 samples::T₂) where {T₁ <: SLFN,
                                     T₂ <: AbstractMatrix}
    X = samples
    H = project(model.hidden_layer, X)
    Β = model.output_weights
    Y = Β * H
    typeof(Y) <: AbstractMatrix && first(size(Y)) == 1 ? vec(Y) : Y
end

function predict(model::T₁,
                 sample::T₂) where {T₁ <: SLFN,
                                    T₂ <: AbstractVector}
    predict(model, reshape(sample, 1, :))
end


text:
@defcomp impactcardiovascularrespiratory begin
    regions = Index()

    basecardvasc = Variable(index=[time,regions])
    baseresp = Variable(index=[time,regions])

    cardheat = Variable(index=[time,regions])
    resp = Variable(index=[time,regions])
    cardcold = Variable(index=[time,regions])

    cardvasc90 = Parameter(index=[regions])
    plus90 = Parameter(index=[regions])
    resp90 = Parameter(index=[regions])
    chplbm = Parameter(index=[regions])
    chmlbm = Parameter(index=[regions])
    chpqbm = Parameter(index=[regions])
    chmqbm = Parameter(index=[regions])
    rlbm = Parameter(index=[regions])
    rqbm = Parameter(index=[regions])
    ccplbm = Parameter(index=[regions])
    ccmlbm = Parameter(index=[regions])
    ccpqbm = Parameter(index=[regions])
    ccmqbm = Parameter(index=[regions])

    plus = Parameter(index=[time,regions])
    temp = Parameter(index=[time,regions])
    urbpop = Parameter(index=[time,regions])
    population = Parameter(index=[time,regions])

    cvlin = Parameter(default = 0.025901)
    rlin = Parameter(default = 0.001583)
    maxcardvasc = Parameter(default = 0.05)

    function run_timestep(p, v, d, t)

        if !is_first(t)
            for r in d.regions
                v.basecardvasc[t, r] = p.cardvasc90[r] + p.cvlin * (p.plus[t, r] - p.plus90[r])
                if v.basecardvasc[t, r] > 1.0
                    v.basecardvasc[t, r] = 1.0
                end

                v.baseresp[t, r] = p.resp90[r] + p.rlin * (p.plus[t, r] - p.plus90[r])
                if v.baseresp[t, r] > 1.0
                    v.baseresp[t, r] = 1.0
                end

                v.cardheat[t, r] = (p.chplbm[r] * p.plus[t, r] + p.chmlbm[r] * (1.0 - p.plus[t, r])) * p.temp[t, r] +
                        (p.chpqbm[r] * p.plus[t, r] + p.chmqbm[r] * (1.0 - p.plus[t, r])) * p.temp[t, r]^2
                v.cardheat[t, r] = v.cardheat[t, r] * p.urbpop[t, r] * p.population[t, r] * 10
                if v.cardheat[t, r] > 1000.0 * p.maxcardvasc * v.basecardvasc[t, r] * p.urbpop[t, r] * p.population[t, r]
                    v.cardheat[t, r] = 1000 * p.maxcardvasc * v.basecardvasc[t, r] * p.urbpop[t, r] * p.population[t, r]
                end
                if v.cardheat[t, r] < 0.0
                    v.cardheat[t, r] = 0
                end

                v.resp[t, r] = p.rlbm[r] * p.temp[t, r] + p.rqbm[r] * p.temp[t, r]^2
                v.resp[t, r] = v.resp[t, r] * p.urbpop[t, r] * p.population[t, r] * 10
                if v.resp[t, r] > 1000 * p.maxcardvasc * v.baseresp[t, r] * p.urbpop[t, r] * p.population[t, r]
                    v.resp[t, r] = 1000 * p.maxcardvasc * v.baseresp[t, r] * p.urbpop[t, r] * p.population[t, r]
                end
                if v.resp[t, r] < 0
                    v.resp[t, r] = 0
                end

                v.cardcold[t, r] = (p.ccplbm[r] * p.plus[t, r] + p.ccmlbm[r] * (1.0 - p.plus[t, r])) * p.temp[t, r] +
                        (p.ccpqbm[r] * p.plus[t, r] + p.ccmqbm[r] * (1.0 - p.plus[t, r])) * p.temp[t, r]^2
                v.cardcold[t, r] = v.cardcold[t, r] * p.population[t, r] * 10
                if v.cardcold[t, r] < -1000 * p.maxcardvasc * v.basecardvasc[t, r] * p.population[t, r]
                    v.cardcold[t, r] = -1000 * p.maxcardvasc * v.basecardvasc[t, r] * p.population[t, r]
                end
                if v.cardcold[t, r] > 0
                    v.cardcold[t, r] = 0
                end
            end
        end
    end
end


text:
using OrdinaryDiffEq, MultiScaleArrays

# Create Structures that define a Cell and Population
struct Cell{B} <: AbstractMultiScaleArrayLeaf{B}
    values::Vector{B}
end

struct Population{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArrayHead{B}
    nodes::Vector{T}
    values::Vector{B}
    end_idxs::Vector{Int}
end

nCells=2
cellConstructor = [Cell([1.0, 2.0]) for i=1:nCells]
cellPop = construct(Population, deepcopy(cellConstructor))

function StateModel(dpop,pop,p,t) #This is arbitrary
    for (cell,dcell) in LevelIter(1,pop,dpop)
        dcell = 1.0
    end
end

prob = ODEProblem(StateModel,cellPop,(0.0,1.0))

integrator = init(prob,Tsit5())

add_node!(integrator,Cell([5.0, 20.0]))

@test integrator.u.nodes[end].values == [5.0, 20.0]


text:
using BenchmarkTools
using Descartes

println("Descartes Benchmarks")

c = Cuboid(3,3,3)

println("Signed Distance Field")
t = @benchmark SignedDistanceField(c)
println(IOContext(stdout, :verbose => true, :compact => false), t)


text:
# Solving VRPTW
# Modified from the example given by Google OR-Tools https://developers.google.com/optimization/routing/vrp

using Routing, OffsetArrays
# include("../src/Routing_include.jl")

using Test

n_customers = 16

# Time windows, only for customer nodes
time_windows = [
    (7, 12),    # customer  1
    (10, 15),   # customer  2
    (16, 18),   # customer  3
    (10, 13),   # customer  4
    (0, 5),     # customer  5
    (5, 10),    # customer  6
    (0, 4),     # customer  7
    (5, 10),    # customer  8
    (0, 3),     # customer  9
    (10, 16),   # customer 10
    (10, 15),   # customer 11
    (0, 5),     # customer 12
    (5, 10),    # customer 13
    (7, 8),     # customer 14
    (10, 15),   # customer 15
    (11, 15),   # customer 16
]
early_time = [time_windows[i][1] for i in 1:n_customers]
late_time = [time_windows[i][2] for i in 1:n_customers]

# Service time is zero in each node in this example
service_time = zeros(n_customers)

# Customer demand, only for customer nodes
load = [1, 1, 2, 4, 2, 4, 8, 8, 1, 2, 1, 2, 4, 4, 8, 8]

# Create a set of requests
requests = Vector{Request}()
for i in 1:n_customers
    request_id = i
    push!(requests, Request(request_id, early_time[i], late_time[i], load[i], service_time[i]))
end

# Create the fleet information
n_vehicles = n_customers    # assume there are enough number of vehicles
capacity = 15               # vehicle capacity
max_travel_time = 35        # the latest time that each vehicle returns to the depot
fleet = Fleet(n_vehicles, capacity, max_travel_time)

# Create a set of nodes
nodes = Vector{Node}()
# if you can supply coordinates (cx, cy) to calculate the distance
# if you don't have coordinates, just put (0.0, 0.0)
push!(nodes, Node(0, 0.0, 0.0))   # (id, cx, cy)
for i in 1:n_customers
    node_id = i
    push!(nodes, Node(node_id, 0.0, 0.0))
end

# Since we didn't supply the coordinates, we need to provide the distance information
dists_data = Float64[
    0 6 9 8 7 3 6 2 3 2 6 6 4 4 5 9 7;
    6 0 8 3 2 6 8 4 8 8 13 7 5 8 12 10 14;
    9 8 0 11 10 6 3 9 5 8 4 15 14 13 9 18 9;
    8 3 11 0 1 7 10 6 10 10 14 6 7 9 14 6 16;
    7 2 10 1 0 6 9 4 8 9 13 4 6 8 12 8 14;
    3 6 6 7 6 0 2 3 2 2 7 9 7 7 6 12 8;
    6 8 3 10 9 2 0 6 2 5 4 12 10 10 6 15 5;
    2 4 9 6 4 3 6 0 4 4 8 5 4 3 7 8 10;
    3 8 5 10 8 2 2 4 0 3 4 9 8 7 3 13 6;
    2 8 8 10 9 2 5 4 3 0 4 6 5 4 3 9 5;
    6 13 4 14 13 7 4 8 4 4 0 10 9 8 4 13 4;
    6 7 15 6 4 9 12 5 9 6 10 0 1 3 7 3 10;
    4 5 14 7 6 7 10 4 8 5 9 1 0 2 6 4 8;
    4 8 13 9 8 7 10 3 7 4 8 3 2 0 4 5 6;
    5 12 9 14 12 6 6 7 3 3 4 7 6 4 0 9 2;
    9 10 18 6 8 12 15 8 13 9 13 3 4 5 9 0 9;
    7 14 9 16 14 8 5 10 6 5 4 10 8 6 2 9 0;
]
# In the above `travel_time_data`, the index '1' means depot and '2' means customer #1, and so on.
# We want the index '0' means the depot, and '1' means customer #1.
dists_oa = OffsetArray(dists_data, 0:n_customers, 0:n_customers)

# create a Solomon
solomon = Solomon(nodes, fleet, requests)

# solve
@time routes, total_distance = solveVRP(solomon, dists_oa, pricing_method="pulse")
# If you have given proper coordinates, don't pass `dists` as follows:
# @time routes, total_distance = solveVRP(solomon, pricing_method="pulse")

@show total_distance
@show routes
@test total_distance == 73.0


text:
abstract type AbstractStorageFormulation <: AbstractDeviceFormulation end
struct BookKeeping <: AbstractStorageFormulation end
struct BookKeepingwReservation <: AbstractStorageFormulation end
#################################################Storage Variables#################################

function active_power_variables!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{St},
) where {St <: PSY.Storage}
    add_variable(
        psi_container,
        devices,
        variable_name(ACTIVE_POWER_IN, St),
        false,
        :nodal_balance_active,
        -1.0;
        lb_value = d -> 0.0,
    )
    add_variable(
        psi_container,
        devices,
        variable_name(ACTIVE_POWER_OUT, St),
        false,
        :nodal_balance_active;
        lb_value = d -> 0.0,
    )
    return
end

function reactive_power_variables!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{St},
) where {St <: PSY.Storage}
    add_variable(
        psi_container,
        devices,
        variable_name(REACTIVE_POWER, St),
        false,
        :nodal_balance_reactive,
    )
    return
end

function energy_storage_variables!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{St},
) where {St <: PSY.Storage}
    add_variable(
        psi_container,
        devices,
        variable_name(ENERGY, St),
        false;
        lb_value = d -> 0.0,
    )
    return
end

function storage_reservation_variables!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{St},
) where {St <: PSY.Storage}
    add_variable(psi_container, devices, variable_name(RESERVE, St), true)
    return
end

###################################################### output power constraints#################################

function active_power_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{St},
    model::DeviceModel{St, BookKeeping},
    ::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {St <: PSY.Storage, S <: PM.AbstractPowerModel}
    constraint_data_in = Vector{DeviceRange}(undef, length(devices))
    constraint_data_out = Vector{DeviceRange}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        name = PSY.get_name(d)
        in_lims = PSY.get_inputactivepowerlimits(d)
        out_lims = PSY.get_outputactivepowerlimits(d)
        constraint_data_in[ix] = DeviceRange(name, in_lims)
        constraint_data_out[ix] = DeviceRange(name, out_lims)
        add_device_services!(constraint_data_in[ix], constraint_data_out[ix], d, model)
    end

    device_range(
        psi_container,
        constraint_data_out,
        constraint_name(OUTPUT_POWER_RANGE, St),
        variable_name(ACTIVE_POWER_OUT, St),
    )

    device_range(
        psi_container,
        constraint_data_in,
        constraint_name(INPUT_POWER_RANGE, St),
        variable_name(ACTIVE_POWER_IN, St),
    )
    return
end

function active_power_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{St},
    model::DeviceModel{St, BookKeepingwReservation},
    ::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {St <: PSY.Storage, S <: PM.AbstractPowerModel}
    constraint_data_in = Vector{DeviceRange}(undef, length(devices))
    constraint_data_out = Vector{DeviceRange}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        name = PSY.get_name(d)
        in_lims = PSY.get_inputactivepowerlimits(d)
        out_lims = PSY.get_outputactivepowerlimits(d)
        constraint_data_in[ix] = DeviceRange(name, in_lims)
        constraint_data_out[ix] = DeviceRange(name, out_lims)
        add_device_services!(constraint_data_in[ix], constraint_data_out[ix], d, model)
    end

    reserve_device_semicontinuousrange(
        psi_container,
        constraint_data_in,
        constraint_name(INPUT_POWER_RANGE, St),
        variable_name(ACTIVE_POWER_IN, St),
        variable_name(RESERVE, St),
    )

    reserve_device_semicontinuousrange(
        psi_container,
        constraint_data_out,
        constraint_name(OUTPUT_POWER_RANGE, St),
        variable_name(ACTIVE_POWER_OUT, St),
        variable_name(RESERVE, St),
    )
    return
end

"""
This function adds the reactive  power limits of generators when there are CommitmentVariables
"""
function reactive_power_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{St},
    model::DeviceModel{St, D},
    ::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {St <: PSY.Storage, D <: AbstractStorageFormulation, S <: PM.AbstractPowerModel}
    constraint_data = Vector{DeviceRange}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        name = PSY.get_name(d)
        limits = PSY.get_reactivepowerlimits(d)
        range_data = DeviceRange(name, limits)
        constraint_data[ix] = range_data
    end

    device_range(
        psi_container,
        constraint_data,
        constraint_name(REACTIVE_RANGE, St),
        variable_name(REACTIVE_POWER, St),
    )
    return
end

########################## Make initial Conditions for a Model #############################
function initial_conditions!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{St},
    ::Type{D},
) where {St <: PSY.Storage, D <: AbstractStorageFormulation}
    storage_energy_init(psi_container.initial_conditions, devices)
    return
end

###################################################### Energy Capacity constraints##########

function energy_capacity_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{St},
    model::DeviceModel{St, D},
    ::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {St <: PSY.Storage, D <: AbstractStorageFormulation, S <: PM.AbstractPowerModel}
    constraint_data = Vector{DeviceRange}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        name = PSY.get_name(d)
        limits = PSY.get_capacity(d)
        range_data = DeviceRange(name, limits)
        add_device_services!(range_data, d, model)
        constraint_data[ix] = range_data
    end

    device_range(
        psi_container,
        constraint_data,
        constraint_name(ENERGY_CAPACITY, St),
        variable_name(ENERGY, St),
    )
    return
end

###################################################### book keeping constraints ############

function make_efficiency_data(
    devices::IS.FlattenIteratorWrapper{St},
) where {St <: PSY.Storage}
    names = Vector{String}(undef, length(devices))
    in_out = Vector{InOut}(undef, length(devices))

    for (ix, d) in enumerate(devices)
        names[ix] = PSY.get_name(d)
        in_out[ix] = PSY.get_efficiency(d)
    end

    return names, in_out
end

function energy_balance_constraint!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{St},
    ::Type{D},
    ::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {St <: PSY.Storage, D <: AbstractStorageFormulation, S <: PM.AbstractPowerModel}
    efficiency_data = make_efficiency_data(devices)
    energy_balance(
        psi_container,
        get_initial_conditions(psi_container, ICKey(EnergyLevel, St)),
        efficiency_data,
        constraint_name(ENERGY_LIMIT, St),
        (
            variable_name(ACTIVE_POWER_OUT, St),
            variable_name(ACTIVE_POWER_IN, St),
            variable_name(ENERGY, St),
        ),
    )
    return
end


text:
using MAT
using StatsBase
std = StatsBase.std
trials = 10
include("../include/HyperLocal/HyperLocal.jl")
include("../include/HyperLocal/SparseCardHyperLocal.jl")
include("../src/hypergraph-clustering-utils.jl")
include("load_stackoverflow.jl")

## Load Delta-linear results
Druns = zeros(45,trials)
Df1s = zeros(45,trials)

for trial = 1:trials
    M = matread("Output/deltalinear_45_$trial.mat")
    stats = M["delta_stats"]
    drun = stats[:,6]
    df1 = stats[:,3]
    Druns[:,trial] = drun
    Df1s[:,trial] = df1
end

d_run_mean = mean(Druns,dims = 2)
d_run_std = std(Druns,dims = 2)
d_f1_mean = mean(Df1s,dims = 2)
d_f1_std = std(Df1s,dims = 2)

## Other methods
run_mean = zeros(45,9)
run_std = zeros(45,9)

f1_mean = zeros(45,9)
f1_std = zeros(45,9)

types = ["clique","x9","sqrt"]

epsis = 3
for tt = 1:3
    type = types[tt]

    for jj = 1:epsis
        Druns = zeros(45,trials)
        Df1s = zeros(45,trials)
        for trial = 1:trials
            M = matread("Output/SCHL_$(type)_45_$trial.mat")
            stats = M["sparse$(jj)_stats"]
            drun = stats[:,6]
            df1 = stats[:,3]
            Druns[:,trial] = drun
            Df1s[:,trial] = df1
        end
        # println("$type \t $jj")
        index = 3*(tt-1)+jj
        run_mean[:,index] = mean(Druns,dims = 2)
        run_std[:,index] = std(Druns,dims = 2)
        f1_mean[:,index] = mean(Df1s,dims = 2)
        f1_std[:,index] = std(Df1s,dims = 2)
    end
end

## Mean results

F1s = [d_f1_mean'; f1_mean[:,1]'; f1_mean[:,4]'; f1_mean[:,7]']
BestF1 = zeros(Int64,size(F1s,1))
for i = 1:45
    a1 = F1s[1,i]
    a2 = F1s[2,i]
    a3 = F1s[3,i]
    a4 = F1s[4,i]
    if a1 == a2 || a1 == a3 || a1 == a4 || a2 == a3 || a2 == a4 || a3 == a4
        println("Same results for two methods")
    end
    mx, mxi = findmax(F1s[:,i])
    BestF1[mxi] += 1
end

## Plot mean f1 scores in each case
s1 = 700
s2 = 400
ms = 6
p = sortperm(vec(d_f1_mean))
pl = scatter(d_f1_mean[p],markersize = ms, color = :blue,markerstrokecolor = :blue,legend = :topleft,label="δ-linear" )
scatter!(pl,f1_mean[p,1],markersize = ms, color = :green,markerstrokecolor = :green,label ="clique")
scatter!(pl,f1_mean[p,4],markersize = ms, color = :orange,markerstrokecolor = :orange, label="x^{0.9}")
scatter!(pl,f1_mean[p,7],markersize = ms, color = :purple, markerstrokecolor = :purple,label="sqrt",size = (s1,s2), grid = false)

ln = "  ".*names[p]
stepy = 1
plot!(pl,xticks = (1:stepy:length(ln), ln[1:stepy:length(ln)]),xrotation = 40,
xtickfont=font(8),
ytickfont=font(11),
guidefont=font(12),
titlefont=font(10),
legendfont=font(9))

cl = f1_mean[:,1]
sum(cl .> d_f1_mean)
mean(cl - d_f1_mean)
pl
# savefig("Figures/MeanPlots.pdf")

max_stars_repo_name: nveldt/SparseCardDSFM

max_issues_repo_name: nveldt/SparseCardDSFM

max_forks_repo_name: nveldt/SparseCardDSFM


text:
# Load some data and rows for the train/test sets:
using Koala
X, y = load_ames();
const train, test = split(eachindex(y), 0.8); # 80:20 split

# Instantiate a tree model:
using KoalaTrees
t = TreeRegressor(regularization=0.5)
showall(t)

# Build and train a machine:
tree = Machine(t, X, y, train)
fit!(tree, train)
showall(tree)

# Compute the RMS error on the test set:
err(tree, test)

# Tune the regularization parameter:
u, v = @curve r logspace(-3,2,100) begin
    t.regularization = r
    fit!(tree, train)
    err(tree, test)
end

t.regularization = u[indmin(v)]
fit!(tree, train)
err(tree, test)


text:
export SymmetricGroup, infinity_norm, infinity_distance, fubini_study,
    PrecisionOption, PRECISION_ADAPTIVE, PRECISION_FIXED_64, PRECISION_FIXED_128

@enum ActiveCoeffs begin
    COEFFS_EVAL
    COEFFS_DT
    COEFFS_UNKNOWN
end

"""
    PrecisionOption

Controls the used precision for computing the residual in Newton's method.
See [[Tisseur01]](https://epubs.siam.org/doi/abs/10.1137/S0895479899359837) for the analysis behind this approach.

## Values
* `PRECISION_FIXED_64`: Only use default 64 bit machine precision
* `PRECISION_FIXED_128`: Always use emulated 128 bit floating point numbers. These are provided by the [DoubleFloats.jl](https://github.com/JuliaMath/DoubleFloats.jl) package.
* `PRECISION_ADAPTIVE`: Adaptively switch between 64 and 128 bit floating point numbers.
"""
@enum PrecisionOption begin
    PRECISION_FIXED_64
    PRECISION_FIXED_128
    PRECISION_ADAPTIVE
end

"""
    unpack(a::Union{Nothing, T}, b::T)

Returns `a` if it is not `nothing`, otherwise `b`.
"""
unpack(a::Union{Nothing, T}, b::T) where {T} = a === nothing ? b : a

"""
    SymmetricGroup(n)

Group action of the symmetric group S(n).
"""
struct SymmetricGroup{N}
    permutations::Vector{SVector{N,Int}}
end
SymmetricGroup(N::Int) = SymmetricGroup(permutations(Val(N)))
function permutations(::Val{N}) where {N}
    s = StaticArrays.MVector{N}(1:N)
    perms = [SVector(s)]
    while true
        i = N - 1
        while i>=1 && s[i] >= s[i+1]
            i -= 1
        end
        if i > 0
            j = N
            while j > i && s[i] >= s[j]
                j -= 1
            end
            s[i], s[j] = s[j], s[i]
            reverse!(s, i+1)
        else
            s[1] = N+1
        end

        s[1] > N && break
        push!(perms, SVector(s))
    end
    perms
end

Base.eltype(::Type{SymmetricGroup{N}}) where {N} = SVector{N, Int}
Base.length(p::SymmetricGroup{N}) where {N} = length(p.permutations)
Base.iterate(p::SymmetricGroup) = iterate(p.permutations)
Base.iterate(p::SymmetricGroup, s) = iterate(p.permutations, s)

"""
    deprecatekwarg(oldkw, newkw)

Logs a deprecation warning and assigns `oldkw` to `newkw` if `oldkw` is not `nothing`.
"""
macro deprecatekwarg(oldkw, newkw)
    quote
        if $(esc(oldkw)) !== nothing
            old = $(Expr(:quote, oldkw))
            new = $(Expr(:quote, newkw))
            @warn("`$(old)=$($(esc(oldkw)))` is deprecated, use `$(new)=$($(esc(oldkw)))` instead.")
            $(esc(newkw)) = $(esc(oldkw))
        end
    end
end

"""
    invalid_keywords(kwargs, allowed)

Look for invalid keywords in the `kwargs`.
"""
function invalid_kwargs(kwargs, allowed)
    invalids = []
    for kwarg in kwargs
        kw = first(kwarg)
        if !any(isequal(kw), allowed)
            push!(invalids, kwarg)
        end
    end
    invalids
end

"""
    nthroot(x, n)

Compute the `n`-th root of `x`.
"""
function nthroot(x::Real, N::Integer)
    if N == 3
        cbrt(x)
    elseif N == 2
        √(x)
    elseif N == 4
        √(√(x))
    elseif N == 1
        x
    elseif N == 0
        one(x)
    else
        x^(1/N)
    end
end

"""
    log₁₀(x)

Same as `log10` but faster.
"""
log₁₀(x) = log(x) / 2.302585092994046

"""
    @modulenum(name, block)

This is a modification of `@enum` and creates an intermediate enum.

## Example

The definition

```julia
@moduleenum Car begin
    audi
    volkswagen
    bmw
end
```

expands into

```julia
module Car
    @enum t begin
        audi
        volkswagen
        bmw
    end
end
import .Car
"""
macro moduleenum(name, content...)
    mname = esc(name)
    quote
        @eval module $name
            @enum t $(content...)
        end
        import .$name
    end
end

"""
     print_fieldnames(io::IO, obj)

 A better default printing for structs.
 """
 function print_fieldnames(io::IO, obj)
     println(io, typeof(obj), ":")
     for name in fieldnames(typeof(obj))
         if getfield(obj, name) !== nothing
             println(io, " • ", name, " → ", getfield(obj, name))
         end
     end
 end

"""
    check_kwargs_empty(kwargs, [allowed_kwargs])

Chack that the list of `kwargs` is empty. If not, print all unsupported keywords
with their arguments.
"""
function check_kwargs_empty(kwargs, allowed_kwargs=[])
    if !isempty(kwargs)
        msg = "Unexpected keyword argument(s): "
        first_el = true
        for kwarg in kwargs
            if !first_el
                msg *= ", "
            end
            msg *= "$(first(kwarg))=$(last(kwarg))"
            first_el = false
        end
        if !isempty(allowed_kwargs)
            msg *= "\nAllowed keywords are\n"
            msg *= join(allowed_kwargs, ", ")
        end
        throw(ArgumentError(msg))
    end
end

"""
    is_real_vector(v::AbstractVector, tol=1e-6)

Check whether the 2-norm of the imaginary part of `v` is at most `tol`.
"""
is_real_vector(z::AbstractVector{<:Real}, tol=1e-6) = true
function is_real_vector(z::AbstractVector{<:Complex}, tol=1e-6)
    total = zero(real(eltype(z)))
    for zᵢ in z
        total += abs2(imag(zᵢ))
    end
    sqrt(total) < tol
end
is_real_vector(z::NTuple{N, T}, tol=1e-6) where {N,T} = is_real_vector(SVector{N}(z), tol)

"""
    randseed(range=1_000:1_000_000)

Return a random seed in the range `range`.
"""
randseed(range=1_000:1_000_000) = rand(range)

"""

    infinity_distance(u, v)

Compute the ∞-norm of `u-v`.
"""
function infinity_distance(z₁::AbstractVector, z₂::AbstractVector)
    m = abs(z₁[1] - z₂[1])
    n₁, n₂ = length(z₁), length(z₂)
    if n₁ ≠ n₂
        return convert(typeof(m), Inf)
    end
    @inbounds for k=2:n₁
        @fastmath m = max(m, abs(z₁[k] - z₂[k]))
    end
    m
end

"""

    infinity_norm(z)

Compute the ∞-norm of `z`. If `z` is a complex vector this is more efficient
than `norm(z, Inf)`.
"""
infinity_norm(z::AbstractVector{<:Complex}) = sqrt(maximum(abs2, z))

"""
    fubini_study(x::PVector, y::PVector)

Computes the Fubini-Study distance between `x` and `y`.
"""
function fubini_study(x::PVector{<:Number, 1}, y::PVector{<:Number, 1})
    acos(min(1.0, abs(first(LinearAlgebra.dot(x,y)))))
end
function fubini_study(x::PVector{<:Number, M}, y::PVector{<:Number, M}) where {M}
    sqrt(sum(abs2.(acos.(min.(1.0, abs.(LinearAlgebra.dot(x,y)))))))
end

function randomish_gamma()
    # Usually values near 1, i, -i, -1 are not good randomization
    # Therefore we artificially constrain the choices
    theta = rand() * 0.30 + 0.075 + (rand(Bool) ? 0.0 : 0.5)
    cis(2π * theta)
end

"""
    splitkwargs(kwargs, supported_keywords)

Split the vector of `kwargs` in two vectors, the first contains all `kwargs`
whose keywords appear in `supported_keywords` and the rest the other one.
"""
function splitkwargs(kwargs, supported_keywords)
    supported = []
    rest = []
    for kwarg in kwargs
        if any(kw -> kw == first(kwarg), supported_keywords)
            push!(supported, kwarg)
        else
            push!(rest, kwarg)
        end
    end
    supported, rest
end

start_solution_sample(xs) = promote_start_solution(first(xs))
start_solution_sample(x::AbstractVector{<:Number}) = promote_start_solution(x)
promote_start_solution(x::AbstractVector{ComplexF64}) = x
promote_start_solution(x::SVector) = map(xᵢ -> first(promote(xᵢ, 0.0im)), x)
function promote_start_solution(x)
    x_new =similar(x, promote_type(eltype(x), ComplexF64), length(x))
    copyto!(x_new, x)
    x_new
end

"""
    ComplexSegment(start, target)

Represents the line segment from `start` to `finish`.
Supports indexing of the values `t ∈ [0, length(target-start)]` in order to get the
corresponding point on the line segment.
"""
struct ComplexSegment
    start::ComplexF64
    target::ComplexF64
    # derived
    Δ_target_start::ComplexF64
    abs_target_start::Float64
end
function ComplexSegment(start, target)
    Δ_target_start = convert(ComplexF64, target) - convert(ComplexF64, start)
    abs_target_start = abs(Δ_target_start)

    ComplexSegment(start, target, Δ_target_start, abs_target_start)
end

function Base.getindex(segment::ComplexSegment, t::Real)
    Δ = t / segment.abs_target_start
    if 1.0 - Δ < 2eps()
        Δ = 1.0
    end
    segment.start + Δ * segment.Δ_target_start
end
Base.length(segment::ComplexSegment) = segment.abs_target_start

function Base.show(io::IO, segment::ComplexSegment)
    print(io, "ComplexSegment($(segment.start), $(segment.target))")
end
Base.show(io::IO, ::MIME"application/prs.juno.inline", opts::ComplexSegment) = opts

# Parallelization

"""
    single_thread_blas()

We set the number of BLAS threads to 1 since we multithread by ourself.
But even if we would not, the overhead of the threading is so large
that the single threaded version is around two times faster (at least on Mac Julia v0.6.2).
Returns the previous number of BLAS threads.
"""
function single_thread_blas()
    nblas_threads = get_num_BLAS_threads()
    set_num_BLAS_threads(1)
    nblas_threads
end

set_num_BLAS_threads(n) = LinearAlgebra.BLAS.set_num_threads(n)
get_num_BLAS_threads() = convert(Int, _get_num_BLAS_threads())
# This is into 0.7 but we need it for 0.6 as well
const _get_num_BLAS_threads = function() # anonymous so it will be serialized when called
    blas = LinearAlgebra.BLAS.vendor()
    # Wrap in a try to catch unsupported blas versions
    try
        if blas == :openblas
            return ccall((:openblas_get_num_threads, Base.libblas_name), Cint, ())
        elseif blas == :openblas64
            return ccall((:openblas_get_num_threads64_, Base.libblas_name), Cint, ())
        elseif blas == :mkl
            return ccall((:MKL_Get_Max_Num_Threads, Base.libblas_name), Cint, ())
        end

        # OSX BLAS looks at an environment variable
        if Sys.isapple()
            return ENV["VECLIB_MAXIMUM_THREADS"]
        end
    catch
    end

    return nothing
end

"""

    partition_work(unit_range, n_chunks)

Partition a unit_range in `n_chunks` equally sized chunks.
"""
partition_work(arg, n_chunks) = partition_work!(Vector{UnitRange{Int}}(undef, n_chunks), arg, n_chunks)
function partition_work!(dst, R::UnitRange, n_chunks)
    max_elems_per_chunk = ceil(Int, length(R) / n_chunks)
    a = R.start
    for i in 1:n_chunks
        if i == 1
            a = R.start
        else
            a = dst[i - 1].stop + 1
        end
        b = min(a + max_elems_per_chunk - 1, R.stop)
        dst[i] = a:b
    end
    dst
end

"""
    BatchIterator(iter, batch_size)

An iterator which returns the elements of `iter` in batches of at most `batch_size`.
"""
struct BatchIterator{Iter, T}
    iter::Iter
    batch_size::Int
	batch::Vector{T}
end

function BatchIterator(iter, size::Int)
	batch = Vector{eltype(iter)}(undef, size)
	BatchIterator(iter, size, batch)
end

Base.length(iter::BatchIterator) = ceil(Int, length(iter.iter) / iter.batch_size)
Base.eltype(iter::BatchIterator) = Vector{eltype(iter.iter)}

function Base.iterate(iter::BatchIterator, state=nothing)
    next = state === nothing ? iterate(iter.iter) : iterate(iter.iter, state)
    if next === nothing
        return nothing
    end
    el, state = next
	k = 1
    iter.batch[k] = el
    while k < iter.batch_size
        next = iterate(iter.iter, state)
		if next === nothing
			resize!(iter.batch, k)
			break
		end
        el, state = next
		k += 1
        iter.batch[k] = el
    end

    iter.batch, state
end


text:
using MomentClosure
using DiffEqJump, JLD2
using Test

JLD2.@load joinpath(@__DIR__, "test_data.jld2") sol

μ_t = Dict([(0, 1, 0) => [30.0, 33.6, 37.6]
            (0, 0, 2) => [3600.0, 3461.4, 3207.0]
            (0, 0, 1) => [60.0, 58.8, 56.6]
            (1, 0, 0) => [70.0, 68.4, 71.0]
            (0, 2, 0) => [900.0, 1134.2, 1420.8]
            (1, 0, 1) => [4200.0, 4021.8, 4018.9]
            (0, 1, 1) => [1800.0, 1973.8, 2124.2]
            (2, 0, 0) => [4900.0, 4687.2, 5058.2]
            (1, 1, 0) => [2100.0, 2299.1, 2668.3]])

μ_1 = get_raw_moments(sol, 2)
μ_2 = get_raw_moments(sol, 2, naive=false, b=3)
x_1 = [μ_t[key] ≈ μ_1[key] for key in keys(μ_t)]
x_2 = [μ_t[key] ≈ μ_2[key] for key in keys(μ_t)]
@test all(x_1) && all(x_2)

M_t = Dict([(0, 0, 2) => [0.0, 3.96, 3.44]
            (0, 1, 1) => [0.0, -1.88, -3.96]
            (0, 2, 0) => [0.0, 5.24, 7.04]
            (1, 0, 1) => [0.0, -0.12, 0.3]
            (2, 0, 0) => [0.0, 8.64, 17.2]
            (1, 1, 0) => [0.0, 0.86, -1.3]])

M_1 = get_central_moments(sol, 2)
M_2 = get_central_moments(sol, 2, naive=false, b=3)
x_1 = [M_t[key] ≈ M_1[key] for key in keys(M_t)]
x_2 = [M_t[key] ≈ M_2[key] for key in keys(M_t)]
@test all(x_1) && all(x_2)

κ_t = Dict([(0, 1, 0) => [30.0, 33.6, 37.6]
            (0, 0, 2) => [0.0, 3.96, 3.44]
            (0, 0, 1) => [60.0, 58.8, 56.6]
            (1, 0, 0) => [70.0, 68.4, 71.0]
            (0, 2, 0) => [0.0, 5.24, 7.04]
            (1, 0, 1) => [0.0, -0.12, 0.3]
            (0, 1, 1) => [0.0, -1.88, -3.96]
            (2, 0, 0) => [0.0, 8.64, 17.2]
            (1, 1, 0) => [0.0, 0.86, -1.3]])

κ_1 = get_cumulants(sol, 2)
κ_2 = get_cumulants(sol, 2, naive=false, b=3)
x_1 = [κ_t[key] ≈ κ_1[key] for key in keys(κ_t)]
x_2 = [κ_t[key] ≈ κ_2[key] for key in keys(κ_t)]
@test all(x_1) && all(x_2)


text:
###
  Name     : c9_34_Kmean_randomNumbers.jl
  Book     : Hands-on Data Science with Anaconda )
  Publisher: Packt Publishing Ltd.
  Author   : Yuxing Yan and James Yan
  Date     : 4/6/2018
  email    : yany@canisius.edu
             paulyxy@hotmail.com
###

using Clustering
srand(1234)
nRow=5
nCol=1000
x = rand(nRow,nCol)
maxInter=200  #max interation
nCluster=20
R = kmeans(x,nCluster;maxiter=maxInter,display=:iter)
@assert nclusters(R) ==nCluster
c = counts(R)
clusters= R.centers


text:
using Documenter, CmdStan

makedocs(
    modules = [CmdStan],
    format = :html,
    sitename = "StanJulia/CmdStan.jl",
    pages = Any[
        "Introduction" => "INTRO.md",
        "Installation" => "INSTALLATION.md",
        "Walkthrough" => "WALKTHROUGH.md",
        "Versions" => "VERSIONS.md",
        "Home" => "index.md",
    ]
)

deploydocs(
    repo = "github.com/StanJulia/CmdStan.jl.git",
    target = "build",
    julia = "1.0",
    osname = "linux",
    deps = nothing,
    make = nothing
 )


text:
# Overloadable method, e.g to insert OffsetPrecalc's precalculated stride multiples
@inline tdot(ptr::AbstractStridedPointer, ::Tuple{}, ::Tuple{}) = (pointer(ptr), Zero())
@inline tdot(ptr::AbstractStridedPointer{T}, a, b) where {T} = tdot(pointer(ptr), a, b)

@inline function tdot(p::Ptr{T}, a::Tuple{A}, b::Tuple{B,Vararg}) where {T,A,B}
  p, lazymul(first(a), first(b))
end
@inline function tdot(
  p::Ptr{T},
  a::Tuple{A},
  b::Tuple{B,Vararg},
  c::Tuple{C,Vararg},
) where {T,A,B,C}
  p, lazymul(first(a), first(b), first(c))
end

@inline function tdot(
  p::Ptr{T},
  a::Tuple{A1,A2,Vararg},
  b::Tuple{B1,B2,Vararg},
) where {T,A1,A2,B1,B2}
  i = lazymul(first(a), first(b))
  p, j = tdot(p, tail(a), tail(b))
  add_indices(p, i, j)
end
@inline function tdot(
  p::Ptr{T},
  a::Tuple{A1,A2,Vararg},
  b::Tuple{B1,B2,Vararg},
  c::Tuple{C1,C2,Vararg},
) where {T,A1,A2,B1,B2,C1,C2}
  i = lazymul(first(a), first(b), first(c))
  p, j = tdot(p, tail(a), tail(b), tail(c))
  add_indices(p, i, j)
end

@inline function tdot(
  p::Ptr{T},
  a::Tuple{A1,A2,Vararg},
  b::Tuple{B1,B2,Vararg},
  c::Tuple{C1},
) where {T,A1,A2,B1,B2,C1}
  i = lazymul(first(a), first(b), first(c))
  p, j = tdot(p, tail(a), tail(b))
  add_indices(p, i, j)
end


text:
@doc raw"""
    rigid_body_velocity(vel_P, ω, r_PQ)

Calculate rigid solid velocity field.

Return velocity of a point Q of a rigid solid given the velocity of a point P (vel_P), the
rotational velocity of the solid (ω) and the relative position of Q with respect to P.

If the reference frame 1 is attached to the solid and the velocity is calculated with
respect to reference frame 0:

``v_{10}^{Q} = v_{10}^{P} + \omega_{10} \times r^{PQ}``

being:
- ``v_{10}^{Q}`` the velocity of point Q, fixed to 1, wrt 0
- ``\omega_{10}`` the angular velocity of the solid 1 wrt 0
- ``r^{PQ}`` the position of Q wrt P (``r^{Q}-r^{P}``)

Every vector needs to be expressed in the same coordinate system.

# References

1. Stevens, B. L., Lewis, F. L., (1992). Aircraft control and simulation: dynamics, controls
   design, and autonomous systems. John Wiley & Sons. (Section 1.3, page 26)
"""
function rigid_body_velocity(vel_P, ω, r_PQ)
    vel_Q = vel_P + ω × r_PQ
    return vel_Q
end

@doc raw"""
    rigid_body_acceleration(acc_P, ω, ω_dot, r_PQ)

Calcualte rigid body acceleration field.

Return the acceleration of a point Q of a rigid solid given the acceleration of a point P
(acc_P), the rotational velocity of the solid (ω), the rotational acceleration of the solid
(ω_dot) and the relative position of Q with respect to P.

``a_{10}^{Q} = a_{10}^{P} + \omega_{10} \times (\omega_{10} \times r^{PQ}) + \dot{\omega}_{10} \times r^{PQ}``

being:
- ``a_{10}^{Q}`` the acceleration of point Q, fixed to 1, wrt 0
- ``\omega_{10}`` the angular velocity of the solid 1 wrt 0
- ``\dot{\omega}_{10}`` the angular acceleration of the solid 1 wrt 0
- ``r^{PQ}`` the position of Q wrt P (``r^{Q}-r^{P}``)

# References

1. Stevens, B. L., Lewis, F. L., (1992). Aircraft control and simulation: dynamics, controls
   design, and autonomous systems. John Wiley & Sons. (Section 1.3, Formaula 1.3-14c, page 26)
"""
function rigid_body_acceleration(acc_P, ω, ω_dot, r_PQ)
    acc_Q = acc_P + ω × (ω × r_PQ) + ω_dot × r_PQ
    return acc_Q
end

"""
    pqr_2_ψθϕ_dot(p, q, r, ψ, θ, ϕ)

Transform body angular velocity (p, q, r) [rad/s] to Euler angles rates (ψ_dot, θ_dot,
ϕ_dot) [rad/s] given the euler angles (θ, ϕ) [rad] using kinematic angular equations.

# See also

[`ψθϕ_dot_2_pqr`](@ref), [`pqr_2_quat_dot`](@ref)

# References

1. Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control and simulation:
   dynamics, controls design, and autonomous systems. John Wiley & Sons. Equation (1.4-4) (page 20)
"""
function pqr_2_ψθϕ_dot(p, q, r, θ, ϕ)

    sθ, cθ = sin(θ), cos(θ)
    sϕ, cϕ = sin(ϕ), cos(ϕ)

    ψ_dot = (q * sϕ + r * cϕ) / cθ
    θ_dot = q * cϕ - r * sϕ
    # ϕ_dot = p + (q * sϕ + r * cϕ) * tan(θ)
    ϕ_dot = p + ψ_dot * sθ

    return [ψ_dot, θ_dot, ϕ_dot]
end

"""
    ψθϕ_dot_2_pqr(ψ_dot, θ_dot, ϕ_dot, ψ, θ, ϕ)

Transform Euler angles rates (ψ_dot, θ_dot, ϕ_dot) [rad/s] to body angular velocity (p, q,
r) [rad/s] given the euler angles (θ, ϕ) [rad] using kinematic angular equations.

# See also

[`pqr_2_ψθϕ_dot`](@ref)

# References

1. Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control and simulation:
   dynamics, controls design, and autonomous systems. John Wiley & Sons. Equation (1.4-3) (page 20)
"""
function ψθϕ_dot_2_pqr(ψ_dot, θ_dot, ϕ_dot, θ, ϕ)

    sθ, cθ = sin(θ), cos(θ)
    sϕ, cϕ = sin(ϕ), cos(ϕ)

    p = ϕ_dot - sθ * ψ_dot
    q = cϕ * θ_dot + sϕ*cθ * ψ_dot
    r = -sϕ * θ_dot + cϕ*cθ * ψ_dot

    return [p, q, r]
end

"""
    pqr_2_quat_dot(p, q, r, q0, q1, q2, q3)

Transform body angular velocity (p, q, r) [rad/s] to quaternion rates [1/s].

# See also

[`pqr_2_ψθϕ_dot`](@ref)

# References

1. Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control and simulation:
   dynamics, controls design, and autonomous systems. John Wiley & Sons. Equation (1.8-15)
   (page 51).
"""
function pqr_2_quat_dot(p, q, r, q0, q1, q2, q3)

    Ω = [
        0 -p -q -r;
        p  0  r -q;
        q -r  0  p;
        r  q -p  0;
    ]
    q = [q0; q1; q2; q3]
    q_dot = 0.5 * Ω * q

    return q_dot
end

"""
    uvw_to_tasαβ(u, v, w)

Calculate true air speed (TAS), angle of attack (α) and angle of side slip (β) from velocity
expressed in body axis.

# Notes

This function assumes that u, v, w are the body components of the aerodynamic speed. This is
not true in genreal (wind speed different from zero), as u, v, w represent velocity with
respect to an inertial reference frame.

# References

1. Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control and simulation:
   dynamics, controls design, and autonomous systems. John Wiley & Sons. Equation (2.3-6b)
   (page 78).
"""
function uvw_to_tasαβ(u, v, w)

    tas = sqrt(u*u + v*v + w*w)
    α = atan(w, u)
    β = asin(v / tas)
    return [tas, α, β]
end

"""
    uvw_dot_to_tasαβ_dot(u, v, w, u_dot, v_dot, w_dot)

Calculate time derivatives of velocity expressed as TAS, AOA, AOS.

# Notes
Note that tas here is not necessarily true air speed. Could also be inertial speed in the
direction of airspeed. It will concide whith TAS for no wind.

# See also

[`tasαβ_dot_to_uvw_dot`](@ref)

# References
1. Morelli, Eugene A., and Vladislav Klein. Aircraft system identification: Theory and
   practice. Williamsburg, VA: Sunflyte Enterprises, 2016. Equation 3.33 (page 44).
2. Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control and simulation:
   dynamics, controls design, and autonomous systems. John Wiley & Sons. Equation (2.3-10)
   (page 81).
"""
function uvw_dot_to_tasαβ_dot(u, v, w, u_dot, v_dot, w_dot)

    # [1] 3.31
    tas = sqrt(u*u + v*v + w*w)
    # [1] 3.33a
    tas_dot = (u * u_dot + v * v_dot + w * w_dot) / tas
    # [1] 3.33b
    β_dot = ((u*u + w*w) * v_dot - v * (u * u_dot + w * w_dot)) / (tas^2 * sqrt(u*u + w*w))
    # [2] 2.3.10b
    # β_dot = (v_dot * tas - v * tas_dot) / (tas * sqrt(u*u + w*w))
    # [1] 3.33c
    α_dot = (w_dot * u - w * u_dot) / (u*u + w*w)

    return [tas_dot, α_dot, β_dot]
end

"""
    tasαβ_dot_to_uvw_dot(tas, α, β, tas_dot, α_dot, β_dot)

Obatain body velocity derivatives given velocity in wind axis and its derivatives.

# See also

[`uvw_dot_to_tasαβ_dot`](@ref)

# References

1. Morelli, Eugene A., and Vladislav Klein. Aircraft system identification: Theory and
   practice. Williamsburg, VA: Sunflyte Enterprises, 2016. Derived from equation 3.32
   (page 44).
"""
function tasαβ_dot_to_uvw_dot(tas, α, β, tas_dot, α_dot, β_dot)
    u_dot = tas_dot * cos(α) * cos(β) - tas * (α_dot * sin(α) * cos(β) + β_dot * cos(α) * sin(β))
    v_dot = tas_dot * sin(β) + tas * β_dot * cos(β)
    w_dot = tas_dot * sin(α) * cos(β) + tas * (α_dot * cos(α) * cos(β) - β_dot * sin(α) * sin(β))
    return [u_dot, v_dot, w_dot]
end

@doc raw"""
    rate_of_climb_constrain_no_wind(γ, α, β, ϕ)

Calculate pitch angle (θ rad) to obtain a flight path angle (γ rad) at certain angle of
attack (α rad), angle of sideslip (β rad) and roll (ϕ rad).

Inertial velocity and aerodynamic velocity are equivalent in the absence of wind:

`` v_{cg-e}^e = R_{eb} R_{bw} v_{cg-air}^w ``

# References

1. Stevens, B. L., Lewis, F. L., (1992). Aircraft control and simulation: dynamics, controls
   design, and autonomous systems. John Wiley & Sons. (Section 3.6, equation 3.6-3,
   page 187)
"""
function rate_of_climb_constrain_no_wind(γ, α, β, ϕ)
    a = cos(α) * cos(β)
    b = sin(ϕ) * sin(β) + cos(ϕ) * sin(α) * cos(β)
    sq = sqrt(a^2 - sin(γ)^2 + b^2)
    θ = (a * b + sin(γ) * sq) / (a^2 - sin(γ)^2)
    θ = atan(θ)
    return θ
end