text:
using Test
using LinearAlgebra, SparseArrays
using SummationByPartsOperators, BandedMatrices

D_test_list = (MattssonNordström2004(), MattssonSvärdNordström2004(),
                MattssonSvärdShoeybi2008(), Mattsson2014(),
                MattssonAlmquistCarpenter2014Extended(),
                MattssonAlmquistCarpenter2014Optimal())
Di_test_list = (MattssonSvärdNordström2004(),)
D2var_test_list = (Mattsson2012(),)

for T in (Float32, Float64), acc_order in (2,4,6,8), diss_order in (2,4,6,8), D_source in D_test_list, Di_source in Di_test_list
    xmin = zero(T)
    xmax = 5*one(T)
    N = 101
    D_serial = try
        derivative_operator(D_source, 1, acc_order, xmin, xmax, N, Val{:serial}())
    catch err
        !isa(err, ArgumentError) && throw(err)
        nothing
    end
    D_serial == nothing && continue

    D_full   = Matrix(D_serial)
    D_sparse = sparse(D_serial)
    D_banded = BandedMatrix(D_serial)
    x = grid(D_serial)
    u = cospi.(x)
    dest1 = fill(zero(eltype(u)), length(u))
    dest2 = fill(zero(eltype(u)), length(u))

    @test BandedMatrices.isbanded(D_serial) == BandedMatrices.isbanded(D_banded)
    @test bandwidth(D_serial, 1) == bandwidth(D_banded, 1)
    @test bandwidth(D_serial, 2) == bandwidth(D_banded, 2)

    mul!(dest1, D_serial, u)
    mul!(dest2, D_full, u)
    @test all(i->isapprox(dest1[i], dest2[i], atol=5000*eps(T)), eachindex(u))
    mul!(dest2, D_sparse, u)
    @test all(i->isapprox(dest1[i], dest2[i], atol=5000*eps(T)), eachindex(u))
    mul!(dest2, D_banded, u)
    @test all(i->isapprox(dest1[i], dest2[i], atol=5000*eps(T)), eachindex(u))

    Di_serial = try
        dissipation_operator(Di_source, D_serial, order=diss_order)
    catch err
        !isa(err, ArgumentError) && throw(err)
        nothing
    end
    Di_serial == nothing && continue

    Di_full   = Matrix(Di_serial)
    Di_sparse = sparse(Di_serial)
    Di_banded = BandedMatrix(Di_serial)

    @test BandedMatrices.isbanded(Di_serial) == BandedMatrices.isbanded(Di_banded)
    @test bandwidth(Di_serial, 1) == bandwidth(Di_banded, 1)
    @test bandwidth(Di_serial, 2) == bandwidth(Di_banded, 2)

    mul!(dest1, Di_serial, u)
    mul!(dest2, Di_full, u)
    @test all(i->isapprox(dest1[i], dest2[i], atol=500*eps(T)), eachindex(u))
    mul!(dest2, Di_sparse, u)
    @test all(i->isapprox(dest1[i], dest2[i], atol=500*eps(T)), eachindex(u))
    mul!(dest2, Di_banded, u)
    @test all(i->isapprox(dest1[i], dest2[i], atol=500*eps(T)), eachindex(u))
end

