text:
using Agents
using Random
using FileIO

const AXIAL_DIRECTIONS = [-1, 0, 1]
const NEIGHBORHOOD =
    [CartesianIndex(a) for a in Iterators.product([-1:1 for φ = 1:2]...) if a != (0, 0)]

@agent Bacterium GridAgent{2} begin
    species::Int
    age::Int        # age, in iterations. They only live so long
    energy::Float64 # energy. Less than zero, it dies. More than threshold, it doubles
    sensory_radius::Float64 # how far it can see food
    reproduction_threshold::Float64 # threshold energy to reproduce by binary fission
    speed::Int      # how far it can move in one iteration
    food_target::Dims{2}    # the food it is moving toward
end

struct FoodData
    neighborhood::Vector{CartesianIndex{2}} # what neighbors can food spread to
    food_cap::Float64   # limit on maximum food in one tile
    regen_rate::Float64 # regeneration rate for food tiles
    spread_multiplier::Float64  # proportionality constant for food spreading
    random_spread_chance::Float64   # probability of random spread
end

### Constructors
FoodData(;
    neighborhood::Vector{CartesianIndex{2}} = NEIGHBORHOOD,
    food_cap::Real = 200.0,
    regen_rate::Real = 10.0,
    spread_multiplier::Float64 = 0.07,
    random_spread_chance::Float64 = 0.0005,
) = FoodData(neighborhood, food_cap, regen_rate, spread_multiplier, random_spread_chance)

Bacterium(;
    energy::Real = 20.0,
    sensory_radius::Real = 4.0,
    reproduction_threshold::Real = 80.0,
    speed::Int = 3,
) = Bacterium(
    -1,
    (-1, -1),
    0,
    0,
    energy,
    sensory_radius,
    reproduction_threshold,
    speed,
    (-1, -1),
)

Bacterium(id::Int, pos::Dims{2}, template::Bacterium) = Bacterium(
    id,
    pos,
    template.species,
    0,
    template.energy,
    template.sensory_radius,
    template.reproduction_threshold,
    template.speed,
    (-1, -1),
)

"""
    inherit(parameter::Float64, rng::MersenneTwister, std::Float64)
    inherit(parameter::Int, rng::MersenneTwister, std::Float64)

Simulates inheritance of `parameter` with given `standard deviation` using normal distribution.
"""
inherit(parameter::Float64, rng::MersenneTwister, std::Float64) =
    max(parameter + randn(rng) * std, 1.0)
inherit(parameter::Int, rng::MersenneTwister, std::Float64) =
    max(parameter + floor(Int, randn(rng) * std), 1)

"""
function initialize_model(food::Array{Float64,2}, n_bacteria::Union{Int,Vector{Int}}; kwargs...)

Initialize a model with given food distribution and number of bacteria

Keyword arguments
-----------------
- `initial_bacterium`: `Bacterium` denoting the initial values of parameters. `Vector{Bacterium}`
  for multiple species of bacteria. In such a case, `n_bacteria` should be a vector corresponding
  to the population of each species.
- `food_data` : `FoodData` parametrising the food spread
- `lifetime` : Maximum lifetime of a bacterium
- `eat_rate_factor` : Multiplied by a bacterium's `speed` to obtain the amount of food it will
  eat in one iteration
- `sensory_radius_cost`: Multipled by a bacterium's `sensory_radius` to obtain the corresponding
  energy cost
- `distance_cost` : Multipled by the distance a bacterium moves in one iteration to obtain the
  corresponding energy cost
- `σ_sensory_radius` : Standard deviation of normal distribution describing variation in
  `sensory_radius` over generations
- `σ_reproduction_threshold` : Similar to above, for `reproduction_threshold`
- `σ_speed` : Similar to above, for `speed`
- `reproduction_energy_cost` : Energy cost of reproduction
- `random_seed` : Seed for RNG used by the model
"""
function initialize_model(
    food::Array{Float64,2};
    n_bacteria::Union{Int,Vector{Int}} = 10,
    initial_bacterium::Union{Bacterium,Vector{Bacterium}} = Bacterium(),
    food_data::FoodData = FoodData(food_cap = maximum(food)),
    lifetime::Int = 200,
    eat_rate_factor::Float64 = 4.0,
    sensory_radius_cost::Float64 = 5.0,
    distance_cost::Float64 = 5.0,
    σ_sensory_radius::Float64 = 0.5,
    σ_reproduction_threshold::Float64 = 0.5,
    σ_speed::Float64 = 0.6,
    reproduction_energy_cost::Float64 = 10.0,
    random_seed::Int = 42,
)
    # sanity checks
    @assert lifetime > 0
    @assert (
        n_bacteria isa Array &&
        initial_bacterium isa Array &&
        length(n_bacteria) == length(initial_bacterium)
    ) || (!(n_bacteria isa Array) && !(initial_bacterium isa Array))

    dims = size(food)
    if initial_bacterium isa Array
        for i = 1:length(initial_bacterium)
            initial_bacterium[i].species = i
        end
    else
        initial_bacterium.species = 1
    end

    # reproducibility matters
    rng = MersenneTwister(random_seed)

    properties = (
        nspecies = initial_bacterium isa Array ? length(initial_bacterium) : 1,
        food = food,
        food_update_order = shuffle!(rng, collect(CartesianIndices(food))),
        food_data = food_data,
        lifetime = lifetime,
        eat_rate_factor = eat_rate_factor,
        sensory_radius_cost = sensory_radius_cost,
        distance_cost = distance_cost,
        reproduction_energy_cost = reproduction_energy_cost,
        σ_sensory_radius = σ_sensory_radius,
        σ_reproduction_threshold = σ_reproduction_threshold,
        σ_speed = σ_speed,
    )

    # non-periodic euclidean space
    space = GridSpace(dims; periodic = false, metric = :euclidean)
    model = ABM(Bacterium, space; properties, rng)

    # add bacteria at random places
    if n_bacteria isa Int
        for _ = 1:n_bacteria
            add_agent!(Bacterium, model, initial_bacterium)
        end
    else
        for species = 1:properties.nspecies
            for _ = 1:n_bacteria[species]
                add_agent!(Bacterium, model, initial_bacterium[species])
            end
        end
    end

    return model
end

"""
    image_to_foodmap(path::String)

Reads image at given path, returns a matrix representing relative food concentrations,
using the red channel of RGB value at each pixel.
"""
image_to_foodmap(path::String) = map(x -> x.r, load(path))

max_stars_repo_name: jiaangou/SimpleEvolution

max_issues_repo_name: jiaangou/SimpleEvolution

max_forks_repo_name: jiaangou/SimpleEvolution


text:
module ModelingToolkit

export Operation, Expression
export calculate_jacobian, generate_jacobian, generate_function
export @register

using DiffEqBase
using StaticArrays, LinearAlgebra

using MacroTools
import MacroTools: splitdef, combinedef

abstract type Expression <: Number end
abstract type AbstractSystem end

Base.promote_rule(::Type{<:Number},::Type{<:Expression}) = Expression
Base.zero(::Type{<:Expression}) = Constant(0)
Base.one(::Type{<:Expression}) = Constant(1)

function calculate_jacobian end
function generate_jacobian end
function generate_function end

@enum FunctionVersion ArrayFunction=1 SArrayFunction=2

include("variables.jl")
include("operations.jl")
include("differentials.jl")
include("equations.jl")
include("systems/diffeqs/diffeqsystem.jl")
include("systems/diffeqs/first_order_transform.jl")
include("systems/nonlinear/nonlinear_system.jl")
include("function_registration.jl")
include("simplify.jl")
include("utils.jl")

end # module

ext: jl


text:
"""
    AbstractAtlas{𝔽}

An abstract class for atlases whith charts that have values in the vector space `𝔽ⁿ`
for some value of `n`. `𝔽` is a number system determined by an [`AbstractNumbers`](@ref)
object.
"""
abstract type AbstractAtlas{𝔽} end

"""
    RetractionAtlas{
        𝔽,
        TRetr<:AbstractRetractionMethod,
        TInvRetr<:AbstractInverseRetractionMethod,
        TBasis<:AbstractBasis,
    } <: AbstractAtlas{𝔽}

An atlas indexed by points on a manifold, such that coordinate transformations are performed
using retractions, inverse retractions, and coordinate calculation for a given basis.

# See also

[`AbstractAtlas`](@ref), [`AbstractInverseRetractionMethod`](@ref),
[`AbstractRetractionMethod`](@ref), [`AbstractBasis`](@ref)
"""
struct RetractionAtlas{
    𝔽,
    TRetr<:AbstractRetractionMethod,
    TInvRetr<:AbstractInverseRetractionMethod,
    TBasis<:AbstractBasis{<:Any,TangentSpaceType},
} <: AbstractAtlas{𝔽}
    retr::TRetr
    invretr::TInvRetr
    basis::TBasis
end

function RetractionAtlas(
    retr::AbstractRetractionMethod,
    invretr::AbstractInverseRetractionMethod,
)
    basis = DefaultOrthonormalBasis()
    return RetractionAtlas{ℝ,typeof(retr),typeof(invretr),typeof(basis)}(
        retr,
        invretr,
        basis,
    )
end
RetractionAtlas() = RetractionAtlas(ExponentialRetraction(), LogarithmicInverseRetraction())

"""
    get_default_atlas(::AbstractManifold)

Determine the default real-valued atlas for the given manifold.
"""
function get_default_atlas(::AbstractManifold)
    return RetractionAtlas()
end

"""
    get_point_coordinates(M::AbstractManifold, A::AbstractAtlas, i, p)

Calculate coordinates of point `p` on manifold `M` in chart from an [`AbstractAtlas`](@ref)
`A` at index `i`. Coordinates are in the number system determined by `A`.

# See also

[`get_point`](@ref), [`get_chart_index`](@ref)
"""
get_point_coordinates(::AbstractManifold, ::AbstractAtlas, ::Any, ::Any)

function get_point_coordinates(M::AbstractManifold, A::AbstractAtlas, i, p)
    a = allocate_result(M, get_point_coordinates, p)
    get_point_coordinates!(M, a, A, i, p)
    return a
end

function allocate_result(M::AbstractManifold, f::typeof(get_point_coordinates), p)
    T = allocate_result_type(M, f, (p,))
    return allocate(p, T, manifold_dimension(M))
end

function get_point_coordinates!(M::AbstractManifold, a, A::RetractionAtlas, i, p)
    return get_coordinates!(M, a, i, inverse_retract(M, i, p, A.invretr), A.basis)
end

function get_point_coordinates(M::AbstractManifold, A::RetractionAtlas, i, p)
    return get_coordinates(M, i, inverse_retract(M, i, p, A.invretr), A.basis)
end

"""
    get_point(M::AbstractManifold, A::AbstractAtlas, i, a)

Calculate point at coordinates `a` on manifold `M` in chart from an [`AbstractAtlas`](@ref)
`A` at index `i`.

# See also

[`get_point_coordinates`](@ref), [`get_chart_index`](@ref)
"""
get_point(::AbstractManifold, ::AbstractAtlas, ::Any, ::Any)

function get_point(M::AbstractManifold, A::AbstractAtlas, i, a)
    p = allocate_result(M, get_point, a)
    get_point!(M, p, A, i, a)
    return p
end

function allocate_result(M::AbstractManifold, f::typeof(get_point), a)
    T = allocate_result_type(M, f, (a,))
    return allocate(a, T, representation_size(M)...)
end

function get_point(M::AbstractManifold, A::RetractionAtlas, i, a)
    return retract(M, i, get_vector(M, i, a, A.basis), A.retr)
end

function get_point!(M::AbstractManifold, p, A::RetractionAtlas, i, a)
    return retract!(M, p, i, get_vector(M, i, a, A.basis), A.retr)
end

"""
    get_chart_index(M::AbstractManifold, A::AbstractAtlas, p)

Select a chart from an [`AbstractAtlas`](@ref) `A` for manifold `M` that is suitable for
representing the neighborhood of point `p`. This selection should be deterministic, although
different charts may be selected for arbitrarily close but distinct points.

# See also

[`get_default_atlas`](@ref)
"""
get_chart_index(::AbstractManifold, ::AbstractAtlas, ::Any)

get_chart_index(::AbstractManifold, ::RetractionAtlas, p) = p

"""
    transition_map(M::AbstractManifold, A_from::AbstractAtlas, i_from, A_to::AbstractAtlas, i_to, a)
    transition_map(M::AbstractManifold, A::AbstractAtlas, i_from, i_to, a)

Given coordinates `a` in chart `(A_from, i_from)` of a point on manifold `M`, returns
coordinates of that point in chart `(A_to, i_to)`. If `A_from` and `A_to` are equal, `A_to`
can be omitted.

# See also

[`AbstractAtlas`](@ref)
"""
function transition_map(
    M::AbstractManifold,
    A_from::AbstractAtlas,
    i_from,
    A_to::AbstractAtlas,
    i_to,
    a,
)
    return get_point_coordinates(M, A_to, i_to, get_point(M, A_from, i_from, a))
end

function transition_map(M::AbstractManifold, A::AbstractAtlas, i_from, i_to, a)
    return transition_map(M, A, i_from, A, i_to, a)
end

function transition_map!(
    M::AbstractManifold,
    y,
    A_from::AbstractAtlas,
    i_from,
    A_to::AbstractAtlas,
    i_to,
    a,
)
    return get_point_coordinates!(M, y, A_to, i_to, get_point(M, A_from, i_from, a))
end

function transition_map!(M::AbstractManifold, y, A::AbstractAtlas, i_from, i_to, a)
    return transition_map!(M, y, A, i_from, A, i_to, a)
end

"""
    induced_basis(M::AbstractManifold, A::AbstractAtlas, i, p, VST::VectorSpaceType)

Basis of vector space of type `VST` at point `p` from manifold `M` induced by
chart (`A`, `i`).

# See also

[`VectorSpaceType`](@ref), [`AbstractAtlas`](@ref)
"""
induced_basis(M::AbstractManifold, A::AbstractAtlas, i, VST::VectorSpaceType)

function induced_basis(
    ::AbstractManifold,
    A::RetractionAtlas{
        <:AbstractRetractionMethod,
        <:AbstractInverseRetractionMethod,
        <:DefaultOrthonormalBasis,
    },
    i,
    p,
    ::TangentSpaceType,
)
    return A.basis
end
function induced_basis(
    M::AbstractManifold,
    A::RetractionAtlas{
        <:AbstractRetractionMethod,
        <:AbstractInverseRetractionMethod,
        <:DefaultOrthonormalBasis,
    },
    i,
    p,
    ::CotangentSpaceType,
)
    return dual_basis(M, p, A.basis)
end

"""
    InducedBasis(vs::VectorSpaceType, A::AbstractAtlas, i)

The basis induced by chart with index `i` from an [`AbstractAtlas`](@ref) `A` of vector
space of type `vs`.

# See also

[`VectorSpaceType`](@ref), [`AbstractBasis`](@ref)
"""
struct InducedBasis{𝔽,VST<:VectorSpaceType,TA<:AbstractAtlas,TI} <: AbstractBasis{𝔽,VST}
    vs::VST
    A::TA
    i::TI
end

"""
    induced_basis(::AbstractManifold, A::AbstractAtlas, i, VST::VectorSpaceType)

Get the basis induced by chart with index `i` from an [`AbstractAtlas`](@ref) `A` of vector
space of type `vs`. Returns an object of type [`InducedBasis`](@ref).
"""
function induced_basis(
    ::AbstractManifold{𝔽},
    A::AbstractAtlas,
    i,
    VST::VectorSpaceType,
) where {𝔽}
    return InducedBasis{𝔽,typeof(VST),typeof(A),typeof(i)}(VST, A, i)
end

function dual_basis(
    M::AbstractManifold{𝔽},
    ::Any,
    B::InducedBasis{𝔽,TangentSpaceType},
) where {𝔽}
    return induced_basis(M, B.A, B.i, CotangentSpace)
end
function dual_basis(
    M::AbstractManifold{𝔽},
    ::Any,
    B::InducedBasis{𝔽,CotangentSpaceType},
) where {𝔽}
    return induced_basis(M, B.A, B.i, TangentSpace)
end

"""
    local_metric(M::AbstractManifold, p, B::InducedBasis)

Compute the local metric tensor for vectors expressed in terms of coordinates
in basis `B` on manifold `M`. The point `p` is not checked.
"""
local_metric(::AbstractManifold, ::Any, ::InducedBasis)


text:
# test draw pose

using Caesar, CloudGraphs, IncrementalInference

using DrakeVisualizer

vc = startdefaultvisualization()

# include(joinpath(dirname(@__FILE__),"..","database","blandauthremote.jl"))
# addrdict["session"] = "SESSTURT21"
# cloudGraph, addrdict = standardcloudgraphsetup(addrdict=addrdict)
cloudGraph, addrdict = standardcloudgraphsetup(addrdict=addrdict)

session = ddrdict["session"]

# fg = Caesar.initfg(sessionname=addrdict["session"], cloudgraph=cloudGraph)

IDs = getPoseExVertexNeoIDs(cloudGraph.neo4j.connection, session=addrdict["session"], reqbackendset=false);

cv = CloudGraphs.get_vertex(cloudGraph, IDs[10][2])
vert = cloudVertex2ExVertex(cv)

drawpose!(vc, vert, session=session, drawtype=:max )
drawposepoints!(vc, vert, session=session )

#


text:
using Test
using Dorothy

@DorothyAll()

@testset "XTC" begin

	@testset "Single" begin
		f1 = infer("$(Dorothy.datapath)/eMHC.xtc")
		@test getformat(f1) == "trajectory/x-xtc"
		f2 = infer(open("$(Dorothy.datapath)/eMHC.xtc"))
		@test getformat(f2) == "trajectory/x-xtc"
		m1 = read(f2)
		@test m1.header.step == 0
		@test m1.header.time == 0.0
		@test length(m1) == 6002
		@test volume(m1.header.cell) ≈ 665416.28
		m2 = read!(f2, MolecularModel())
		@test m2.header.step == 500e3
		@test m2.header.time == 2e3
		@test isapprox(m2.R[1], [88.65, 82.36, 47.93], atol=0.01)
		@test isapprox(m2.R[end], [72.15, 56.0, 41.75], atol=0.01)
	end

	@testset "Trajectory" begin
		streamf("$(Dorothy.datapath)/eMHC.xtc") do s
			frame = MolecularModel()
			i = 0
			while !eof(s)
				read!(s, frame)
				i += 1
				if frame.header.step == 500e3
					@test frame.header.time == 2e3
					@test isapprox(frame.R[1], [88.65, 82.36, 47.93], atol=0.01)
					@test isapprox(frame.R[end], [72.15, 56.0, 41.75],
							atol=0.01)
				end
			end
			@test i == 501
		end
	end

end # @testset


text:
### Concrete Types

type RODESolution{uType,uType2,DType,tType,IType,randType,P,A} <: AbstractRODESolution
  u::uType
  u_analytic::uType2
  errors::DType
  t::tType
  W::randType
  prob::P
  alg::A
  maxstacksize::Int
  maxstacksize2::Int
  interp::IType
  dense::Bool
  tslocation::Int
  retcode::Symbol
end
(sol::RODESolution)(t,deriv::Type=Val{0};idxs=nothing) = sol.interp(t,idxs,deriv)
(sol::RODESolution)(v,t,deriv::Type=Val{0};idxs=nothing) = sol.interp(v,t,idxs,deriv)

function build_solution{uType,tType,isinplace,ND}(
        prob::AbstractRODEProblem{uType,tType,isinplace,ND},
        alg,t,u;W=[],timeseries_errors=true,dense_errors=false,calculate_error=true,
        maxstacksize=0,maxstacksize2=0,
        interp = (tvals) -> nothing,
        dense = false, retcode = :Default, kwargs...)

  if has_analytic(prob.f)
    u_analytic = Vector{uType}(0)
    errors = Dict{Symbol,eltype(prob.u0)}()
    sol = RODESolution(u,u_analytic,errors,t,W,prob,alg,maxstacksize,maxstacksize2,interp,dense,0,retcode)
    if calculate_error
      calculate_solution_errors!(sol;timeseries_errors=timeseries_errors,dense_errors=dense_errors)
    end
    return sol
  else
    return RODESolution(u,nothing,nothing,t,W,prob,alg,maxstacksize,maxstacksize2,interp,dense,0,retcode)
  end
end

function calculate_solution_errors!(sol::AbstractRODESolution;fill_uanalytic=true,timeseries_errors=true,dense_errors=true)
  if fill_uanalytic
    for i in 1:size(sol.u,1)
      push!(sol.u_analytic,sol.prob.f(Val{:analytic},sol.t[i],sol.prob.u0,sol.W[i]))
    end
  end

  save_everystep = length(sol.u) > 2

  if !isempty(sol.u_analytic)
    sol.errors[:final] = mean(abs.(sol.u[end]-sol.u_analytic[end]))
    if save_everystep && timeseries_errors
      sol.errors[:l∞] = maximum(vecvecapply((x)->abs.(x),sol.u-sol.u_analytic))
      sol.errors[:l2] = sqrt(mean(vecvecapply((x)->float.(x).^2,sol.u-sol.u_analytic)))
    end
  end
end

function build_solution(sol::AbstractRODESolution,u_analytic,errors)
  RODESolution(sol.u,u_analytic,errors,sol.t,sol.W,sol.prob,sol.alg,sol.maxstacksize,sol.dense,sol.tslocation,sol.retcode)
end


text:
struct Variable{Name} <: AbstractVariable
end

MP.name(::Type{Variable{N}}) where {N} = N
MP.name(v::Variable) = name(typeof(v))
MP.name_base_indices(v::Variable) = name_base_indices(typeof(v))
function MP.name_base_indices(v::Type{Variable{N}}) where N
    name = string(N)
    splits = split(string(N), r"[\[,\]]\s*", keepempty=false)
    if length(splits) == 1
        return name, Int[]
    else
        return splits[1], parse.(Int, splits[2:end])
    end
end

MP.variables(v::Variable) = (v,)
MP.variables(::Type{V}) where {V <: Variable} = (V(),)
Base.hash(v::Variable, u::UInt) = hash(name(v), u)

checksorted(x::Tuple{Any}, cmp) = true
checksorted(x::Tuple{}, cmp) = true
checksorted(x::Tuple, cmp) = cmp(x[1], x[2]) && checksorted(Base.tail(x), cmp)

struct Monomial{V, N} <: AbstractMonomial
    exponents::NTuple{N, Int}

    function Monomial{V, N}(exponents::NTuple{N, Int}=ntuple(_ -> 0, Val{N}())) where {V, N}
        @assert checksorted(V, >)
        new{V, N}(exponents)
    end
    Monomial{V}(exponents::NTuple{N, Integer}=()) where {V, N} = Monomial{V, N}(exponents)
    Monomial{V}(exponents::AbstractVector{<:Integer}) where {V} = Monomial{V}(NTuple{length(V), Int}(exponents))
end

Monomial(v::Variable) = monomialtype(v)((1,))
MP.monomial(v::Variable) = Monomial(v)

MP.variables(::Type{<:Monomial{V}}) where {V} = V
MP.variables(m::Monomial) = variables(typeof(m))
MP.nvariables(::Type{<:Monomial{V}}) where {V} = length(V)
MP.nvariables(m::Monomial) = nvariables(typeof(m))
MP.monomialtype(::Type{V}) where V<:Variable = monomialtype(V())
MP.monomialtype(v::Variable) = Monomial{(v,), 1}

MP.exponents(m::Monomial) = m.exponents
MP.exponent(m::Monomial, i::Integer) = m.exponents[i]
_exponent(v::V, p1::Tuple{V, Integer}, p2...) where {V <: Variable} = p1[2]
_exponent(v::Variable, p1::Tuple{Variable, Integer}, p2...) = _exponent(v, p2...)
_exponent(v::Variable) = 0
MP.exponent(m::Monomial, v::Variable) = _exponent(v, powers(m)...)

struct Term{CoeffType, M <: Monomial} <: AbstractTerm{CoeffType}
    coefficient::CoeffType
    monomial::M
end
Term(v::Variable) = Term(Monomial(v))
Term(x, v::Variable) = Term(x, Monomial(v))
Term(m::Monomial) = Term(1, m)

MP.coefficient(t::Term) = t.coefficient
MP.monomial(t::Term) = t.monomial
coefftype(::Type{<:Term{C}}) where {C} = C
MP.termtype(::Union{Term{C, M}, Type{<:Term{C, M}}}, ::Type{T}) where {C, M, T} = Term{T, M}
MP.termtype(::Union{M, Type{M}}, ::Type{T}) where {M<:Monomial, T} = Term{T, M}
MP.monomialtype(::Type{<:Term{C, M}}) where {C, M} = M
MP.monomialtype(t::Term) = monomialtype(typeof(t))
MP.variables(t::Union{Term, Type{<:Term}}) = variables(monomialtype(t))
MP.nvariables(t::Union{Term, Type{<:Term}}) = nvariables(monomialtype(t))

struct Polynomial{CoeffType, T <: Term{CoeffType}, V <: AbstractVector{T}} <: AbstractPolynomial{CoeffType}
    terms::V

    Polynomial{C, T, V}(terms::AbstractVector{T}) where {C, T, V} = new{C, T, V}(terms)
    Polynomial{C, T, V}(p::Polynomial{C, T, V}) where {C, T, V} = p
end
Polynomial(terms::AbstractVector{T}) where {C, T <: Term{C}} = Polynomial{C, T, typeof(terms)}(terms)
Polynomial(t::AbstractVector) = Polynomial(Term.(t))
Polynomial(x) = Polynomial(Term(x))
MP.termtype(::Type{<:Polynomial{C, T}}) where {C, T} = T
changeeltype(::Type{<:Vector}, ::Type{T}) where T = Vector{T}
function MP.polynomialtype(::Type{<:Polynomial{C, T, V}}, ::Type{NewC}) where {C, T, V, NewC}
    NewT = termtype(T, NewC)
    Polynomial{NewC, NewT, changeeltype(V, NewT)}
end
MP.polynomialtype(::Type{Term{C, M}}) where {C, M} = Polynomial{C, Term{C, M}, Vector{Term{C, M}}}
Polynomial(term::TT) where TT<:Term = Polynomial(iszero(term) ? TT[] : [term])

MP.terms(p::Polynomial) = p.terms
MP.variables(::Union{Polynomial{C, T}, AbstractArray{<:Polynomial{C, T}}, Type{<:Polynomial{C, T}}}) where {C, T} = variables(T)
MP.nvariables(::Union{Polynomial{C, T}, AbstractArray{<:Polynomial{C, T}}, Type{<:Polynomial{C, T}}}) where {V, N, C, M<:Monomial{V, N}, T<:Term{C, M}} = N
const MonomialLike = Union{Variable, Monomial}
const TermLike = Union{MonomialLike, Term}
const PolynomialLike = Union{TermLike, Polynomial}

MP.variables(::Union{AbstractVector{T}, Type{<:AbstractVector{T}}}) where {T <: PolynomialLike} = variables(T)
MP.nvariables(::Union{AbstractVector{T}, Type{<:AbstractVector{T}}}) where {T <: PolynomialLike} = nvariables(T)
MP.variable_union_type(::Union{PolynomialLike, Type{<:PolynomialLike}}) = Variable
MP.constantmonomial(p::PolynomialLike) = Monomial{variables(p), nvariables(p)}()
MP.constantmonomial(::Type{TT}) where {TT<:PolynomialLike} = Monomial{variables(TT), nvariables(TT)}()