for T in (Float32, Float64), acc_order in (2,4,6), D2var_source in D2var_test_list
    xmin = zero(T)
    xmax = 5*one(T)
    N = 101
    D2var_serial = try
        var_coef_derivative_operator(D2var_source, 2, acc_order, xmin, xmax, N, one, Val{:serial}())
    catch err
        !isa(err, ArgumentError) && throw(err)
        nothing
    end
    D2var_serial == nothing && continue

    D2var_full   = Matrix(D2var_serial)
    D2var_sparse = sparse(D2var_serial)
    D2var_banded = BandedMatrix(D2var_serial)
    x = grid(D2var_serial)
    u = cospi.(x)
    dest1 = fill(zero(eltype(u)), length(u))
    dest2 = fill(zero(eltype(u)), length(u))

    @test BandedMatrices.isbanded(D2var_serial) == BandedMatrices.isbanded(D2var_banded)
    @test bandwidth(D2var_serial, 1) == bandwidth(D2var_banded, 1)
    @test bandwidth(D2var_serial, 2) == bandwidth(D2var_banded, 2)

    mul!(dest1, D2var_serial, u)
    mul!(dest2, D2var_full, u)
    @test all(i->isapprox(dest1[i], dest2[i], atol=5000*eps(T)), eachindex(u))
    mul!(dest2, D2var_sparse, u)
    @test all(i->isapprox(dest1[i], dest2[i], atol=5000*eps(T)), eachindex(u))
    mul!(dest2, D2var_banded, u)
    @test all(i->isapprox(dest1[i], dest2[i], atol=5000*eps(T)), eachindex(u))
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""PolicyRule holds information that describes a policy rule, but does not contain information about who the rule applies to or which namespace the rule applies to.

    IoK8sApiRbacV1beta1PolicyRule(;
        apiGroups=nothing,
        nonResourceURLs=nothing,
        resourceNames=nothing,
        resources=nothing,
        verbs=nothing,
    )

    - apiGroups::Vector{String} : APIGroups is the name of the APIGroup that contains the resources.  If multiple API groups are specified, any action requested against one of the enumerated resources in any API group will be allowed.
    - nonResourceURLs::Vector{String} : NonResourceURLs is a set of partial urls that a user should have access to.  *s are allowed, but only as the full, final step in the path Since non-resource URLs are not namespaced, this field is only applicable for ClusterRoles referenced from a ClusterRoleBinding. Rules can either apply to API resources (such as \"pods\" or \"secrets\") or non-resource URL paths (such as \"/api\"),  but not both.
    - resourceNames::Vector{String} : ResourceNames is an optional white list of names that the rule applies to.  An empty set means that everything is allowed.
    - resources::Vector{String} : Resources is a list of resources this rule applies to.  '*' represents all resources in the specified apiGroups. '*/foo' represents the subresource 'foo' for all resources in the specified apiGroups.
    - verbs::Vector{String} : Verbs is a list of Verbs that apply to ALL the ResourceKinds and AttributeRestrictions contained in this rule.  VerbAll represents all kinds.
"""
mutable struct IoK8sApiRbacV1beta1PolicyRule <: SwaggerModel
    apiGroups::Any # spec type: Union{ Nothing, Vector{String} } # spec name: apiGroups
    nonResourceURLs::Any # spec type: Union{ Nothing, Vector{String} } # spec name: nonResourceURLs
    resourceNames::Any # spec type: Union{ Nothing, Vector{String} } # spec name: resourceNames
    resources::Any # spec type: Union{ Nothing, Vector{String} } # spec name: resources
    verbs::Any # spec type: Union{ Nothing, Vector{String} } # spec name: verbs

    function IoK8sApiRbacV1beta1PolicyRule(;apiGroups=nothing, nonResourceURLs=nothing, resourceNames=nothing, resources=nothing, verbs=nothing)
        o = new()
        validate_property(IoK8sApiRbacV1beta1PolicyRule, Symbol("apiGroups"), apiGroups)
        setfield!(o, Symbol("apiGroups"), apiGroups)
        validate_property(IoK8sApiRbacV1beta1PolicyRule, Symbol("nonResourceURLs"), nonResourceURLs)
        setfield!(o, Symbol("nonResourceURLs"), nonResourceURLs)
        validate_property(IoK8sApiRbacV1beta1PolicyRule, Symbol("resourceNames"), resourceNames)
        setfield!(o, Symbol("resourceNames"), resourceNames)
        validate_property(IoK8sApiRbacV1beta1PolicyRule, Symbol("resources"), resources)
        setfield!(o, Symbol("resources"), resources)
        validate_property(IoK8sApiRbacV1beta1PolicyRule, Symbol("verbs"), verbs)
        setfield!(o, Symbol("verbs"), verbs)
        o
    end
end # type IoK8sApiRbacV1beta1PolicyRule

const _property_map_IoK8sApiRbacV1beta1PolicyRule = Dict{Symbol,Symbol}(Symbol("apiGroups")=>Symbol("apiGroups"), Symbol("nonResourceURLs")=>Symbol("nonResourceURLs"), Symbol("resourceNames")=>Symbol("resourceNames"), Symbol("resources")=>Symbol("resources"), Symbol("verbs")=>Symbol("verbs"))
const _property_types_IoK8sApiRbacV1beta1PolicyRule = Dict{Symbol,String}(Symbol("apiGroups")=>"Vector{String}", Symbol("nonResourceURLs")=>"Vector{String}", Symbol("resourceNames")=>"Vector{String}", Symbol("resources")=>"Vector{String}", Symbol("verbs")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiRbacV1beta1PolicyRule }) = collect(keys(_property_map_IoK8sApiRbacV1beta1PolicyRule))
Swagger.property_type(::Type{ IoK8sApiRbacV1beta1PolicyRule }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiRbacV1beta1PolicyRule[name]))}
Swagger.field_name(::Type{ IoK8sApiRbacV1beta1PolicyRule }, property_name::Symbol) =  _property_map_IoK8sApiRbacV1beta1PolicyRule[property_name]

function check_required(o::IoK8sApiRbacV1beta1PolicyRule)
    (getproperty(o, Symbol("verbs")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiRbacV1beta1PolicyRule }, name::Symbol, val)
end

ext: jl


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modeling language for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# src/norms.jl
# Defines all types relating to norms
# - GenericNorm                 Lp norm of vector of GenericAffExprs
#   - Norm                      Lp norm of vector of AffExprs
# - GenericNormExpr             Expression of norms and affine expressions
#   - GenericSOCExpr            Alias for P=2
#       - SOCExpr               Alias for P=2 and AffExprs
# - GenericSOCConstraint        Constraint of form α||Ax-b||₂ + cᵀx + d ≤ 0
#   - SOCConstraint             Alias for AffExprs
# Operator overloads in src/operators.jl
#############################################################################

#############################################################################
# GenericNorm
# Container for Lp norms, ‖Ax‖p
type GenericNorm{P,C,V}
    terms::Vector{GenericAffExpr{C,V}}
end
# Preferred constructor, validates the norm type based on the expression
function GenericNorm{C,V}(P, terms::Vector{GenericAffExpr{C,V}})
    if C == Float64 && V == Variable
        # JuMP doesn't support anything else than the L2 Norm
        # Throw an error now before going any further
        P ==   1 && error("JuMP doesn't support L₁ norms.")
        P == Inf && error("JuMP doesn't support L∞ norms.")
        P !=   2 && error("JuMP only supports L₂ norms.")
    end
    GenericNorm{P,C,V}(terms)
end
Base.copy{P,C,V}(x::GenericNorm{P,C,V}) = GenericNorm{P,C,V}(copy(x.terms))

# Handle the norm() function by flattening arguments into a vector
Base.norm{V<:AbstractJuMPScalar}(x::V,           p::Real=2) = vecnorm(x,p)
Base.norm{V<:AbstractJuMPScalar}(x::Array{V},    p::Real=2) = vecnorm(x,p)
Base.norm{V<:AbstractJuMPScalar}(x::JuMPArray{V},p::Real=2) = vecnorm(x,p)
Base.norm{V<:AbstractJuMPScalar}(x::JuMPDict{V}, p::Real=2) = vecnorm(x,p)
Base.norm{C,V}(x::GenericAffExpr{C,V},           p::Real=2) = vecnorm(x,p)
Base.norm{C,V}(x::Array{GenericAffExpr{C,V}},    p::Real=2) = vecnorm(x,p)
Base.norm{C,V}(x::JuMPArray{GenericAffExpr{C,V}},p::Real=2) = vecnorm(x,p)
Base.norm{C,V}(x::JuMPDict{GenericAffExpr{C,V}}, p::Real=2) = vecnorm(x,p)

_vecaff(C,V,x) = map(GenericAffExpr{C,V},vec(x))
Base.vecnorm{V<:AbstractJuMPScalar}(x::V,           p::Real=2) = GenericNorm(p, [GenericAffExpr{Float64,V}(x)] )
Base.vecnorm{V<:AbstractJuMPScalar}(x::Array{V},    p::Real=2) = GenericNorm(p, _vecaff(Float64,V,x) )
Base.vecnorm{V<:AbstractJuMPScalar}(x::JuMPArray{V},p::Real=2) = GenericNorm(p, _vecaff(Float64,V,x.innerArray) )
Base.vecnorm{V<:AbstractJuMPScalar}(x::JuMPDict{V}, p::Real=2) = GenericNorm(p, _vecaff(Float64,V,collect(values(x))) )
Base.vecnorm{C,V}(x::GenericAffExpr{C,V},           p::Real=2) = GenericNorm(p, [x])
Base.vecnorm{C,V}(x::Array{GenericAffExpr{C,V}},    p::Real=2) = GenericNorm(p, vec(x))
Base.vecnorm{C,V}(x::JuMPArray{GenericAffExpr{C,V}},p::Real=2) = GenericNorm(p, vec(x.innerArray))
Base.vecnorm{C,V}(x::JuMPDict{GenericAffExpr{C,V}}, p::Real=2) = GenericNorm(p, collect(values(x)))

# Called by the parseNorm macro for e.g. norm2{...}
# If the arguments are tightly typed, just pass to the constructor
_build_norm{C,V}(P,terms::Vector{GenericAffExpr{C,V}}) = GenericNorm(P,terms)
# The terms vector produced by the macro may not be tightly typed,
# so we need to repackage in a tighter typed vector ourselves.
# This function is needed for performance reasons on only 0.3, as
# the following works just as well on 0.4:
_build_norm(Lp, terms::Vector{GenericAffExpr}) = _build_norm(Lp, [terms...])

# Alias for AffExprs. Short-hand used in operator overloads, etc.
typealias Norm{P} GenericNorm{P,Float64,Variable}

getvalue{P,C,V}(n::GenericNorm{P,C,V}) = norm(getvalue(n.terms),P)

##########################################################################
# GenericNormExpr
# Container for expressions containing GenericNorm and GenericAffExpr
type GenericNormExpr{P,C,V}
    norm::GenericNorm{P,C,V}
    coeff::C
    aff::GenericAffExpr{C,V}
end

GenericNormExpr{P,C,V}(norm::GenericNorm{P,C,V}) =
    GenericNormExpr{C,V}(norm, one(C), zero(GenericAffExpr{C,V}))
Base.copy{P,C,V}(x::GenericNormExpr{P,C,V}) =
    GenericNormExpr{P,C,V}(copy(x.norm), copy(x.coeff), copy(x.aff))
Base.convert{P,C,V}(::Type{GenericNormExpr{P,C,V}}, x::GenericNorm{P,C,V}) =
    GenericNormExpr{P,C,V}(x, one(C), zero(GenericAffExpr{C,V}))

# Alias for ‖Ax‖₂ case
typealias GenericSOCExpr{C,V} GenericNormExpr{2,C,V}

# Alias for ‖Ax‖₂ and AffExpr case
typealias SOCExpr GenericSOCExpr{Float64,Variable}

getvalue(n::GenericNormExpr) = n.coeff * getvalue(n.norm) + getvalue(n.aff)

##########################################################################
# GenericSOCConstraint
# Second-order cone constraint of form
# α||Ax-b||₂ + cᵀx + d ≤ 0
type GenericSOCConstraint{T<:GenericSOCExpr} <: AbstractConstraint
    normexpr::T
    function GenericSOCConstraint{T}(normexpr::T)
        if normexpr.coeff < 0
            # The coefficient in front of the norm is negative, which
            # means we have `norm >= c`, which is not convex.
            error("Invalid second-order cone constraint $(normexpr) ≤ 0")
        end
        new(normexpr)
    end
end

# Alias for the AffExpr case
typealias SOCConstraint GenericSOCConstraint{SOCExpr}

function addconstraint(m::Model, c::SOCConstraint)
    push!(m.socconstr,c)
    m.internalModelLoaded = false
    ConstraintRef{Model,SOCConstraint}(m,length(m.socconstr))
end


text:
function recpow(X, n)
	#%UNTITLED2 Summary of this function goes here
	#%   Detailed explanation goes here

	if round(n)!= n
	    error(string("The power has to be an integer", n,"\n"));
	end

	if n==Inf || isnan(n)
		Y = X.^n;
		return Y
	end

	if n==0
		Y = ones(size(X));
	else
		if n==1
			Y = X;
		else
			if mod(n,2) == 0
			Z = recpow(X, n/2);
			Y = Z .* Z;
			else
				Z = recpow(X, (n-1)/2);
				Y = Z .* Z .* X;
			end
		end
	end
	Y
end

max_stars_repo_name: barzan/dbseer

max_issues_repo_name: barzan/dbseer

max_forks_repo_name: barzan/dbseer


text:
const _show_pad = 12
const _show_signif = 5

function Base.show(io::IO, bs::BootstrapSample)
    r = _show_pad
    s = _show_signif
    df = estimate_summary(bs)
    df = split(string(df), '\n')[2:end]
    df[1] = replace(df[1], "Row", "Var")
    println(io, "Bootstrap Sampling")
    println(io, rpad("  Estimates:", r))
    for i in 1:length(df)
        println(io, string("    ", df[i]))
    end
    println(io, rpad("  Sampling:", r),
            replace(string(typeof(sampling(bs))), "Bootstrap.", ""))
    println(io, rpad("  Samples:", r), nrun(bs))
    println(io, rpad("  Data:", r), data_summary(bs.data))
end

"""
Summarizes an object in the form of 'TYPE: { SIZE }'
"""
function data_summary(x)
    t = string(typeof(x))
    n = length(size(x)) > 0 ? size(x) : length(x)
    ns = join(n, " \u00D7 ") # Unicode 'times'
    s = "$t: { $ns }"
    return s
end

function estimate_summary(bs::BootstrapSample)
    n = nvar(bs)
    df = DataFrame(Estimate = [original(bs)...],
                   Bias = bias(bs),
                   StdError = se(bs))
    return df
end


text:
"""
The normal form of transcritical bifurcations.

https://en.wikipedia.org/wiki/Transcritical_bifurcation
"""
module Transcritical

using DiffEqBase: ODEProblem
using Setfield: @lens

using ...Bifurcations: BifurcationProblem

f(u::Real, p, t) = p * u - u^2

u0 = 0.0
tspan = (0.0, 1.0)
p = -1.0
ode = ODEProblem(f, u0, tspan, p)

param_axis = @lens _
prob = BifurcationProblem(ode, param_axis, (-1.0, 1.0))

using StaticArrays: @SVector

function closest_analytic(u, p)
    analytic = @SVector [0, p]
    _, i = findmin(abs.(u .- analytic))
    return [analytic[i], p]
end

end  # module


text:
using Unitful

@testset "utils" begin
    r = toHz(12.3)
    @test typeof(r) == Float64
    r = toHz(12.3u"kHz")
    @test typeof(r) == Float64
    @test r == 12300

    r = tos(12.3)
    @test typeof(r) == Float64
    r = tos(12.3u"s")
    @test typeof(r) == Float64
    @test r == 12.3
end


text:
export one_step_size

"""
Returns 1.0 step size
"""
function one_step_size(h  :: LineModel,  stop_ls :: LS_Stopping;
                       f_meta = LS_Function_Meta(),
                       φ_dφ   :: Function = (x, y) -> phi_dphi(x, y),
                       verboseLS :: Bool = false, kwargs...)
    state = stop_ls.current_state
    update!(state, x = 1.0)

    return (state, true)
end # function


text:
# FileIO compatible load function
function load(rnxfile::File{format"RNX"})
    open(rnxfile) do rnxstream
        iostream = stream(rnxstream) #convert to regular ::IOStream
        rinex_header = rinex_read(iostream)
    end
end


text:
@testset "https://github.com/aviatesk/JET.jl/issues/71" begin
    stats = @benchmark_freshexec ntimes = 1 @report_call println(QuoteNode(nothing))
    @test stats.time ≤ (get(ENV, "CI", nothing) == "true" ? 20 : 10) # a CI runner might be slow
end


text:
"""
```
eval(slm::SLM{T}, x::AbstractVector{T}, evalmode::Int = 0) where {T <: Real}
```

evaluates the spline prescribed by `slm` at the values `x`. The optional input
`evalmode` specifies the evaluation of the spline according to

0  -> spline's value at each point in `x`
1  -> first derivative of the spline
2  -> second derivative of the spline
3  -> third derivative of the spline

Note that the inverse of points which fall above the maximum
or below the minimum value of the function will be returned as a NaN.

Currently, no constant extrapolation is supported.
"""
function eval(slm::SLM{T}, x::AbstractVector{T}, evalmode::Int = 0) where {T <: Real}

    ## Set up
    nₓ     = length(x)
    knots  = get_knots(slm)
    nk     = length(knots)
    dknots = diff(knots)
    coef   = get_coef(slm)

    # Extrapolation
    if slm[:extrap] != :none
        error("Extrapolation method $(slm[:extrap]) is not supported.")
    end

    ## Evaluate spline

    # Sort points into bins. Inverse case is handled later -> just zeros for now.
    x_bins = evalmode >= 0 ? bin_sort(x, knots) : zeros(Int, nₓ)

    if get_type(slm) == :cubic
        y = if evalmode == 0
            eval_cubic(x, knots, dknots, coef, x_bins)
        elseif evalmode == 1
            eval_cubic_derivative_1(x, knots, dknots, coef, x_bins)
        elseif evalmode == 2
            eval_cubic_derivative_2(x, knots, dknots, coef, x_bins)
        elseif evalmode == 3
            eval_cubic_derivative_3(x, knots, dknots, coef, x_bins)
        elseif evalmode == -1
            error("The function inverse has not been implemented yet.")
        else
            error("The input `evalmode` cannot be $(evalmode). It must be one of [0, 1, 2, 3].")
        end
    else
        error("Cannot use the spline type $(get_type(slm))")
    end

    return y
end

function eval_cubic(x::AbstractVector{T}, knots::AbstractVector{T}, dknots::AbstractVector{T},
                    coef::AbstractMatrix{T}, x_bins::AbstractVector{Int}) where {T <: Real}

    # Set up
    t  = (x - knots[x_bins]) ./ dknots[x_bins]
    t² = t .^ 2
    t³ = t .^ 3
    s² = (1. .- t) .^ 2
    s³ = (1. .- t) .^ 3

    # No extrapolation yet so directly return
    x_bins_up = x_bins .+ 1

    return (coef[x_bins, 2] .* (s² - s³) +
            coef[x_bins_up, 2] .* (t³ - t²)) .* dknots[x_bins] +
            coef[x_bins, 1] .* (3 .* s² - 2. .* s³) +
            coef[x_bins_up, 1] .* (3. .* t² - 2. .* t³)

end

function eval_cubic_derivative_1(x::AbstractVector{T}, knots::AbstractVector{T}, dknots::AbstractVector{T},
                                 coef::AbstractMatrix{T}, x_bins::AbstractVector{Int}) where {T <: Real}

    # Set up
    t  = (x - knots[x_bins]) ./ dknots[x_bins]
    t² = t .^ 2
    s  = 1. .- t
    s² = (1. .- t) .^ 2

    # No extrapolation yet so directly return
    x_bins_up = x_bins .+ 1

    return coef[x_bins, 2] .* (3. .* s² - 2. .* s) +
        coef[x_bins_up, 2] .* (3. .* t² - 2. .* t) +
        6 .* (coef[x_bins, 1] .* (-s + s²) +
              coef[x_bins_up, 1] .* (t - t²)) ./ dknots[x_bins]
end

function eval_cubic_derivative_2(x::AbstractVector{T}, knots::AbstractVector{T}, dknots::AbstractVector{T},
                                 coef::AbstractMatrix{T}, x_bins::AbstractVector{Int}) where {T <: Real}

    # Set up
    t  = (x - knots[x_bins]) ./ dknots[x_bins]
    s  = 1. .- t

    # No extrapolation yet so directly return
    x_bins_up = x_bins .+ 1

    return (coef[x_bins, 2] .* (2. .- 6. .* s) +
            coef[x_bins_up, 2] .* (6 .* t .- 2.)) ./ dknots[x_bins] +
            (coef[x_bins, 1] .* (6. .- 12. .* s) +
             coef[x_bins_up, 1] .* (6. .- 12. .* t)) ./ (dknots[x_bins] .^ 2)

end

function eval_cubic_derivative_3(x::AbstractVector{T}, knots::AbstractVector{T}, dknots::AbstractVector{T},
                                 coef::AbstractMatrix{T}, x_bins::AbstractVector{Int}) where {T <: Real}

    # No extrapolation yet so directly return
    x_bins_up = x_bins .+ 1

    return 6. .* (coef[x_bins, 2] .+ coef[x_bins_up, 2]) ./ (dknots[x_bins] .^ 2) +
        12. .* (coef[x_bins, 1] .- coef[x_bins_up, 1]) ./ (dknots[x_bins] .^ 3)
end


text:
"""
    @pretty @cast A[...] := B[...]

Prints an approximately equivalent expression with the macro expanded.
Compared to `@macroexpand1`, generated symbols are replaced with animal names (from MacroTools),
comments are deleted, module names are removed from functions,
and the final expression is fed to `println()`.

To copy and run the printed expression, you may need various functions which aren't exported.
Try something like `using TensorCast: orient, star, rview, @assert_, red_glue, sliceview`
"""
macro pretty(ex)
    if @capture(ex, @cast_str str_)
        full = "@cast " * cast_string(str)
        println("# " * full)
        return :(@pretty $(Meta.parse(full)) )

    elseif @capture(ex, @reduce_str str_)
        full = reduce_string(str)
        println("# " * full)
        return :(@pretty $(Meta.parse(full)) )

    elseif @capture(ex, @matmul_str str_)
        full = matmul_string(str)
        println("# " * full)
        return :(@pretty $(Meta.parse(full)) )

    else
        :( macroexpand($(__module__), $(ex,)[1], recursive=false) |> pretty |> println )
    end
end

function pretty(ex::Union{Expr,Symbol})
    ex = MacroTools.alias_gensyms(ex) # animal names
    ex = MacroTools.striplines(ex)    # remove line references
    ex = MacroTools.postwalk(ex) do x # remove esc()
        (x isa Expr && x.head == :escape && length(x.args) == 1) ? x.args[1] : x
    end
    pretty(string(ex))
end

function pretty(str::String)
    str = replace(str, r"\n(\s+\n)" => "\n")      # remove empty lines

    str = replace(str, r"\w+\.var\"(@\w+)\"" => s"\1")   # Julia 1.3, deal with e.g. Strided.var"@strided"

    str = replace(str, r"\w+\.(\w+)\(" => s"\1(") # remove module names on functions
    str = replace(str, r"\w+\.(\w+)!\(" => s"\1!(")
    str = replace(str, r"\w+\.(@\w+)" => s"\1")   # and on macros
    str = replace(str, r"\w+\.(\w+){" => s"\1{")  # and on structs...

    str = replace(str, "Colon()" => ":")
    str = replace(str, r"\$\(QuoteNode\((\d+)\)\)" => s":\1")
    str = replace(str, r"for (\w+) = " => s"for \1 in ")
end

pretty(tup::Tuple) = pretty(string(tup))

pretty(vec::Vector) = pretty(join(something.(vec, "nothing"), ", ")) # for indices, [i,j,k]
pretty(many...) = pretty(join(something.(many, "nothing"), ""))      # for unparse, no commas

pretty(x) = string(x)

function unparse(str::String, exs...)
    @capture(exs[1], left_ = right_ ) && return pretty(str, " ", left, " = ", right, "  ", join(exs[2:end],"  "))
    @capture(exs[1], left_ := right_ ) && return pretty(str, " ", left, " := ", right, "  ", join(exs[2:end],"  "))
    @capture(exs[1], left_ |= right_ ) && return pretty(str, " ", left, " |= ", right, "  ", join(exs[2:end],"  "))

    @capture(exs[1], left_ += right_ ) && return pretty(str, " ", left, " += ", right, "  ", join(exs[2:end],"  "))
    @capture(exs[1], left_ -= right_ ) && return pretty(str, " ", left, " -= ", right, "  ", join(exs[2:end],"  "))
    @capture(exs[1], left_ *= right_ ) && return pretty(str, " ", left, " *= ", right, "  ", join(exs[2:end],"  "))

    @capture(exs[1], red_(ind__) ) && return pretty(str, " ", join(exs,"  "))

    return pretty(exs)
end


text:
using CANalyze.Signals
using CANalyze.Messages
using CANalyze.Databases

signal1 = Signals.NamedSignal("ABC", nothing, nothing, Signals.Float32Signal(start=0, byte_order=:little_endian))
signal2 = Signals.NamedSignal("ABCD", nothing, nothing, Signals.Unsigned(start=40, length=17, factor=2, offset=20, byte_order=:big_endian))
signal3 = Signals.NamedSignal("ABCDE", nothing, nothing, Signals.Unsigned(start=32, length=8, factor=2, offset=20, byte_order=:little_endian))

m1 = Messages.Message(0x1FE, 8, "ABC", signal1; strict=true)
m2 = Messages.Message(0x1FF, 8, "ABD", signal1, signal2, signal3; strict=true)

d = Databases.Database(m1, m2)
# println(d)

m = d["ABC"]
println(m)

m = d[0x1FF]
println(m)

m = get(d, 0x1FA)
println(m)


text:
using LinearAlgebra: norm

using ..Bodies: CelestialBody, sun, mean_radius, ellipsoid
using ..Ephemerides: AbstractEphemeris, position
using ..Frames: Rotation, icrf, iau_earth
using ..TwoBody: isretrograde
using ..Util: angle, point_on_limb

abstract type Detector end
abstract type DiscreteDetector <: Detector end
abstract type IntervalDetector <: Detector end

struct Event{Scale, TType, T<:Detector}
    type::Symbol
    epoch::Epoch{Scale, TType}
    detector::T
end

struct Apocenter <: DiscreteDetector end

function detect(::Apocenter, t, tra)
    y = tra(t)
    el = keplerian(y[1:3], y[4:6], grav_param(centralbody(tra)))
    ano = el[6]
    if ano > pi/2
        ano = abs(ano - pi)
    elseif ano < -pi/2
        ano = -abs(ano + pi)
    end
    return isretrograde(el[3]) ? ano : -ano
end

struct Pericenter <: DiscreteDetector end

function detect(::Pericenter, t, tra)
    y = tra(t)
    el = keplerian(y[1:3], y[4:6], grav_param(centralbody(tra)))
    return isretrograde(el[3]) ? -el[6] : el[6]
end

struct Eclipse{T<:AbstractEphemeris, C1<:CelestialBody, C2<:CelestialBody} <: IntervalDetector
    eph::T
    occulting::C1
    occulted::C2
    total::Bool
end

function detect(ecl::Eclipse, t, tra)
    find_eclipse(ecl, epoch(tra) + t, tra(t)[1:3], centralbody(tra))
end

function find_eclipse(ecl::Eclipse, ep, pos, center)
    # FIXME: Use proper rotation instead of hard-coded values
    rot = Rotation(icrf, iau_earth, ep)
    pocc = position(ecl.eph, ep, center, ecl.occulting)
    pted, _ = rot(position(ecl.eph, ep, ecl.occulting, ecl.occulted), zeros(3))
    psat, _ = rot(pos .- pocc, zeros(3))
    plimb = point_on_limb(ellipsoid(ecl.occulting), psat, pted)
    ps = psat .- pted
    pi = psat .- plimb
    ang = angle(ps, psat)
    rs = asin(mean_radius(ecl.occulted) / norm(ps))
    isnan(rs) && return π
    ro = angle(pi, psat)
    val = ecl.total ? (ang - ro + rs) : (ang - ro - rs)
    return -val
end


text:
"""
# LoggingTestSets

`LoggingTestSet` is an [`AbstractTestSet`](https://docs.julialang.org/en/v1/stdlib/Test/#Creating-Custom-AbstractTestSet-Types)
that logs test results using `@info` and `@error` from the
[Logging][https://docs.julialang.org/en/v1/stdlib/Logging/] module.

```julia
julia> using LoggingTestSets
julia> using Test
julia> @testset LoggingTestSet "Some tests" begin
           @test 1 == 1
           @test 1 == 2
       end

┌ Error: Some tests: Test Failed at REPL[6]:3
│   Expression: 1 == 2
│    Evaluated: 1 == 2
└ @ LoggingTestSets LoggingTestSets/src/LoggingTestSets.jl:49
┌ Info: Test Summary: Some tests
│   passes = 1
│   fails = 1
│   errors = 0
└   broken = 0
Test Summary: | Pass  Fail  Total
Some tests    |    1     1      2
┌ Error: Some tests did not pass: 1 passed, 1 failed, 0 errored, 0 broken.
└ @ LoggingTestSets LoggingTestSets/src/LoggingTestSets.jl:69
ERROR: Some tests did not pass: 1 passed, 1 failed, 0 errored, 0 broken.
```

`TestFileLogger` is an [`AbstractLogger`](https://docs.julialang.org/en/v1/stdlib/Logging/#AbstractLogger-interface)
that writes log messages to a file.

```julia
julia> using LoggingTestSets
julia> using Logging
julia> using LoggingExtras
julia> tee = TeeLogger(global_logger(), TestFileLogger("test.log"))

julia> global_logger(tee)
```
"""
module LoggingTestSets
export LoggingTestSet
export TestFileLogger

using Test
using Test: DefaultTestSet, AbstractTestSet, Fail, Error, scrub_backtrace

using Logging
using LoggingExtras

using Dates: now

# Logging test results.

struct LoggingTestSet <: AbstractTestSet
    ts::DefaultTestSet
    LoggingTestSet(args...; kw...) = new(DefaultTestSet(args...; kw...))
end

function log_failure(ts::LoggingTestSet, t::Fail)
    str = "Test failed in testset $(ts.ts.description)"

    @error(str,
           Expression=t.orig_expr,
           Evaluated=t.data,
           exception=(DomainError(-1.0, "blah"), scrub_backtrace(backtrace())),
           _file=string(t.source.file),
           _line=t.source.line)

    return
end

function log_failure(ts::LoggingTestSet, t::Error)
    str = "Error during test in testset $(ts.ts.description)\n$(t.backtrace)"

    @error(str,
           _file=string(t.source.file),
           _line=t.source.line)

    return
end

function Test.record(ts::LoggingTestSet, t::Union{Fail, Error})

    io = IOBuffer()

    begin # copied from: stdlib/Test/src/Test.jl https://git.io/JqZQk
        print(io, ts.ts.description, ": ")
        # don't print for interrupted tests
        if !(t isa Error) || t.test_type !== :test_interrupted
            print(io, t)
            if !isa(t, Error) # if not gets printed in the show method
                Base.show_backtrace(io, scrub_backtrace(backtrace()))
            end
            println(io)
        end
    end

    @error String(take!(io)) _file=string(t.source.file) _line=t.source.line

    # log_failure(ts, t)

    push!(ts.ts.results, t)
end

Test.record(ts::LoggingTestSet, args...) = Test.record(ts.ts, args...)

Test.finish(ts::LoggingTestSet) = Test.finish(ts.ts)

# Logging to file.

struct TestFileLogger <: AbstractLogger
    io::IO
    TestFileLogger(filename) = new(open(filename, append=true))
end

Logging.shouldlog(::TestFileLogger, args...) = true
Logging.min_enabled_level(::TestFileLogger) = Logging.Debug
Logging.catch_exceptions(::TestFileLogger) = false

function Logging.handle_message(l::TestFileLogger,
                                level, message, _module, group, id, file, line;
                                kwargs...)

    # Copied from: base/logging.jl https://git.io/JqZ5Q
    msglines = split(chomp(string(message)), '\n')

    # Copied from: ConsoleLogger.jl https://git.io/JqZdY
    valbuf=IOBuffer()
    for (key,val) in pairs(kwargs)
        Logging.showvalue(valbuf, val)
        vallines = split(String(take!(valbuf)), '\n')
        if length(vallines) == 1
            push!(msglines, "$key = $(vallines[1])")
        else
            push!(msglines, "$key =")
            append!(msglines, vallines)
        end
    end

    # Copied from: base/logging.jl https://git.io/JqZ5Q
    println(l.io, "┌ ", string(level), " ", now(),  ": ", msglines[1])
    for i in 2:length(msglines)
        println(l.io, "│ ", msglines[i])
    end
    println(l.io, "└ @ ", _module, " ", basename(file), ":", line)
    flush(l.io)
    nothing
end

# Documentation.

readme() = Docs.doc(@__MODULE__)

end # module


text:
function ReceiverAntenna(arg0::JString, arg1::JString, arg2::Map, arg3::JString)
    return ReceiverAntenna((JString, JString, Map, JString), arg0, arg1, arg2, arg3)
end

function get_serial_number(obj::ReceiverAntenna)
    return jcall(obj, "getSerialNumber", JString, ())
end


text:
n = 7 # rand(1:10)

r = mdopen("deriv2", n, false)
@test r !== nothing && r.A !== nothing
rf32 = mdopen("deriv2", Float32, n, false)
@test rf32 !== nothing && rf32.A !== nothing

r2 = mdopen("deriv2", Float32, n, 2, false)
@test r2 !== nothing && r2.b !== nothing
@test mdopen("deriv2", Float32, (n ÷ 2) * 2, 3, false).A !== nothing

@test_throws ArgumentError matrixdepot("deriv2", Float64, n, 4, false)

A = matrixdepot("deriv2", n)
@test A == r.A

@test r.A*r.x ≈ r.b
@test issymmetric(r.A)

r = mdopen("shaw", 2*n, false)
@test r !== nothing && r.A !== nothing
A = matrixdepot("shaw", 2*n)
@test A !== nothing

# print(r)

@test issymmetric(r.A)

rf32 = matrixdepot("shaw", Float32, 2*n, false)
@test rf32 !== nothing

r = mdopen("wing", n, false)
@test r !== nothing && r.A !== nothing
A = matrixdepot("wing", n)
@test A !== nothing && A == r.A

@test r.A == matrixdepot("wing", n, 1/3, 2/3, false)

r = mdopen("foxgood", n, false)
@test r !== nothing && r.A !== nothing
rf32 = mdopen("foxgood", Float32, n, false)
@test rf32 !== nothing && rf32.A !== nothing
A = matrixdepot("foxgood", n)
@test A == r.A

r = mdopen("heat", Float32, 2*n, false)
@test r !== nothing && r.A !== nothing
rf32 = matrixdepot("heat", Float32, 2*n)
@test rf32 == r.A

r = mdopen("baart", 2*n, false)
@test r !== nothing && r.b !== nothing
rf32 = mdopen("baart", Float32, 2*n, false)
@test rf32 !== nothing && rf32.A !== nothing
A = matrixdepot("baart", 2*n)
@test A == r.A

n = 8 # rand(3:10)
r = mdopen("phillips", 4*n, false)
@test r !== nothing && r.A !== nothing
rf32 = matrixdepot("phillips", Float32, 4*n)
@test rf32 !== nothing

A = matrixdepot("gravity", n)
@test A !== nothing

@test matrixdepot("gravity", n, 1, 0, 1, 0.25) == A

@test mdopen("gravity", Float32, n, 1, false).A !== nothing
@test mdopen("gravity", Float32, n, 2, false).A !== nothing
@test mdopen("gravity", Float32, n, 3, false).A !== nothing
@test_throws ErrorException matrixdepot("gravity", Float32, n, 4, false)

A = matrixdepot("blur", n)
@test A !== nothing

@test matrixdepot("blur", n, 3, 0.7) == A

r1 = mdopen("blur", Float32, n, false)
@test r1 !== nothing && r1.A !== nothing

A = matrixdepot("spikes", n)
@test A !== nothing

@test matrixdepot("spikes", n, 5) == A

n = 7 # rand(5:10)
r1 = mdopen("spikes", Float32, n, false)
@test r1 !== nothing && r1.A !== nothing

A = matrixdepot("ursell", n)
@test A !== nothing
@test matrixdepot("ursell", Float64, n) == A
r1 = mdopen("ursell", Float32, n, false)
@test r1 !== nothing && r1.A !== nothing

A = matrixdepot("parallax", n)
@test A !== nothing

m, n = size(A)
@test m == 26
@test matrixdepot("parallax", Float64, n) == A
r1 = mdopen("parallax", Float32, n, false)
@test r1 !== nothing && r1.A !== nothing


text:
using Revise
using BenchmarkTools
using LarSurf
using Statistics
using SparseArrays
# x = 1

# paper
# Calculation of lar model (141 ms on HP) is faster than reading LARmodel from
# file (923 ms on HP)

block_size = [16,16,16]
# block_size = [64,64,64]
# pure computation
LarSurf.set_param(larmodel_allow_memory=false, larmodel_allow_write_files=false; larmodel_allow_read_files=false)
bm = @benchmark LarSurf.get_larmodel(block_size)
println("Mean time: $(mean(bm.times) * 10^-9) [s]")
# mean time [16,16,16]
# 141 ms with JLD
# 124 with JLD2
# mean time [64,64,64]
# 20.67 s with JLD2

# calculate and write
LarSurf.set_param(larmodel_allow_memory=false, larmodel_allow_write_files=true; larmodel_allow_read_files=false)
bm = @benchmark LarSurf.get_larmodel(block_size)
println("Mean time: $(mean(bm.times) * 10^-9) [s]")
# mean time [16,16,16]
# 923 ms with JLD
# 130 with JLD2
# mean time [64,64,64]
# 21.92 s with JLD2

# just read
LarSurf.set_param(larmodel_allow_memory=false, larmodel_allow_write_files=false; larmodel_allow_read_files=true)
bm = @benchmark LarSurf.get_larmodel(block_size)
println("Mean time: $(mean(bm.times) * 10^-9) [s]")
# mean time [16,16,16]
# 983 ms with JLD
# 134 with JLD2
# mean time [64,64,64]
# 20.86 s with JLD2

# Pure computation
LarSurf.set_param(larmodel_allow_memory=false, larmodel_allow_write_files=false; larmodel_allow_read_files=false)
LarSurf.set_param(boundary_allow_memory=false, boundary_allow_write_files=false; boundary_allow_read_files=false)
bm = @benchmark LarSurf.get_boundary3(block_size, false)
println("Mean time: $(mean(bm.times) * 10^-9) [s]")
# mean 16
# 214 ms

# Computation and write
LarSurf.set_param(larmodel_allow_memory=false, larmodel_allow_write_files=false; larmodel_allow_read_files=false)
LarSurf.set_param(boundary_allow_memory=false, boundary_allow_write_files=true; boundary_allow_read_files=false)
bm = @benchmark LarSurf.get_boundary3(block_size, false)
println("Mean time: $(mean(bm.times) * 10^-9) [s]")
# mean 16
# 270 ms

# Computation and read
LarSurf.set_param(larmodel_allow_memory=false, larmodel_allow_write_files=false; larmodel_allow_read_files=false)
LarSurf.set_param(boundary_allow_memory=false, boundary_allow_write_files=false; boundary_allow_read_files=true)
bm = @benchmark LarSurf.get_boundary3(block_size, false)
println("Mean time: $(mean(bm.times) * 10^-9) [s]")
# mean 16
# 139 ms
#


text:
module CanonicalSymplecticForms

import Base
import LinearAlgebra

export canonical_one_form, CanonicalSymplecticMatrix, canonical_two_form

function checkn(T, n::Integer)
    n > 0 || throw(ArgumentError("$T must have positive size"))
    iseven(n) || throw(ArgumentError("$T must have even size"))
end

struct CanonicalSymplecticVector{T, VT <: AbstractVector{T}} <: AbstractVector{T}
    p::VT

    function CanonicalSymplecticVector{T, VT}(p) where {T, VT}
        Base.require_one_based_indexing(p)
        new{T, VT}(p)
    end
end

CanonicalSymplecticVector{T}(p::VT) where {T, VT} = CanonicalSymplecticVector{T, VT}(p)
CanonicalSymplecticVector(p::VT) where VT = CanonicalSymplecticVector{eltype(p), VT}(p)

Base.size(C::CanonicalSymplecticVector) = (length(C.p) * 2,)

function Base.getindex(C::CanonicalSymplecticVector{T}, i::Int) where T
    mid = length(C.p)
    @boundscheck let n = mid * 2
        if !(1 ≤ i ≤ n)
            msg = "attempt to access $n-element CanonicalSymplecticVector at index [$i]"
            throw(BoundsError(msg))
        end
    end

    if i ≤ mid
        return @inbounds C.p[i]
    else
        return zero(T)
    end
end

function LinearAlgebra.dot(C::CanonicalSymplecticVector, x::AbstractVector)
    m = length(C.p)
    axes(x, 1) == 1:2m || throw(DimensionMismatch())

    s = zero(promote_type(eltype(x), eltype(C)))
    @inbounds for i in 1:m
        s += C.p[i] * x[i]
    end

    s
end

# LinearAlgebra.dot(x::AbstractVector, C::CanonicalSymplecticVector) = LinearAlgebra.dot(C, x)

function canonical_one_form(z, t, p)
    n = length(z)
    iseven(n) || throw(ArgumentError("z must have even length"))
    mid = n ÷ 2
    CanonicalSymplecticVector(view(z, mid+1:n))
end

"""
    CanonicalSymplecticMatrix{T}(n::Integer)

constructs a canonical symplectic matrix of size `(n, n)` with eltype `T`.
`n` must be even and positive. See the examples to see the form of the canonical symplectic
matrix as defined here.

# Examples

```jldoctest
julia> CanonicalSymplecticMatrix(4)
4×4 CanonicalSymplecticMatrix{Int64}:
 0  0  -1   0
 0  0   0  -1
 1  0   0   0
 0  1   0   0

julia> CanonicalSymplecticMatrix{Int32}(6)
6×6 CanonicalSymplecticMatrix{Int32}:
 0  0  0  -1   0   0
 0  0  0   0  -1   0
 0  0  0   0   0  -1
 1  0  0   0   0   0
 0  1  0   0   0   0
 0  0  1   0   0   0
```
"""
struct CanonicalSymplecticMatrix{T} <: AbstractMatrix{T}
    mid::Int
    function CanonicalSymplecticMatrix{T}(n::Integer) where T
        checkn(CanonicalSymplecticMatrix, n)
        new{T}(n ÷ 2)
    end
end

CanonicalSymplecticMatrix(n::Integer) = CanonicalSymplecticMatrix{Int}(n)

Base.size(C::CanonicalSymplecticMatrix) = (sz = C.mid * 2; (sz, sz))

function Base.getindex(C::CanonicalSymplecticMatrix{T}, i1::Int, i2::Int) where T
    @boundscheck let n = C.mid * 2
        if !(1 ≤ i1 ≤ n && 1 ≤ i2 ≤ n)
            msg = "attempt to access $n-element CanonicalSymplecticMatrix{$T} at index [$i1, $i2]"
            throw(BoundsError(msg))
        end
    end

    if i2 - C.mid == i1
        return T(-1)
    elseif i1 - C.mid == i2
        return T(1)
    else
        return T(0)
    end
end

function LinearAlgebra.dot(x::AbstractVector, C::CanonicalSymplecticMatrix, y::AbstractVector)
    m = C.mid
    axes(x, 1) == axes(y, 1) == 1:2m || throw(DimensionMismatch())
    s = zero(promote_type(eltype(x), eltype(y)))
    @inbounds for i in 1:m
        s += y[i] * x[m + i]
        s -= y[m + i] * x[i]
    end

    s
end

canonical_two_form(z, t, p) = CanonicalSymplecticMatrix{eltype(z)}(length(z))

end  # module


text:
"""
#
# PEJROOT calculates (multiple) roots of polynomial using an iterative
# method developed by Z. Zeng. For details, you can request the paper
#
#     "Computing Multiple roots of polynomials: pejorative condition and accurate
#     computation", by Zhonggang Zeng
#
# This code is dated Dec. 4, 2002. Please sent bug/failure report to
#
#             Prof. Zhonggang Zeng
#             Department of Mathematics
#             Northeastern Illinois University
#             Chicago, IL 60625
#
#             email: zzeng@neiu.edu
#
# This code is freely released for research exchange only. The author
# is not responsible for any demage caused by using this code.
#
# Calling syntax: The simplest way to call is
#
#     pejroot(f,y,m)
#
# where the input items are
#          f --- (row vector) the target polynomial
#          y --- (row vector) the initial iterate (of the roots)
#          m --- (row vector) the multiplicity structure
#
# For more advanced usage:
#
#	 [z, e, c] = pejroot(f, y, l; noi, tol, style, prtsty)
#
# The output
#          z --- (matrix) distinct roots (1st column) and corresponding
#                         multiplicities (2nd column)
#          e --- (scalar) the backward error of the roots
#          c --- (scalar) the pejorative condition number
#
#
#  Additional input parameters
#
#        noi --- (integer)    the number of iteration allowed
#                               (default: 10)
#        tol --- (numeric)    the error tolerance
#                               (default: 1.0d-8)
#      style --- (integer)    backward error minimization style:
#                               1 : overall (absolute)
#                               2 : coefficientwise (relative,default)
#     prtsty --- (integer)    intermediate results showing style
#                               0: minimal (default)
#                               1: plus intermediate backward error and
#                                     root corrections
#                               2: plus intemediate root approximations
#  Example:
#
#  f = poly([ones(1,20) 2*ones(1,10) 3*ones(1,5)]); # construct the test
#                					# polynomial with multiple roots
#  y = [0.95,2.05,2.95];           # prepare initial iterate
#  m = [20, 10, 5];                # prepare the multiplicity structure
#  pejroot(f,y,m)                  # running the program
# ans =
#
#   3.00000000000000   5.00000000000000
#   2.00000000000000  10.00000000000000
#   1.00000000000000  20.00000000000000
#
"""
function pejroot(f::AbstractVector{T}, z::AbstractVector{S}, mult; noi::Int = 10, tol::U = 1e-8, style::Int = 2, prtsty::Int = 0) where {T<:Number,S<:Number,U<:AbstractFloat}

    E = one(U)
    Z = zero(U)

    m = length(mult)	# number of variables (different roots)
    n = sum(mult)		# number of equations (degree of polynomial f)

    length(f) == n+1 || throw(ArgumentError("degree of polynomial $(length(f)-1) does not match multiplicies $(n)"))
    f[1] != Z || throw(ArgumentError("polynomial has not full degree compared to its length"))

    f = f / f[1]		# make the polynomial monic

    # sort initial values to enhance accuracy
    # It is interesting, sorting really improves accuracy
    jj = sortperm(abs.(z), rev=false)	# greatest first
    y0 = z[jj]
    ll = mult[jj]
    CT, zs = zerostructure(y0, ll)
    if !(CT<:T)
        f = CT.(f)
    end

    job = 0						# initialize job success return value
    y = copy(y0)				# pass the initial iterate
    h = copy(f[2:n+1])	        # make the RHS
    delta = zeros(U, 1, noi)	# space for sizes of the correction
    bcker = zeros(U, 1, noi)	# space for backward errors

    w = ones(T, n)				# set weight
    if style == 2
       for j = 1:n
           afj = abs(f[j+1])
           if afj >= E
               w[j] = E / afj
           end
       end
    end

    if prtsty >= 1
       println("  step   bckerr        delta")
    end

    Df = zeros(T, n ,m)			# open the space for Df
    bkerr = Z
    pjcnd = Z
    kk = 0

    for k = 1:noi
        kk = k

        if prtsty >= 2
            println(y)
        end

        # evaluate the coefficient operator and its Jacobian
        pr = products(y, ll)
        b = polyG(pr) - h
        Df = polyDG(y, ll, pr, zs)

        if style == 2   # scale
            b = w .* b
            lmul!(Diagonal(w), Df)
         end

         d = realtocoeff(Df \ b, S, zs) # least squares

         delta[k] = norm(d, 2)
         bcker[k] = norm(b, 2)

         if prtsty >= 1
             @printf(" %2.0f    %10.2e    %10.2e\n", k, bcker[k], delta[k])
         end

         if delta[k] < tol
            job = 1
         end     							# convergence criterion 1
         if k > 1
             if delta[k] >= delta[k-1] && bcker[k] >= bcker[k-1]
                 if job == 1
                     bkerr = bcker[k-1]
                     # accept result of previous iteration
                     break
                 end
             elseif delta[k] < delta[k-1]	# criterion 2
                 if delta[k]^2 / (delta[k-1] - delta[k]) < tol
                     job = 1
                 end
            end
        end

         yy = y
         step = 1.0
         stepmin = 1/16
         while step >= stepmin
             y = yy - d * step
             pr = products(y, ll)
             b = (polyG(pr) - h)
             if style == 2
                 b = w .* b
             end
             bks = norm(b, 2)
             if isnan(bks)
                 stepmin *= 0.5
             end
             if bks < bcker[k]
                 break
             end
             step *= 0.5
         end
         bkerr = bcker[k]
      end # k-loop

      if job == 1
          s = svdvals(Df)
          pjcnd = 1 / s[m]  # get pej. cond. number
          jj = sortperm(ll)
          y = y[jj]
          ll = ll[jj]   # sort by multiplicities
      end

      PolyZeros(y, ll), bkerr, pjcnd, job, kk
end


text:
function mle_for_θ(y::AbstractVector, μ::AbstractVector, wts::AbstractVector;
                   maxiter=30, tol=1.e-6)
    function first_derivative(θ::Real)
        tmp(yi, μi) = (yi+θ)/(μi+θ) + log(μi+θ) - 1 - log(θ) - digamma(θ+yi) + digamma(θ)
        unit_weights ? sum(tmp(yi, μi) for (yi, μi) in zip(y, μ)) :
                       sum(wti * tmp(yi, μi) for (wti, yi, μi) in zip(wts, y, μ))
    end
    function second_derivative(θ::Real)
        tmp(yi, μi) = -(yi+θ)/(μi+θ)^2 + 2/(μi+θ) - 1/θ - trigamma(θ+yi) + trigamma(θ)
        unit_weights ? sum(tmp(yi, μi) for (yi, μi) in zip(y, μ)) :
                       sum(wti * tmp(yi, μi) for (wti, yi, μi) in zip(wts, y, μ))
    end

    unit_weights = length(wts) == 0
    if unit_weights
        n = length(y)
        θ = n / sum((yi/μi - 1)^2 for (yi, μi) in zip(y, μ))
    else
        n = sum(wts)
        θ = n / sum(wti * (yi/μi - 1)^2 for (wti, yi, μi) in zip(wts, y, μ))
    end
    δ, converged = one(θ), false

    for t = 1:maxiter
        θ = abs(θ)
        δ = first_derivative(θ) / second_derivative(θ)
        if abs(δ) <= tol
            converged = true
            break
        end
        θ = θ - δ
    end
    if !converged
        info_msg = "Estimating dispersion parameter failed, which may " *
            "indicate Poisson distributed data."
        throw(ConvergenceException(maxiter, NaN, NaN, info_msg))
    end
    θ
end

"""
    negbin(formula,
           data,
           link;
           initialθ::Real=Inf,
           maxiter::Integer=30,
           atol::Real=1e-6,
           rtol::Real=1.e-6,
           verbose::Bool=false,
           kwargs...)

Fit a negative binomial generalized linear model to data, while simultaneously
estimating the shape parameter θ. Extra arguments and keyword arguments will be
passed to [`glm`](@ref).

# Keyword Arguments
- `initialθ::Real=Inf`: Starting value for shape parameter θ. If it is `Inf`
  then the initial value will be estimated by fitting a Poisson distribution.
- `maxiter::Integer=30`: See `maxiter` for [`glm`](@ref)
- `atol::Real=1.0e-6`: See `atol` for [`glm`](@ref)
- `rtol::Real=1.0e-6`: See `rtol` for [`glm`](@ref)
- `verbose::Bool=false`: See `verbose` for [`glm`](@ref)
"""
function negbin(F,
                D,
                args...;
                initialθ::Real=Inf,
                maxiter::Integer=30,
                atol::Real=1e-6,
                rtol::Real=1.e-6,
                verbose::Bool=false,
                kwargs...)
    if haskey(kwargs, :maxIter)
        Base.depwarn("'maxIter' argument is deprecated, use 'maxiter' instead", :negbin)
        maxiter = kwargs[:maxIter]
    end
    if haskey(kwargs, :minStepFac)
        Base.depwarn("'minStepFac' argument is deprecated, use 'minstepfac' instead", :negbin)
        minstepfac = kwargs[:minStepFac]
    end
    if haskey(kwargs, :convTol)
        Base.depwarn("`convTol` argument is deprecated, use `atol` and `rtol` instead", :negbin)
        rtol = kwargs[:convTol]
    end
    if !issubset(keys(kwargs), (:maxIter, :minStepFac, :convTol))
        throw(ArgumentError("unsupported keyword argument"))
    end
    if haskey(kwargs, :tol)
        Base.depwarn("`tol` argument is deprecated, use `atol` and `rtol` instead", :negbin)
        rtol = kwargs[:tol]
    end

    maxiter >= 1 || throw(ArgumentError("maxiter must be positive"))
    atol > 0  || throw(ArgumentError("atol must be positive"))
    rtol > 0  || throw(ArgumentError("rtol must be positive"))
    initialθ > 0 || throw(ArgumentError("initialθ must be positive"))

    # fit a Poisson regression model if the user does not specify an initial θ
    if isinf(initialθ)
        regmodel = glm(F, D, Poisson(), args...;
                       maxiter=maxiter, atol=atol, rtol=rtol, verbose=verbose, kwargs...)
    else
        regmodel = glm(F, D, NegativeBinomial(initialθ), args...;
                       maxiter=maxiter, atol=atol, rtol=rtol, verbose=verbose, kwargs...)
    end

    μ = regmodel.model.rr.mu
    y = regmodel.model.rr.y
    wts = regmodel.model.rr.wts
    lw, ly = length(wts), length(y)
    if lw != ly && lw != 0
        throw(ArgumentError("length of wts must be either $ly or 0 but was $lw"))
    end

    θ = mle_for_θ(y, μ, wts; maxiter=maxiter, tol=rtol)
    d = sqrt(2 * max(1, deviance(regmodel)))
    δ = one(θ)
    ll = loglikelihood(regmodel)
    ll0 = ll + 2 * d

    converged = false
    for i = 1:maxiter
        if abs(ll0 - ll)/d + abs(δ) <= rtol
            converged = true
            break
        end
        verbose && println("[ Alternating iteration ", i, ", θ = ", θ, " ]")
        regmodel = glm(F, D, NegativeBinomial(θ), args...;
                       maxiter=maxiter, atol=atol, rtol=rtol, verbose=verbose, kwargs...)
        μ = regmodel.model.rr.mu
        prevθ = θ
        θ = mle_for_θ(y, μ, wts; maxiter=maxiter, tol=rtol)
        δ = prevθ - θ
        ll0 = ll
        ll = loglikelihood(regmodel)
    end
    converged || throw(ConvergenceException(maxiter))
    regmodel
end


text:
module oneL0

using ..APIUtils

using CEnum

using Printf

using NEO_jll
using oneAPI_Level_Zero_Loader_jll

include("utils.jl")
include("pointer.jl")

# core API
macro check(ex)
    quote
        res = $(esc(ex))
        if res != RESULT_SUCCESS
            throw_api_error(res)
        end

        return
    end
end
include("libze.jl")

# level zero's structure types are often assumed to be zero-initialized (`= {}` in C).
# Julia's memory is not, so define default constructors that ensure everything is zero.
#
# at the same time, our structs are immutable, so we can't just memset to 0 and set fields.
# so instead the constructor we generate has keyword arguments for every fields, and the
# default value for every field is set to 0 (can be overridden by defining `zeroinit`).
#
# TODO: is it really required to (1) zero-initialize memory and (2) set the stype field?
#       none of these are documented...
for (structure_type_enum, _) in CEnum.name_value_pairs(ze_structure_type_t)
    structure_type_name = string(structure_type_enum)
    @assert startswith(structure_type_name, "ZE_STRUCTURE_TYPE")

    T = Symbol("ze_" * lowercase(structure_type_name[19:end]) * "_t")
    if isdefined(oneL0, T)
        struct_typ = getfield(oneL0, T)

        args = Expr[]
        for field in fieldnames(struct_typ)
            field_type = fieldtype(struct_typ, field)
            field_value = if field_type == ze_structure_type_t
                :($structure_type_enum)
            else
                :(zeroinit($field_type))
            end
            push!(args, Expr(:kw, field, field_value))
        end

        @eval begin
            @inline function $T(;$(args...))
                $(T)($(fieldnames(struct_typ)...))
            end
        end
    end
end

# alternative approach: make the structs mutable, and do memset + setfield
zeroinit(::Type{T}) where {T} = convert(T, 0)
zeroinit(::Type{T}) where {T<:CEnum.Cenum} = T(0)
zeroinit(::Type{T}) where {T<:NTuple} = ntuple(_->zeroinit(T.parameters[1]), length(T.parameters))
zeroinit(::Type{ze_driver_uuid_t}) = ze_driver_uuid_t(ntuple(_->zero(UInt8), 16))
zeroinit(::Type{ze_device_uuid_t}) = ze_device_uuid_t(ntuple(_->zero(UInt8), 16))
zeroinit(::Type{_ze_native_kernel_uuid_t}) =
    _ze_native_kernel_uuid_t(ntuple(_->zero(UInt8), 16))
zeroinit(::Type{ze_kernel_uuid_t}) =
    ze_kernel_uuid_t(ntuple(_->zero(UInt8), 16), ntuple(_->zero(UInt8), 16))

# core wrappers
include("error.jl")
include("common.jl")
include("driver.jl")
include("device.jl")
include("context.jl")
include("cmdqueue.jl")
include("cmdlist.jl")
include("fence.jl")
include("event.jl")
include("barrier.jl")
include("module.jl")
include("memory.jl")
include("copy.jl")
include("residency.jl")

const functional = Ref{Bool}(false)

function __init__()
    res = unsafe_zeInit(0)
    if res == RESULT_ERROR_UNINITIALIZED
        @error """No compatible oneAPI driver implementation found.
                  Your hardware probably is not supported by any oneAPI driver.

                  oneAPI.jl currently only supports the Intel Compute runtime,
                  consult their README for a list of compatible hardware:
                  https://github.com/intel/compute-runtime#supported-platforms"""
    elseif res !== RESULT_SUCCESS
        throw_api_error(res)
    else
        functional[] = true
    end
end

end


text:
using ViewerGL
GL = ViewerGL

points,lines = GL.text("PLaSM")
GL.VIEW([
      GL.GLFrame2
      GL.GLLines(points,lines)
]);


text:
# This file is a part of UpROOT.jl, licensed under the MIT License (MIT).

__precompile__(true)

module UpROOT

using Statistics

using ArraysOfArrays
using PyCall
using StatsBase
using Tables
using TypedTables

awkward = PyNULL()
uproot = PyNULL()
numpy = PyNULL()

include("testdata.jl")
include("pywrappers.jl")
include("rootio.jl")
include("opaque.jl")
include("pyjlconv.jl")
include("tdirectory.jl")
include("ttree.jl")

function __init__()
    copy!(awkward, pyimport_conda("awkward0", "awkward0", "conda-forge"))
    copy!(uproot, pyimport_conda("uproot3", "uproot3", "conda-forge"))
    copy!(numpy, pyimport("numpy"))

    for k in _testfile_keys
        testfiles[k] = _testfilename(k)
    end
end

end # module


text:
module RomanNumerals

# Defines the exception InvalidRomanError and the immutable type RomanNumeral
include("types.jl")

# Defines the functions toroman and fromroman, to convert from Roman-formatted
# string to int and vice versa
include("roman_conversion.jl")

# Defines new method dispatching specific to Roman numerals
include("arithmetic.jl")

export RomanNumeral, @rn_str, InvalidRomanError

end # module RomanNumerals


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""
    ObjsPaging(;
        count=nothing,
        page=nothing,
        pages=nothing,
        per_page=nothing,
        spill=nothing,
        total=nothing,
    )

    - count::Int32
    - page::Int32
    - pages::Int32
    - per_page::Int32
    - spill::Int32
    - total::Int32
"""
mutable struct ObjsPaging <: SwaggerModel
    count::Any # spec type: Union{ Nothing, Int32 } # spec name: count
    page::Any # spec type: Union{ Nothing, Int32 } # spec name: page
    pages::Any # spec type: Union{ Nothing, Int32 } # spec name: pages
    per_page::Any # spec type: Union{ Nothing, Int32 } # spec name: per_page
    spill::Any # spec type: Union{ Nothing, Int32 } # spec name: spill
    total::Any # spec type: Union{ Nothing, Int32 } # spec name: total

    function ObjsPaging(;count=nothing, page=nothing, pages=nothing, per_page=nothing, spill=nothing, total=nothing)
        o = new()
        validate_property(ObjsPaging, Symbol("count"), count)
        setfield!(o, Symbol("count"), count)
        validate_property(ObjsPaging, Symbol("page"), page)
        setfield!(o, Symbol("page"), page)
        validate_property(ObjsPaging, Symbol("pages"), pages)
        setfield!(o, Symbol("pages"), pages)
        validate_property(ObjsPaging, Symbol("per_page"), per_page)
        setfield!(o, Symbol("per_page"), per_page)
        validate_property(ObjsPaging, Symbol("spill"), spill)
        setfield!(o, Symbol("spill"), spill)
        validate_property(ObjsPaging, Symbol("total"), total)
        setfield!(o, Symbol("total"), total)
        o
    end
end # type ObjsPaging

const _property_map_ObjsPaging = Dict{Symbol,Symbol}(Symbol("count")=>Symbol("count"), Symbol("page")=>Symbol("page"), Symbol("pages")=>Symbol("pages"), Symbol("per_page")=>Symbol("per_page"), Symbol("spill")=>Symbol("spill"), Symbol("total")=>Symbol("total"))
const _property_types_ObjsPaging = Dict{Symbol,String}(Symbol("count")=>"Int32", Symbol("page")=>"Int32", Symbol("pages")=>"Int32", Symbol("per_page")=>"Int32", Symbol("spill")=>"Int32", Symbol("total")=>"Int32")
Base.propertynames(::Type{ ObjsPaging }) = collect(keys(_property_map_ObjsPaging))
Swagger.property_type(::Type{ ObjsPaging }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ObjsPaging[name]))}
Swagger.field_name(::Type{ ObjsPaging }, property_name::Symbol) =  _property_map_ObjsPaging[property_name]

function check_required(o::ObjsPaging)
    (getproperty(o, Symbol("page")) === nothing) && (return false)
    (getproperty(o, Symbol("total")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ObjsPaging }, name::Symbol, val)
end


text:
# Opencolor : https://yeun.github.io/open-color/
const GRAY4 = colorant"#CED4dA"
const GRAY7 = colorant"#495057"
const RED4 = colorant"#FF8787"
const RED7 = colorant"#F03E3E"
const PINK4 = colorant"#F783AC"
const GRAPE4 = colorant"#DA77F2"
const VIOLET4 = colorant"#9775FA"
const INDIGO4 = colorant"#748FFC"
const BLUE4 = colorant"#4DABF7"
const CYAN4 = colorant"#3BC9DB"
const TEAL4 = colorant"#38D9A9"
const GREEN4 = colorant"#69DB7C"
const LIME4 = colorant"#A9E34b"
const YELLOW4 = colorant"#FFD43B"
const ORANGE4 = colorant"#FFA94D"

function plot_ARIMA_mean(means, tdir::Symbol, ycol::Symbol,
    title_string::String, output_dir::String, output_prefix::String)

    ENV["GKSwstype"] = "100"

    if tdir == :day
        x_means = 1:366*24
        _xlabel = "day * time"
    elseif tdir == :hour
        x_means = 0:23
        _xlabel = "time"
    else
        error("Time directive must be between :day, :hour")
    end

    arima_mean_path = output_dir * "$(output_prefix).png"
    pl = Plots.plot(x_means, means,
        size = (1920, 1080),
        title = title_string,
        guidefontsize = 32, titlefontsize = 48, tickfontsize = 24, legendfontsize = 24, margin=36PlotMeasures.px,
        guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
        background_color = BG_COLOR, color=:black,
        xlabel=_xlabel, ylabel=String(ycol), legend=false)
    Plots.png(pl, arima_mean_path)

    arima_mean_csvpath = output_dir * "$(output_prefix).csv"
    df = DataFrame(tdir => x_means, :mean => means)
    CSV.write(arima_mean_csvpath, df, writeheader = true)

    nothing
end

function plot_ARIMA_mean_smoothing(means, tdir::Symbol, ycol::Symbol, us, vs,
    title_string::String, output_dir::String, output_prefix::String)

    ENV["GKSwstype"] = "100"

    if tdir == :day
        x_means = 1:366*24
        _xlabel = "day * time"
    elseif tdir == :hour
        x_means = 0:23
        _xlabel = "time"
    else
        error("Time directive must be between :day, :hour")
    end

    arima_fluc_path = output_dir * "$(output_prefix).png"
    pl = Plots.plot(x_means, means,
        size = (1920, 1080),
        title = title_string,
        guidefontsize = 32, titlefontsize = 48, tickfontsize = 24, legendfontsize = 24, margin=36PlotMeasures.px,
        guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
        background_color = BG_COLOR, color=:black,
        xlabel=_xlabel, ylabel=String(ycol), legend=false)
    Plots.plot!(us, vs, linewidth=2, color=:red)
    Plots.png(pl,arima_fluc_path)

    arima_fluc_csvpath = output_dir * "$(output_prefix).csv"
    df = DataFrame(tdir => x_means, :mean => means)
    CSV.write(arima_fluc_csvpath, df, writeheader = true)

    nothing
end

function plot_ARIMA_fluc(df::DataFrame, tdir::Symbol, ycol::Symbol,
    title_string::String, output_dir::String, output_prefix::String)

    ENV["GKSwstype"] = "100"

    arima_fluc_path = output_dir * "$(output_prefix).png"
    dates = DateTime.(df[!, :date], Local)
    pl = Plots.plot(dates, df[!, Symbol(tdir, "_fluc")],
        size = (1920, 1080),
        title = title_string,
        guidefontsize = 32, titlefontsize = 48, tickfontsize = 24, legendfontsize = 24, margin=36PlotMeasures.px,
        guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
        background_color = BG_COLOR, color=:black,
        xlabel="dates", ylabel=String(ycol), legend=false)
    Plots.png(pl,arima_fluc_path)

    arima_fluc_csvpath = output_dir * "$(output_prefix).csv"
    df = DataFrame(:date => dates, :fluc => df[!, Symbol(tdir, "_fluc")])
    CSV.write(arima_fluc_csvpath, df, writeheader = true)

    nothing
end

function plot_ARIMA_acf(df::DataFrame, output_dir::String, output_prefix::String, sim_name::String)
    ENV["GKSwstype"] = "100"

    # ACF
    plot_path = output_dir * "$(output_prefix)_acf_$(sim_name).png"

    pl = Plots.plot(df[:, :lags], df[:, :acf],
        size = (2560, 1080),
        line=:solid, linewidth=5,
        ylim=(min(0.0, minimum(df[:, :acf])), 1.0),
        guidefontsize = 18, titlefontsize = 24, tickfontsize = 18, legendfontsize = 18, margin=30PlotMeasures.px,
        guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
        background_color = BG_COLOR, color= :black,
        xlabel = "lags", ylabel = "acf", legend=false)

    Plots.png(pl, plot_path)

    # PACF
    plot_path = output_dir * "$(output_prefix)_pacf_$(sim_name).png"

    pl = Plots.plot(df[:, :lags], df[:, :pacf],
        size = (2560, 1080),
        line=:solid, linewidth=5,
        ylim=(min(0.0, minimum(df[:, :pacf])), 1.0),
        guidefontsize = 18, titlefontsize = 24, tickfontsize = 18, legendfontsize = 18, margin=30PlotMeasures.px,
        guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
        background_color = BG_COLOR, color= :black,
        xlabel = "lags", ylabel = "pacf", legend=false)

    Plots.png(pl, plot_path)

    nothing
end

function plot_ARIMA_corr(df_corr::DataFrame, output_size::Integer, output_dir::String, output_prefix::String, sym_name::Symbol)

    ENV["GKSwstype"] = "100"

    plot_path = output_dir * "$(output_prefix)_corr_hourly_$(String(sym_name))"

    pl = Plots.plot(float.(df_corr[!, :hour]), float.(df_corr[!, :corr]),
        size = (2560, 1080),
        title="Correlation on hourly prediction", xlabel="hour", ylabel="corr",
        line=:solid, linewidth=5, label="OBS",
        guidefontsize = 18, titlefontsize = 24, tickfontsize = 18, legendfontsize = 18, margin=15PlotMeasures.px,
        guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
        background_color = BG_COLOR, linecolor = LN01_COLOR, legend=false)

    png(pl, plot_path)

    nothing
end

function plot_ARIMA_scatter(dfs::Array{DataFrame, 1}, ycol::Symbol,
    output_size::Integer, output_dir::String, output_prefix::String, sym_name::Symbol)

    ENV["GKSwstype"] = "100"

    for i = 1:output_size
        i_pad = lpad(i, 2, '0')
        sc_path = output_dir * "$(i_pad)/" * "$(output_prefix)_scatter_$(String(sym_name))_$(i_pad)h.png"
        df = dfs[i]
        lim = max(maximum(float.(df[!, :y])), maximum(float.(df[!, :yhat])))

        sc = Plots.scatter(float.(df[!, :y]), float.(df[!, :yhat]),
            size = (1080, 1080),
            xlim = (0, lim), ylim = (0, lim), legend=false,
            guidefontsize = 18, titlefontsize = 24, tickfontsize = 18, legendfontsize = 18, margin=15PlotMeasures.px,
            guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
            title="DNN/OBS ($(i_pad)h)", xlabel="OBS", ylabel="DNN",
            background_color = BG_COLOR, markercolor = MK_COLOR)

        # diagonal line (corr = 1)
        Plots.plot!(collect(0:0.1:lim), collect(0:0.1:lim),
            xlim = (0, lim), ylim = (0, lim), legend=false,
            background_color = BG_COLOR, linecolor = LN_COLOR)
        Plots.png(sc, sc_path)

        sc_csvpath = output_dir * "$(i_pad)/" * "$(output_prefix)_scatter_$(i_pad)h.csv"
        outdf = DataFrame(obs = float.(df[!, :y]), model = float.(df[!, :yhat]))
        CSV.write(sc_csvpath, outdf, writeheader = true)
    end

    nothing
end

function plot_ARIMA_lineplot(dfs::Array{DataFrame, 1}, ycol::Symbol,
    output_size::Integer, test_fdate::DateTime, test_tdate::DateTime,
    output_dir::String, output_prefix::String, sym_name::Symbol)

    ENV["GKSwstype"] = "100"

    for i = 1:output_size
        i_pad = lpad(i, 2, '0')
        line_plotpath = output_dir * "$(i_pad)/" * "$(output_prefix)_line_$(String(sym_name))_$(i_pad)h.png"
        # filter by date because batched input makes `undef` DateTime values and throw Overflow
        df = filter(row -> test_fdate <= row[:date] <= test_tdate, dfs[i])
        dates_h = df[!, :date] .+ Dates.Hour(i)

        pl = Plots.plot(dates_h, float.(df[!, :y]),
            size = (2560, 1080),
            ylim = (0.0,
                max(maximum(float.(df[!, :y])), maximum(float.(df[!, :yhat])))),
            line=:solid, linewidth=5, label="OBS",
            guidefontsize = 18, titlefontsize = 24, tickfontsize = 18, legendfontsize = 18, margin=15PlotMeasures.px,
            guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
            background_color = BG_COLOR, color=LN01_COLOR,
            title=String(ycol) * " by dates ($(i_pad)h)",
            xlabel="date", ylabel=String(ycol), legend=:best)

        pl = Plots.plot!(dates_h, float.(df[!, :yhat]),
            line=:solid, linewidth=5, color=LN02_COLOR, label="DNN")
        Plots.png(pl, line_plotpath)

        line_csvpath = output_dir * "$(i_pad)/" * "$(output_prefix)_line_$(i_pad)h.csv"
        outdf = DataFrame(dates = dates_h, obs = float.(df[!, :y]), model = float.(df[!, :yhat]))
        CSV.write(line_csvpath, outdf, writeheader = true)
    end

    nothing
end


text:
type FloatVectorRepresentation <: Representation
  length::Int
  min::Float
  max::Float
  range::FloatRange

  FloatVectorRepresentation(length::Int, min::Float, max::Float) =
    new(length, min, max, min:eps():max)
end

"""
Float vectors are implemented as a fixed-length vector of floating point
numbers.

**Parameters:**
* `length::Int`, the length of this vector. Defaults to 80, if no value is
supplied.
* `min::Float`, the minimum value that a component of this vector may assume.
Defaults to 0.0, if no value is supplied.
* `max::Float`, the maximum value that a component of this vector may assume.
Defaults to 1.0, if no value is supplied.
"""
float_vector{S <: AbstractString}(args::Dict{S, Any}) =
  float_vector( get(args, "length", 80),
                get(args, "min", 0.0),
                get(args, "max", 1.0))
float_vector() = float_vector(80, 0.0, 1.0)
float_vector(length::Int) = float_vector(length, 0.0, 1.0)
float_vector(length::Int, min::Float, max::Float) =
  FloatVectorRepresentation(length, min, max)

minimum_value(r::FloatVectorRepresentation) = r.min
maximum_value(r::FloatVectorRepresentation) = r.max
codon_type(::FloatVectorRepresentation) = Float
chromosome(r::FloatVectorRepresentation) = Vector{Float}
describe(i::Vector{Float}) = "($(join(i, ", ")))"

"""
Generates a float vector at random using its representation.
"""
generate(r::FloatVectorRepresentation) = Base.rand(r.range, r.length)


text:
using MAT

include("../algorithms/FlowSeed-0.6.jl")

# Set resolution parameters
epsi = 0.1

# Choose 2,3, or 5 percent of the target community to use.
# Relative results of algorithms are comparable for all choices of seed set size.
# Choosing 5% leads to larger seed sets and better results for all methods.
percentofset = 5

datasets = ["DBLP","Amazon","LiveJournal","Orkut"]

for graph = datasets

     # "nosoft" refers to the fact that we include no soft penalties on
     # excluding seed nodes from the output set
     outputstring = "Output"*string(percentofset)*"/"*graph*"slfs_"*string(epsi)*"epsnosoft_truecond.txt"
     outputmatrix = "Output"*string(percentofset)*"/"*graph*"slfs_"*string(epsi)*"epsnosotf_truecond.mat"

    open(outputstring,"w") do f
         write(f,"Flowseed. Strict penalties on excluding starter nodes
         (which are known to be in the target cluster) and soft penalties of 1 on excluding immediate neighbors.
         We assume we know $percentofset percent of the target cluster,
         and we set epsilon to $epsi.\n")
    end

    mat = matread("datasets/"*graph*"-top10.mat")
    A = mat["A"]

    matS = matread("datasets/"*graph*"-seed-starter.mat")
    if percentofset == 2
         S = matS["S2"]
    elseif percentofset == 3
         S = matS["S3"]
    elseif percentofset == 5
         S = matS["S5"]
    else
         println("Choose 2,3, or 5 percent of the target set. Defaulting to 5.")
         S = mat["S5"]
    end

    comm = mat["C"]
    Conds = mat["Conds"]
    d = sum(A,1)'
    volA = sum(nonzeros(A))

    n = size(A,1)
    d = sum(A,1)'
    println("nodes = $n")
    numcom = size(comm,2)

    fs_sets = spzeros(n,numcom)
    fs_stats = zeros(6,numcom)

    for commID = collect(1:numcom)

        ## Choose one community to be the target cluster
        Target = find(comm[:,commID])
        TarSize = length(Target)
        TargetStats = set_stats(A,vec(Target),volA)
        Tcond = round(TargetStats[4],3)
        Tsize = length(Target)
        println("Community $commID has $Tsize nodes and a conductance of $Tcond")

        ## Get the starter for the seed set for this community
        Rstart = find(S[:,commID])
        R = neighborhood(A,Rstart,1)      # Grow by a one-hope neighborhood
        Rn = neighborhood(A,R,1)    # get the immediate neighbors of R...
        Rn = setdiff(Rn,R)          # ...but we exclude R itself
        inRc = ones(n)
        inRc[R] = 0
        Rc = find(inRc)             # complement of R
        numR = length(R)
        volR = sum(d[R])
        fR = volR/(volA - volR)
        epsilon = max(epsi,2*fR)
        epsilon = epsi
        STATS = set_stats(A,R,volA)
        condR = round(STATS[4],4)
        pr,re,f1 = round.(PRF(Target,R),4)
        println("Seed Set: |R| = $numR, condR = $condR, PR = $pr, RE = $re, F1 = $f1")

        open(outputstring,"a") do f
             write(f,"\nCommunity $commID has $Tsize nodes and a conductance of $Tcond \n")
             write(f, "Seed Set: |R| = $numR, condR = $condR, PR = $pr, RE = $re, F1 = $f1 \n")
        end

        ## Run some community detection experiments
        RinS = zeros(numR,1)
        pR = zeros(numR,1)
        # For standard flow seed procedure, include
        # Hard constraint on the starter seed nodes.
        for r = 1:numR
           if in(R[r],Rstart)
                RinS[r] = 1
           end
        end

        tic()
        S1, cond1 = FlowSeed(A,R,Rn,Rc,epsilon,pR,RinS,d,volA,volR,true)
        t1 = toc()
        c1 = round(cond1,4)
        size1 = length(S1)
        prf = PRF(Target,S1)
        pr1 = round(prf[1],3)
        re1 = round(prf[2],3)
        f11 = round(prf[3],3)
        Stats = set_stats(A,S1,volA)
        truecond = round(Stats[4],4)
        println("FlowSeed: \tPR = $pr1 \t RE = $re1 \t F1 = $f11 \t Size = $size1 \t Cond = $c1 \t TrueCond = $truecond\t Time = $t1")

        fs_sets[S1,commID] = 1
        fs_stats[:,commID] = [t1; length(S1); c1; pr1; re1; f11]

        open(outputstring,"a") do f

          write(f, "Flow-Seed: \tPR = $pr1 \t RE = $re1 \t F1 = $f11 \t Size = $size1 \t Cond = $c1 \t TrueCond = $truecond\t Time = $t1 \n")
        end

    end

    matwrite(outputmatrix, Dict( "fs_stats" => fs_stats,"fs_sets" => fs_sets))

end

max_stars_repo_name: nveldt/FlowSeed

max_issues_repo_name: nveldt/FlowSeed

max_forks_repo_name: nveldt/FlowSeed


text:
iy = ix
incx = 1
incy = 1


text:
const MM = 2
const LM = 3
const MP = 4
const LP = 5

struct Event
	id_node1::Int
	id_node2::Int
	position_node1::Array
	position_node2::Array
	type::Int
	time::DateTime
end

struct EventWriter
	file
	queue::Queue{Event}

	function EventWriter(filename::String)
		file = open(filename,"w+")
		queue = Queue{Event}()
		new(file,queue)
	end
end

function add_event(event_writer::EventWriter,event::Event)
	enqueue!(event_writer.queue,event)
end

function write_events(event_writer::EventWriter)
	for event in event_writer.queue
		write(event_writer.file,"$(event.time);$(event.id_node1);$(event.id_node2);$(event.type);$(event.position_node1);$(event.position_node2)\n")
		dequeue!(event_writer.queue)
	end
end

function close_writer(event_writer::EventWriter)
	write_events(event_writer)
	close(event_writer.file)
end

function define_event_and_add(event_writer,node,near_node,time)
	if near_node in node.node_seen
		return false
	end
	type = near_node.status + node.status
	add_event(event_writer, Event(node.id,near_node.id,node.position,near_node.position,type,time))
	return true
end


text:
"""
   A K-Means algorithm for data clustering

# Brief:
    The algorithm takes n-dimensional points represented as a vector where each element of that vector
    corresponds to one dimension of the point. In the beginning, the algorithm select k points.
    Every of the k points is used as a centre of the cluster. In each iteration, every of these k points
    is connected to the nearest points and moved, until it stops changing its position through iterations.

# For more informations:
    https://en.wikipedia.org/wiki/K-means_clustering

# Contributed by:- [Nikola Mircic](https://github.com/Nikola-Mircic)
"""

module KMeans
# A function that takes two points and returns the distance
# between them
function distance(a::Vector{T}, b::Vector{T}) where {T<:Number}
    if length(a) != length(b)
        return 0.0
    end
    sum = 0.0
    for i in 1:length(a)
        sum = sum + (a[i] - b[i]) * (a[i] - b[i])
    end
    return sqrt(sum)
end

# A function that returns the middle point based on the passed
# matrix in which every row contains coordinates of the single point
function findMiddle(points::Vector{Vector{T}}) where {T<:Number}
    middle::Vector{T} = []
    # Calculate arithmetic mean of coordinates for each dimension
    for dim in 1:length(points[1])
        sum = 0.0
        for val in 1:length(points[:, 1])
            sum = sum + points[val][dim]
        end
        push!(middle, sum / length(points[:, 1]))
    end
    return middle
end

# Function that assigns each point to the nearest cluster and then creates a list
# of centres of those clusters
# Returns coordinates of the clusters at the end of the iteration
function moveClusters(
    k::Int,
    dim::Int,
    clusters::Vector{Vector{T}},
    points::Vector{Vector{T}},
) where {T<:Number}
    assign::Vector{Vector{Int}} = [] # A matrix to write point-cluster pairs
    for i in 1:k
        push!(assign, [])
    end
    # For each point we are finding the nearest cluster
    for p in 1:length(points[:, 1])
        min_dis = distance(points[p], clusters[1])
        idx = 1
        for c in 2:k
            if distance(points[p], clusters[c]) < min_dis
                min_dis = distance(points[p], clusters[c])
                idx = c
            end
        end
        push!(assign[idx], p) # Assign a point p to the cluster with index idx
    end

    # After assigning each element to the cluster,
    # we are using those points to calculate the centers of our clusters
    temp_clust = copy(clusters)
    for i in 1:length(assign)
        temp_clust[i] = findMiddle(points[assign[i]])
    end

    return temp_clust
end

# MAIN FUNCTION
# Call this function to find k clusters for passed points
# If you want to have a limited number of iterations, set iter to any positive integer,
# else set it to 0 or any negative number
function cluster(k::Int, iter::Int, points::Vector{Vector{T}}) where {T<:Number}
    points_num = length(points[:, 1]) # Number of points
    clust_dim = length(points[1]) # dimension of the cluster we need to find
    #println("dim: ",points_num,"x",clust_dim)

    if points_num < k
        # If there is more requested clusters than points,
        # set coordinates of each cluster to 0
        return [zeros(clust_dim) for i in 1:k]
    elseif points_num == k
        # Number of points and requested clusters are equal, so just assign every
        # point to its own cluster
        return points
    else
        clusters::Vector{Vector{T}} = [] # Create a matrix to store centre points of the clusters
        step = points_num / k
        idx = 1.0

        # To avoid generating empty clusters, we will assign each of clusters
        # to a different point

        for clust in 1:k
            push!(clusters, points[Int(idx)])
            idx = floor(idx + step)
        end

        #println(clusters)

        if iter <= 0
            # There is no specific number of iteration
            # Program will run untill it finds the perfect solution
            new_clust = moveClusters(k, clust_dim, clusters, points)
            while new_clust != clusters # Check if something has changed in the last iteration
                clusters = new_clust
                new_clust = moveClusters(k, clust_dim, clusters, points)
            end
        else
            new_clust = moveClusters(k, clust_dim, clusters, points)
            for i in 2:iter
                if new_clust == clusters
                    # There is no need for further iterations, return the clusters
                    return clusters
                end
                clusters = new_clust
                new_clust = moveClusters(k, clust_dim, clusters, points)
            end
        end

        return clusters
    end
end
end

max_stars_repo_name: Nikola-Mircic/Julia

max_issues_repo_name: Nikola-Mircic/Julia

max_forks_repo_name: Nikola-Mircic/Julia


text:
using UtilsGJL
using Test

include("utiltests.jl")


text:
#=##############################################################################################
Filename: lift_moment_distribution.jl
Author: Ryan Anderson
Contact: rymanderson@gmail.com
README: define an `Action` object to extract the lifting line lift and moment distributions
=###############################################################################################
using StaticArrays
using LinearAlgebra

"""
    lift_moment_distribution(aircraft, parameters, freestream, environment, step_range, stepi, step_symbol) <: Action

Solves for the aerodynamic force distribution at each step.

NOTE: THIS ACTION DOES NOT SOLVE THE VORTEX LATTICE. Call `solve_wing_CF_CM` prior to calling this action.

# Arguments:

* `aircraft::Aircraft` : `Aircraft` system object
* `parameters<:Parameters` `Parameters` struct
* `freestream::Freestream` : `Freestream` object
* `environment::Environment` `Environment` object
* `step_range::AbstractArray` : array of steps for which the simulation is run
* `stepi::Int` : index of the current step
* `step_symbol::String` : defines the step, e.g. `alpha` or `time`

`parameters <: Parameters` requires the following elements:

* `cfs::Vector{Vector{Array{Float64,2}}}`: [i][j]th element is an array of size (3,nspanwisepanels) of force coefficients cd, cy, cl corresponding to the ith step, jth lifting surface
* `cms::Vector{Vector{Array{Float64,2}}}`: [i][j]th element is an array of size (3,nspanwisepanels) of moment coefficients cmx, cmy, cmz corresponding to the ith step, jth lifting surface

"""
function lift_moment_distribution(aircraft, parameters, freestream, environment, step_range, stepi, step_symbol)

    # extract lift and moment distribution
    cfs, cms = VL.lifting_line_coefficients(aircraft.wing_system.system, aircraft.wing_system.lifting_line_rs, aircraft.wing_system.lifting_line_chords; frame = VL.Wind())
    # store to `parameters`
    nwings = length(aircraft.wing_system.system.surfaces)
    parameters.cfs[stepi] .= cfs
    parameters.cms[stepi] .= cms
    # for iwing in 1:nwings

    # end

    return false
end

"""
    lift_moment_distribution_blownwing(aircraft, parameters, freestream, environment, step_range, stepi, step_symbol) <: Action

Solves for the aerodynamic force distribution at each step. This method solves for the lifting line coefficients without using rotor induced velocity in a blown wing setup, as used in Epema (2017).

NOTE: THIS ACTION DOES NOT SOLVE THE VORTEX LATTICE. Call `solve_CF` prior to calling this action.

# Arguments:

* `aircraft::Aircraft` : `Aircraft` system object
* `parameters<:Parameters` `Parameters` struct
* `freestream::Freestream` : `Freestream` object
* `environment::Environment` `Environment` object
* `step_range::AbstractArray` : array of steps for which the simulation is run
* `stepi::Int` : index of the current step
* `step_symbol::String` : defines the step, e.g. `alpha` or `time`

`parameters <: Parameters` requires the following elements:

* `cfs::Vector{Vector{Array{Float64,2}}}`: [i][j]th element is an array of size (3,nspanwisepanels) of force coefficients cd, cy, cl corresponding to the ith step, jth lifting surface
* `cms::Vector{Vector{Array{Float64,2}}}`: [i][j]th element is an array of size (3,nspanwisepanels) of moment coefficients cmx, cmy, cmz corresponding to the ith step, jth lifting surface

"""
function lift_moment_distribution_blownwing(aircraft, parameters, freestream, environment, step_range, stepi, step_symbol)

    # extract lift and moment distribution
    cfs, cms = VL.lifting_line_coefficients(aircraft.wing_system.system, aircraft.wing_system.lifting_line_rs, aircraft.wing_system.lifting_line_chords; frame = VL.Wind())

    # store to `parameters`
    parameters.cfs[stepi] .= cfs
    parameters.cms[stepi] .= cms

    # correct lift coefficients as done by Epema
    nwings = length(aircraft.wing_system.system.surfaces)
    Vref = aircraft.wing_system.system.reference[1].V
    Vinf = aircraft.wing_system.system.freestream[1].Vinf
    for iwing in 1:nwings
        # first determine cl w/o rotor induced velocities
        panels = aircraft.wing_system.system.surfaces[iwing]
        n_spanwise_sections = size(aircraft.wing_system.system.properties[iwing])[2]
        gammas = [Vref * sum([panel.gamma for panel in aircraft.wing_system.system.properties[iwing][:,i_spanwise_section]])  for i_spanwise_section in 1:n_spanwise_sections]
        velocities = [panel.velocity * Vref for panel in aircraft.wing_system.system.properties[iwing]]
        vwakes = parameters.wake_function[stepi].(VL.top_center.(panels))
        # vwakes_minus_y = [[v[1]; 0.0; v[3]] for v in vwakes]
        # velocities_minus_rotor = velocities .- vwakes_minus_y #w/o rotor-induced velocities
        velocities_minus_rotor = velocities .- vwakes #w/o rotor-induced velocities
        Δs = VL.top_vector.(panels)
        v_perps = LA.cross.(velocities_minus_rotor, Δs) ./ LA.norm.(Δs)
        v_perp_norm = LA.norm.([[v[1]; 0.0; v[3]] for v in v_perps])
        # TODO: interpote to get 1/4 chord velocity
        v_perp = [Statistics.mean(v_perp_norm[:,i]) for i in 1:n_spanwise_sections]
        # this step normalizes by the mac
        cls = 2 * gammas .* v_perp ./ (aircraft.wing_system.cmacs[iwing] * Vinf^2)
        # store to parameters
        parameters.cfs[stepi][iwing][3,:] .= cls
    end

    return false
end

"""
    lift_moment_distribution(aircraft, step_range)

Method returns initialized elements required for the `parameters <: Parameters` struct during simulation.

# Arguments:

* `aircraft::Aircraft` : system to be simulated
* `step_range::AbstractArray` : defines each step of the simulation

# Returns:

* `cfs::Vector{Vector{Array{Float64,2}}}`: [i][j]th element is an array of size (3,nspanwisepanels) of force coefficients cd, cy, cl corresponding to the ith step, jth lifting surface
* `cms::Vector{Vector{Array{Float64,2}}}`: [i][j]th element is an array of size (3,nspanwisepanels) of moment coefficients cmx, cmy, cmz corresponding to the ith step, jth lifting surface

"""
function lift_moment_distribution(aircraft, step_range)

    nwings = length(aircraft.wing_system.system.surfaces)
    nspanwisepanels = [size(surface)[2] for surface in aircraft.wing_system.system.surfaces]
    cfs = [[zeros(3, length(aircraft.wing_system.system.surfaces[i])) for i in 1:nwings] for j in 1:length(step_range)]
    cms = deepcopy(cfs)

    return cfs, cms
end

"""
    lifting_line_coefficients_blownwing(system, r, c, wakefunction)

This version solves for the lifting line coefficients without the rotor induced velocities in the equation. For simplicity this neglects the left and right vortex filaments.

Modeled after `lifting_line_coefficients` in `VortexLattice.jl`.
"""
function lifting_line_coefficients_blownwing(system, r, c, wakefunction)

    TF = promote_type(eltype(system), eltype(eltype(r)), eltype(eltype(c)))
    nsurf = length(system.surfaces)
    cf = Vector{Matrix{TF}}(undef, nsurf)
    cm = Vector{Matrix{TF}}(undef, nsurf)
    for isurf = 1:nsurf
        ns = size(system.surfaces[isurf], 2)
        cf[isurf] = Matrix{TF}(undef, 3, ns)
        cm[isurf] = Matrix{TF}(undef, 3, ns)
    end

    return lifting_line_coefficients_blownwing!(cf, cm, system, r, c, wakefunction)
end

"""
    lifting_line_coefficients_blownwing!(cf, cm, system, r, c)

For simplicity, this neglects effects of the left and right vortex filaments.

Modeled after `lifting_line_coefficients!` in `VortexLattice.jl`.
"""
function lifting_line_coefficients_blownwing!(cf, cm, system, r, c, wakefunction)

    # number of surfaces
    nsurf = length(system.surfaces)

    # check that a near field analysis has been performed
    @assert system.near_field_analysis[] "Near field analysis required"

    # extract reference properties
    ref = system.reference[]

    # iterate through each lifting surface
    for isurf = 1:nsurf

        nc, ns = size(system.surfaces[isurf])

        # extract current surface panels and panel properties
        panels = system.surfaces[isurf]
        properties = system.properties[isurf]
        cr = CartesianIndices(panels)

        # loop through each chordwise set of panels
        for j = 1:ns

            I = cr[j]

            # calculate segment length
            rls = SVector(r[isurf][1,j], r[isurf][2,j], r[isurf][3,j])
            rrs = SVector(r[isurf][1,j+1], r[isurf][2,j+1], r[isurf][3,j+1])
            ds = norm(rrs - rls)

            # calculate reference location
            rs = (rls + rrs)/2

            # calculate reference chord
            cs = (c[isurf][j] + c[isurf][j+1])/2

            # calculate section force and moment coefficients
            cf[isurf][:,j] .= 0.0
            cm[isurf][:,j] .= 0.0
            for i = 1:nc
                # add influence of bound vortex
                rb = VL.top_center(panels[i,j])
                Γ = properties[i,j].gamma
                Vlocal = properties[i,j].velocity
                Vlocal -= wakefunction(rb)

                Δs = VL.top_vector(panels[I])
                cfb = VL.RHO * Γ * cross(Vlocal, Δs)

                # cfb = properties[i,j].cfb
                cf[isurf][:,j] .+= cfb
                cm[isurf][:,j] .+= cross(rb - rs, cfb)
            end

            # update normalization
            # cf[isurf][:,j] .*= 2 ./ (cs * (VL.freestream_velocity(system.freestream[1])) .^ 2)
            cf[isurf][:,j] .*= ref.S/(ds*cs)
            cm[isurf][:,j] .*= ref.S/(ds*cs^2)
        end
    end

    return cf, cm
end

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
using Gnuplot
@gp 2:3 "notit"  xrange = (0,1) yrange =(0,1)
@gp :- "set obj 1 rect from 0.31,0.31 to 0.5,0.5 fc lt 1"
@gp :- "set obj 2 rect from 0.1,0.1 to 0.3,0.3 fc 'red'"
@gp :- "set obj 3 rect from 0.5,0.2 to 0.65,0.3 fs empty border rgb 'blue'"
@gp :- "set obj 4 rect at 0.6,0.62 size 0.3,0.2 fc 'orange'"
@gp :- "set obj 5 rect at 0.5,0.85 size 0.7,0.2 fc 'black'"
@gp :- "set obj 5 fill transparent solid 0.2"
save(term ="pngcairo size 600,400", output ="rectangles.png")

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
#Here will will run a similar file to Bayesian_Opt_RNN_XOR but instead we will only take the means of the deterministic values as the sequence develops

#Bayesian_Opt_RNN_XOR.jl

#Set the random seed:

srand(1234)

include("Reccurent_XOR.jl")
include("Kernals.jl")
include("gaussian_process.jl")

# Here we will create the reccurent neural network:

Layer_1=uniform(0,1,1,2)
Layer_2=uniform(0,1,2,1)
recurrent_layer=uniform(0,1,2,2)

Seq_Len=33
epochs=100 #1000 doesnt demonstrate much as run time is longer

function hyper_curry(h)
    return (x->sigmoid(x,h))
end

function hyper_curry_deriv(h)
    return (x->sigmoid_deriv(x,h))
end

#initialise learning rates and hyper-parameters

a=0.001  #Change to 0.001
b=1
c=0.001
d=1
N=10 #10 does well

#15 gives Minimum Average Square Error for Random Selection = 0.11429806093691719
#Minimum Average Square Error for Bayesian Optimization = 0.1226313941360362

Random_Learning_Rates=uniform(a,b,N,1)

Random_Hyperparameters=uniform(c,d,N,1)

Random_Mat=cat(2,Random_Learning_Rates,Random_Hyperparameters)
Random_Average_of_Deterministic_Terms=zeros(N)

# Random_Varience=zeros(N)

inloop_rand_t=0

for i=1:length(Random_Learning_Rates)
    tic()
    node_function=hyper_curry(Random_Mat[i,2])
    node_deriv=hyper_curry_deriv(Random_Mat[i,2])
    learning_rate=Random_Mat[i,1]

    #Here P is the first epoch of the temporal and the last epoch of the temporal SE Values in a 1x2 matrix
    P=Train_Reccurent_Net_Loop(epochs,Layer_1,Layer_2,recurrent_layer,learning_rate,node_function,node_deriv,Seq_Len)

 #   Random_Varience[i]=maximum(P[2])-minimum(P[2])

    #Here we will only take the MSE average of the values that can be predicted,i.e every third term:
    determ_vals=P[2][3:3:end]

    Random_Average_of_Deterministic_Terms[i]=(1/length(determ_vals))*sum(determ_vals)

    #Random_MSE[i]=Train_Neural_Net_Loop(epochs,Layer_1,Layer_2,learning_rate,node_function,node_deriv)[3]
    println("Epoch Complete")
    inloop_rand_t+=toc()

    #plotting for below- may want to remove this

    if i==1
        y0=P[1]
        yend=P[2]
        xvals=[i for i=1:length(y0)]

        println("The Varience for the first epoch is ", maximum(y0)-minimum(y0))
        println("The Varience for the last epoch is ",maximum(yend)-minimum(yend))

        using PyPlot

        plot(xvals,y0,label="MSE of First Epoch",alpha=0.3,color="#aa231f")
        plot(xvals,yend,label="MSE of Last Epoch",alpha=0.8,color="#40d5bb")
        title("MSE Plot For Each Term in XOR Sequence ")
        # xlabel(L"$n$")
        # ylabel(L"${XOR}_n$")
        legend(loc="upper right",fancybox="true")
        axis("tight")
        grid("off")
        show()
    end

end

println("Random Learning Rates Training Completed")

# Plot the random hyper-parameter and Avergae SE values:

using PyPlot
# fig = figure("pyplot_subplot_mixed",figsize=(7,7))
ax=axes()

surf(reshape(Random_Learning_Rates,size(Random_Average_of_Deterministic_Terms)),reshape(Random_Hyperparameters,size(Random_Average_of_Deterministic_Terms)),Random_Average_of_Deterministic_Terms,alpha=0.7,color="#40d5bb")
title("Average Temporal SE Plot for varied Sigmoid Hyper-Parameters and Learning Rates")
xlabel("Learning Rates")
ylabel("Hyper-Parameters")
zlabel("Mean Square Error")
grid("on")
show()

#Bayesian Optimization Example===================================================

#Here are the points we can pick from in the Optimization
number_of_points=75 # at 75 was 0.12035 vs 0.11 #try

LR_Test=linspace(a,b,number_of_points)
HP_Test=linspace(c,d,number_of_points)

#Here is the carteisan product of these written as a vector
Test=gen_points([LR_Test,HP_Test])[1]

#We first have to pick a random point to begin bayesian optimization:

#currently starts with the midpoint, possibly randomise this:
Bayesian_Points=[Test[Int(round(length(Test)/2))]]

#Bayesian_Points is an vector of arrays where in each array first entry is LR second entry is Hyper-Parameters:

#Define hyperparemeter functions:
node_function=hyper_curry(Bayesian_Points[1][2])
node_deriv=hyper_curry_deriv(Bayesian_Points[1][2])

#Define Learning Rate:
learning_rate=Bayesian_Points[1][1]

#Run first train before Bayesian Optimization:
P=Train_Reccurent_Net_Loop(epochs,Layer_1,Layer_2,recurrent_layer,learning_rate,node_function,node_deriv,Seq_Len)
determ_vals=P[2][3:3:end]
Bayes_Average_of_Deterministic_Terms=[(1/length(determ_vals))*sum(determ_vals)]

mu=zeros(number_of_points)
sigma=zeros(number_of_points)
#Begin Bayesian Optimization:

# sigma_control=linspace(1,0.1,N) #Here we have added a parameter that fades away with time

inloop_bayes_t=0

#Bayes_Varience=zeros(N)

for k=2:N
    tic()
    D=[(Bayesian_Points[i],Bayes_Average_of_Deterministic_Terms[i]) for i=1:length(Bayesian_Points)]
    mu, sigma, D=gaussian_process_chol(std_exp_square_ker,D,1e-6,Test)  #This line here has been changed to chol version of gaussian_process
    println("Gaussian Process", k, " Complete","\r")
    mu=reshape(mu,length(mu));
    sigma=reshape(sigma,length(sigma))

    new_point=findmin(mu-sigma)[2]

    #Here we will need to change the number 2 to k
    Bayesian_Points=cat(1,Bayesian_Points,[Test[new_point]])

    learning_rate=Bayesian_Points[k][1]

    node_function=hyper_curry(Bayesian_Points[k][2])

    node_deriv=hyper_curry_deriv(Bayesian_Points[k][2])

    P=Train_Reccurent_Net_Loop(epochs,Layer_1,Layer_2,recurrent_layer,learning_rate,node_function,node_deriv,Seq_Len)

 #   Bayes_Varience[k]=maximum(P[2])-minimum(P[2])

    determ_vals=P[2][3:3:end]
    value_to_be_appended=(1/length(determ_vals))*sum(determ_vals)

    inloop_bayes_t+=toc()

    if k==N

        y0=P[1]
        yend=P[2]
        xvals=[i for i=1:length(y0)]

        println("The Varience for the first epoch is ", maximum(y0)-minimum(y0))
        println("The Varience for the last epoch is ",maximum(yend)-minimum(yend))

        using PyPlot

        plot(xvals,y0,label="MSE of First Epoch",alpha=0.3,color="#aa231f")
        plot(xvals,yend,label="MSE of Last Epoch",alpha=0.8,color="#40d5bb")
        title("MSE Plot For Each Term in XOR Sequence ")
        xlabel(L"$n$")
        ylabel(L"${XOR}_n$")
        legend(loc="upper right",fancybox="true")
        axis("tight")
        grid("off")
        show()
    end

    if value_to_be_appended !=Bayes_Average_of_Deterministic_Terms[k-1]
        Bayes_Average_of_Deterministic_Terms=cat(1,Bayes_Average_of_Deterministic_Terms,[value_to_be_appended])
        println("Epoch ", k, " Complete")

    else
        println("Found Optimum on the ", k-1, " iteration of ", N, " iterations")
        Bayesian_Points=Bayesian_Points[1:length(Bayesian_Points)-1]
        break
    end

end

println("Bayesian_Learning_Rates Training Complete")
println("total time rand = ",inloop_rand_t)
println("total time Bayes = ",inloop_bayes_t)

# Bayesian Plotting =========================================================

LR=[Bayesian_Points[i][1] for i=1:length(Bayesian_Points)]
HP=[Bayesian_Points[i][2] for i=1:length(Bayesian_Points)]

#Move this to the bottom
println("Minimum Average Square Error for Random Selection = ", minimum(Random_Average_of_Deterministic_Terms))
println("Minimum Average Square Error for Bayesian Optimization = ", minimum(Bayes_Average_of_Deterministic_Terms))

#println("The varience for the error of Random Selection is = ",Random_Varience[  findmin(Random_Average_of_Deterministic_Terms)[2]]  )
#println("The varience for the error by Bayesian Optimization is = ",Bayes_Varience[ findmin(Bayes_Average_of_Deterministic_Terms)[2]] )
println("The Minimum Average Square Error for Random Selection is found at epoch ",findmin(Random_Average_of_Deterministic_Terms)[2])
println("The Minimum Average Square Error for Bayesian Optimization is found at epoch ",findmin(Bayes_Average_of_Deterministic_Terms)[2])

using PyPlot
# fig = figure("pyplot_subplot_mixed",figsize=(7,7))
# ax=axes()

surf(LR,HP,Bayes_Average_of_Deterministic_Terms,alpha=0.65,color="#40d5bb")
surf(LR,HP,mu+2*sigma,alpha=0.3) #This should be just mu not 2*mu
surf(LR,HP,mu-2*sigma,alpha=0.3)

title("Average Temporal SE Plot for Optimized Sigmoid Hyper-Parameters and Learning Rates")
xlabel("Learning Rates")
ylabel("Hyper-Parameters")
zlabel("Average Square Error")

grid("off")
show()