# Based on fillZfordeg!() from MultivariatePolynomials.jl by Benoit Legat
# https://github.com/blegat/MultivariatePolynomials.jl/blob/d85ad85de413afa20fc8f5354c980387218ced2c/src/mono.jl#L186-L259
function monomial_powers(::Val{N}, degree) where N
    result = Vector{NTuple{N, Int}}()
    powers = zeros(Int, N)
    powers[1] = degree
    while true
        push!(result, NTuple{N, Int}(powers))
        if powers[end] == degree
            break
        end
        total = 1
        for j in (N - 1):-1:1
            if powers[j] != 0
                powers[j] -= 1
                powers[j+1] += total
                break
            else
                total += powers[j+1]
                powers[j+1] = 0
            end
        end
    end
    result
end

function MP.monomials(vars::Tuple{Vararg{<:Variable}}, degree::Integer, filter::Function=m->true)
    checksorted(vars, >) || throw(ArgumentError("Variables must be in order"))
    Monomial{vars, length(vars)}[Monomial{vars}(p) for p in monomial_powers(Val{length(vars)}(), degree) if filter(Monomial{vars}(p))]
end

function MP.monomials(vars::Tuple{Vararg{<:Variable}}, degrees::AbstractArray, filter::Function=m->true)
    checksorted(vars, >) || throw(ArgumentError("Variables must be in order"))
    if isempty(degrees)
        # Otherwise, the following error is thrown: "ArgumentError: argument to Flatten must contain at least one iterator"
        Monomial{vars, length(vars)}[]
    else
        Monomial{vars, length(vars)}[Monomial{vars}(p) for d in sort(degrees, rev=true)
            for p in monomial_powers(Val{length(vars)}(), d) if filter(Monomial{vars}(p))]
    end
end

MP.similarvariable(::Union{PolynomialLike, Type{<:PolynomialLike}}, ::Type{Val{N}}) where N = Variable{N}()


text:
# v.2.2

# CHANGE LOG
# 2.0: negativesense() function depreciated and merged into parseIO
# 2.1: Sequences won't parse

# REQUIRED FILE STRUCTURE
# .
# ├── /reference       	Directory containing blast search subject sequences in fasta format
# ├── /blast_out       	Directory containing *_out.txt blast output in outfmt6 format
# └── blast_extract.jl

# INFO
# Updated README to come..

using Dates
function logIO(str::AbstractString)
	if @isdefined log
		println(log, str, now())
		flush(log)
	else
		global log = open(ARGS[1]*"/blast_extract.log", "w")
		println(log, "log file for blast_extract.jl - ", now())
		flush(log)
	end
end

using JuliaDB
logIO("Initiating referenceIO function - ")
function referenceIO(dir::AbstractString = ".") #DONE
	global reference = Dict{String, Any}()
	for file in readdir("$dir/reference")
		temp_counter = 0
		file_name = SubString(file, 1, Integer(findfirst('.', file))-1)
		logIO("Starting referenceIO for $file_name - ")
		temp_dict = Dict{String, String}()
		temp_key = ""
		temp_value = ""
		for line in eachline("$dir/reference/$file")
			if temp_counter <= 6 #ONLY LOADS FIRST 7 SCAFFOLDS
				if startswith(line, ">") === true
					if findfirst(' ', line) === nothing
						temp_key = SubString(line, 2:length(line))
					else
						temp_index = findfirst(' ', line)
						temp_key = SubString(line, 2:(temp_index-1))
					end
					temp_counter = temp_counter + 1
				else
					temp_value = temp_value * line
				end
			else
				break
			end
			merge!(temp_dict, Dict{String, String}(temp_key => temp_value))
		end
		merge!(reference, Dict(file_name => temp_dict))
	end
	return reference
end
logIO("referenceIO function initiated! - ")

logIO("Initiating blastIO function - ")
function blastIO(dir::AbstractString = ".") #DONE
	global blast_out = Dict{String, IndexedTable}()
	for file in readdir("$dir/blast_out")
		file_name = SubString(file, 1, Integer(findfirst('.', file))-1)
		logIO("Starting blastIO for $file_name - ")
		merge!(blast_out, Dict(file_name => loadtable("$dir/blast_out/$file", spacedelim = true, header_exists = false,
			colnames = ["qseqid", "sseqid", "pident", "length", "mismatch", "gapopen", "qstart", "qend", "sstart", "send", "evalue", "bitscore"])))
		blast_out[file_name] = insertcols(blast_out[file_name], 1, :index => range(1, length(blast_out[file_name])))
	end
	return blast_out
end
logIO("blastIO function initiated! - ")

logIO("Initiating sseqid_merge function - ")
function sseqid_merge(base::Integer) #DONE
	global sseqid_parse = Dict{String, Any}()
	temp_counter = 0
	if length(keys(reference)) == 1 #FOR ONE REFERENCE AND MULTIPLE BLAST_OUTS
		for key in keys(reference)
			for n in 1:length(keys(blast_out))
				key_bridge = SubString(key, 1, Integer(findfirst('_', key))-1)
				temp_dict = Dict{String, Vector}()
				for sid in Set(select(blast_out["$key_bridge"*"$n"*"_out"], :sseqid))
					temp_table_s = filter(val -> (val.sseqid == sid), blast_out["$key_bridge"*"$n"*"_out"])
					for qid in Set(select(blast_out["$key_bridge"*"$n"*"_out"], :qseqid))
						temp_table_sq = filter(val -> (val.qseqid == qid), temp_table_s)
						temp_table_sq = sort(temp_table_sq, :send)
						if length(temp_table_sq) == 1
							if !haskey(temp_dict, sid)
								merge!(temp_dict, Dict{String, Array}(sid => [temp_table_sq[1][:sstart], temp_table_sq[1][:send]]))
							else
								append!(temp_dict[sid], [temp_table_sq[1][:sstart], temp_table_sq[1][:send]])
							end
						elseif length(temp_table_sq) == 2
							if temp_table_sq[1][:send]+base >= temp_table_sq[2][:sstart]
								if !haskey(temp_dict, sid)
									merge!(temp_dict, Dict{String, Array}(sid => [temp_table_sq[1][:sstart], temp_table_sq[2][:send]]))
								else
									append!(temp_dict[sid], [temp_table_sq[1][:sstart], temp_table_sq[2][:send]])
								end
							else
								if !haskey(temp_dict, sid)
									merge!(temp_dict, Dict{String, Array}(sid => [temp_table_sq[1][:sstart], temp_table_sq[1][:send]]))
								else
									append!(temp_dict[sid], [temp_table_sq[1][:sstart], temp_table_sq[1][:send]])
								end
							end
						elseif length(temp_table_sq) > 2
							temp_send = nothing
							for i in range(1, length(temp_table_sq))
								if i+1 > length(temp_table_sq)
									if temp_send === nothing
										if !haskey(temp_dict, sid)
											merge!(temp_dict, Dict{String, Array}(sid => [temp_table_sq[i][:sstart], temp_table_sq[i][:send]]))
										else
											append!(temp_dict[sid], [temp_table_sq[i][:sstart], temp_table_sq[i][:send]])
										end
									else
										if !haskey(temp_dict, sid)
											try
												merge!(temp_dict, Dict{String, Array}(sid => [temp_table_sq[i-temp_counter][:sstart], temp_table_sq[i][:send]]))
												temp_send = nothing
												temp_counter = 0
											catch
												merge!(temp_dict, Dict{String, Array}(sid => [temp_table_sq[i-temp_counter][:sstart], temp_table_sq[i][:send]]))
												temp_send = nothing
												temp_counter = 0
											end
										else
											try
												append!(temp_dict[sid], [temp_table_sq[i-temp_counter][:sstart], temp_send])
												temp_send = nothing
												temp_counter = 0
											catch
												append!(temp_dict[sid], [temp_table_sq[1][:sstart], temp_send])
												temp_send = nothing
												temp_counter = 0
											end
										end
									end
								elseif temp_table_sq[i][:send]+base >= temp_table_sq[i+1][:sstart]
									temp_send = temp_table_sq[i+1][:send]
									temp_counter += 1
								elseif temp_table_sq[i][:send]+base < temp_table_sq[i+1][:sstart]
									if temp_send === nothing
										if !haskey(temp_dict, sid)
											merge!(temp_dict, Dict{String, Array}(sid => [temp_table_sq[i][:sstart], temp_table_sq[i][:send]]))
										else
											append!(temp_dict[sid], [temp_table_sq[i][:sstart], temp_table_sq[i][:send]])
										end
									else
										try
											append!(temp_dict[sid], [temp_table_sq[i-temp_counter][:sstart], temp_send])
											temp_send = nothing
											temp_counter = 0
										catch
											append!(temp_dict[sid], [temp_table_sq[1][:sstart], temp_send])
											temp_send = nothing
											temp_counter = 0
										end
									end
								end
							end
						end
					end
				end
				merge!(sseqid_parse, Dict(key => temp_dict))
			end
		end
	else #MULTIPLE REFERENCE AND MULTIPLE BLAST_OUTS
		for key in keys(reference)
			temp_dict = Dict{String, Vector}()
			key_bridge = SubString(key, 1, Integer(findfirst('_', key))-1)
			for sid in Set(select(blast_out["$key_bridge"*"$n"*"_out"], :sseqid))
				temp_table_s = filter(val -> (val.sseqid == sid), blast_out["$key_bridge"*"$n"*"_out"])
				for qid in Set(select(blast_out["$key_bridge"*"$n"*"_out"], :qseqid))
					temp_table_sq = filter(val -> (val.qseqid == qid), temp_table_s)
					temp_table_sq = sort(temp_table_sq, :send)
					if length(temp_table_sq) == 1
						if !haskey(temp_dict, sid)
							merge!(temp_dict, Dict{String, Array}(sid => [temp_table_sq[1][:sstart], temp_table_sq[1][:send]]))
						else
							append!(temp_dict[sid], [temp_table_sq[1][:sstart], temp_table_sq[1][:send]])
						end
					elseif length(temp_table_sq) == 2
						if temp_table_sq[1][:send]+base >= temp_table_sq[2][:sstart]
							if !haskey(temp_dict, sid)
								merge!(temp_dict, Dict{String, Array}(sid => [temp_table_sq[1][:sstart], temp_table_sq[2][:send]]))
							else
								append!(temp_dict[sid], [temp_table_sq[1][:sstart], temp_table_sq[2][:send]])
							end
						else
							if !haskey(temp_dict, sid)
								merge!(temp_dict, Dict{String, Array}(sid => [temp_table_sq[1][:sstart], temp_table_sq[1][:send]]))
							else
								append!(temp_dict[sid], [temp_table_sq[1][:sstart], temp_table_sq[1][:send]])
							end
						end
					elseif length(temp_table_sq) > 2
						println(base)
						temp_send = nothing
						for i in range(1, length(temp_table_sq))
							if i+1 > length(temp_table_sq)
								if temp_send === nothing
									if !haskey(temp_dict, sid)
										merge!(temp_dict, Dict{String, Array}(sid => [temp_table_sq[i][:sstart], temp_table_sq[i][:send]]))
									else
										append!(temp_dict[sid], [temp_table_sq[i][:sstart], temp_table_sq[i][:send]])
									end
								else
									if !haskey(temp_dict, sid)
										try
											merge!(temp_dict, Dict{String, Array}(sid => [temp_table_sq[i-temp_counter][:sstart], temp_table_sq[i][:send]]))
											temp_send = nothing
											temp_counter = 0
										catch
											merge!(temp_dict, Dict{String, Array}(sid => [temp_table_sq[i-temp_counter][:sstart], temp_table_sq[i][:send]]))
											temp_send = nothing
											temp_counter = 0
										end
									else
										try
											append!(temp_dict[sid], [temp_table_sq[i-temp_counter][:sstart], temp_send])
											temp_send = nothing
											temp_counter = 0
										catch
											append!(temp_dict[sid], [temp_table_sq[1][:sstart], temp_send])
											temp_send = nothing
											temp_counter = 0
										end
									end
								end
							elseif temp_table_sq[i][:send]+base >= temp_table_sq[i+1][:sstart]
								temp_send = temp_table_sq[i+1][:send]
								temp_counter += 1
							elseif temp_table_sq[i][:send]+base < temp_table_sq[i+1][:sstart]
								if temp_send === nothing
									if !haskey(temp_dict, sid)
										merge!(temp_dict, Dict{String, Array}(sid => [temp_table_sq[i][:sstart], temp_table_sq[i][:send]]))
									else
										append!(temp_dict[sid], [temp_table_sq[i][:sstart], temp_table_sq[i][:send]])
									end
								else
									try
										append!(temp_dict[sid], [temp_table_sq[i-temp_counter][:sstart], temp_send])
										temp_send = nothing
										temp_counter = 0
									catch
										append!(temp_dict[sid], [temp_table_sq[1][:sstart], temp_send])
										temp_send = nothing
										temp_counter = 0
									end
								end
							end
						end
					end
				end
			end
			merge!(sseqid_parse, Dict(key => temp_dict))
		end
	end
	return sseqid_parse
end
logIO("sseqid_merge function initiated! - ")

logIO("Initiating parseIO function - ")
using IterTools
function parseIO() #DONE
	for file in keys(sseqid_parse)
		logIO("Parsing sequences for hits found in $file - ")
		file_bridge = SubString(file, 1, Integer(findfirst('_', file))-1)
		io = open("$file_bridge"*"_parsed.fa", "w")
		for id in keys(sseqid_parse[file])
			if match(r"^Chromosome", id) !== nothing
				for (i, (sstart, send)) in enumerate(partition(sseqid_parse[file][id], 2, 2))
					if sstart > send
						temp_substring = SubString(reference[file][id], send:sstart)
						write(io, ">$id"*"_$i\n") + write(io, "$temp_substring\n")
					else
						temp_substring = SubString(reference[file][id], sstart:send)
						write(io, ">$id"*"_$i\n") + write(io, "$temp_substring\n")
					end
				end
			else
				continue
			end
		end
		close(io)
	end
end
logIO("parseIO function initiated! - ")

referenceIO(ARGS[1]) #upload subject sequences
logIO("referenceIO completed! - ")
println(log, reference["test_subject"]["Chromosome1"]) #sanity check # works
blastIO(ARGS[1]) #upload blast output
logIO("blastIO completed! - ")
println(log, blast_out["test1_out"][1]) #sanity check #works
sseqid_merge(10) #merge hits within n basepairs and extract hits
logIO("Successfully merged hits within the specified number of basepairs! - ")
println(log, sseqid_parse["test_subject"]) #sanity check #Need nested dictionary for each blast out
parseIO() #output extracted hits
close(log)

max_stars_repo_name: RDRidenbaugh/odds-n-ends

max_issues_repo_name: RDRidenbaugh/odds-n-ends

max_forks_repo_name: RDRidenbaugh/odds-n-ends


text:
using Test, TestSetExtensions
using AssociatedLegendrePolynomials # For doctests, include AssociatedLegendrePolynomials as a binding in Main

include("testsuite.jl")
using .TestSuite: LMAX, NumTypes
const Norms = (LegendreUnitNorm(), LegendreUnitCoeff{BigFloat}(LMAX))

function prettytime(t)
    v, u = t < 1e3 ? (t, "ns") :
           t < 1e6 ? (t/1e3, "μs") :
           t < 1e9 ? (t/1e6, "ms") :
                     (t/1e9, "s")
    return string(round(v, digits=3), " ", u)
end
macro include(file, desc)
    mod = gensym(first(splitext(file)))
    quote
        print($desc, ": ")
        t0 = time_ns()
        @testset $desc begin
            @eval module $mod
                using Test, AssociatedLegendrePolynomials
                import ..NumTypes
                include($file)
            end
        end
        printstyled("  ", prettytime(time_ns() - t0), "\n", color=:light_black)
    end
end

@testset ExtendedTestSet "AssociatedLegendrePolynomials" begin
    @include "norms.jl" "Normalizations"
    @include "coeffs.jl" "Precomputed coefficients"
    @include "implementation.jl" "Implementation details"
    @include "analytic.jl" "Analytic checks"
    @include "errors.jl" "Error checking"
    @include "scalar.jl" "Utility: scalar container"
    @include "doctests.jl" "Doctests"
end


text:
using DuckDB
using Documenter

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

makedocs(;
    modules = [DuckDB],
    authors = "Kimmo Linna <kimmo.linna@gmail.com> and contributors",
    repo = "https://github.com/kimmolinna/DuckDB.jl/blob/{commit}{path}#{line}",
    sitename = "DuckDB.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://kimmolinna.github.io/DuckDB.jl",
        assets = String[],
    ),
    pages = ["Home" => "index.md"],
)

deploydocs(; repo = "github.com/kimmolinna/DuckDB.jl", devbranch = "main")


text:
import Pkg

function package_list()::Vector{String}
    project_toml::String = package_directory("Project.toml",)
    pkg_list::Vector{String} = sort(
        unique(
            strip.(
                collect(
                    keys(
                        Pkg.TOML.parsefile(project_toml)["deps"]
                        )
                    )
                )
            )
        )
    return pkg_list
end

function print_list_of_package_imports(io::IO = stdout)::Nothing
    pkg_list::Vector{String} = sort(unique(package_list()))
    println("##### Beginning of file")
    println()
    println("# Automatically generated by: ")
    println("# PredictMD.print_list_of_package_imports()")
    println()
    for i = 1:length(pkg_list)
        println(io, "import $(pkg_list[i])",)
    end
    println()
    println("##### End of file")
    println()
    return nothing
end


text:
# MIT License

#Copyright (c) 2021 hondoRandale <jules.rasetaharison@tutanota.com> and contributors

#Permission is hereby granted, free of charge, to any person obtaining a copy
#of this software and associated documentation files (the "Software"), to deal
#in the Software without restriction, including without limitation the rights
#to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
#copies of the Software, and to permit persons to whom the Software is
#furnished to do so, subject to the following conditions:

#The above copyright notice and this permission notice shall be included in all
#copies or substantial portions of the Software.

#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
#OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
#SOFTWARE.

#######################
## 2D-test functions ##
#######################

## Beale function
Beale_2D(x;kwargs) = ( 1.5f0 - x[1] + x[1]*x[2] )^2 + ( 2.25f0 - x[1] + x[1]*x[2]^2 )^2 + ( 2.625f0 - x[1] + x[1]*x[2]^3 )^2

## Himmelblau's function
Himmelblau_2D(x;kwargs) = ( x[1]^2 + x[2] - 11.0f0 )^2 + ( x[1] + x[2]^2 - 7.0f0  )^2

## Lévi function N.13
Levi_no13_2D(x;kwargs) = sin( 3.0f0 * π * x[1] )^2 + ( x[1] - 1.0f0 )^2 * ( 1.0f0 + sin( 1.0f0+3.0f0*π*x[2] )^2 ) +
                                 ( x[2] - 1.0f0 )^2 * ( 1.0f0 + sin( 2.0f0*π*x[2] )^2 )

## Rosenbrock_function
Rosenbrock_2D(x;kwargs) = ( 1.0f0 - x[1] )^2 + 100.0f0 * (x[2] - x[1]^2)^2

## Booth_function
Booth_2D(x;kwargs) = ( x[1] + 2.0f0*x[2] - 7.0f0 )^2 +
                     ( 2.0f0*x[1] + x[2] - 5.0f0 )^2

## Ackley_function
Ackley_2D(x;kwargs) = -20.0f0 * exp( -0.2f0 * sqrt( 0.5f0 * ( x[1]^2 + x[2]^2 ) )  ) -
                      exp( 0.5f0 * ( cos( 2.0f0 * π * x[1] ) + cos( 2.0f0 * π * x[2] ) ) ) +
                      exp( 1.0f0 ) +
                      20.0f0

## Matyas_function
Matyas_2D(x;kwargs) = 0.26f0 * ( x[1]^2 + x[2]^2  ) - 0.48f0 * x[1] * x[2]

## Three_hump_camel_function
Three_hump_camel_2D(x;kwargs) = 2.0f0 * x[1]^2 - 1.05f0 * x[1]^4 +
                                x[1]^6 / 6 + x[1] * x[2] + x[2]^2

## Easom_function
Easom_2D(x;kwargs) = -1.0f0 * cos( x[1] ) * cos( x[2] ) * exp( -1.0f0 * ( (x[1]-π)^2 + (x[2]-π)^2 ) )

## Sphere_function
Sphere_2D(x;kwargs) = x[1]^2 + x[2]^2

## Schaffer_function_no2
Schaffer_function_no2_2D(x;kwargs) = 0.5f0 + ( sin( x[1]^2 - x[2]^2 )^2 - 0.5f0 ) /
                                             ( 1.0f0 + 0.001f0*( x[1]^2 + x[2]^2 ) )^2

## Bukin_function_no6
Bukin_function_no6_2D(x;kwargs) = 100.0f0 * sqrt( abs( x[2] - 0.01f0 * x[1]^2 ) ) + 0.01f0*abs( x[1] + 10.0f0 )

obj_functions_2D = [ Beale_2D,
                     Himmelblau_2D, Levi_no13_2D,
                     Ackley_2D, Matyas_2D, Booth_2D,
                     Three_hump_camel_2D, Easom_2D,
                     Sphere_2D, Rosenbrock_2D,
                     Schaffer_function_no2_2D,
                     Bukin_function_no6_2D ];

#######################
## ND-test functions ##
#######################

## Rastrigin function
Rastrigin_ND(x;kwargs) = 10.0f0 * Float32( length( x ) ) + sum( ( x.^2 ) .- 10.0f0 * cos.( 2.0f0*π*x )  )

## Sphere function
Sphere_ND(x;kwargs) = sum( x.^2 )

## Styblinski–Tang function
Styblinski_Tang_ND(x;kwargs) = 0.5f0 * sum( (x).^4 .- 16.0f0*(x).^2 .+ 5.0f0.*x )

## Axis Parallel Hyper-Ellipsoid function
AxisParallelHyperEll_ND(x;kwargs) = sum( Float32.( 1:1:length( x ) ) .* ( x.^2 ) )

## Zakharov Function
Zakharov_ND(x;kwargs) = sum( x.^2 ) + sum( 0.5f0 * Float32.( 1:1:length( x ) ) .* x )^2 + sum( 0.5f0 * Float32.( 1:1:length( x ) ) .* x )^4

## Michalewicz function
Michalewicz_ND(x;kwargs) = -sum( sin.( x ) .* sin.( ( Float32.( 1:1:length( x ) ) .* ( x.^2 ) ) / π ) )

## Exponential function
Exponential_ND(x;kwargs) = -exp( -0.5f0 * sum( x.^2 ) )

## Griewank function
Griewank_ND(x;kwargs) = 1.0f0 + sum( (x).^2 / 4000.0f0 ) - prod( cos.( x ./ sqrt.( Float32.( 1:1:length( x ) ) ) ) )

## Qing function
Qing_ND(x;kwargs) = sum( ( (x).^2 .- Float32.( 1:1:length( x ) ) ).^2 )

## Salomon function
Salomon_ND(x;kwargs) = 1.0f0 - cos( 2*π*sqrt( sum( (x).^2 ) ) ) + 0.1f0 * sqrt( sum( (x).^2 ) )

## Schwefel function
Schwefel_ND(x;kwargs) = 418.9829f0 * length( x ) - sum( x .* sin.( sqrt.( abs.( x ) ) ) )

## Schwefel 2.20 function
Schwefel_2_2_0_ND(x;kwargs) = sum( abs.( x ) )

## Schwefel 2.21 function
Schwefel_2_2_1_ND(x;kwargs) = maximum( abs.( x ) )

## Schwefel 2.22 function
Schwefel_2_2_2_ND(x;kwargs) = sum( abs.( x ) ) + prod( abs.( x ) )

## Schwefel 2.23 function
Schwefel_2_2_3_ND(x;kwargs) = sum( (x).^10 )

obj_functions_Nd = [ Rastrigin_ND, Sphere_ND,
                     Styblinski_Tang_ND, AxisParallelHyperEll_ND,
                     Zakharov_ND, Michalewicz_ND,
                     Exponential_ND, Griewank_ND,
                     Qing_ND, Salomon_ND,
                     Schwefel_ND, Schwefel_2_2_0_ND,
                     Schwefel_2_2_1_ND, Schwefel_2_2_2_ND,
                     Schwefel_2_2_3_ND ]


text:
# benchmarking
# [extract_rotor_info, build_operating_points, select_rotor, solve_ccblade, nondimensionalize]
labels = ["(x10^5) extract rotor info", "(x10^3) build_operating_points", "(x10^4) select_rotor", "(x10^0) solve_ccblade", "(x10^3) nondimensionalize"]
kronos_epema_rotor_sweep = [8.53e-2, 4.2694e-2, 1.6446e-2, 0.034703221, 9.54845e-2]
kronos_epema_blown_wing = [9.96e-2, 4.9103e-2, 6.97e-3, 0.245689797, 7.3299e-2]
kronos_epema_blown_wing_simplified = [1.136e-1, 4.1276e-2, 1.371e-2, 0.068328487, 7.6982e-2]
# kronos_epema_rotor_sweep = [8.53e-7, 4.2694e-5, 1.644e-6, 0.034703221, 9.5484e-5]
# kronos_epema_blown_wing = [9.96e-7, 4.9103e-5, 6.9e-7, 0.245689797, 7.3299e-5]
# kronos_epema_blown_wing_simplified = [1.136e-6, 4.1276e-5, 1.371e-6, 0.068328487, 7.6982e-5]
kronos = [kronos_epema_rotor_sweep, kronos_epema_blown_wing, kronos_epema_blown_wing_simplified]
tick_labels = ["epema_rotor_sweep", "epema_blown_wing", "epema_blown_wing_simplified"]

benchmark_fig = plt.figure("benchmark_ccblade")
benchmark_fig.clear()
benchmark_ax = benchmark_fig.add_subplot(1,1,1)
towers = []
for (i,time) in enumerate(kronos)
    local xs = ones(length(time)) .* i
    local heights = time
    # heights = deepcopy(time)
    local widths = ones(length(time)) .* 0.4
    local ys = [vcat(0,heights)[i] for i = 1:length(heights)]
    for j in 1:length(xs)
        local cratio = (j-1)/(length(xs)-1)
        if j == 1
            benchmark_ax.bar(xs[j], heights[j], widths[j], ys[j], tick_label=tick_labels[i], label = i == 1 ? labels[j] : nothing, color=(0.05, 0.85-cratio*0.7, 0.15 + 0.75 * cratio))
        else
            benchmark_ax.bar(xs[j], heights[j], widths[j], ys[j], label = i == 1 ? labels[j] : nothing, color=(0.05, 0.85-cratio*0.7, 0.15 + 0.75 * cratio))
        end
    end
end

benchmark_ax.set_ylabel("time [s]")
benchmark_ax.set_xticks(1:length(kronos))
benchmark_ax.set_xticklabels(["epema rotor sweep", "epema blown wing", "epema simplified blown wing"])
# benchmark_ax.set_xticklabels(vcat([tick_labels[i] for i in 1:length(kronos) for j in 1:length(kronos_epema_blown_wing)]))
benchmark_ax.legend()
benchmark_fig.tight_layout()

benchmark_fig.savefig(joinpath(plot_directory, "benchmark_ccblade.png"), bbox_inches="tight")

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
struct Integrator{T}
    t::Array{T, 1}
    q::Array{Array{T, 1}, 1}
    u::Array{Array{T, 1}, 1}
    Λ::Array{Array{T, 1}, 1}
    m::Moreau
    Δt::T
    # extradynamics::Function
    extra_state::Array{Array{T, 1}, 1}
end

function Integrator(gap::Function, dynamics::Function, q0::AbstractArray,
        u0::AbstractArray; Δt::S=1e-3, ε::S=0.5) where {S <: Real}
    t = Array{S, 1}()
    q = Array{Array{S, 1},1}()
    u = Array{Array{S, 1},1}()
    Λ = Array{Array{S, 1},1}()
    extra_state = Array{Array{S, 1},1}()
    push!(t, 0.0)
    push!(q, q0)
    push!(u, u0)
    m = Moreau(gap, dynamics, q[1], u[1], Δt, ε)
    set_state(m, q0, u0)

    Integrator{S}(t, q, u, Λ, m, Δt, extra_state)
end

function Integrator(gap::Function, dynamics::Function, hcon::Function,
        jac::Function, jacdot::Function, q0::AbstractArray, u0::AbstractArray;
        Δt::S=1e-3, ε::S=0.5) where {S <: Real}
    t = Array{S, 1}()
    q = Array{Array{S, 1},1}()
    u = Array{Array{S, 1},1}()
    Λ = Array{Array{S, 1},1}()
    extra_state = Array{Array{S, 1},1}()
    push!(t, 0.0)
    push!(q, q0)
    push!(u, u0)
    m = Moreau(gap, dynamics, hcon, jac, jacdot, q[1], u[1], Δt, ε)
    set_state(m, q0, u0)

    Integrator{S}(t, q, u, Λ, m, Δt, extra_state)
end

function Integrator(gap::Function, dynamics::Function, q0::AbstractArray,
    u0::AbstractArray, extra_state0::AbstractArray; Δt::S=1e-3, ε::S=0.1) where {S <: Real}
    t = Array{S, 1}()
    q = Array{Array{S, 1},1}()
    u = Array{Array{S, 1},1}()
    Λ = Array{Array{S, 1},1}()
    extra_state = Array{Array{S, 1},1}()
    push!(t, 0.0)
    push!(q, q0)
    push!(u, u0)
    push!(extra_state, extra_state0)
    m = Moreau(gap, dynamics, q[1], u[1], extra_state[1], Δt, ε)
    set_state(m, q0, u0, extra_state0)

    Integrator{S}(t, q, u, Λ, m, Δt, extra_state)
end

function integrate(system::Integrator, final_time::S) where {S <: Real}
    for time in range(system.t[1]; step=system.Δt, stop=final_time)
        push!(system.t, time)
        step(system.m)
        push!(system.q, system.m.qE)
        push!(system.u, system.m.uE)
        push!(system.Λ, system.m.Λ)
        set_state(system.m, system.q[end], system.u[end])
    end
end

function integrate(system::Integrator, extradynamics::Function, final_time::S) where {S <: Real}
    for time in range(system.t[1]; step=system.Δt, stop=final_time)
        push!(system.t, time)
        step(system.m)

        next_extra_state = system.extra_state[end] +
            system.Δt * extradynamics(system.extra_state[end], system.q[end], system.u[end])

        # The extra state is assumed to be current. If not, remove this saturation.
        current_limit = 10.0
        push!(system.extra_state, clamp.(next_extra_state, -current_limit, current_limit))

        push!(system.q, system.m.qE)
        push!(system.u, system.m.uE)
        push!(system.Λ, system.m.Λ)
        set_state(system.m, system.q[end], system.u[end], system.extra_state[end])
    end
end


text:
# Benchmarks to compare the performance of
#
#   - devec macro (in Devectorize)
#   - hand-coded devectorized for-loop
#   - vectorized Julia code
#

using Devectorize

# benchmark helpers

repeat = 20

macro my_bench(Name)
    quote
        # create task
        println("bench: ", $(string(Name)))
        task = ($Name)()

        vec_eval(task, a, b, c)
        t1 = @elapsed for i = 1 : repeat
            vec_eval(task, a, b, c)
        end
        @printf "    vec_eval   : %8.4f sec  |  gain = %7.4f\n" t1 t1 / t1

        devec_eval(task, a, b, c)
        t2 = @elapsed for i = 1 : repeat
            devec_eval(task, a, b, c)
        end
        @printf "    devec_eval : %8.4f sec  |  gain = %7.4f\n" t2 t1 / t2

        hand_loop(task, a, b, c)
        t3 = @elapsed for i = 1 : repeat
            hand_loop(task, a, b, c)
        end
        @printf "    hand_loop  : %8.4f sec  |  gain = %7.4f\n" t3 t1 / t3

        println()
    end
end

# prepare data

m = 1000
n = 1000
a = zeros(m, n)
b = zeros(m, n)
c = zeros(m, n)

# simple element-wise

type simple_ewise end
formula(::simple_ewise) = "(a - b).^2 + c"

function vec_eval{T<:Real}(::simple_ewise, a::Array{T}, b::Array{T}, c::Array{T})
    r = sqr(a - b) + c   # Note: Devectorize defines a sqr on arrays
end

function devec_eval{T<:Real}(::simple_ewise, a::Array{T}, b::Array{T}, c::Array{T})
    @devec r = sqr(a - b) + c
end

function hand_loop{T<:Real}(::simple_ewise, a::Array{T}, b::Array{T}, c::Array{T})
    r = similar(a)
    for i = 1 : length(a)
        v = (a[i] - b[i])
        r[i] = v * v + c[i]
    end
end

# complex element-wise

type complex_ewise end
formula(::complex_ewise) = "log(exp((a - b).^2) + exp(a + b)) - c .* log(c)"

function vec_eval{T<:Real}(::complex_ewise, a::Array{T}, b::Array{T}, c::Array{T})
    r = log(exp((a - b).^2) + exp(a + b)) - c .* log(c)
end

function devec_eval{T<:Real}(::complex_ewise, a::Array{T}, b::Array{T}, c::Array{T})
    @devec r = log(exp((a - b).^2) + exp(a + b)) - c .* log(c)
end

function hand_loop{T<:Real}(::complex_ewise, a::Array{T}, b::Array{T}, c::Array{T})
    r = similar(a)
    for i = 1 : length(a)
        r[i] = log(exp((a[i] - b[i]).^2) + exp(a[i] + b[i])) - c[i] .* log(c[i])
    end
end

# shift dot

type shift_dot end
formula(::shift_dot) = "sum( (a - mean(a)) .* (b - mean(b)) )"

function vec_eval{T<:Real}(::shift_dot, a::Array{T}, b::Array{T}, c::Array{T})
    r = sum( (a - mean(a)) .* (b - mean(b)) )
end

function devec_eval{T<:Real}(::shift_dot, a::Array{T}, b::Array{T}, c::Array{T})
    @devec r = sum( (a - mean(a)) .* (b - mean(b)) )
end

function hand_loop{T<:Real}(::shift_dot, a::Array{T}, b::Array{T}, c::Array{T})
    n = length(a)

    # calculate mean(a)
    sa = 0.
    for i = 1 : n
        sa += a[i]
    end
    ma = sa / n

    # calculate mean(b)
    sb = 0.
    for i = 1 : n
        sb += b[i]
    end
    mb = sb / n

    # calculate shift dot
    r = 0.
    for i = 1 : n
        r += (a[i] - ma) * (b[i] - mb)
    end
end

# column wise sum

type colwise_sum end
formula(::colwise_sum) = "sum(a, 1)"

function vec_eval{T<:Real}(::colwise_sum, a::Array{T}, b::Array{T}, c::Array{T})
    r = sum(a, 1)
end

function devec_eval{T<:Real}(::colwise_sum, a::Array{T}, b::Array{T}, c::Array{T})
    @devec r = sum(a, 1)
end

function hand_loop{T<:Real}(::colwise_sum, a::Array{T}, b::Array{T}, c::Array{T})
    m, n = size(a)
    r = zeros(1, n)
    for j = 1 : n
        s = 0.
        for i = 1 : m
            s += a[i, j]
        end
        r[j] = s
    end
end

# row-wise sum

type rowwise_sum end
formula(::rowwise_sum) = "sum(a, 2)"

function vec_eval{T<:Real}(::rowwise_sum, a::Array{T}, b::Array{T}, c::Array{T})
    r = sum(a, 2)
end

function devec_eval{T<:Real}(::rowwise_sum, a::Array{T}, b::Array{T}, c::Array{T})
    @devec r = sum(a, 2)
end

function hand_loop{T<:Real}(::rowwise_sum, a::Array{T}, b::Array{T}, c::Array{T})
    m, n = size(a)
    r = zeros(1, n)
    for j = 1 : n
        for i = 1 : m
            r[i] += a[i,j]
        end
    end
end

# col-wise euclidean

type colwise_eucdist end
formula(::colwise_eucdist) = "sqrt(sum((a - b).^2, 1))"

function vec_eval{T<:Real}(::colwise_eucdist, a::Array{T}, b::Array{T}, c::Array{T})
    r = sqrt(sum(sqr(a - b), 1))
end

function devec_eval{T<:Real}(::colwise_eucdist, a::Array{T}, b::Array{T}, c::Array{T})
    @devec r = sqrt(sum(sqr(a - b), 1))
end

function hand_loop{T<:Real}(::colwise_eucdist, a::Array{T}, b::Array{T}, c::Array{T})
    m, n = size(a)
    r = zeros(1, n)
    for j = 1 : n
        s = 0.
        for i = 1 : m
            v = a[i, j] - b[i, j]
            s += v * v
        end
        r[j] = sqrt(s)
    end
end

#  Benchmarks

@my_bench simple_ewise
@my_bench complex_ewise
@my_bench shift_dot
@my_bench colwise_sum
@my_bench rowwise_sum
@my_bench colwise_eucdist


text:
using NoteSequences: temporalstretch!, transpose

@testset "Temporal stretch" begin
    @testset "Speed up" begin
        ns = NoteSequence()
        addnotes!(ns, 0, [(60, 100, 0, 220), (64, 100, 220, 440), (67, 100, 440, 660)])
        temporalstretch!(ns, 3)
        @test notes(ns, 0) == [(60, 100, 0, 660), (64, 100, 660, 1320), (67, 100, 1320, 1980)]
        @test ns.total_time == 1980
    end

    @testset "Slow down" begin
        ns = NoteSequence()
        addnotes!(ns, 0, [(60, 100, 0, 220), (64, 100, 220, 440), (67, 100, 440, 660)])
        temporalstretch!(ns, 0.25)
        @test notes(ns, 0) == [(60, 100, 0, 55), (64, 100, 55, 110), (67, 100, 110, 165)]
        @test ns.total_time == 165
    end
end

@testset "Transpose" begin
    @testset "Higher and lower transpositions" begin
        ns = NoteSequence()
        addnotes!(ns, 1, [(116, 4, 0, 1), (42, 75, 103, 635), (67, 100, 216, 773)])

        higher_ns, num_deleted = transpose(ns, 1, 0, 127)
        @test notes(higher_ns, 1) == [(117, 4, 0, 1), (43, 75, 103, 635), (68, 100, 216, 773)]
        @test num_deleted == 0

        lower_ns, num_deleted = transpose(ns, -3, 0, 127)
        @test notes(lower_ns, 1) == [(113, 4, 0, 1), (39, 75, 103, 635), (64, 100, 216, 773)]
        @test num_deleted == 0
    end

    @testset "Out of range transpositions" begin
        ns = NoteSequence()
        addnotes!(ns, 1, [(116, 4, 0, 1), (42, 75, 103, 635), (67, 100, 216, 773)])

        removed_ns, num_deleted = transpose(ns, -12, 20, 60)
        @test notes(removed_ns, 1) == [(30, 75, 103, 635), (55, 100, 216, 773)]
        @test num_deleted == 1

        removed_ns, num_deleted = transpose(ns, 4, 100, 120)
        @test notes(removed_ns, 1) == [(120, 4, 0, 1)]
        @test num_deleted == 2
    end
end


text:
"""
Plot the performance statistics of greedy-multi-path vs the number of end-user pairs
"""

using QuNet
using LightGraphs, SimpleWeightedGraphs, GraphPlot, MetaGraphs, Plots, Colors, Statistics
using LaTeXStrings
using JLD
using Parameters

datafile = "data/userpairs"

# Params
max_pairs = 50::Int64
num_trials = 5000::Int64

generate_new_data = false
if generate_new_data == true

    # The average routing costs between end-users sampled over num_trials for different numbers of end-users
    perf_data = []
    # The associated errors of the costs sampled over num_trials
    perf_err = []
    # Average numbers of paths used, sampled over num_trials for different numbers of end-users
    # e.g. [3,4,5]: 3 end-users found no path on average, 4 end-users found 1 path on average etc.
    path_data = []
    # Associated errors of path_data
    path_err = []
    # Numer of paths found by the last userpair
    last_count = []
    # Rate at which last user finds a path
    prob_last_nopath = []

    grid_size = 10

    for i in 1:max_pairs
        println("Collecting for pairsize: $i")

        net = GridNetwork(grid_size, grid_size)

        # Collect performance statistics
        p, p_e, pat, pat_e, lc, num_zeropath = net_performance(net, num_trials, i, max_paths=4)
        push!(perf_data, p)
        push!(perf_err, p_e)
        push!(path_data, pat)
        push!(path_err, pat_e)
        push!(last_count, lc)
        push!(prob_last_nopath, num_zeropath/num_trials)
    end

    # Collect data for conditional probability of purification: (N2+N3)/∑N_i
    cpp = []
    cpp_err = []
    for i in 1:max_pairs
        tableau = path_data[i]
        errors = path_err[i]
        off = 1

        # Numerator and denominator
        num = tableau[2 + off] + tableau[3 + off] + tableau[4 + off]
        denom = tableau[1 + off] + tableau[2 + off] + tableau[3 + off] + tableau[4 + off]

        # Percentage errors (Ignore factor 100. Not needed)
        num_perr = (errors[2 + off] + errors[3 + off] + errors[4 + off])/num
        denom_perr = (errors[1 + off] + errors[2 + off] + errors[3 + off] + errors[4 + off])/denom

        data = num/denom
        data_err = (num_perr + denom_perr) * data
        push!(cpp, data)
        push!(cpp_err, data_err)
        i += 1
    end

    # Collect data for average number of paths used
    avepath = []
    for i in 1:max_pairs
        tableau = convert(Vector{Float64}, path_data[i])
        data = QuNet.ave_paths_used(tableau)

        # TODO: Include errors
        push!(avepath, data)
        i += 1
    end

    # Save data
    d = Dict{Symbol, Any}()
    @pack! d = max_pairs, num_trials, grid_size, perf_data, perf_err, path_data, last_count, prob_last_nopath, path_err, cpp, cpp_err, avepath
    save("$datafile.jld", "data", d)
else
    # Load data
    if !isfile("$datafile.jld")
        error("$datafile.jld not found in working directory")
    end
    d = load("$datafile.jld")["data"]
    @unpack max_pairs, num_trials, grid_size, perf_data, perf_err, path_data, last_count, prob_last_nopath, path_err, cpp, cpp_err, avepath = d
end

# Get values for x axis
x = collect(1:max_pairs)

# Extract data from performance
loss = collect(map(x->x["loss"], perf_data))
z = collect(map(x->x["Z"], perf_data))
loss_err = collect(map(x->x["loss"], perf_err))
z_err = collect(map(x->x["Z"], perf_err))

# Extract data from path: PX is the rate of using X paths
P0 = [path_data[i][1]/i for i in 1:max_pairs]
P1 = [path_data[i][2]/i for i in 1:max_pairs]
P2 = [path_data[i][3]/i for i in 1:max_pairs]
P3 = [path_data[i][4]/i for i in 1:max_pairs]
P4 = [path_data[i][5]/i for i in 1:max_pairs]

P0e = [path_err[i][1]/i for i in 1:max_pairs]
P1e = [path_err[i][2]/i for i in 1:max_pairs]
P2e = [path_err[i][3]/i for i in 1:max_pairs]
P3e = [path_err[i][4]/i for i in 1:max_pairs]
P4e = [path_err[i][5]/i for i in 1:max_pairs]

##### Plot cost_userpair #####
# color_palette = palette(:lightrainbow, 4)
plot(x, loss, ylims=(0,0.12), seriestype = :scatter, yerror = loss_err, label=L"$\eta$",
legend=:bottomright, guidefontsize=14, tickfontsize=12, legendfontsize=10, fontfamily="computer modern",
color=:Dodgerblue, markersize=5, markershape=:utriangle)
plot!(x, x./9 .* .093 .*(1 .-P0), label="average number that get path")
# xaxis!(L"$\textrm{Number of End User Pairs}$")
yaxis!("Costs")

# Collect and plot data for efficiency and fidelity per user pair
loss_per_user = []
for (n, e) in enumerate(loss)
    push!(loss_per_user, (1-P0[n])*e)
end
z_per_user = []
for (n, f) in enumerate(z)
    push!(z_per_user, (1-P0[n])*f)
end
plot!(x, loss_per_user, seriestype = :scatter, label=L"$\eta \textrm{ per user}$",
markersize=5, markershape=:utriangle, color=:LightSkyBlue)

plot!(x, z, seriestype = :scatter, yerror = z_err, label=L"$F$", markersize=5,
color = :Crimson)

# Try combining two plots:
combo_per_user = loss_per_user .* z
plot!(x, combo_per_user, seriestype = :scatter, label=L"$F \times \eta \textrm{ per user}$",
markersize=4, markershape=:square, color=:mediumpurple)

savefig("plots/cost_userpair.png")
savefig("plots/cost_userpair.pdf")

##### Plot path_userpair #####
plot(x, P0, ylims=(0,1), seriestype=:scatter, yerr = P0e, label=L"$P_0$", legend= :right,
guidefontsize=14, tickfontsize=12, legendfontsize=10, fontfamily="computer modern",
color_palette = palette(:plasma, 5), markersize=5)
plot!(x, P1, seriestype=:scatter, yerr = P1e, label=L"$P_1$", markersize=5)
plot!(x, P2, seriestype=:scatter, yerr = P2e, label=L"$P_2$", markersize=5)
plot!(x, P3, seriestype=:scatter, yerr = P3e, label=L"$P_3$", markersize=5)
plot!(x, P4, seriestype=:scatter, yerr = P4e, label=L"$P_4$", markersize=5)
plot!(x, cpp, seriestype=:scatter, yerr = cpp_err, label=L"$P_{P}$", color=:MediumSpringGreen,
markersize=5)
yaxis!("Path probability")
xaxis!(L"$\textrm{Number of End User Pairs}$")
# xaxis!(L"$\textrm{Number of End User Pairs}$")
savefig("plots/path_userpair.png")
savefig("plots/path_userpair.pdf")

plot(x, avepath, linewidth=2, legend = false,
xguidefontsize=14, tickfontsize=12, legendfontsize=8, fontfamily="computer modern")
xaxis!(L"$\textrm{Number of End User Pairs}$")
yaxis!(L"$\textrm{Average Number of Paths Used}$")
savefig("plots/avepath_userpair.png")
savefig("plots/avepath_userpair.pdf")

##### Plot nopath_rate #####
# yaxis=:log not working because zero -> -Inf log scale
plot(x[1:30], last_count[1:30], seriestype = :scatter, legend = false, yaxis=(:log10, [10^(-5), :auto]))
# println(last_count[1:30])
# plot(x, last_count, seriestype = :scatter, legend = false)
# xaxis!(L"$\textrm{Number of End User Pairs}$")
yaxis!(L"$\textrm{Number of Paths found by last pair}$")
savefig("plots/lastpair.png")
savefig("plots/lastpair.pdf")

function get_firstzero_loc(arr)
    for (idx, val) in enumerate(prob_last_nopath)
        if val == 1
            return idx
        end
    end
    return idx
end

first_zero = get_firstzero_loc(prob_last_nopath)

plot(x[1:first_zero], 1 .- prob_last_nopath[1:first_zero], seriestype = :scatter, legend = false,
guidefontsize=14, tickfontsize=12, legendfontsize=8, fontfamily="computer modern", markersize=5)
plot!(x[first_zero+1:max_pairs], 1 .- prob_last_nopath[first_zero+1:max_pairs], seriestype =:scatter,
legend = false, color=:white, markersize=5)
# xaxis!(L"$\textrm{Number of End User Pairs}$")
yaxis!(L"$\textrm{Rate that last pair finds path}$")

plot!(x[1:20], 1 .- prob_last_nopath[1:20], seriestype = :scatter, inset = (1, bbox(0.05, 0.05, 0.33, 0.33, :top, :right)),
subplot=2, legend=false, yaxis=(:log10, [10^(-5), :auto]))
savefig("plots/nopath_rate.png")
savefig("plots/nopath_rate.pdf")
# println(last_count)

# Simple graph saturation explanation?
plot(x, x.*(1 .-P0), label="average number of users that get a path")

# plot!(x[1:20], 1 .- prob_last_nopath[1:20], seriestype = :scatter, inset = (1, bbox(0.10, 0.10, 0.5, 0.5, :bottom, :left)))

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
"""
# SMM.jl
This is the documentation for SMM.jl

"""

module SMM

  # Dependencies
  #-------------
  using BlackBoxOptim
  using Optim
  using JLD2
  using Plots
  using CSV
  using ProgressMeter
  using Calculus
  using CompEcon
  using Distributions
  using DataFrames

  @static if VERSION < v"0.7.0"
      using DataStructures
  else
      using OrderedCollections
  end

  # Exports from BlackBoxOptim
  #---------------------------
  export best_candidate

  # Includes
  #---------
  # Types
  #------
  include("types.jl");

  # API
  #----
  include("api.jl")

  # General functions, useful at several places
  #---------------------------------------------
  include("generic.jl")

  # Functions to load and save
  #--------------------------
  include("save_load.jl")

  # Functions to minimize the objective function
  #---------------------------------------------
  include("optimize.jl")

  # Functions to do inference
  #---------------------------------------------
  include("econometrics.jl")

  # Functions to do plots
  #---------------------------------------------
  include("plots.jl")

  # Exports
  #--------
  # Functions and types in types.jl
  #----------------------------------
  export SMMOptions, SMMProblem
  export default_function, rosenbrock2d
  export is_global_optimizer, is_local_optimizer
  export convert_to_optim_algo, convert_to_fminbox
  export is_bb_optimizer, is_optim_optimizer

  # Functions and types in api.jl
  #-------------------------------

  # Functions and types in generic.jl
  #----------------------------------
  export set_simulate_empirical_moments!, construct_objective_function!
  export set_priors!, set_empirical_moments!, set_Sigma0!
  export set_bbSetup!, generate_bbSearchRange
  export create_lower_bound, create_upper_bound
  export set_global_optimizer!
  export create_grid, create_grid_stochastic, generate_std, latin_hypercube_sampling
  export get_now

  # Functions and types in save_load.jl
  #------------------------------------
  export read_priors, read_empirical_moments
  export saveSMMOptim, loadSMMOptim

  # Functions and types in optimize.jl
  #-----------------------------------
  export smmoptimize!, smm_minimizer
  export smm_refine_globalmin!, smm_local_minimizer
  export smm_local_minimum
  export smm_localmin, local_to_global!

  # Functions in econometrics.jl
  #-----------------------------
  export calculate_D, calculate_Avar!, calculate_se, calculate_t, calculate_pvalue, calculate_CI
  export summary_table

  # Functions in plots.jl
  #----------------------
  export smm_slices

end # module


text:
"""
    applyrule(rule::Rule, data, state, index)

Updates cell values based on their current state and the state of other cells
as defined in the Rule.

### Arguments:
- `rule` : [`Rule`](@ref)
- `data` : [`FrameData`](@ref)
- `state`: the value of the current cell
- `index`: a (row, column) tuple of Int for the current cell coordinates - `t`: the current time step

Returns a value to be written to the current cell.
"""
function applyrule end

"""
    applyrule!(rule::PartialRule, data, state, index)

A rule that manually writes to the dest array, used in rules inheriting
from [`PartialRule`](@ref).

### Arguments:
see [`applyrule`](@ref)
"""
function applyrule! end

"""
    applyinteraction(interacttion::PartialRule, data, state, index)

Applay an interation that returns a tuple of values.
### Arguments:
see [`applyrule`](@ref)
"""
function applyinteraction end

"""
    applyinteraction!(interacttion::PartialRule, data, state, index)

Applay an interation that manually writes to the passed in dest arrays.
### Arguments:
see [`applyrule`](@ref)
"""
function applyinteraction! end

"""
    precalcrule!(rule, data)

Run any precalculations needed to run a rule for a particular frame.

It may be better to do this in a functional way with an external precalc object
passed into a rule via the `data` object, but it's done statefully for now for simplicity.
"""
function precalcrule! end

"""
    neighbors(hood::Neighborhood, state, indices, t, source, args...)

Checks all cells in neighborhood and combines them according
to the particular neighborhood type.
"""
function neighbors end

"""
    mapreduceneighbors(f, data, neighborhood, rule, state, index)

Run `f` over all cells in the neighborhood and sums its return values.
`f` is a function or functor with the form:
`f(data, neighborhood, rule, state, hood_index, dest_index)`.
"""
function mapreduceneighbors end

"""
Set value of a cell in the neighborhood.
Usually called in `mapreduceneighbors`.
"""
function setneighbor! end

"""
Return the radius of a rule or ruleset if it has one, otherwise zero.
"""
function radius end


text:
using LinearAlgebra: I, qr
using Printf: @sprintf

push!(LOAD_PATH, @__DIR__)
using MultiprecisionFloats
using RKTK2
using DZMisc

const AccurateReal = Float64x3
const BASE_EVALUATOR = RKOCEvaluator{AccurateReal}(10, 16)
const ERROR_EVALUATOR = RKOCEvaluator{AccurateReal}(11, 16)

const NUM_VARS = BASE_EVALUATOR.num_vars
const NUM_CONSTRS = BASE_EVALUATOR.num_constrs
const NUM_ERR_COEFFS = ERROR_EVALUATOR.num_constrs

const CONSTR_IDXS = [
    1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 18, 19, 20, 21, 22, 23,
    24, 25, 26, 27, 28, 31, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 51, 52,
    58, 59, 62, 75, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99,
    100, 101, 106, 107, 110, 115, 116, 123, 134, 135, 201, 202, 203, 205, 206,
    207, 208, 210, 211, 214, 215, 221, 222, 225, 226, 227, 230, 231, 238, 249,
    250, 269, 270, 316, 404, 406, 487, 488, 489, 490, 491, 492, 493, 496, 497,
    498, 499, 500, 501, 502, 507, 508, 511, 516, 517, 524, 535, 536, 539, 544,
    555, 556, 602, 603, 650, 690, 984]

################################################################################

function log(args...)
    println(args...)
    flush(stdout)
end

function constraint_residual(x)
    result = Vector{AccurateReal}(undef, NUM_CONSTRS)
    evaluate_residual!(result, x, BASE_EVALUATOR)
    result[CONSTR_IDXS]
end

function constraint_jacobian(x)
    result = Matrix{AccurateReal}(undef, NUM_CONSTRS, NUM_VARS)
    evaluate_jacobian!(result, x, BASE_EVALUATOR)
    result[CONSTR_IDXS, :]
end

function constraint_objective(x)
    norm2(constraint_residual(x))
end

function error_coefficients(x)
    result = Vector{AccurateReal}(undef, NUM_ERR_COEFFS)
    evaluate_error_coefficients!(result, x, ERROR_EVALUATOR)
    result[NUM_CONSTRS+1:end]
end

function error_jacobian(x)
    result = Matrix{AccurateReal}(undef, NUM_ERR_COEFFS, NUM_VARS)
    evaluate_error_jacobian!(result, x, ERROR_EVALUATOR)
    result[NUM_CONSTRS+1:end, :]
end

function error_objective(x)
    norm2(error_coefficients(x))
end

function error_gradient(x)
    dbl.(transpose(error_jacobian(x)) * error_coefficients(x))
end

function constrained_error_gradient(x)
    constraint_jac = copy(transpose(constraint_jacobian(x)))
    orthonormalize_columns!(constraint_jac)
    grad = error_gradient(x)
    grad - constraint_jac * (transpose(constraint_jac) * grad)
end

function constrain_step(x, step)
    constraint_jac = copy(transpose(constraint_jacobian(x)))
    orthonormalize_columns!(constraint_jac)
    step - constraint_jac * (transpose(constraint_jac) * step)
end

################################################################################

function constrain(x)
    x_old, obj_old = x, constraint_objective(x)
    while true
        direction = qr(constraint_jacobian(x_old)) \ constraint_residual(x_old)
        x_new = x_old - direction
        obj_new = constraint_objective(x_new)
        if obj_new < obj_old
            x_old, obj_old = x_new, obj_new
        else
            break
        end
    end
    x_old, obj_old
end

function constrained_step_value(step_size, x, step_direction, step_norm)
    x_new, obj_new = constrain(x - (step_size / step_norm) * step_direction)
    if obj_new < AccurateReal(1e-75)
        error_objective(x_new)
    else
        AccurateReal(NaN)
    end
end

struct ConstrainedBFGSOptimizer
    objective_value::Ref{AccurateReal}
    last_step_size::Ref{AccurateReal}
end

const INPUT_FILE_NAME = readdir()[end]
log("Opening ", INPUT_FILE_NAME)
const INPUT_POINT = AccurateReal.(BigFloat.(split(read(INPUT_FILE_NAME, String))))
@assert(length(INPUT_POINT) == NUM_VARS)
log("Successfully read input file.")

const opt = ConstrainedBFGSOptimizer(
    error_objective(INPUT_POINT),
    AccurateReal(0.00001))

function main()

    x = copy(INPUT_POINT)
    inv_hess = Matrix{AccurateReal}(I, NUM_VARS, NUM_VARS)

    cons_grad = constrained_error_gradient(x)
    cons_grad_norm = norm(cons_grad)
    while true

        # log("Performing gradient descent step...")
        grad_step_size, obj_grad = quadratic_line_search(constrained_step_value,
            opt.objective_value[], opt.last_step_size[], x, cons_grad, cons_grad_norm)

        # log("Performing BFGS step...")
        bfgs_step = constrain_step(x, inv_hess * cons_grad)
        bfgs_step_norm = norm(bfgs_step)
        bfgs_step_size, obj_bfgs = quadratic_line_search(constrained_step_value,
            opt.objective_value[], opt.last_step_size[], x, bfgs_step, bfgs_step_norm)

        if obj_bfgs < opt.objective_value[] && obj_bfgs <= obj_grad
            x, _ = constrain(x - (bfgs_step_size / bfgs_step_norm) * bfgs_step)
            cons_grad_new = constrained_error_gradient(x)
            cons_grad_norm_new = norm(cons_grad_new)
            update_inverse_hessian!(inv_hess,
                -bfgs_step_size / bfgs_step_norm,
                bfgs_step,
                cons_grad_new - cons_grad,
                Vector{AccurateReal}(undef, NUM_VARS))
                opt.last_step_size[] = bfgs_step_size
                opt.objective_value[] = obj_bfgs
            cons_grad = cons_grad_new
            cons_grad_norm = cons_grad_norm_new
            log(Float64(opt.objective_value[]), " | ",
                Float64(cons_grad_norm), " | ",
                Float64(opt.last_step_size[]), " | BFGS")
        elseif obj_grad < opt.objective_value[]
            x, _ = constrain(x - (grad_step_size / cons_grad_norm) * cons_grad)
            inv_hess = Matrix{AccurateReal}(I, NUM_VARS, NUM_VARS)
            opt.last_step_size[] = grad_step_size
            opt.objective_value[] = obj_grad
            cons_grad = constrained_error_gradient(x)
            cons_grad_norm = norm(cons_grad)
            log(Float64(opt.objective_value[]), " | ",
                Float64(cons_grad_norm), " | ",
                Float64(opt.last_step_size[]), " | GRAD")
        else
            log("Neither BFGS step (", Float64(obj_bfgs),
                ") or gradient step (", Float64(obj_grad),
                ") are better than previous point (", Float64(opt.objective_value[]), ")")
            break
        end
    end
end

main()


text:
using DataStructures
import GlobalSensitivityAnalysis

mutable struct DeltaData <: AbstractSimulationData

    function DeltaData()
        return new()
    end
end

const DeltaSimulationDef = SimulationDef{DeltaData}
const DeltaSimulationInstance = SimulationInstance{DeltaData}

function create_GSA_payload(sim_inst::DeltaSimulationInstance)
    rv_info = OrderedDict{Symbol, Any}([name => _get_dist(rv) for (name, rv) in sim_inst.sim_def.rvdict])

    # create payload
    return GlobalSensitivityAnalysis.DeltaData(params = rv_info, N = sim_inst.trials)
end

function sample!(sim_inst::DeltaSimulationInstance, samplesize::Int)
    rvdict = sim_inst.sim_def.rvdict
    sim_inst.trials = samplesize

    # get the samples to plug in to trials
    payload = create_GSA_payload(sim_inst)
    samples = GlobalSensitivityAnalysis.sample(payload)

    for (i, rv) in enumerate(values(rvdict))
        # use underlying distribution, if known
        orig_dist = _get_dist(rv)

        name = rv.name
        values = samples[:, i]
        rvdict[name] = RandomVariable(name, SampleStore(values, orig_dist))
    end
end

function analyze(sim_inst::DeltaSimulationInstance, model_input::AbstractArray{<:Number, N1}, model_output::AbstractArray{<:Number, N2}; num_resamples::Int = 1_000, conf_level::Number = 0.95, N_override::Union{Nothing, Int}=nothing, progress_meter::Bool = true) where {N1, N2}

    if sim_inst.trials == 0
        error("Cannot analyze simulation with 0 trials.")
    end

    payload = create_GSA_payload(sim_inst)

    return GlobalSensitivityAnalysis.analyze(payload, model_input, model_output; num_resamples = num_resamples, conf_level = conf_level, N_override = N_override, progress_meter = progress_meter)
end


text:
module RTLSDR

using DSP: welch_pgram
import Base.open, Base.close

# I'm not sure I want this in here
using PyPlot

export RtlSdr, open, close, read_samples, get_strength, get_strength2
export get_rate, set_rate, get_freq, set_freq

include("c_interface.jl")

mutable struct RtlSdr
	valid_ptr::Bool
	dongle_ptr::Ptr{rtlsdr_dev}

	function RtlSdr()
		dp = rtlsdr_open()

		r = new(true, dp)

		# set sample rate and center freq
		set_rate(r, 2.0e6)
		set_freq(r, 88.5e6)

		return r
	end
end

function open(r::RtlSdr, device_index::Int)
end

#close(r::RtlSdr) = rtlsdr_close(r.dongle_ptr)
function close(r::RtlSdr)
	rtlsdr_close(r.dongle_ptr)
	r.valid_ptr = false
end

# Sample rate in MHz
function get_rate(r::RtlSdr)
	@assert r.valid_ptr
	rate = rtlsdr_get_sample_rate(r.dongle_ptr)
	return Int(rate)
end
function set_rate(r::RtlSdr, sample_rate)
	@assert r.valid_ptr
	rtlsdr_set_sample_rate(r.dongle_ptr, sample_rate)
end

# center frequency in MHz
function set_freq(r::RtlSdr, freq)
	@assert r.valid_ptr
	rtlsdr_set_center_freq(r.dongle_ptr, freq)
end
function get_freq(r::RtlSdr)
	@assert r.valid_ptr
	freq = rtlsdr_get_center_freq(r.dongle_ptr)
	return Int(freq)
end

"""
`read_samples(r::RtlSdr, num_samples)`

Returns a vector of length `num_samples` with complex numbers.
"""
function read_samples(r::RtlSdr, num_samples)
	@assert r.valid_ptr
	num_bytes = 2num_samples
	raw_data = read_bytes(r.dongle_ptr, num_bytes)
	return packed_bytes_to_iq(raw_data)
end
function packed_bytes_to_iq(bytes)
	num_bytes = length(bytes)
	num_iq = round(Int, num_bytes/2.0, RoundDown)
	iq_vals = zeros(Complex{Float64}, num_iq)
	den = 255.0/2.0

	for i = 1:num_iq
		iq_vals[i] = bytes[2i-1]/den - 1.0 + im*(bytes[2i]/den - 1.0)
	end

	return iq_vals
end

# maximum
function get_strength(r::RtlSdr, n=10; plot_max::Bool=false)
	max_sample = -Inf
	max_p = 0
	sample_rate = get_rate(r)
	for i = 1:n
		#samples = read_samples(r,256*10240)
		#samples = read_samples(r, 256 * 1024 * 8)
		#samples = read_samples(r,256*1024)
		#samples = read_samples(r,256*940)
		#samples = read_samples(r,256*920)
		#samples = read_samples(r,256*900)
		#samples = read_samples(r,256*800)
		samples = read_samples(r, 256*500)
		#samples = read_samples(r, 256*50)
		p = welch_pgram(samples, fs=sample_rate)
		temp_max = maximum(p.power)
		if temp_max > max_sample
			max_sample = temp_max
			max_p = deepcopy(p)
		end
	end
	if plot_max
		center_freq = get_freq(r)
		plot(max_p.freq + center_freq, 10log10(max_p.power))
	end
	return 10log10(max_sample)
end

function get_strength2(r::RtlSdr, n=10; plot_max::Bool=false)
	max_sample = -Inf
	max_p = 0
	sample_rate = get_rate(r)
	max_freqs = 0
	max_pows = 0
	for i = 1:n
		samples = read_samples(r, 256*500)
		p = welch_pgram(samples, fs=sample_rate)
		n_plot = round(Int, length(p.power) / 2, RoundDown)-3
		pows = p.power[2:n_plot]
		freqs = p.freq[2:n_plot]

		# taking maximum...
		temp_max = maximum(pows)

		# taking spectrum max...
		temp_max = sum(pows)

		if temp_max > max_sample
			max_sample = temp_max
			max_p = deepcopy(p)
			max_pows = deepcopy(pows)
			max_freqs = deepcopy(freqs)
		end
	end
	if plot_max
		center_freq = get_freq(r)
		#plot(max_p.freq + center_freq, 10log10(max_p.power))
		plot(max_freqs + center_freq, 10log10(max_pows))
	end
	return 10log10(max_sample*146.)
end

end # module


text:
include("./aug_lagrangian.jl")

function aug_test()
    f(x, ::Val{false}) = x[1] ^ 2 + x[2] ^ 2 - 2
    f(x, ::Val{true}) = x[1] ^ 2 + x[2] ^ 2 - 2, 2 .* [x[1], x[2]]
    c(x, ::Val{false}) = [x[1] - 2]
    c(x, ::Val{true}) = [x[1] - 2], [1 0;]
    x, ok = augmented_lagrangian(f, c, [10.0, 10.0], [5.0], 10.0; rtol=1e-10, maxiter=20, μpenalty=10)
    display(x)
    display(ok)
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
# test_helpers.jl - Funciones de ayuda para probar los tipos de este paquete

"""
    getrandomweights(T=Float32, G=218)
Función para generar pesos aleatorios
"""
function getrandomweights(T=Float32, G=218)
    w = rand(T, G)
    w = 100 * w / sum(w)
    w
end

"""
    getbasedates(vmat, startdate=Date(2000, 12))
Función para generar fechas a partir de matriz de variaciones intermensuales
"""
function getbasedates(vmat, startdate=Date(2000, 12))
    T = size(vmat, 1)
    startdate:Month(1):(startdate + Month(T-1))
end

"""
    getzerobase(;
        [T_type = Float32,
        G = 218,
        T_periods = 120,
        startdate = Date(2001,1),
        baseindex = 100*one(T_type)])
Función para obtener base de tipo `VarCPIBase` con variaciones intermensuales
iguales a cero.
"""
function getzerobase(;
    T_type = Float32,
    G = 218,
    T_periods = 120,
    startdate = Date(2001,1),
    baseindex = 100*one(T_type))

    getzerobase(T_type, G, T_periods, startdate, baseindex)
end

function getzerobase(T_type=Float32, G=218, T_periods=120, startdate=Date(2001,1), baseindex=100*one(T_type))
    vmat = zeros(T_type, T_periods, G)
    w = getrandomweights(T_type, G)
    dates = getbasedates(vmat, startdate)
    VarCPIBase(vmat, w, dates, baseindex)
end

"""
    getzerocountryst(T_type = Float32)
Obtiene un `UniformCountryStructure` cuyas variaciones intermensuales son todas
iguales a cero en la configuración de períodos del IPC de Guatemala.
"""
function getzerocountryst(T_type = Float32)
    gt00 = getzerobase(T_type, 218, 120, Date(2001, 1))
    gt10 = getzerobase(T_type, 279, 120, Date(2011, 1))
    UniformCountryStructure(gt00, gt10)
end

"""
    getrandomcountryst(T_type = Float32)
Obtiene un `UniformCountryStructure` cuyas variaciones intermensuales son todas
aleatorias en la configuración de períodos del IPC de Guatemala.
"""
function getrandomcountryst(T_type = Float32)
    gt00 = getrandombase(T_type, 218, 120, Date(2001, 1))
    gt10 = getrandombase(T_type, 279, 120, Date(2011, 1))
    UniformCountryStructure(gt00, gt10)
end

"""
    getrandombase(;
        [T_type = Float32,
        G = 218,
        T_periods = 120,
        startdate = Date(2001,1),
        baseindex = 100*one(T_type)])
Función para obtener una base de tipo `VarCPIBase` con variaciones
intermensuales aleatorias.
"""
function getrandombase(;
    T_type = Float32,
    G = 218,
    T_periods = 120,
    startdate = Date(2001,1),
    baseindex = 100*one(T_type))

    getrandombase(T_type, G, T_periods, startdate, baseindex)
end

function getrandombase(T_type=Float32, G=218, T_periods=120, startdate=Date(2001,1), baseindex=100*one(T_type))
    vmat = rand(T_type, T_periods, G) / 2
    w = getrandomweights(T_type, G)
    dates = getbasedates(vmat, startdate)
    VarCPIBase(vmat, w, dates, baseindex)
end


text:
using NCDatasets
using Test

filename = tempname()
# The mode "c" stands for creating a new file (clobber)
ds = NCDatasets.NCDataset(filename,"c")

# define the dimension "lon" and "lat" with the size 10 and 11 resp.
ds.dim["lon"] = 10
ds.dim["lat"] = 11

fv = fillvalue(Float32)

v = NCDatasets.defVar(ds,"var_with_missing_data",Float32,("lon","lat"), fillvalue = fv)

data = [Float32(i+j) for i = 1:10, j = 1:11]
# mask the frist element
datam = Array{Union{Float32,Missing}}(data)
datam[1] = missing

v[:,:] = datam
@test ismissing(v[1,1])
@test isequal(v[:,:],datam)

# load without transformation
@test v.var[1,1] == fv

# write/read without transformation
v.var[:,:] = data
@test v.var[:,:] ≈ data

NCDatasets.close(ds)

sz = (4,5)
filename = tempname()
#filename = "/tmp/test-6.nc"
# The mode "c" stands for creating a new file (clobber)

NCDatasets.NCDataset(filename,"c") do ds

    # define the dimension "lon" and "lat"
    ds.dim["lon"] = sz[1]
    ds.dim["lat"] = sz[2]

    # variables
    for T in [UInt8,Int8,UInt16,Int16,UInt32,Int32,UInt64,Int64,Float32,Float64]
    #for T in [Float32]
        v = NCDatasets.defVar(ds,"var-$T",T,("lon","lat"); fillvalue = 124)
        v[:,:] = fill(T(123),size(v))
        @test all(v[:,:][:] .== 123)

        @test NCDatasets.fillvalue(v) == 124
        @test NCDatasets.fillmode(v) == (false,124)

    end
end

# NaN as fillvalue

filename = tempname()
# The mode "c" stands for creating a new file (clobber)
ds = NCDatasets.NCDataset(filename,"c")

# define the dimension "lon" and "lat" with the size 10 and 11 resp.
ds.dim["lon"] = 10
ds.dim["lat"] = 11

fv = NaN32
v = NCDatasets.defVar(ds,"var_with_missing_data",Float32,("lon","lat"), fillvalue = fv)

data = [Float32(i+j) for i = 1:10, j = 1:11]
# mask the frist element
datam = Array{Union{Float32,Missing}}(data)
datam[1] = missing

v[:,:] = datam
@test ismissing(v[1,1])
@test isequal(v[:,:],datam)

# load without transformation
@test isnan(v.var[1,1])

NCDatasets.close(ds)

# all fill-values
filename = tempname()
ds = NCDatasets.NCDataset(filename,"c")
ds.dim["lon"] = 3

v = NCDatasets.defVar(ds,"var_with_all_missing_data",Float32,("lon",), fillvalue = fv)
data = [missing, missing, missing]

v[:] = data
@test all(ismissing.(v[:]))

v[1] = 1234.
@test !ismissing(v[1])

v[1] = missing
@test ismissing(v[1])

v[1] = 1234.
v[1:1] .= missing
@test ismissing(v[1])

NCDatasets.close(ds)

# test nomissing

data = [missing, Float64(1.), Float64(2.)]
@test_throws ErrorException NCDatasets.nomissing(data)

dataf = NCDatasets.nomissing(data,-9999.)
@test eltype(dataf) == Float64
@test dataf == [-9999., 1., 2.]

data = Union{Float64,Missing}[1., 2.]
dataf = NCDatasets.nomissing(data)
@test eltype(dataf) == Float64
@test dataf == [1., 2.]

@test nomissing(Union{Int64,Missing}[]) == []

@test nomissing([1,2,3,4],-9999) == [1,2,3,4]
@test nomissing([missing,2,3,4],-9999) == [-9999,2,3,4]

# issue 39
using NCDatasets

filename = tempname()

NCDataset(filename, "c") do ds
    defDim(ds, "lon", 2)
    defDim(ds, "lat", 2)
    v = defVar(ds, "Char variable", Char, ("lon","lat"), fillvalue = ' ')
    v[:,:] = ['a' 'b'; 'c' 'd']
end

NCDataset(filename, "r") do ds
    @test ds["Char variable"][:,:] == ['a' 'b'; 'c' 'd']
end

rm(filename)


text:
"""
**AnalyticVI**

Variational Inference solver for conjugate or conditionally conjugate likelihoods (non-gaussian are made conjugate via augmentation)
All data is used at each iteration (use AnalyticSVI for Stochastic updates)

```julia
AnalyticVI(;ϵ::T=1e-5)
```
**Keywords arguments**

    - `ϵ::T` : convergence criteria
"""
mutable struct AnalyticVI{T<:Real} <: Inference{T}
    ϵ::T #Convergence criteria
    nIter::Integer #Number of steps performed
    optimizer::LatentArray{Optimizer} #Learning rate for stochastic updates
    Stochastic::Bool #Use of mini-batches
    nSamples::Int64
    nSamplesUsed::Int64 #Size of mini-batches
    MBIndices::Vector{Int64} #Indices of the minibatch
    ρ::T #Stochastic Coefficient
    HyperParametersUpdated::Bool #To know if the inverse kernel matrix must updated
    ∇η₁::LatentArray{Vector{T}}
    ∇η₂::LatentArray{Matrix{T}} #Stored as a matrix since symmetric sums do not help for the moment WARNING
    x::SubArray{T,2,Matrix{T}}#,Tuple{Base.Slice{Base.OneTo{Int64}},Base.Slice{Base.OneTo{Int64}}},true}
    y::LatentArray{SubArray}
    function AnalyticVI{T}(ϵ::T,nIter::Integer,optimizer::AbstractVector{<:Optimizer},Stochastic::Bool,nSamples::Integer,nSamplesUsed::Integer,MBIndices::AbstractVector,ρ::T,flag::Bool) where T
        return new{T}(ϵ,nIter,optimizer,Stochastic,nSamples,nSamplesUsed,MBIndices,ρ,flag)
    end
    function AnalyticVI{T}(ϵ::T,nIter::Integer,optimizer::AbstractVector{<:Optimizer},Stochastic::Bool,nSamples::Integer,nSamplesUsed::Integer,MBIndices::AbstractVector,ρ::T,flag::Bool,∇η₁::AbstractVector{<:AbstractVector},
    ∇η₂::AbstractVector{<:AbstractMatrix}) where T
        return new{T}(ϵ,nIter,optimizer,Stochastic,nSamples,nSamplesUsed,MBIndices,ρ,flag,∇η₁,∇η₂)
    end
end

function AnalyticVI(;ϵ::T=1e-5) where {T<:Real}
    AnalyticVI{Float64}(ϵ,0,[VanillaGradDescent(η=1.0)],false,1,1,[1],1.0,true)
end

"""
**AnalyticSVI**
Stochastic Variational Inference solver for conjugate or conditionally conjugate likelihoods (non-gaussian are made conjugate via augmentation)

```julia
AnalyticSVI(nMinibatch::Integer;ϵ::T=1e-5,optimizer::Optimizer=InverseDecay())
```
    - `nMinibatch::Integer` : Number of samples per mini-batches

**Keywords arguments**

    - `ϵ::T` : convergence criteria
    - `optimizer::Optimizer` : Optimizer used for the variational updates. Should be an Optimizer object from the [GradDescent.jl](https://github.com/jacobcvt12/GradDescent.jl) package. Default is `InverseDecay()` (ρ=(τ+iter)^-κ)
"""
function AnalyticSVI(nMinibatch::Integer;ϵ::T=1e-5,optimizer::Optimizer=InverseDecay()) where {T<:Real}
    AnalyticVI{T}(ϵ,0,[optimizer],true,1,nMinibatch,1:nMinibatch,1.0,true)
end

function Base.show(io::IO,inference::AnalyticVI{T}) where T
    print(io,"Analytic$(inference.Stochastic ? " Stochastic" : "") Variational Inference")
end

"""Initialize the final version of the inference object"""
function init_inference(inference::AnalyticVI{T},nLatent::Integer,nFeatures::Integer,nSamples::Integer,nSamplesUsed::Integer) where {T<:Real}
    inference.nSamples = nSamples
    inference.nSamplesUsed = nSamplesUsed
    inference.MBIndices = 1:nSamplesUsed
    inference.ρ = nSamples/nSamplesUsed
    inference.optimizer = [copy(inference.optimizer[1]) for _ in 1:nLatent]
    inference.∇η₁ = [zeros(T,nFeatures) for _ in 1:nLatent];
    inference.∇η₂ = [Matrix(Diagonal(ones(T,nFeatures))) for _ in 1:nLatent]
    return inference
end

"""Generic method for variational updates using analytical formulas"""
function variational_updates!(model::AbstractGP{T,L,AnalyticVI{T}}) where {T,L}
    local_updates!(model)
    natural_gradient!(model)
    global_update!(model)
end

#Coordinate ascent updates on the natural parameters
function natural_gradient!(model::Union{VGP{T,L,AnalyticVI{T}},VStP{T,L,AnalyticVI{T}}}) where {T,L}
    model.η₁ .= ∇E_μ(model) .+ invK(model).*model.μ₀
    model.η₂ .= -Symmetric.(Diagonal{T}.(∇E_Σ(model)).+0.5.*invK(model))
end

#Computation of the natural gradient for the natural parameters
function natural_gradient!(model::SVGP{T,L,AnalyticVI{T}}) where {T,L}
    model.inference.∇η₁ .= ∇η₁.(∇E_μ(model),fill(model.inference.ρ,model.nLatent),model.κ,model.invKmm,model.μ₀,model.η₁)
    model.inference.∇η₂ .= ∇η₂.(∇E_Σ(model),fill(model.inference.ρ,model.nLatent),model.κ,model.invKmm,model.η₂)
end

function ∇η₁(∇μ::AbstractVector{T},ρ::Real,κ::AbstractMatrix{T},invKmm::Symmetric{T,Matrix{T}},μ₀::PriorMean,η₁::AbstractVector{T}) where {T <: Real}
    transpose(κ)*(ρ*∇μ) + invKmm*μ₀ - η₁
end

function ∇η₂(θ::AbstractVector{T},ρ::Real,κ::AbstractMatrix{<:Real},invKmm::Symmetric{T,Matrix{T}},η₂::Symmetric{T,Matrix{T}}) where {T<:Real}
    -(ρκdiagθκ(ρ,κ,θ)+0.5.*invKmm) - η₂
end

#Conversion from natural to standard distribution parameters
function global_update!(model::Union{VGP{T,L,AnalyticVI{T}},VStP{T,L,AnalyticVI{T}}}) where {T,L}
    model.Σ .= -0.5.*inv.(model.η₂)
    model.μ .= model.Σ.*model.η₁
end

#Update of the natural parameters and conversion from natural to standard distribution parameters
function global_update!(model::SVGP{T,L,AnalyticVI{T}}) where {T,L}
    if model.inference.Stochastic
        for k in 1:model.nLatent
            Δ = GradDescent.update(model.inference.optimizer[k],vcat(model.inference.∇η₁[k],model.inference.∇η₂[k][:]))
            model.η₁[k] .= model.η₁[k] + Δ[1:model.nFeatures]
            model.η₂[k] .= Symmetric(model.η₂[k] + reshape(Δ[(model.nFeatures+1):end],model.nFeatures,model.nFeatures))
            # model.η₁ .= model.η₁ .+ GradDescent.update.(model.inference.optimizer_η₁,model.inference.∇η₁)
            # model.η₂ .= Symmetric.(model.η₂ .+ GradDescent.update.(model.inference.optimizer_η₂,model.inference.∇η₂))
        end
    else
        model.η₁ .+= model.inference.∇η₁
        model.η₂ .= Symmetric.(model.inference.∇η₂ .+ model.η₂)
    end
    model.Σ .= -0.5.*inv.(model.η₂)
    model.μ .= model.Σ.*model.η₁
end