text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    get_bounds(model::MOI.ModelLike, ::Type{T}, x::MOI.VariableIndex)

Return a tuple `(lb, ub)` of type `Tuple{T, T}`, where `lb` and `ub` are lower
 and upper bounds, respectively, imposed on `x` in `model`.
"""
function get_bounds(
    model::MOI.ModelLike,
    ::Type{T},
    x::MOI.VariableIndex,
) where {T}
    xval = x.value
    c_lt = MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{T}}(xval)
    c_gt = MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{T}}(xval)
    c_int = MOI.ConstraintIndex{MOI.VariableIndex,MOI.Interval{T}}(xval)
    c_eq = MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{T}}(xval)
    c_sc = MOI.ConstraintIndex{MOI.VariableIndex,MOI.Semicontinuous{T}}(xval)
    c_si = MOI.ConstraintIndex{MOI.VariableIndex,MOI.Semiinteger{T}}(xval)
    if MOI.is_valid(model, c_int)
        # It is assumed that none of the other ConstraintIndexs are valid
        int::MOI.Interval{T} = MOI.get(model, MOI.ConstraintSet(), c_int)
        return int.lower, int.upper
    elseif MOI.is_valid(model, c_eq)
        # It is assumed that none of the other ConstraintIndexs are valid
        eq::MOI.EqualTo{T} = MOI.get(model, MOI.ConstraintSet(), c_eq)
        return eq.value, eq.value
    elseif MOI.is_valid(model, c_sc)
        # It is assumed that none of the other ConstraintIndexs are valid
        sc::MOI.Semicontinuous{T} = MOI.get(model, MOI.ConstraintSet(), c_sc)
        return min(zero(T), sc.lower), max(zero(T), sc.upper)
    elseif MOI.is_valid(model, c_si)
        # It is assumed that none of the other ConstraintIndexs are valid
        si::MOI.Semiinteger{T} = MOI.get(model, MOI.ConstraintSet(), c_si)
        return min(zero(T), si.lower), max(zero(T), si.upper)
    elseif MOI.is_valid(model, c_lt)
        lt::MOI.LessThan{T} = MOI.get(model, MOI.ConstraintSet(), c_lt)
        # It is valid to have both LessThan and GreaterThan constraints on the
        # same variable.
        if MOI.is_valid(model, c_gt)
            gt_1::MOI.GreaterThan{T} = MOI.get(model, MOI.ConstraintSet(), c_gt)
            return gt_1.lower, lt.upper
        else
            return typemin(T), lt.upper
        end
    elseif MOI.is_valid(model, c_gt)
        gt_2::MOI.GreaterThan{T} = MOI.get(model, MOI.ConstraintSet(), c_gt)
        return gt_2.lower, typemax(T)
    else
        return typemin(T), typemax(T)
    end
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
@testset "NamedMetaAxisArray" begin
    nma1 = NamedMetaAxisArray{(:dim_1,:dim_2)}(ones(Int, 2, 2), ["a", "b"], [:one, :two]);
    nma2 = NamedMetaAxisArray(ones(Int, 2, 2), (dim_1 = ["a", "b"], dim_2 = [:one, :two]));
    nma3 = NamedMetaAxisArray{(:dim_1,:dim_2),Int}(undef, (["a", "b"], [:one, :two]));
    nma4 = NamedMetaAxisArray{(:dim_1,:dim_2),Int}(undef, ["a", "b"], [:one, :two]);
    nma5 = NamedMetaAxisArray{(:dim_1,:dim_2),Int,2}(undef, (["a", "b"], [:one, :two]));
    nma6 = NamedMetaAxisArray{(:dim_1,:dim_2),Int,2}(undef, ["a", "b"], [:one, :two]);

    @test typeof(nma1) <: typeof(nma2)
    @test typeof(nma3) <: typeof(nma4) <: typeof(nma5) <: typeof(nma6)
    @test axes(nma1) == axes(nma2) == axes(nma3) == axes(nma4)
    @test dimnames(nma1) == dimnames(nma2) == dimnames(nma3) == dimnames(nma4)

    @test dimnames(NamedMetaAxisArray(ones(2,2))) == (:dim_1, :dim_2)
    @test axes_keys(NamedMetaAxisArray(ones(Int, 2, 2), (["a", "b"], [:one, :two]))) == (["a", "b"], [:one, :two])
    @test axes_keys(NamedMetaAxisArray(ones(Int, 2, 2), ["a", "b"], [:one, :two])) == (["a", "b"], [:one, :two])
end

@testset "NamedDims" begin
    A = NamedAxisArray(reshape(1:24, 2, 3, 4), x=["a", "b"], y =["one", "two", "three"], z=2:5);
    @test @inferred(ArrayInterface.parent_type(typeof(A))) <: typeof(parent(A))
    @test @inferred(A[1,1,1]) == @inferred(A["a", "one", ==(2)]) == 1
    @test @inferred(A[CartesianIndex(1,1,1)]) == 1
    # FIXME
    @test @inferred(A[["a", "b"], 1:2, 1:2]) ==
          @inferred(A[["a", "b"], 1:2, 1:2, 1]) ==
          @inferred(parent(A)[1:2, 1:2, 1:2]) ==
          @inferred(parent(A)[["a", "b"], 1:2, 1:2, 1])

    @test @inferred(A[1:10]) == 1:10

    @test keys(@inferred(AxisIndices.named_axes(A))) == (:x,:y,:z)
    @test keys(@inferred(AxisIndices.named_axes(parent(A)))) == (:dim_1, :dim_2, :dim_3)

    @test dimnames(@inferred(NamedAxisArray{(:x, :y),Int}(undef, 1:2, 1:2))) == (:x, :y)
    @test dimnames(@inferred(NamedAxisArray{(:x, :y),Int}(undef, (2,2)))) == (:x, :y)
    @test dimnames(@inferred(NamedAxisArray{(:x, :y),Int,2}(undef, 1:2, 1:2))) == (:x, :y)
end

@testset "MetaAxisArray" begin
    meta_axis_array = MetaAxisArray(ones(2,2));
    @test typeof(parent(meta_axis_array)) <: typeof(AxisArray(ones(2, 2)))

    axis_array = AxisArray(ones(2, 2), ["a", "b"], [:one, :two])
    meta_axis_array = MetaAxisArray(ones(2, 2), ["a", "b"], [:one, :two])
    @test typeof(parent(meta_axis_array)) <: typeof(axis_array)

    axis_array = AxisArray{Int}(undef, ["a", "b"], [:one, :two])
    meta_axis_array = MetaAxisArray{Int}(undef, ["a", "b"], [:one, :two])
    @test typeof(parent(meta_axis_array)) <: typeof(axis_array)

    axis_array = AxisArray{Int,2}(undef, ["a", "b"], [:one, :two])
    meta_axis_array = MetaAxisArray{Int,2}(undef, ["a", "b"], [:one, :two])
    @test typeof(parent(meta_axis_array)) <: typeof(axis_array)
end

ext: jl


text:
type Space
    ptr::Ptr{Void}
end

type Body
    ptr::Ptr{Void}
end

type Arbiter
    ptr::Ptr{Void}
end

abstract Shape
type CircleShape <: Shape
    ptr::Ptr{Void}
end

type SegmentShape <: Shape
    ptr::Ptr{Void}
end

type PolyShape <: Shape
    ptr::Ptr{Void}
end

abstract Constraint
type SimpleMotor <: Constraint
    ptr::Ptr{Void}
end

type PivotJoint <: Constraint
    ptr::Ptr{Void}
end

type PinJoint <: Constraint
    ptr::Ptr{Void}
end

type GearJoint <: Constraint
    ptr::Ptr{Void}
end

type SlideJoint <: Constraint
    ptr::Ptr{Void}
end

type DampedSpring <: Constraint
    ptr::Ptr{Void}
end

type CollisionHandler
    ptr::Ptr{Void}
end

type BB
    l::Cdouble
    b::Cdouble
    r::Cdouble
    t::Cdouble
end

type Vect
    x::Cdouble
    y::Cdouble
end

type Transform
    a::Cdouble
    b::Cdouble
    c::Cdouble
    d::Cdouble
    tx::Cdouble
    ty::Cdouble
end

type Mat2x2
    a::Cdouble
    b::Cdouble
    c::Cdouble
    d::Cdouble
end

export Vect, Body, Space, Shape, CircleShape, SegmentShape,
PolyShape, Transform, Mat2x2, BB, CollisionHandler, Constraint,
SimpleMotor, PivotJoint, PinJoint, GearJoint, SlideJoint,
DampedSpring


text:
# This file is part of TaylorSeries.jl, MIT licensed
#
# Tests for TaylorSeries

testfiles = (
    "onevariable.jl",
    "manyvariables.jl",
    "mixtures.jl",
    "mutatingfuncts.jl",
    "intervals.jl",
    "broadcasting.jl",
    "identities_Euler.jl",
    "fateman40.jl"
    )

for file in testfiles
    # Skipping tests with intervals, since IntervalArithmetic.jl requires Julia v1.1+
    VERSION < v"1.1" && file == "intervals.jl" && continue
    include(file)
end


text:
function step(x, s, n)
	x_trj = zeros(2, n+1)
	x_trj[:,1] = x
	for i = 2:n+1
		x, y = x_trj[1,i-1],x_trj[2,i-1]
		sx, sy = sin(x), sin(2*y)/2
		x1 = (2*x +
			(s[1] + s[2]*sy)*sx)
		y1 = (0.5*y +
			(s[4] + s[3]*sx)*sy)

		x_trj[1,i] = x < pi ? x1 : x1 - 2*pi
		x_trj[2,i] = x < pi ? y1 : y1 + pi

		x_trj[1,i] = x_trj[1,i] % (2*pi)
		x_trj[2,i] = x_trj[2,i] % (2*pi)

	end
	return x_trj
end
function next(u, s)
	x, y = u[1], u[2]
	sx, sy = sin(x), sin(2*y)/2
	x1 = (2*x +
			(s[1] + s[2]*sy)*sx)
	y1 = (0.5*y +
			(s[4] + s[3]*sx)*sy)

	x_next = x < pi ? x1 : x1 - 2*pi
	y_next = x < pi ? y1 : y1 + pi

	x_next = x_next % (2*pi)
	y_next = y_next % (2*pi)

	return [x_next, y_next]
end

function dstep(u::Array{Float64,1},s::Array{Float64,1})
	du = zeros(2,2)
	x, y = u[1], u[2]
	sx, sy = sin(x), sin(2*y)/2
	dsx, dsy = cos(x), cos(2*y)
	du[1,1] = 2 + s[1]*dsx + s[2]*sy*dsx
	du[1,2] = s[2]*sx*dsy
	du[2,1] = s[3]*dsx*sy
	du[2,2] = 0.5 + s[4]*dsy + s[3]*sx*dsy
	return du
end

function dstep(u::Array{Float64,2},s::Array{Float64,1})
	n = size(u)[2]
	du = zeros(2,2,n)
	x, y = view(u,1,:),view(u,2,:)
	sx, sy = sin.(x), sin.(2*y)/2
	dsx, dsy = cos.(x), cos.(2*y)
	du[1,1,:] = 2 .+ s[1]*dsx .+
				s[2]*sy.*dsx
	du[1,2,:] = s[2]*sx.*dsy
	du[2,1,:] = s[3]*dsx.*sy
	du[2,2,:] = 0.5 .+ s[4]*dsy .+
				s[3]*sx.*dsy
	return du
end
function pert(u::Array{Float64,2}, p::Int64)
	n = size(u)[2]
	x, y = view(u,1,:), view(u,2,:)
	sx, sy = sin.(x), sin.(2*y)./2
	if p==1
		return reshape([sx zeros(n)],n,2)'
	elseif p==2
		return reshape([sy.*sx zeros(n)],n,2)'
	elseif p==3
		return reshape([zeros(n) sy.*sx],n,2)'
	elseif p==4
		return reshape([zeros(n) sy],n,2)'
	else
		println("parameter perturbation is only defined
				for indices 1 through 4")
		return 0
	end
end
function pert(u::Array{Float64,1}, p::Int64)
	x, y = u[1], u[2]
	sx, sy = sin(x), sin(2*y)/2
	if p==1
		return [sx, 0]
	elseif p==2
		return [sy.*sx, 0]
	elseif p==3
		return [0, sy.*sx]
	elseif p==4
		return [0, sy]
	else
		println("parameter perturbation is only defined
				for indices 1 through 4")
		return 0
	end
end
function dpert_x(u::Array{Float64,1},s::Array{Float64,1})
	# here we asssume s = [a, a, 0, a]
	x, y = u[1], u[2]
	sx, sy = sin(x), sin(2*y)/2
	cx = cos(x)
	dxu = 2 .+ s[1]*cx .+ s[2]*sy.*cx
	dxpx1 = cx + cx*sy
	return [dxpx1/dxu, 0.]
end
function d2step(u::Array{Float64,2}, s::Array{Float64,1})
	n = size(u)[2]
	ddu = zeros(2,4,n)
	for i = 1:n
		x, y = u[1,i], u[2,i]
		sx, sy = sin(x), sin(2*y)/2
		dsx, dsy = cos(x), cos(2*y)
		dxdsx, dydsy = -sin(x), -2*sin(2*y)
		ddu[:,1,i] = [s[1]*dxdsx + s[2]*sy*dxdsx, s[2]*dsy*dsx]
		ddu[:,2,i] = [s[3]*dxdsx*sy, s[3]*dsx*dsy]
		ddu[:,3,i] = [s[2]*dsx*dsy, s[2]*sx*dydsy]
		ddu[:,4,i] = [s[3]*dsx*dsy, s[4]*dydsy + s[3]*sx*dydsy]
	end
	return ddu
end
function d2step(u::Array{Float64,1}, s::Array{Float64,1})
	ddu = zeros(2,4)
	x, y = u[1], u[2]
	sx, sy = sin(x), sin(2*y)/2
	dsx, dsy = cos(x), cos(2*y)
	dxdsx, dydsy = -sin(x), -2*sin(2*y)
	ddu[:,1] = [s[1]*dxdsx + s[2]*sy*dxdsx, s[2]*dsy*dsx]
	ddu[:,2] = [s[3]*dxdsx*sy, s[3]*dsx*dsy]
	ddu[:,3] = [s[2]*dsx*dsy, s[2]*sx*dydsy]
	ddu[:,4] = [s[3]*dsx*dsy, s[4]*dydsy + s[3]*sx*dydsy]
	return ddu
end

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
myread(fn) = readSBML(fn, doc -> begin
    set_level_and_version(3, 2)(doc)
    convert_simplify_math(doc)
end)

# assignment rule
fn = "data/00038-sbml-l3v2.xml" # this case is for observable eqs
m = myread(fn)
rs = ReactionSystem(m)

sys = convert(ODESystem, rs; include_zero_odes = false)
@test length(equations(sys)) == 3
ssys = structural_simplify(sys)
@test length(observed(ssys)) == 1
prob = ODEProblem(ssys, [], (0, 10.0))
sol = solve(prob, Tsit5())

@variables t S3(t)
obsvar_sol = sol[S3]
@test !all(isequal(first(obsvar_sol)), obsvar_sol) # check its nonconstant

# rate rule
fn = "data/00031-sbml-l3v2.xml" # this case is for observable eqs
m = myread(fn)
rs = ReactionSystem(m)
sys = convert(ODESystem, rs; include_zero_odes = false)
@variables t S1(t)
D = Differential(t)
@test equations(sys) == [D(S1) ~ 7]

# algebraic rule
fn = "data/00039-sbml-l3v2.xml" # this case is for observable eqs
m = myread(fn)
rs = ReactionSystem(m)
sys = convert(ODESystem, rs; include_zero_odes = false)
ssys = structural_simplify(sys)
prob = ODEProblem(ssys, [], (0, 10.0))
sol = solve(prob, Tsit5())
@variables t S2(t)
obsvar_sol = sol[S2]


text:
using Test
import sbp
using sbp.LSRK4: lsrk4

y = zeros(1)

@testset "Integrate cubic monomial" begin

        e = ones(1)

        function g(y::AbstractArray, t::Float64)
                return e * t ^ 3
        end

        nt = 10
        dt = 0.1
        Y = zeros(nt+1)
        u = zeros(1)
        dy = zeros(1)

        ans = zeros(nt+1)
        for k in 1:nt
                ans[k+1] = 0.25 * (k * dt) ^ 4
        end

        tk=0.0
        for k in 1:nt
                tk = (k - 1) * dt
                @time global y = lsrk4(g, y, dy, tk, dt)
                Y[k+1] = y[1]
        end
        @test isapprox(ans, Y)
end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
immutable Geometric <: DiscreteUnivariateDistribution
    prob::Float64
    function Geometric(p::Real)
        zero(p) < p < one(p) || error("prob must be in (0, 1)")
    	new(float64(p))
    end
end

Geometric() = Geometric(0.5) # Flips of a fair coin

@_jl_dist_1p Geometric geom

function cdf(d::Geometric, q::Real)
    q < zero(q) ? 0.0 : -expm1(log1p(-d.prob) * (floor(q) + 1.0))
end

function ccdf(d::Geometric, q::Real)
    q < zero(q) ? 1.0 : exp(log1p(-d.prob) * (floor(q + 1e-7) + 1.0))
end

entropy(d::Geometric) = (-xlogx(1.0 - d.prob) - xlogx(d.prob)) / d.prob

insupport(::Geometric, x::Real) = isinteger(x) && zero(x) <= x
insupport(::Type{Geometric}, x::Real) = isinteger(x) && zero(x) <= x

kurtosis(d::Geometric) = 6.0 + d.prob^2 / (1.0 - d.prob)

mean(d::Geometric) = (1.0 - d.prob) / d.prob

function median(d::Geometric)
    iceil(-1.0 / log(2.0, 1.0 - d.prob)) - 1
end

mode(d::Geometric) = 0
modes(d::Geometric) = [0]

function mgf(d::Geometric, t::Real)
    p = d.prob
    if t >= -log(1.0 - p)
        error("MGF does not exist for all t")
    end
    (p * exp(t)) / (1.0 - (1.0 - p) * exp(t))
end

function cf(d::Geometric, t::Real)
    p = d.prob
    (p * exp(im * t)) / (1.0 - (1.0 - p) * exp(im * t))
end

skewness(d::Geometric) = (2.0 - d.prob) / sqrt(1.0 - d.prob)

var(d::Geometric) = (1.0 - d.prob) / d.prob^2

## Fit model

immutable GeometricStats
    sx::Float64
    tw::Float64

    GeometricStats(sx::Real, tw::Real) = new(float64(sx), float64(tw))
end

suffstats{T<:Integer}(::Type{Geometric}, x::Array{T}) = GeometricStats(sum(x), length(x))

function suffstats{T<:Integer}(::Type{Geometric}, x::Array{T}, w::Array{Float64})
    n = length(x)
    if length(w) != n
        throw(ArgumentError("Inconsistent argument dimensions."))
    end
    sx = 0.
    tw = 0.
    for i = 1:n
        wi = w[i]
        sx += wi * x[i]
        tw += wi
    end
    GeometricStats(sx, tw)
end

fit_mle(::Type{Geometric}, ss::GeometricStats) = Geometric(1.0 / (ss.sx / ss.tw + 1.0))


text:
abstract type AbstractHashlessTileCoder <: AbstractFeatureConstructor
end

"""
    HashlessTileCoder(tiles_per_dim, bounds_per_dim, num_tilings, wrap, offset)

This is a struct for tile-coding raw states to a list of active indices. The constructor described here uses the helper function `hashlesstilecoder_args` to convert human-readable variables into the quantities used to compute the tiled indices.

# Arguments
- `tiles_per_dim::Vector{<:Integer}`: A vector describing the number of tiles per dimension of the state input. To tile two dimensions of state together with 2 tiles on the first dimension and 3 tiles on the second, set `tiles_per_dim = [2, 3]`.
- `bounds_per_dim::Matrix{<:Real}`: A matrix of real values denoting the upper and lower bounds of each state dimension. If both state dimensions are bounded in `[0, 1]`, pass `bounds_per_dim = [0 0 ; 1 1]`.
- `num_tilings::Integer`: The tilecoder will map state to `num_tilings` indices.
- `wrap::Union{AbstractVector{Bool}, Nothing}`: Nothing, indicating no wrapping, or a vector of booleans indicating whether that index of the state input should be wrapped. Setting `wrap[k] = true` means that the two ends of `bounds_per_dim[:,k]` are equivalent. For example, if `state[k]` describes an angle between 0 and 2π, `wrap[k]` should be set to `true`.
- `offset`: A function that describes how each tiling should be offset from the others. The recommended setting is to use the default odd-numbered offsets, called asymmetric displacement by [Parks and Militzer](https://doi.org/10.1016/S1474-6670(17)54222-6).

"""
mutable struct HashlessTileCoder <: AbstractHashlessTileCoder
    limits::Matrix{Float64}
    norm_dims::Vector{Float64}
    tiling_dims::Vector{Int32}
    wrap_any_dims::Bool
    offsets::Matrix{Float64}
    tiling_loc::Vector{Int}
    tile_loc::Vector{Int}
    num_features::Int
    num_active_features::Int

    function HashlessTileCoder(a...;k...)
        l, n, t, w, o, tgl, til, nf, naf = hashlesstilecoder_args(a...;k...)
        new(l, n, t, w, o, tgl, til, nf, naf)
    end
end

feature_size(fc::T) where {T<:HashlessTileCoder} = fc.num_features
Base.size(fc::HashlessTileCoder) = fc.num_features

"""
Helper function to make constructing each feature vector later easier
"""
function hashlesstilecoder_args(tiles_per_dim::Vector{<:Integer},
                                bounds_per_dim::Matrix{<:Real},
                                num_tilings::Integer;
                                wrap::Union{AbstractVector{Bool}, Nothing}=nothing,
                                offset= n -> collect(1:2:2*n-1) # 1st n odd nums
                                )
    n = length(tiles_per_dim)

    # these normalize the ith input float to be between 0 and dim[i] + 1
    limits = bounds_per_dim
    norm_dims = tiles_per_dim ./ (limits[2, :] .- limits[1, :])

    # wrapping means not adding 1 to the ith dim
    if wrap == nothing
        bonus = ones(Bool, n)
    else
        bonus = .!wrap
    end
    wrap_any_dims = any(.!bonus)
    tiling_dims = tiles_per_dim .+ bonus

    # displacement matrix; default is assymetric displacement a la Parks
    # and Militzer https://doi.org/10.1016/S1474-6670(17)54222-6
    offset_vec = offset(n)
    offsets = (offset_vec
               .* hcat([collect(0:num_tilings-1) for _ in 1:n]...)'
               ./ num_tilings
               .% 1)

    # these send each displaced float to the proper index
    tiling_loc = collect(0:num_tilings-1) .* prod(tiling_dims)
    tile_loc = [prod(tiling_dims[1:i-1]) for i in 1:n]

    # the total number of indices needed
    num_features = num_tilings * prod(tiling_dims)
    num_active_features = num_tilings

    return (limits, norm_dims, tiling_dims, wrap_any_dims, offsets, tiling_loc,
            tile_loc, num_features, num_active_features)
end

"""
Constructs a list of active indices from the quantities defined by
`hashlesstilecoder_args`.
"""
function _create_features(fc::T, s) where {T<:AbstractHashlessTileCoder}
    if fc.wrap_any_dims
        # wrapping means modding by dim[i] instead of dim[i] + 1
        off_coords = map(x -> floor(Int, x),
                         ((s .- fc.limits[1, :])
                          .* fc.norm_dims
                          .+ fc.offsets)
                         .% fc.tiling_dims)
    else
        # don't need to mod here, because dim[i] + 1 is bigger than the
        # displaced floats
        off_coords = Int.(floor.((s .- fc.limits[1, :]) .* fc.norm_dims .+ fc.offsets))
    end

    return fc.tiling_loc .+ off_coords' * fc.tile_loc .+ 1
end
create_features(fc::HashlessTileCoder, s) = _create_features(fc, s)

"""
    HashlessTileCoderFull(tiles_per_dim, bounds_per_dim, num_tilings, wrap, offset)

This is a struct for tile-coding raw states. Unlike HashlessTileCoder, HashlessTileCoderFull returns a full binary array where ones denote that an index is active. The constructor described here uses the helper function `hashlesstilecoder_args` to convert human-readable variables into the quantities used to compute the tiled indices.

# Arguments
- `tiles_per_dim::Vector{<:Integer}`: A vector describing the number of tiles per dimension of the state input. To tile two dimensions of state together with 2 tiles on the first dimension and 3 tiles on the second, set `tiles_per_dim = [2, 3]`.
- `bounds_per_dim::Matrix{<:Real}`: A matrix of real values denoting the upper and lower bounds of each state dimension. If both state dimensions are bounded in `[0, 1]`, pass `bounds_per_dim = [0 0 ; 1 1]`.
- `num_tilings::Integer`: The tilecoder will map state to `num_tilings` indices.
- `wrap::Union{AbstractVector{Bool}, Nothing}`: Nothing, indicating no wrapping, or a vector of booleans indicating whether that index of the state input should be wrapped. Setting `wrap[k] = true` means that the two ends of `bounds_per_dim[:,k]` are equivalent. For example, if `state[k]` describes an angle between 0 and 2π, `wrap[k]` should be set to `true`.
- `offset`: A function that describes how each tiling should be offset from the others. The recommended setting is to use the default odd-numbered offsets, called asymmetric displacement by [Parks and Militzer](https://doi.org/10.1016/S1474-6670(17)54222-6).

"""
mutable struct HashlessTileCoderFull <: AbstractHashlessTileCoder
    limits::Matrix{Float64}
    norm_dims::Vector{Float64}
    tiling_dims::Vector{Int32}
    wrap_any_dims::Bool
    offsets::Matrix{Float64}
    tiling_loc::Vector{Int}
    tile_loc::Vector{Int}
    num_features::Int

    function HashlessTileCoderFull(a...;k...)
        l, n, t, w, o, tgl, til, nf = hashlesstilecoder_args(a...;k...)
        new(l, n, t, w, o, tgl, til, nf)
    end
end
feature_size(fc::T) where {T<:HashlessTileCoderFull} = fc.num_features

function create_features(fc::HashlessTileCoderFull, s)
    x = zeros(Int, fc.num_features)

    x[_create_features(fc, s)] .= 1

    return x
end


text:
using PrettyTables
using ProgressMeter

using RollingFunctions
using Clustering

if VERSION < v"1.1-DEV"
    isnothing(::Nothing) = true
    isnothing(::Any) = false
end

function estimate_convergence_rate(ρ, ϵ)
    i = findall(!iszero, ϵ)
    x = log10.(ρ[i])
    y = log10.(abs.(ϵ[i]))

    f = (x,y) -> ([x ones(length(x))] \ y)[1]

    # Find all slopes for consecutive groups of three errors
    slopes = rolling(f, x, y, 3)
    e = findfirst(s -> s < 0, slopes)
    if e !== nothing && e > 5
        slopes = slopes[1:max(1,e-1)]
    end

    k = 3 # Number of clusters requested
    fit = kmeans(reshape(slopes, 1, :), k)
    maximum(fit.centers)
end

# https://stackoverflow.com/a/52507859
function squeeze(A::AbstractArray)
    singleton_dims = tuple((d for d in 1:ndims(A) if size(A, d) == 1)...)
    dropdims(A, dims=singleton_dims)
end

function test_convergence_rates(fun::Function,
                                Rfun::Function, rₘₐₓ::Number,
                                Ns::AbstractVector{<:Integer},
                                ρcol::Int, # Which column contains the grid spacing?
                                errors::Vector{<:Tuple{Int,String,<:Real,<:Real}},
                                headers::Vector{String},
                                other_vectors::AbstractVector...;
                                convergence_rate_atol=1e-1,
                                plot_fun::Function = (i,data) -> nothing,
                                verbosity=1)
    ovs = Iterators.product(other_vectors...)

    progress = Progress(length(Ns))
    data = map(Ns) do N
        R,ρ = Rfun(rₘₐₓ, N)
        r = axes(R,1)
        cache = (R = R,
                 ρ = ρ,
                 V = vandermonde(R),
                 S = R'R,
                 r⁻¹ = R'*QuasiDiagonal(1 ./ r)*R,
                 poissons = Dict{Int,PoissonProblem}())
        d = map(o -> fun(cache, o...), ovs)
        verbosity > 1 || ProgressMeter.next!(progress)
        d
    end
    data = permutedims(squeeze(hcat(data...)))
    convergence_rates = map(enumerate(ovs)) do (i,o)
        case_data = vcat(filter(!isnothing,data[i,:])...)
        if isempty(case_data)
            ProgressMeter.next!(progress)
            return nothing
        end
        verbosity > 1 && pretty_table(case_data, headers)

        ρ = case_data[:,ρcol]
        label = join(string.(filter(c -> !(c isa Function), [o...])), ", ")

        plot_fun(i, case_data)

        e = map(errors) do (ecol,elabel,eexpected,rateexpected)
            curerr = case_data[:,ecol]
            minerr = minimum(abs, curerr)
            γ = estimate_convergence_rate(ρ, curerr)
            # We deem that the integration has passed if the minimal
            # error is below the requested threshold and the
            # convergence rate is approximately equal to or larger
            # than the requested order.
            pass = minerr < eexpected &&
                (isapprox(γ, rateexpected, atol=convergence_rate_atol) || γ > rateexpected)
            Any[minerr γ pass] # Without Any, pass will be coerced to Float64
        end |> e -> hcat(label, e...)

        verbosity > 1 || ProgressMeter.next!(progress)

        e
    end |> e -> vcat(filter(!isnothing, e)...)
    println()

    passjs = 1 .+ 3*(1:length(errors))

    if verbosity > 0
        pass = Highlighter((data,i,j) -> j>1 && mod(j-1,3)==0 && data[i,j],
                           bold = true, foreground = :green)
        fail = Highlighter((data,i,j) -> j>1 && mod(j-1,3)==0 && !data[i,j],
                           bold = true, foreground = :red)
        # formatter = merge(ft_printf("%2.3e", passjs .- 2),
        #                   ft_printf("%2.3f", passjs .- 1),
        #                   Dict(j => (v,i) -> v ? "✓" : "⨯" for j in passjs))

        fails = vec(.!reduce(&, convergence_rates[:,passjs], dims=2))
        sel = verbosity == 1 ? findall(fails) : (1:size(convergence_rates,1))
        if !isempty(sel)
            pretty_table(convergence_rates[sel, :],
                         vcat("Case", [["$elabel, minδ [$eexpected]", "$elabel, rate [$rateexpected]", "Pass"]
                                       for (_,elabel,eexpected,rateexpected) in errors]...),
                         highlighters=(pass,fail),
                         # formatter=formatter,
                         crop=:none)
        end
        fails = sum([count(.!convergence_rates[:,j]) for j in passjs])
        num_cases = length(errors)*size(convergence_rates,1)
        println("Number of fails: $fails ≡ $(100.0fails/num_cases)% fail rate")
    end

    all([all(convergence_rates[:,j]) for j in passjs])
end


text:
# integration with LinearAlgebra.jl

CuMatOrAdj{T} = Union{CuMatrix, LinearAlgebra.Adjoint{T, <:CuMatrix{T}}, LinearAlgebra.Transpose{T, <:CuMatrix{T}}}
CuOrAdj{T} = Union{CuVecOrMat, LinearAlgebra.Adjoint{T, <:CuVecOrMat{T}}, LinearAlgebra.Transpose{T, <:CuVecOrMat{T}}}

# matrix division

function Base.:\(_A::CuMatOrAdj, _B::CuOrAdj)
    A, B = copy(_A), copy(_B)
    A, ipiv = CUSOLVER.getrf!(A)
    return CUSOLVER.getrs!('N', A, ipiv, B)
end

# patch JuliaLang/julia#40899 to create a CuArray
# (see https://github.com/JuliaLang/julia/pull/41331#issuecomment-868374522)
if VERSION >= v"1.7-"
_zeros(::Type{T}, b::AbstractVector, n::Integer) where {T} = CUDA.zeros(T, max(length(b), n))
_zeros(::Type{T}, B::AbstractMatrix, n::Integer) where {T} = CUDA.zeros(T, max(size(B, 1), n), size(B, 2))
function Base.:\(F::Union{LinearAlgebra.LAPACKFactorizations{<:Any,<:CuArray},
                          Adjoint{<:Any,<:LinearAlgebra.LAPACKFactorizations{<:Any,<:CuArray}}},
                 B::AbstractVecOrMat)
    m, n = size(F)
    if m != size(B, 1)
        throw(DimensionMismatch("arguments must have the same number of rows"))
    end

    TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F)))
    FF = Factorization{TFB}(F)

    # For wide problem we (often) compute a minimum norm solution. The solution
    # is larger than the right hand side so we use size(F, 2).
    BB = _zeros(TFB, B, n)

    if n > size(B, 1)
        # Underdetermined
        copyto!(view(BB, 1:m, :), B)
    else
        copyto!(BB, B)
    end

    ldiv!(FF, BB)

    # For tall problems, we compute a least squares solution so only part
    # of the rhs should be returned from \ while ldiv! uses (and returns)
    # the complete rhs
    return LinearAlgebra._cut_B(BB, 1:n)
end
end

# qr

using LinearAlgebra: AbstractQ

# AbstractQ's `size` is the size of the full matrix,
# while `Matrix(Q)` only gives the compact Q.
# See JuliaLang/julia#26591 and JuliaGPU/CUDA.jl#969.
CuMatrix{T}(Q::AbstractQ{S}) where {T,S} = convert(CuArray, Matrix{T}(Q))
CuMatrix(Q::AbstractQ{T}) where {T} = CuMatrix{T}(Q)
CuArray{T}(Q::AbstractQ) where {T} = CuMatrix{T}(Q)
CuArray(Q::AbstractQ) = CuMatrix(Q)

# dot

function LinearAlgebra.dot(x::AnyCuArray{T1}, y::AnyCuArray{T2}) where {T1,T2}
    n = length(x)
    n==length(y) || throw(DimensionMismatch("dot product arguments have lengths $(length(x)) and $(length(y))"))

    # custom kernel using simple linear indexing and atomic additions,
    # resulting in about 10% speed-up compared to a simple mapreduce.
    function kernel(x, y, res::AbstractArray{T}, shuffle) where {T}
        local_val = zero(T)

        # grid-stride loop
        i = threadIdx().x + (blockIdx().x - 1i32)*blockDim().x
        while i <= length(x)
            @inbounds local_val += dot(x[i], y[i])
            i += blockDim().x * gridDim().x
        end

        val = reduce_block(+, local_val, zero(T), shuffle)
        if threadIdx().x == 1i32
            # NOTE: introduces nondeterminism
            @inbounds @atomic res[] += val
        end

        return
    end

    dev = device()
    let T = promote_type(T1, T2)
        # only use the above kernel if we don't care about determinism
        # and if atomic operations are supported on these inputs
        atomic = T <: Union{Int16, Int32, Int64, Float16, Float32, Float64}
        if math_mode() == PEDANTIC_MATH || !atomic
            return mapreduce((x,y)->dot(x, y), +, x, y)
        end

        res = zeros(T, 1)

        # be conservative about using shuffle instructions
        shuffle = T <: Union{Bool,
                             UInt8, UInt16, UInt32, UInt64, UInt128,
                             Int8, Int16, Int32, Int64, Int128,
                             Float16, Float32, Float64,
                             ComplexF16, ComplexF32, ComplexF64}

        # how many threads do we want?
        # reduce_block(shuffle=true) requires the block to consist of full warps.
        wanted_threads = shuffle ? nextwarp(dev, n) : n
        function compute_threads(max_threads)
            if wanted_threads > max_threads
                shuffle ? prevwarp(dev, max_threads) : max_threads
            else
                wanted_threads
            end
        end

        # how many threads can we launch?
        kernel = @cuda launch=false kernel(x, y, res, Val(shuffle))
        compute_shmem(threads) = shuffle ? 0 : threads*sizeof(T)
        config = launch_configuration(kernel.fun; shmem=compute_shmem∘compute_threads)
        threads = compute_threads(config.threads)
        blocks = min(config.blocks, cld(n, config.blocks))
        shmem = compute_shmem(threads)
        kernel(x, y, res, Val(shuffle); threads, blocks, shmem)

        @allowscalar res[]
    end
end


text:
using Mathieu
using Test, Printf, DelimitedFiles
using Statistics: mean, median
using Base.MathConstants: φ

const refdir = joinpath(@__DIR__,"ref")

## Test vs Mathematica
n1 = 0:24
q1 = 0:2φ:80
z1 = LinRange(0,4φ,10)
cases = ((q->Mathieu.ap(n1,q),           "MathieuCharA_Periodic"),
         (q->Mathieu.aa(n1,q),           "MathieuCharA_Antiperiodic"),
         (q->Mathieu.ba(n1,q),           "MathieuCharB_Antiperiodic"),
         (q->Mathieu.bp(n1,q),           "MathieuCharB_Periodic"),
         (q->Mathieu.cep(n1,q,z1),       "MathieuC_Periodic"),
         (q->Mathieu.cea(n1,q,z1),       "MathieuC_Antiperiodic"),
         (q->Mathieu.sea(n1,q,z1),       "MathieuS_Antiperiodic"),
         (q->Mathieu.sep(n1,q,z1),       "MathieuS_Periodic"),
         (q->Mathieu.cep_prime(n1,q,z1), "MathieuCPrime_Periodic"),
         (q->Mathieu.cea_prime(n1,q,z1), "MathieuCPrime_Antiperiodic"),
         (q->Mathieu.sea_prime(n1,q,z1), "MathieuSPrime_Antiperiodic"),
         (q->Mathieu.sep_prime(n1,q,z1), "MathieuSPrime_Periodic"),
         )
for (fjl,file) in cases
    println(string("Comparing ",file," with Mathematica:"))
    # Read Mathematica reference values
    arr_ref = readdlm(joinpath(refdir,"$file.csv"),',',Float64)
    # Calculate corresponding values with Mathieu.jl
    arr_jul = mapreduce(fjl,hcat,q1)
    # Make sure the values are approximately equal
    @test arr_jul ≈ arr_ref
    a = abs.(arr_jul-arr_ref)./eps.(arr_jul)
    println(@sprintf("  mean error: %.2f eps, median: %d eps, max: %d eps.",mean(a),median(a),maximum(a)))
end

println(string("Comparing ","Int_Per_PerPrime"," with Mathematica:"))
# Read Mathematica reference values
arr_ref = readdlm(joinpath(refdir,"Int_Per_PerPrime.csv"),',',Float64)
# Calculate corresponding values with Mathieu.jl
arr_jul = Mathieu.int_p_pprime(n1,n1,223/100,[0.5,3])
# Make sure the values are approximately equal
@test arr_jul ≈ arr_ref[1:25,1:25]
a = abs.(arr_jul-arr_ref[1:25,1:25])./eps.(arr_jul)
println(@sprintf("  mean error: %.2f eps, median: %d eps, max: %d eps.",mean(a),median(a),maximum(a)))

# Test q change of sign
n2 = [0:20,7:21,10:35,12:46,50:100,100:150,150:200]
q2 = LinRange(0,100,202) #0:2φ:30
z2 = LinRange(0,4φ,10)
for nn in n2
    println("Testing change of sign for n = $nn:")
    N2 = maximum(nn)+12 # Make N different to hit different eigendecompositions
    print("  ap")
    A1 = mapreduce(q->Mathieu.ap(nn,q),hcat,-q2)
    A2 = mapreduce(q->Mathieu.ap(nn,q,N2),hcat,q2)
    @test A1 ≈ A2
    print(", aa, ba")
    A1 = mapreduce(q->Mathieu.aa(nn,q),hcat,-q2)
    A2 = mapreduce(q->Mathieu.ba(nn,q,N2),hcat,q2)
    @test A1 ≈ A2
    print(", bp")
    A1 = mapreduce(q->Mathieu.bp(nn,q),hcat,-q2)
    A2 = mapreduce(q->Mathieu.bp(nn,q,N2),hcat,q2)
    @test A1 ≈ A2
    print(", cep")
    A1 = mapreduce(q->Mathieu.cep(nn,q,z2),vcat,-q2)
    A2 = mapreduce(q->Mathieu.cep(nn,q,π/2 .- z2,N2).*(-1).^nn',vcat,q2)
    @test A1 ≈ A2
    print(", cea")
    A1 = mapreduce(q->Mathieu.cea(nn,q,z2),vcat,-q2)
    A2 = mapreduce(q->Mathieu.sea(nn,q,π/2 .- z2,N2).*(-1).^nn',vcat,q2)
    @test A1 ≈ A2
    print(", sea")
    A1 = mapreduce(q->Mathieu.sea(nn,q,z2),vcat,-q2)
    A2 = mapreduce(q->Mathieu.cea(nn,q,π/2 .- z2,N2).*(-1).^nn',vcat,q2)
    @test A1 ≈ A2
    println(", sep.")
    A1 = mapreduce(q->Mathieu.sep(nn,q,z2),vcat,-q2)
    A2 = mapreduce(q->Mathieu.sep(nn,q,π/2 .- z2,N2).*(-1).^nn',vcat,q2)
    @test A1 ≈ A2
end


text:
## This page is a work in progress. In particular, need to implement switching function that is typically
## employed within LAMMPS.

struct ZBL <: EmpiricalPotential
    Z_1
    Z_2
    e
    ϵ0
    rcutoff
    species :: AbstractVector
end

get_parameters(zbl::ZBL) = Parameter{}(())
set_parameters(p::Parameter{}, zbl::ZBL) = zbl

deserialize_parameters(p::Parameter{()}, zbl::ZBL) = []
serialize_parameters(p::Vector, zbl::ZBL) = Parameter{()}( () )

get_hyperparameters(zbl::ZBL) = Parameter{:rcutoff}((zbl.rcutoff, ))
set_hyperparameters(p::Parameter{(:rcutoff,)}, zbl::ZBL) = ZBL(zbl.Z_1, zbl.Z_2, zbl.e, zbl.ϵ0, p.rcutoff, zbl.species)

deserialize_hyperparameters(p::Parameter{(:rcutoff,)}, zbl::ZBL) = [p.rcutoff]
serialize_hyperparameters(p::Vector, zbl::ZBL) = Parameter{(:rzblutoff, )}( (p[1],) )

phi(d::AbstractFloat) =  0.18175 * e^(-3.19980 * d) + 0.50986 * e^(-0.94229 * d) + 0.28022 * e^(-0.40290 * d) + 0.02817 * e^(-0.20162 * d)
dphi_dr(d::AbstractFloat) = -3.19980 * 0.18175 * e^(-3.19980 * d) - 0.94229 * 0.50986 * e^(-0.94229 * d) - 0.40290 * 0.28022 * e^(-0.40290 * d) - 0.20162 * 0.02817 * e^(-0.20162 * d)
############################# Energies ##########################################

function potential_energy(R::AbstractFloat, zbl::ZBL)
    a = 0.8854 * 0.529 / (Z_1^(0.23) + Z_2^(0.23))
    d = R / a
    zbl.Z_1 * zbl.Z_2 * e^2 * phi(d) / (4.0 * π * zbl.ϵ0 * R)
end

############################### Forces ##########################################

force(r::SVector{3,<:AbstractFloat}, zbl::ZBL) = force(r, norm(r), zbl)

function force(R::AbstractFloat, r::SVector{3,<:AbstractFloat}, zbl::ZBL)
    a = 0.8854 * 0.529 / (Z_1^(0.23) + Z_2^(0.23))
    d = R / a
    f = zbl.Z_1 * zbl.Z_2 * e^2 * phi(d) / (4.0 * π * zbl.ϵ0 * R^2)
    f += zbl.Z_1 * zbl.Z_2 * e^2 * dphi_dr(d) / (a * 4.0 * π * zbl.ϵ0 * R)
    SVector(f .* r / R)
end

# ##############################   Gradients  ###################################
# TO DO


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# 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.
# *****************************************************************************

# According to file format specified in Matlab script
# written by Mykel Kochenderfer, mykel@stanford.edu
#
# This script converts a log file produced by RLES and converts it to a "waypoints" file
# of the following format.
#
#   WAYPOINTS FILE:
#   The waypoints file contains a set of encounters. Each encounter is
#   defined by a set of waypoints associated with a fixed number of
#   aircraft. The waypoints are positions in space according to a fixed,
#   global coordinate system. All distances are in feet. Time is specified
#   in seconds since the beginning of the encounter. The file is organized
#   as follows:
#
#   [Header]
#   uint32 (number of encounters)
#   uint32 (number of aircraft)
#       [Encounter 1]
#           [Initial positions]
#               [Aircraft 1]
#               double (north position in feet)
#               double (east position in feet)
#               double (altitude in feet)
#               ...
#               [Aircraft n]
#               double (north position in feet)
#               double (east position in feet)
#               double (altitude in feet)
#           [Updates]
#               [Aircraft 1]
#               uint16 (number of updates)
#                   [Update 1]
#                   double (time in seconds)
#                   double (north position in feet)
#                   double (east position in feet)
#                   double (altitude in feet)
#                   ...
#                   [Update m]
#                   double (time in seconds)
#                   double (north position in feet)
#                   double (east position in feet)
#                   double (altitude in feet)
#               ...
#               [Aircraft n]
#                   ...
#       ...
#       [Encounter k]

module Log_To_Waypoints

export log_to_waypoints

using ..DefineSave
using ..SaveHelpers

include("corr_aem_save_scripts.jl")

using JSON
using DataFrames

function log_to_waypoints{T<:AbstractString}(filenames::Vector{T};
    outfile::AbstractString="waypoints.dat")

    d = trajLoad(filenames[1]) #use the first one as a reference
    num_aircraft = get_num_aircraft(d)
    num_encounters = length(filenames) #one encounter per log file
    encounters = Array(Dict{String, Array{Float64, 2}}, num_aircraft, num_encounters)

    #encounter i
    for (i, file) in enumerate(filenames)
        d = trajLoad(file)
        #make sure all of them have the same number of aircraft
        @assert num_aircraft == get_num_aircraft(d)

        #aircraft j
        for j = 1:num_aircraft
            encounters[j, i] = Dict{String, Array{Float64, 2}}()
            encounters[j, i]["initial"] = j2w_initial(d, j)
            encounters[j, i]["update"] = j2w_update(d, j)'
        end
    end
    save_waypoints(outfile, encounters, numupdatetype=UInt16)
    encounters
end

function j2w_initial(d::TrajLog, aircraft_number::Int64)
  #d is the loaded encounter
  out = Array(Float64, 1, 3)
  wm = get_log(d, "WorldModel", aircraft_number)

  pos_north     = wm[1, :y]
  pos_east      = wm[1, :x]
  altitude      = wm[1, :z]
  out[1, :] = [pos_north, pos_east, altitude]
  out
end

function j2w_update(d::TrajLog, aircraft_number::Int64)
  #d is the loaded encounter,
  wm = get_log(d, "WorldModel", aircraft_number)
  t_end = nrow(wm)
  out = Array(Float64, t_end-1, 4)
  for t = 2:t_end #ignore init values
    t_         = wm[t, :t]
    pos_north  = wm[t, :y]
    pos_east   = wm[t, :x]
    altitude   = wm[t, :z]
    out[t - 1, :] = Float64[t_, pos_north, pos_east, altitude]
  end
  out
end

function log_to_waypoints(filename::AbstractString)
  log_to_waypoints([filename], outfile = string(getSaveFileRoot(filename), "_waypoints.dat"))
end

end #module

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
# Autogenerated wrapper script for libsharp2_jll for i686-linux-gnu
export libsharp2

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("libsharp2")
JLLWrappers.@declare_library_product(libsharp2, "libsharp2.so.0")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libsharp2,
        "lib/libsharp2.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using YAML

load_meta(name) = YAML.load_file(
            joinpath(artifact"dash_resources", "$(name).yaml")
            )

deps_path(name) = joinpath(artifact"dash_resources", "$(name)_deps")

dash_dependency_resource(meta) = Resource(
        relative_package_path = meta["relative_package_path"],
        external_url = meta["external_url"]
    )

nothing_if_empty(v) = isempty(v) ? nothing : v

dash_module_resource(meta) = Resource(
        relative_package_path = nothing_if_empty(get(meta, "relative_package_path", "")),
        external_url = nothing_if_empty(get(meta, "external_url", "")),
        dev_package_path = nothing_if_empty(get(meta, "dev_package_path", "")),
        dynamic = get(meta, "dynamic", nothing),
        type = Symbol(meta["type"]),
        async = haskey(meta, "async") ? string(meta["async"]) : nothing
    )

dash_module_resource_pkg(meta; resource_path, version) = ResourcePkg(
    meta["namespace"],
    resource_path, version = version,
    dash_module_resource.(meta["resources"])
)

function setup_renderer_resources()
    renderer_meta = _metadata.dash_renderer
    renderer_resource_path = joinpath(artifact"dash_resources", "dash_renderer_deps")
    DashBase.main_registry().dash_dependency = (
        dev = ResourcePkg(
            "dash_renderer",
            renderer_resource_path, version = renderer_meta["version"],
            dash_dependency_resource.(renderer_meta["js_dist_dependencies"]["dev"])
        ),
        prod = ResourcePkg(
            "dash_renderer",
            renderer_resource_path, version = renderer_meta["version"],
            dash_dependency_resource.(renderer_meta["js_dist_dependencies"]["prod"])
        )
    )

    DashBase.main_registry().dash_renderer = dash_module_resource_pkg(
            renderer_meta["deps"][1],
            resource_path = renderer_resource_path,
            version = renderer_meta["version"]
            )
end

function load_all_metadata()
    dash_meta = load_meta("dash")
    renderer_meta = load_meta("dash_renderer")
    components = Dict{Symbol, Any}()
    for comp in dash_meta["embedded_components"]
        components[Symbol(comp)] = filter(v->v.first!="components", load_meta(comp))
    end
    return (
        dash = dash_meta,
        dash_renderer = renderer_meta,
        embedded_components = (;components...)
    )
end

function setup_dash_resources()
    meta = _metadata.dash
    path = deps_path("dash")
    version = meta["version"]
    for dep in meta["deps"]
        DashBase.register_package(
            dash_module_resource_pkg(
                dep,
                resource_path = path,
                version = version
            )
        )
    end
end


text:
export LEFT, RIGHT, move, add_tile, LEFT_REWARD, RIGHT_REWARD, show, initbboard, display, move_reward

import Base.show, Base.display

# A gameboard is stored on a UInt (64bit)
# where each cell is
const ROWMASK=UInt16(2^16-1)
const CELLMASK=UInt16(2^4-1)

const LEFT, LEFT_REWARD, RIGHT, RIGHT_REWARD = make_row_lookup()

struct Bitboard
    board::UInt64
end

function initbboard()
    zero(UInt64) |> Bitboard |> add_tile |> add_tile
end

"""
    board2cols(board::UInt64)::Vector{Vector{UInt16}}

Return the columns from a board
"""
function board2cols(bitboard::Bitboard)::Vector{UInt16}
    board = bitboard.board
    # initialise the columns
    cols = zeros(UInt16, 4)

    # take one row at a time
    for row_shift in 48:-16:0
        row = (board >> row_shift) & ROWMASK

        # populate the right cell for each column
        for i in 4:-1:1
            cols[i] <<= 4 # doing this here effectively means each col only gets shifted 3 times
            cols[i] |= row & CELLMASK
            row >>= 4
        end
    end
    cols
end

"""
    move(board::UInt64, dir::Dirs)::UInt64
    move(board::UInt64, LOOKUP::Vector{UInt16}::UInt64

Bitboard game
"""
function move(bitboard::Bitboard, LOOKUP::Vector{UInt16})::Bitboard
    board = bitboard.board

    idx = board >> 48
    new_board = zero(UInt64)
    new_board |= LOOKUP[idx+1] #can I add reward lookup to this so I get the new board but also the reward?

    for i in 32:-16:0
        new_board <<= 16
        idx = (board >> i) & ROWMASK
        new_board |= LOOKUP[idx+1]
    end
    Bitboard(new_board)
end

function move(bitboard::Bitboard, LOOKUP::Vector{UInt16}, LOOKUP_REWARD::Vector{Int32})::Tuple{Bitboard, Int32}
    board = bitboard.board

    idx = board >> 48
    reward = LOOKUP_REWARD[idx + 1]
    new_board = zero(UInt64)
    new_board |= LOOKUP[idx+1] #can I add reward lookup to this so I get the new board but also the reward?

    for i in 32:-16:0
        new_board <<= 16
        idx = (board >> i) & ROWMASK
        new_board |= LOOKUP[idx+1]
        reward += LOOKUP_REWARD[idx + 1]
    end
    (Bitboard(new_board), reward)
end

function move_updown(bitboard::Bitboard, LOOKUP::Vector{UInt16})::Bitboard
    cols = board2cols(bitboard)
    cols_moved = getindex.(Ref(LOOKUP), cols .+ 1)

    new_board = zero(UInt64)

    # construct the new board by putting the columns in to rows
    for rowid in 1:4
        for colid in 1:4
            new_board <<= 4
            new_board |= (cols_moved[colid] >> 4(4-rowid)) & CELLMASK
        end
    end

    Bitboard(new_board)
end

function move_updown(bitboard::Bitboard, LOOKUP::Vector{UInt16}, LOOKUP_REWARD::Vector{Int32})::Tuple{Bitboard, Int32}
    cols = board2cols(bitboard)
    cols_moved = getindex.(Ref(LOOKUP), cols .+ 1)

    new_board = zero(UInt64)

    # construct the new board by putting the columns in to rows
    for rowid in 1:4
        for colid in 1:4
            new_board <<= 4
            new_board |= (cols_moved[colid] >> 4(4-rowid)) & CELLMASK
        end
    end

    (Bitboard(new_board), LOOKUP_REWARD[idx+1])
end

function move(bitboard::Bitboard, dir::Dirs)::Bitboard
    if dir == left
        new_bitboard = move(bitboard, LEFT)
    elseif dir == right
        new_bitboard = move(bitboard, RIGHT)
    elseif dir == up
        new_bitboard = move_updown(bitboard, LEFT)
    else
        new_bitboard = move_updown(bitboard, RIGHT)
    end

    new_bitboard
end

function move_reward(bitboard::Bitboard, dir::Dirs)::Tuple{Bitboard, Int32}
    if dir == left
        move(bitboard, LEFT, LEFT_REWARD)
    elseif dir == right
        move(bitboard, RIGHT, RIGHT_REWARD)
    elseif dir == up
        move_updown(bitboard, LEFT, LEFT_REWARD)
    else
        move_updown(bitboard, RIGHT, RIGHT_REWARD)
    end
end

function count0(bitboard::Bitboard)
    board = bitboard.board
    # firstly count how many empty spots there are
    cnt_empty = 0

    for shift in 0:4:60
        cnt_empty += ((board >> shift) & CELLMASK) == 0
    end

    cnt_empty
end

"""this function should randomly add a one or a 2"""
function add_tile(bitboard::Bitboard)::Bitboard
    cnt_empty = count0(bitboard)
    if cnt_empty == 0
        return bitboard
    end

    board = bitboard.board
    for shift in 0:4:60
        if ((board >> shift) & CELLMASK) == 0
            if rand() <= 1/cnt_empty
                board += rand() < 0.1 ? 2*2^shift : 2^shift
                return Bitboard(board)
            end
            cnt_empty -= 1
        end
    end
    Bitboard(board)
end

function add_tile(bitboard::Bitboard, selected::Integer, two_or_four)
    board = bitboard.board
    i = 0
    for shift in 0:4:60
        if ((board >> shift) & CELLMASK) == 0
            i += 1
            if i == selected
                board += two_or_four*2^shift
                return Bitboard(board)
            end
        end
    end
    Bitboard(board)
end

function bitboard_to_array(bitboard::Bitboard)::Array{Int8, 2}
    board = bitboard.board

    outboard = Array{Int8, 2}(undef, 4, 4)

    rowid = 1
    # take one row at a time
    for row_shift in 48:-16:0
        row = (board >> row_shift) & ROWMASK

        # populate the right cell for each column
        for colid in 4:-1:1
            outboard[rowid, colid] = row & CELLMASK
            row >>= 4
        end
        rowid += 1
    end
    outboard
end

function Base.show(io::IO, bitboard::Bitboard)
    show(io, bitboard_to_array(bitboard))
end

function Base.display(bitboard::Bitboard)
    display(bitboard_to_array(bitboard))
end


text:
using GLWindow, ModernGL, GLAbstraction, Romeo, Color, Images, ImmutableArrays

begin
local const BOUNDINGBOX_FRAMEBUFFER = glGenFramebuffers()
glBindFramebuffer(GL_FRAMEBUFFER, BOUNDINGBOX_FRAMEBUFFER)
framebuffsize = [1,1]

minbuffer = Texture(Vec4, framebuffsize)
maxbuffer = Texture(Vec4, framebuffsize)
glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_TEXTURE_2D, minbuffer.id, 0)
glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT1, GL_TEXTURE_2D, maxbuffer.id, 0)
println(GLENUM(glCheckFramebufferStatus(GL_FRAMEBUFFER)).name)
glBindFramebuffer(GL_FRAMEBUFFER, Romeo.RENDER_FRAMEBUFFER)

function boundingboxshader(name::String; view::Dict{ASCIIString, ASCIIString} = Dict{ASCIIString, ASCIIString}(), attributes::Dict{Symbol, Any}=Dict{Symbol, Any}())
    TemplateProgram(
    joinpath(shaderdir, "boundingbox", name), joinpath(shaderdir, "boundingbox", "boundingbox.frag"),
    view=view, attributes=attributes, fragdatalocation=[(0, "minbuffer"),(1, "maxbuffer")]
  )
end
function boundingboxshader(vert::String, frag::String; view::Dict{ASCIIString, ASCIIString} = Dict{ASCIIString, ASCIIString}(), attributes::Dict{Symbol, Any}=Dict{Symbol, Any}())
    TemplateProgram(
    joinpath(shaderdir, "boundingbox", vert), joinpath(shaderdir, "boundingbox", frag),
    view=view, attributes=attributes, fragdatalocation=[(0, "minbuffer"),(1, "maxbuffer")]
  )
end

function boundingbox(renderobject)
    renderobject.boundingbox(renderobject)
end

const FLOAT32MIN = -999999999999999999999999f0
function boundingbox_gpu(renderobject::RenderObject)
	glClearColor(FLOAT32MIN, FLOAT32MIN, FLOAT32MIN, FLOAT32MIN)
    glBindFramebuffer(GL_FRAMEBUFFER, BOUNDINGBOX_FRAMEBUFFER)
    glViewport(0,0, framebuffsize...)
    glClampColor(GL_CLAMP_VERTEX_COLOR, GL_FALSE)
    glClampColor(GL_CLAMP_READ_COLOR, GL_FALSE)
    glClampColor(GL_CLAMP_FRAGMENT_COLOR, GL_FALSE)
    glDrawBuffers(2, [GL_COLOR_ATTACHMENT0, GL_COLOR_ATTACHMENT1])
    glClear(GL_COLOR_BUFFER_BIT)

    glDisable(GL_DEPTH_TEST)
    glDisable(GL_CULL_FACE)
    glDisable(GL_ALPHA_TEST)
    glEnable(GL_BLEND)
    glBlendFunc(GL_ONE, GL_ONE)
    glBlendEquation(GL_MAX)

    program = renderobject.vertexarraybb.program
    glUseProgram(program.id)
    for (key,value) in program.uniformloc
        gluniform(value..., renderobject.uniforms[key])
    end
    if haskey(renderobject.postrenderfunctions, render)
        render(renderobject.vertexarraybb, GL_POINTS)
    elseif haskey(renderobject.postrenderfunctions, renderinstanced)
        renderinstanced(renderobject.vertexarraybb, renderobject[:postrender, renderinstanced][2], GL_POINTS)
    else
        error("RenderObject doesn't have a renderfunction. RenderObject: \n", renderobject)
    end

    println("min: ", -data(minbuffer))
    println("max: ", data(maxbuffer))
    glBindFramebuffer(GL_FRAMEBUFFER, Romeo.RENDER_FRAMEBUFFER)
    glClampColor(GL_CLAMP_VERTEX_COLOR, GL_TRUE)
    glClampColor(GL_CLAMP_READ_COLOR, GL_TRUE)
    glClampColor(GL_CLAMP_FRAGMENT_COLOR, GL_TRUE)
end
end

const testdata = [rgba(0f0,1f0,0f0,1f0) for i=1:9, j=1:11]
obj = visualize(rgba(0,0,0,1))
obj1 = visualize("alksjdlaskjdlkasjd\n"^15)

println(boundingbox(obj))
println(boundingbox(obj1))

GLFW.Terminate()


text:
import Plots
using LinearAlgebra
using DiffEqBase
using OrdinaryDiffEq: SplitODEProblem, solve, IMEXEuler
import SciMLBase

a,b, n = 0, 1, 10               # zmin, zmax, number of cells
n̂_min, n̂_max = -1, 1            # Outward facing unit vectors
α = 100;                        # thermal diffusivity, larger means more stiff
β, γ = 10000, π;                # source term coefficients
Δt = 1000;                      # timestep size
N_t = 10;                       # number of timesteps to take
FT = Float64;                   # float type
Δz = FT(b-a)/FT(n)
Δz² = Δz^2;
∇²_op = [1/Δz², -2/Δz², 1/Δz²]; # interior Laplacian operator
∇T_bottom = 10;                 # Temperature gradient at the top
T_top = 1;                      # Temperature at the bottom
S(z) = β*sin(γ*z)               # source term, (sin for easy integration)
zf = range(a, b, length=n+1);   # coordinates on cell faces

function T_analytic(z) # Analytic steady state solution
    c1 = ∇T_bottom-β*cos(γ*a)/(γ*α)
    c2 = T_top-(β*sin(γ*b)/(γ^2*α)+c1*b)
    return β*sin(γ*z)/(γ^2*α)+c1*z+c2
end

# Initialize interior and boundary stencils:
∇² = Tridiagonal(
    ones(FT, n) .* ∇²_op[1],
    ones(FT, n+1)   .* ∇²_op[2],
    ones(FT, n) .* ∇²_op[3]
);

# Modify boundary stencil to account for BCs

∇².d[1] = -2/Δz²
∇².du[1] = +2/Δz²

# Modify boundary stencil to account for BCs
∇².du[n] = 0  # modified stencil
∇².d[n+1] = 0 # to ensure `∂_t T = 0` at `z=zmax`
∇².dl[n] = 0  # to ensure `∂_t T = 0` at `z=zmax`
D = α .* ∇²

AT_b = zeros(FT, n+1);
AT_b[1] = α*2/Δz*∇T_bottom*n̂_min;
AT_b[end-1] = α*T_top/Δz²;

T = zeros(FT, n+1);
T .= 1;
T[n+1] = T_top; # set top BC

function rhs!(dT, T, params, t)
    n = params.n
    i = 1:n # interior domain
    dT[i] .= S.(zf[i]) .+ AT_b[i]
    return dT
end;

params = (;n)

tspan = (FT(0), N_t*FT(Δt))

prob = SplitODEProblem(
    SciMLBase.DiffEqArrayOperator(
        D,
    ),
    rhs!,
    T,
    tspan,
    params
)
alg = IMEXEuler(linsolve=LinSolveFactorize(lu!))
println("Solving...")
sol = solve(
    prob,
    alg,
    dt = Δt,
    saveat = range(FT(0), N_t*FT(Δt), length=5),
    progress = true,
    progress_message = (dt, u, p, t) -> t,
);

T_end = sol.u[end]

p1 = Plots.plot(zf, T_analytic.(zf), label="analytic", markershape=:circle, markersize=6)
p1 = Plots.plot!(p1, zf, T_end, label="numerical", markershape=:diamond)
p1 = Plots.plot!(p1, title="T ∈ cell faces")

p2 = Plots.plot(zf, abs.(T_end .- T_analytic.(zf)), label="error", markershape=:circle, markersize=6)
p2 = Plots.plot!(p2, title="T ∈ cell faces")

Plots.plot(p1, p2)

max_stars_repo_name: SciML/SciMLTutorialsOutput

max_issues_repo_name: SciML/SciMLTutorialsOutput

max_forks_repo_name: SciML/SciMLTutorialsOutput


text:
# coding:utf-8

doc=
"""
  Dialogue Act Classification Based on Memory Networks
  修正
  2017/02/05 ~
"""

include("switch_depwarn.jl")
switch_depwarn!(false)

using DataFrames
using DataStructures
using PyCall
using ScikitLearn

# python module
unshift!(PyVector(pyimport("sys")["path"]), "")
@pyimport prep
@pyimport memnet

# sklearn module
@sk_import metrics: classification_report

# 単語->単語インデックス
function convert_sequence(word_sequence, word_d, n_count)
  rank =  Dict(w[1] => i for (i, w) in enumerate(sort(collect(word_d), by=x->x[2], rev = true)))
  for (i, seq) in enumerate(word_sequence)
      for (j,w) in enumerate(seq)
        if word_d[w] > n_count
          word_sequence[i][j] = rank[w]
        else
          word_sequence[i][j] = 0
        end
      end
    end
  return word_sequence
end

# 学習用データ作成
function make_feature(d, session_ids)
  stories, questions = [], []
  da_function, da_domain = [], []
  @printf "push ...\n"
  for sid in session_ids
    dialogue_history = []
    for i in collect(1:length(d[sid][:utterance]))
      if i == 1
        dialogue_history = vcat(dialogue_history, ["<S>", "</S>"])
      else
        dialogue_history = vcat(dialogue_history, d[sid][:unigram][i-1])
      end
      push!(stories, dialogue_history)
      push!(questions, d[sid][:unigram][i])
      push!(da_function, d[sid][:function_da_label][i][1])
      push!(da_domain, d[sid][:domain_da_label][i][1])
    end
  end
  @printf "to sequence ...\n"
  stories = convert_sequence(stories, d[:unigram], 0)
  questions = convert_sequence(questions, d[:unigram], 0)
  @printf "pad sequence ...\n"
  story_maxlen = 200 #maximum([length(x) for x in stories])
  query_maxlen = 200 #maximum([length(x) for x in questions])
  stories = prep.pad_sequence(stories, story_maxlen)
  questions = prep.pad_sequence(questions, query_maxlen)
  return Dict(:stories => stories, :questions =>  questions, :da_domain => da_domain, :da_function => da_function)
end

function main(tgt)
  @printf "read corpus ...\n"
  d = open(deserialize, "./data/annotated_corpus.dic")
  session_ids = filter(x -> !(in(x, [:unigram, :bigram])), collect(keys(d)))
  @printf "make feature ...\n"
  feature =  make_feature(d, session_ids)
  @printf "dummy variable...\n"
  feature[tgt], ans_label = prep.get_dummies(feature[tgt])
  @printf "10 fold cross validation ...\n"
  data = prep.cross_validation(collect(1:length(feature[tgt][:, 1])), 10)
  for (i,(train_id, test_id)) in enumerate(zip(data[:,1], data[:, 2]))
    println("---------- $i ----------")
    @printf "split feature ...\n"
    train_feature, test_feature = Dict(), Dict()
    for k in keys(feature)
      train_feature[k]  = feature[k][train_id, :]
      test_feature[k] =  feature[k][test_id, :]
    end
    @printf "read model ...\n"
    vocab_size = length(collect(keys(d[:unigram])))+1
    story_maxlen, query_maxlen = length(feature[:stories][1, :]), length(feature[:questions][1,:])
    answer_dim, embedding_dim = length(feature[tgt][1, :]), 128
    batch_size, nb_epoch = 64, 15 # 256 30 64 15
    @show vocab_size; @show story_maxlen; @show query_maxlen; @show answer_dim
    @show batch_size; @show nb_epoch
    model = memnet.memnet(vocab_size, story_maxlen, query_maxlen, answer_dim, embedding_dim)
    @printf "model fit ...\n"
    model = memnet.best_fit(model, batch_size, nb_epoch, train_feature, test_feature, tgt)
    @printf "predict ...\n"
    memnet.predict(model, test_feature, tgt)
  end
end

main(:da_domain)


text:
export nonmissing, nonmissings

## Utilities for Missing values ##

"""
    nonmissing(vec::T) where T<:AbstractArray
Convenience function to count the number of non-`missing` values
in a vector.
"""
@inline function nonmissing(vec::T) where T<:AbstractArray
    mapreduce(!ismissing, +, vec)
end

"""
    nonmissing(data::PopData, locus::String)
Convenience function to count the number of non-`missing` samples
at a locus.
"""
@inline function nonmissing(data::PopData, locus::String)
    data.genodata[data.genodata.locus .== locus, :genotype] |> nonmissing
end

"""
    nonmissings(vec1::AbstractVector, vec2::AbstractVector)
Return a vector of indices where neither input vectors have a `missing` value, i.e. an
intersection of the indices of their non-missing elements.
"""
@inline function nonmissings(vec1::T, vec2::T) where T <: AbstractVector
    mapreduce(i -> findall(!ismissing, i), intersect, (vec1, vec2))
end


text:
@reexport module Systems

using Distributions
using StatsBase
using DynamicalSystems

include("discretemaps/ar1.jl")
include("discretemaps/ar1_bidir.jl")
include("discretemaps/anishchenko1.jl")
include("discretemaps/henon2.jl")
include("discretemaps/henontriple.jl")
include("discretemaps/ikeda.jl")
include("discretemaps/linearmap.jl")
include("discretemaps/nonlinear3D_linear_and_nonlinear_coupling.jl")
include("discretemaps/nontrivial_pegiun.jl")
include("discretemaps/logistic2_unidir.jl")
include("discretemaps/logistic2_bidir.jl")
include("discretemaps/logistic3.jl")
include("discretemaps/logistic4.jl")
include("discretemaps/var1.jl")
include("discretemaps/verdes.jl")

include("continuous_systems/chuacircuits_driven.jl")
include("continuous_systems/chuacircuit_nscroll_sine.jl")
include("continuous_systems/hindmarsh_rose.jl")
include("continuous_systems/mediated_link.jl")

include("continuous_systems/lorenz_lorenz_bidir.jl")
include("continuous_systems/lorenz_lorenz_lorenz_bidir_forced.jl")
include("continuous_systems/lorenz_lorenz_lorenz_transitive.jl")

include("continuous_systems/rossler_rossler_bidir.jl")
include("continuous_systems/rossler_rossler_rossler_bidir_forced.jl")

include("continuous_systems/rossler_lorenz.jl")

include("noise.jl")
# We export both the equations of motion and function to generate
# the systems for every example. The equations of motion functions
# are always prepended with `eom_`.
# Some systems (involving multiple lags) do not return instances
# of DiscreteDynamicalSystem, but a Dataset with the finished
# iterated map.

export
###############
# Discrete maps
###############
eom_ar1_unidir, ar1_unidir,
eom_anishchenko1, anishchenko1,
eom_henon2, henon2,
eom_henon_triple, henon_triple,
eom_linearmap1, linearmap1,
eom_nonlinear3d, nonlinear3d,
eom_logistic2_unidir, logistic2_unidir,
eom_logistic2_bidir, logistic2_bidir,
eom_logistic3, logistic3,
eom_logistic4, logistic4,
eom_nontrivial_pegiun, nontrivial_pegiun,
eom_var1, var1,
eom_ar1_bidir, ar1_bidir,
eom_verdes, verdes

export
####################
# Continuous systems
###################
eom_chuacircuit_nscroll_sine, chuacircuit_nscroll_sine,
eom_chuacircuits_driven, chuacircuits_driven,
eom_hindmarsh_rose, hindmarsh_rose,
eom_lorenz_lorenz_lorenz_transitive, lorenz_lorenz_lorenz_transitive,
eom_mediated_link, mediated_link,
eom_rossler_rossler_bidir, rossler_rossler_bidir,
eom_rossler_rossler_rossler_bidir_forced, rossler_rossler_rossler_bidir_forced,
eom_rossler_lorenz, rossler_lorenz,
eom_lorenz_lorenz_bidir, lorenz_lorenz_bidir,
eom_lorenz_lorenz_lorenz_bidir_forced, lorenz_lorenz_lorenz_bidir_forced

# Initialise all the systems once, generating a trajectory.
#trajectory(anishchenko1(), 10)
#trajectory(henon2(), 10)
#trajectory(henon4(), 10)
#trajectory(logistic2(), 10)
#trajectory(logistic3(), 10)
#trajectory(logistic4(), 10)
#trajectory(var1(), 10)
#trajectory(ar1_bidirA(), 10)
#trajectory(verdes(), 10)

#trajectory(rossler_rossler(), 10)
#trajectory(rossler_lorenz(), 10)
#trajectory(lorenz_triple(), 10)
#trajectory(mediated_link(), 10)
#trajectory(chuacircuits_driven(), 10)
#trajectory(chuacircuit_nscroll_sine(), 10)

end


text:
using Documenter, ClosedLoopReachability

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

# generate models
include("generate.jl")

makedocs(
    sitename = "ClosedLoopReachability.jl",
    modules = [ClosedLoopReachability],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true",
                             collapselevel = 1,
                             assets = ["assets/aligned.css"]),
    pages = [
        "Home" => "index.md",
        "Examples" => Any["Adaptive Cruise Controller (ACC)" => "models/ACC.md",
                          "Translational Oscillations (TORA) with ReLU controller" => "models/Sherlock-Benchmark-9-TORA.md",
                          "Translational Oscillations (TORA) with sigmoid controller" => "models/Sherlock-Benchmark-9-TORA-Sigmoid.md",
                          "Translational Oscillations (TORA) with ReLU/tanh controller" => "models/Sherlock-Benchmark-9-TORA-ReluTanh.md",
                          "Unicycle Car Model" => "models/Sherlock-Benchmark-10-Unicycle.md",
                          "Vertical Collision Avoidance (VCAS)" => "models/VertCAS.md",
                          "Single Inverted Pendulum" => "models/Single-Pendulum.md",
                          "Double Inverted Pendulum" => "models/Double-Pendulum.md",
                          "Airplane" => "models/Airplane.md"],
                          #"Sliding controller" => "models/Sherlock-Benchmark-7.md",
                          #"Nonlinear Cart-Pole" => "models/Cart-Pole.md"],
        "API Reference" => Any["Problem types"=>"lib/problems.md",
                               "Solvers"=>"lib/solvers.md",
                               "Utilities"=>"lib/utils.md"]
    ],
    strict = false
)

deploydocs(
    repo = "github.com/JuliaReach/ClosedLoopReachability.jl.git",
    push_preview=true
)


text:
@doc raw"""
    NonlocalMean(λ [, r_p=2, r_s=2r_p+1])
    NonlocalMean(λ, img)

Non-local mean denoising filter for additive white Gaussian noise using the
following formula:

```math
    y[p] = \sum_{q \in \mathscr{N}(r\_s, p)} w[p, q]x[q]
```

where `w[p, q]` represents the similarity of two patches centered in `p` and
`q`, calculated by

```math
    w[p, q] = \frac{1}{Z[p]}e^{-\frac{\lVert\mathscr{N}(r\_w, p) - \mathscr{N}(r\_w, q) \rVert^2_{2, a}}{λ^2}}
```

`Z[p]` is the normalizing constant so that `∑w` over `q` equals to `1`.

# Arguments

* `λ::Float32` is the degree of filtering, larger `λ` produces a smoother result.

* `r_p::Int` is the radius of image patch size. By default it's 2.

* `r_s::Int` is the radius of search window, large `r_s` would slow down
the filtering significantly. By default it's `2r_p + 1`.

* If you pass `img` to `NonlocalMean`, it will be used to estimate `r_p` and `r_s`.

# Examples

```julia
img = testimage("lena_color_256")

n = AdditiveWhiteGaussianNoise(0.1)
noisy_img = apply_noise(img, n)

# use default filter arguments
f_denoise = NonlocalMean(0.1)
denoised_img = reduce_noise(noisy_img, f_denoise)

# estimate filter arguments with noisy_img
f_denoise = NonlocalMean(0.1, noisy_img)
denoised_img = reduce_noise(noisy_img, f_denoise)
```

See also: [`reduce_noise`](@ref), [`reduce_noise!`](@ref)

# References

[1] Buades, A., Coll, B., & Morel, J. M. (2005, June). A non-local algorithm for image denoising. In _2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR'05)_ (Vol. 2, pp. 60-65). IEEE.
"""
struct NonlocalMean <: AbstractImageDenoiseAlgorithm
    """degree of filtering"""
    λ::Float32
    """radius of image patch"""
    r_p::Int
    """radius of search window"""
    r_s::Int
    function NonlocalMean(λ, r_p, r_s)
        r_p >= 1 || throw(ArgumentError("radius of image patch r_p should >= 1, instead it's $(r_p)"))
        r_s >= 1 || throw(ArgumentError("radius of image patch r_p should >= 1, instead it's $(r_p)"))
        λ > 0 || @warn "λ is supposed to be positive"
        new(λ, r_p, r_s)
    end
end
NonlocalMean(λ, r_p=2) = NonlocalMean(λ, r_p, 2r_p+1)
NonlocalMean(λ, img::GenericImage) = NonlocalMean(λ, max(round.(size(img)./128)..., 1))

function (f::NonlocalMean)(out::AbstractArray{<:NumberLike, 2},
                           img::AbstractArray{<:NumberLike, 2})
    axes(out) == axes(img) || ArgumentError("Images should have the same axes.")
    r_p, r_s, λ = f.r_p, f.r_s, f.λ^2

    T = floattype(eltype(img))
    kernel = make_kernel(T, f.r_p)
    img = of_eltype(T, img)

    R = CartesianIndices(img) # indices without padding

    oₚ = ntuple(_->r_p, ndims(img))
    padded_axes = map(axes(img), oₚ) do ax, o
        first(ax)-o:last(ax)+o
    end
    img = PaddedView(zero(eltype(img)), img, padded_axes)
    # PaddedView's `getindex` has some overhead so we collect it into a dense array
    # and then pass to OffsetArray. Becuase we're doing a heavy amount of `getindex`
    # operation, this extra memory allocation worths.
    img = OffsetArray(collect(img), padded_axes)
    Δₚ = CartesianIndex(oₚ)
    Δₛ = CartesianIndex(ntuple(_->r_s, ndims(img)))

    patch_p = zeros(T, size(kernel))
    for p in R
        # patch_p will be indexed for many times, thus preallocating it into a contiguous memeory layout
        # helps improve the performance
        patch_p .= @view img[_colon(p-Δₚ, p+Δₚ)]

        wmax = zero(T) # set w[p, p] as wmax instead of 1
        ∑wq  = zero(T)
        ∑w   = zero(T) # Z[p] in the docstring
        @inbounds @simd for q in _colon(max(first(R), p-Δₛ), min(p+Δₛ, last(R)))
            if p != q
                patch_q = @view img[_colon(q-Δₚ, q+Δₚ)]

                # calculate weight
                w = T(exp(-wsqeucliean(kernel, patch_p, patch_q)/λ))

                # weighted sum over q
                ∑wq += w*img[q]
                ∑w  += w;

                w > wmax && (wmax = w)
            end
        end
        # add w[p, p] in the end
        ∑wq += wmax*img[p]
        ∑w  += wmax

        out[p] = ∑wq / ∑w
    end
    return out
end

function (f::NonlocalMean)(out::AbstractArray{T, 2},
                           img::AbstractArray{T, 2}) where T<: AbstractRGB
    cv_out = channelview(out)
    cv_img = channelview(img)
    for i in 1:size(cv_img, 1)
        f(view(cv_out, i, :, :), view(cv_img, i, :, :))
    end
    return out
end

function (f::NonlocalMean)(out::AbstractArray{T, 2},
                  img::AbstractArray{T, 2}) where T<:Color3
    img = RGB.(img)
    tmp = similar(img)
    out .= f(tmp, img)
    # f(of_eltype(RGB, out), of_eltype(RGB, img))
end

""" gaussian-like kernel """
function make_kernel(T, r)
    kernel = OffsetArrays.centered(zeros(2r+1, 2r+1))
    R = CartesianIndices(kernel)
    for d in 1:r
        v = 1/(2d+1)^2
        for i in R[-d:d, -d:d]
            kernel[i] += v
        end
    end
    OffsetArrays.centered(T.(kernel ./ sum(kernel)))
end

""" weighted squared euclidean """
function wsqeucliean(W, X, Y)
    rst = zero(eltype(W))
    # use linear indexing
    @inbounds @simd for i = 1:length(W)
        rst += W[i] * _abs2(X[i] - Y[i])
    end
    return rst
end


text:
using ITensors
using Random

include(joinpath(ITensors.examples_dir(), "src", "electronk.jl"))
include(joinpath(ITensors.examples_dir(), "src", "hubbard.jl"))

function main(; Nx::Int = 6, Ny::Int = 3, U::Float64 = 4.0,
                t::Float64 = 1.0, maxdim::Int = 3000,
                conserve_ky = true, use_splitblocks = true,
                seed = 1234)
  Random.seed!(seed)
  @show Threads.nthreads()
  @show ITensors.using_threaded_blocksparse()

  N = Nx * Ny

  sweeps = Sweeps(10)
  maxdims = min.([100, 200, 400, 800, 2000, 3000, maxdim], maxdim)
  setmaxdim!(sweeps, maxdims...)
  setcutoff!(sweeps, 1e-6)
  setnoise!(sweeps, 1e-6, 1e-7, 1e-8, 0.0)
  @show sweeps

  sites = siteinds("ElecK", N; conserve_qns = true,
                   conserve_ky = conserve_ky, modulus_ky = Ny)

  ampo = hubbard(Nx = Nx, Ny = Ny, t = t, U = U, ky = true)
  H = MPO(ampo, sites)

  # This step makes the MPO more sparse.
  # It generally improves DMRG performance
  # at large bond dimensions but makes DMRG slower at
  # small bond dimensions.
  if use_splitblocks
    H = splitblocks(linkinds, H)
  end

  # Number of structural nonzero elements in a bulk
  # Hamiltonian MPO tensor
  @show nnz(H[end÷2])
  @show nnzblocks(H[end÷2])

  # Create start state
  state = Vector{String}(undef, N)
  for i in 1:N
    x = (i - 1) ÷ Ny
    y = (i - 1) % Ny
    if x % 2 == 0
      if y % 2 == 0
        state[i] = "Up"
      else
        state[i] = "Dn"
      end
    else
      if y % 2 == 0
        state[i] = "Dn"
      else
        state[i] = "Up"
      end
    end
  end

  psi0 = randomMPS(sites, state, 10)

  energy, psi = @time dmrg(H, psi0, sweeps; svd_alg = "divide_and_conquer")
  @show Nx, Ny
  @show t, U
  @show flux(psi)
  @show maxlinkdim(psi)
  @show energy
  return energy, H, psi
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Serd"
version = v"0.30.10"

# Collection of sources required to complete build
sources = [
    ArchiveSource("http://download.drobilla.net/serd-0.30.10.tar.bz2", "affa80deec78921f86335e6fc3f18b80aefecf424f6a5755e9f2fa0eb0710edf")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd serd-0.30.10/
install_license COPYING
./waf configure --prefix=$prefix
./waf
./waf install
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("i686", "linux"; libc = "glibc"),
    Platform("x86_64", "linux"; libc = "glibc"),
    Platform("aarch64", "linux"; libc = "glibc"),
    Platform("armv7l", "linux"; call_abi = "eabihf", libc = "glibc"),
    Platform("powerpc64le", "linux"; libc = "glibc"),
    Platform("i686", "linux"; libc = "musl"),
    Platform("x86_64", "linux"; libc = "musl"),
    Platform("aarch64", "linux"; libc = "musl"),
    Platform("armv7l", "linux"; call_abi = "eabihf", libc = "musl"),
    Platform("x86_64", "macos"; ),
    Platform("x86_64", "freebsd"; )
]

# The products that we will ensure are always built
products = [
    LibraryProduct("libserd-0", :libserd)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: jameskermode/Yggdrasil

max_issues_repo_name: jameskermode/Yggdrasil

max_forks_repo_name: jameskermode/Yggdrasil


text:
struct Path
  folder::String
  file_name::String
end

struct ScalarData
  scalarArr::Vector{Vector{Float64}}
end

struct ScalarNames
  nameArr::Vector{String}
end

struct VectorData
  vectorArr::Vector{Vector{Vector{Float64}}}
end

struct VectorNames
  nameArr::Vector{String}
end

# -----------
# constructors
# -----------

# single entry
ScalarData(u::Vector{Float64}) = ScalarData([u])
ScalarNames(u::String) = ScalarNames([u])
VectorData(u::Vector{Vector{Float64}}) = VectorData([u])
VectorNames(u::String) = VectorNames([u])

# multiple entries
ScalarData(u::Vector{Float64}...) = ScalarData([u[i] for i=1:length(u)])
ScalarNames(u::String...) = ScalarNames([u[i] for i=1:length(u)])
VectorData(u::Vector{Vector{Float64}}...) = VectorData([u[i] for i=1:length(u)])
VectorNames(u::String...) = VectorNames([u[i] for i=1:length(u)])

Path(file_name::String) = Path("solutions",file_name)

##########################
#### vtsave functions ####
##########################

"""
  New scalar solutions
"""
function vtksave(mesh::T,sd::ScalarData,sn::ScalarNames,path::Path;time=0.0) where T<:AbstractMesh
  Nscalars = length(sd.scalarArr)

  # get xy & cm into correct form
  xy = NodesToArray(mesh.xy)
  cm = ElementsToArray(mesh.cm)

  newsol = zeros(Float64,length(mesh.xy),Nscalars)
  for i=1:length(mesh.xy), j=1:Nscalars
    newsol[i,j] = sd.scalarArr[j][i]
  end

  # creates folder if it doesn't already exist
  run(`mkdir -p $(path.folder)`);
  file_name = string(path.folder,"/",path.file_name,".vtk")
  TWOD_UNSTRUCTURED_TO_VTK(xy,cm,newsol,sn.nameArr,file_name;time=time)
end
"""
  New vector solutions
"""
function vtksave(mesh::T,vd::VectorData,vn::VectorNames,path::Path;time=0.0) where T<:AbstractMesh
  Nvectors = length(vd.vectorArr)

  # get xy & cm into correct form
  xy = NodesToArray(mesh.xy)
  cm = ElementsToArray(mesh.cm)

  newsol = zeros(Float64,length(mesh.xy),2,Nvectors)
  for i=1:length(mesh.xy), j=1:2, k=1:Nvectors
    newsol[i,j,k] = vd.vectorArr[k][j][i]
  end

  # creates folder if it doesn't already exist
  run(`mkdir -p $(path.folder)`);
  file_name = string(path.folder,"/",path.file_name,".vtk")
  TWOD_UNSTRUCTURED_TO_VTK(xy,cm,[],[],newsol,vn.nameArr,file_name;time=time)
end
"""
  New scalar + vector solutions -- exports linear element visualization
"""
function vtksave(mesh::T,sd::ScalarData,sn::ScalarNames,vd::VectorData,vn::VectorNames,path::Path;time=0.0) where T<:AbstractMesh
  #if true #mesh is Q2!!!! (or fluid)
  #end
  Nscalars = length(sd.scalarArr)
  Nvectors = length(vd.vectorArr)

  # get xy & cm into correct form
  if typeof(mesh) == FluidMesh
    xy = NodesToArray(mesh.xy)
    cm = ElementsToArray(mesh.cm)
    xyp = NodesToArray(mesh.xyp)
    cmp = ElementsToArray(mesh.cmp)
  else
    xy = NodesToArray(mesh.xy)
    cm = ElementsToArray(mesh.cm)
    xyp = NodesToArray(mesh.xy)
    cmp = ElementsToArray(mesh.cm)
  end

  # if fluid mesh, account for pressure
  if typeof(mesh) == FluidMesh
    newscalar = zeros(Float64,size(xyp,1),Nscalars)
    for j=1:Nscalars
      if (sn.nameArr[j] != "pressure") && (sn.nameArr[j] != "head")
        newscalar[:,j] = lowerScalar(xy,cm,sd.scalarArr[j])
      else
        newscalar[:,j] = sd.scalarArr[j]
      end
    end

    newvector = zeros(Float64,size(xyp,1),2,Nvectors)
    for k=1:Nvectors
      uL,vL = lowerUV(xy,cm,vd.vectorArr[k][1],vd.vectorArr[k][2])
      for ti=1:size(xyp,1)
        newvector[ti,1,k] = uL[ti]
        newvector[ti,2,k] = vL[ti]
      end
    end
  else
    newscalar = zeros(Float64,size(xyp,1),Nscalars)
    for j=1:Nscalars
      newscalar[:,j] = sd.scalarArr[j]
    end

    newvector = zeros(Float64,size(xyp,1),2,Nvectors)
    for k=1:Nvectors
      for ti=1:size(xyp,1)
        newvector[ti,1,k] = vd.vectorArr[k][1][ti]
        newvector[ti,2,k] = vd.vectorArr[k][2][ti]
      end
    end
  end

  # creates folder if it doesn't already exist
  run(`mkdir -p $(path.folder)`);
  file_name = string(path.folder,"/",path.file_name,".vtk")
  TWOD_UNSTRUCTURED_TO_VTK(xyp,cmp,newscalar,sn.nameArr,newvector,vn.nameArr,file_name;time=time)
end

"""
  vtksave for particle output
"""
function vtksave(pList::Vector,radius::Float64,fn::String,timeval::Float64)
  Nparticles = length(pList)
  xy = zeros(Float64,Nparticles,2)
  for i=1:Nparticles
    xy[i,1] = pList[i].x
    xy[i,2] = pList[i].y
  end

  PARTICLES_TO_VTK(xy,radius,fn;time=timeval)
end

"""
  TWOD_UNSTRUCTURED_TO_VTK(xy,cm,fval,fname,vecVal,vecName,file_name[;time=0.0])

  Outputs visualization file in vtk legacy format from FEM data

Arguments:
`xy::Matrix{Float64}`
`cm::Matrix{Float64}`
`fval`
`fname`
`vecVal`
`vecName`
`file_name::String`
`time::Float64`
"""
function TWOD_UNSTRUCTURED_TO_VTK(xy::Matrix{Float64},cm::Matrix{Int},
                                  fval,fname,vecVal,vecName,
                                  file_name::String;time=0.0)
  # prepare
  fvalPresent   = false
  vecValPresent = false
  Npts = size(xy,1)
  Nels = size(cm,1)
  NPE = size(cm,2)      # nodes per element
  if NPE == 4
    NPEcode = 9
  elseif NPE == 9
    NPE = 8
    NPEcode = 23
  else
    error("Invalid cell type. Must be linear or biquadratic quadrilateral")
  end

  # check that fnames is scalar, and turns it into array of scalars
  if typeof(fname) == String
    fname = [fname]
  end

  # check that fvals or vecVals are nonempty
  if size(fval,1) == Npts
    Numfvals    = size(fval,2)
    fvalPresent = true
  end
  if size(vecVal,1) == Npts
    NumVecVals = size(vecVal,3)

    vecValPresent = true
  end

  # open file to write
  f = open(file_name,"w")

  # header
  println(f,"# vtk DataFile Version 3.0")
  println(f,"")
  println(f,"ASCII")

  printBlankLine(f)

  # print domain points
  println(f,"DATASET UNSTRUCTURED_GRID")

  # print time
  println(f,"FIELD FieldData 1")
  println(f,"TIME 1 1 double")
  println(f,time)

  # print points
  println(f,"POINTS ",Npts," float")
  for i=1:Npts
    println(f,xy[i,1]," ",xy[i,2]," ",0.0)
  end

  printBlankLine(f)

  # print cell information
  println(f,"CELLS ",Nels," ",Nels*(NPE+1))
  for i=1:Nels
    print(f,NPE)
    for elNode = 1:NPE
      print(f," ",cm[i,elNode]-1)
    end
    printBlankLine(f)
  end

  printBlankLine(f)

  # print cell type information
  println(f,"CELL_TYPES ",Nels)
  for i=1:Nels
    println(f,NPEcode)
  end

  printBlankLine(f)

  println(f,"POINT_DATA ",Npts)
  # print scalar values
  if fvalPresent
    for j=1:Numfvals
      println(f,"SCALARS ",fname[j]," double")
      println(f,"LOOKUP_TABLE default")
      for i=1:Npts
        println(f,fval[i,j])
      end
    end
  end

  # print vector values -- 2D vectors
  if vecValPresent
    for j=1:NumVecVals
      println(f,"VECTORS ",vecName[j]," double")
      for i=1:Npts
        println(f,vecVal[i,1,j]," ",vecVal[i,2,j]," ",0.0)
      end
    end
  end

  printBlankLine(f)

	# close file
	close(f)
end
function TWOD_UNSTRUCTURED_TO_VTK(xy::Matrix{Float64},cm::Matrix{Int},
                                  fval,fname,file_name::String;time=0.0)
  vName = []; vVal  = []
  return TWOD_UNSTRUCTURED_TO_VTK(xy,cm,fval,fname,vVal,vName,file_name;time=time)
end

# exports surface in 3D, where connectivity matrix has been generated by mesh function
function SURF_UNSTRUCTURED_TO_VTK(xyz::Matrix{Float64},cm::Matrix{Int},
                                  fval,fname,vecVal,vecName,
                                  file_name::String;time=0.0)
  # prepare
  fvalPresent   = false
  vecValPresent = false
  Npts = size(xyz,1)
  Nels = size(cm,1)
  NPE = size(cm,2)      # nodes per element
  if NPE == 4
    NPEcode = 9
  elseif NPE == 9
    NPE = 8
    NPEcode = 23
  else
    error("Invalid cell type. Must be linear or biquadratic quadrilateral")
  end

  # check that fnames is scalar, and turns it into array of scalars
  if typeof(fname) == String
    fname = [fname]
  end

  # check that fvals or vecVals are nonempty
  if size(fval,1) == Npts
    Numfvals    = size(fval,2)
    fvalPresent = true
  end
  if size(vecVal,1) == Npts
    NumVecVals = size(vecVal,3)

    vecValPresent = true
  end

  # open file to write
  f = open(file_name,"w")

  # header
  println(f,"# vtk DataFile Version 3.0")
  println(f,"")
  println(f,"ASCII")

  printBlankLine(f)

  # print domain points
  println(f,"DATASET UNSTRUCTURED_GRID")

  # print time
  println(f,"FIELD FieldData 1")
  println(f,"TIME 1 1 double")
  println(f,time)

  # print points
  println(f,"POINTS ",Npts," float")
  for i=1:Npts
    println(f,xyz[i,1]," ",xyz[i,2]," ",xyz[i,3])
  end

  printBlankLine(f)

  # print cell information
  println(f,"CELLS ",Nels," ",Nels*(NPE+1))
  for i=1:Nels
    print(f,NPE)
    for elNode = 1:NPE
      print(f," ",cm[i,elNode]-1)
    end
    printBlankLine(f)
  end

  printBlankLine(f)

  # print cell type information
  println(f,"CELL_TYPES ",Nels)
  for i=1:Nels
    println(f,NPEcode)
  end

  printBlankLine(f)

  println(f,"POINT_DATA ",Npts)
  # print scalar values
  if fvalPresent
    for j=1:Numfvals
      println(f,"SCALARS ",fname[j]," double")
      println(f,"LOOKUP_TABLE default")
      for i=1:Npts
        println(f,fval[i,j])
      end
    end
  end

  # print vector values -- 2D vectors
  if vecValPresent
    for j=1:NumVecVals
      println(f,"VECTORS ",vecName[j]," double")
      for i=1:Npts
        println(f,vecVal[i,1,j]," ",vecVal[i,2,j]," ",vecVal[i,3,j])
      end
    end
  end

  printBlankLine(f)

  # close file
  close(f)
end
function SURF_UNSTRUCTURED_TO_VTK(xyz::Matrix{Float64},cm::Matrix{Int},
                                  fval,fname,file_name::String;time=0.0)
  vName = []; vVal  = []
  return SURF_UNSTRUCTURED_TO_VTK(xyz,cm,fval,fname,vVal,vName,file_name;time=time)
end
function SURF_UNSTRUCTURED_TO_VTK(xyz::Matrix{Float64},cm::Matrix{Int},file_name::String;time=0.0)
  fName = []; fVal  = []
  vName = []; vVal  = []
  return SURF_UNSTRUCTURED_TO_VTK(xyz,cm,fVal,fName,vVal,vName,file_name;time=time)
end

"""
    PARTICLES_TO_VTK(xy,radius,file_name)

Saves particle data in legacy *.vtk file format

# Arguments
- `xy::Matrix{Float64}`: Nx2 matrix of particle positions
- `radius::Float64`: radius of particles
- `file_name::String`: file name for export
"""
function PARTICLES_TO_VTK(xy::Matrix{Float64},radius::Float64,file_name::String;time=0.0)
  # prepare
  Npts = size(xy,1)

  # open file to write
  fnvtk = string(file_name,".vtk")
  f = open(fnvtk,"w")

  # header
  println(f,"# vtk DataFile Version 3.0")
  printBlankLine(f)
  println(f,"ASCII")
  printBlankLine(f)

  # print domain points
  println(f,"DATASET UNSTRUCTURED_GRID")

  # print time
  println(f,"FIELD FieldData 1")
  println(f,"TIME 1 1 double")
  println(f,time)

  # print points
  println(f,"POINTS ",Npts," float")
  for i=1:Npts
    println(f,xy[i,1]," ",xy[i,2]," ",0.0)
  end

  printBlankLine(f)

  # print cell type information
  println(f,"CELL_TYPES ",Npts)
  for i=1:Npts
    println(f,"1")
  end

  printBlankLine(f)

  println(f,"POINT_DATA ",Npts)
  # print scalar values
  println(f,"SCALARS ","radius"," double")
  println(f,"LOOKUP_TABLE default")
  for i=1:Npts
    println(f,radius)
  end

	# close file
	close(f)
end


text:
export collision, next_collision

#####################################################################################
# Accuracy & Convenience Functions
#####################################################################################
"""
Approximate arccos(1 - x) for x very close to 0.
"""
@inline (acos1mx(x::T)::T) where {T} = sqrt(2x) + sqrt(x)^3/sixsqrt
const sixsqrt = 6sqrt(2)

@inline cross2D(a, b) = a[1]*b[2] - a[2]*b[1]

@inline accuracy(::Type{T}) where {T} = sqrt(eps(T))
@inline accuracy(::Type{BigFloat}) = BigFloat(1e-32)
@inline nocollision(::Type{T}) where {T} = (T(Inf), SV{T}(0.0, 0.0))

#######################################################################################
## Particle
#######################################################################################
"""
    collision(p::AbstractParticle, o::Obstacle) → t, cp
Find the collision (if any) between given particle and obstacle.
Return the time until collision and the estimated collision point `cp`.

Return `Inf, SV(0, 0)` if the collision is not possible *or* if the
collision happens backwards in time.

**It is the duty of `collision` to avoid incorrect collisions when the particle is
on top of the obstacle (or very close).**
"""
function collision(p::Particle{T}, w::Wall{T}) where {T}
    n = normalvec(w, p.pos)
    denom = dot(p.vel, n)
    if denom ≥ 0.0
        return nocollision(T)
    else
        t = dot(w.sp - p.pos, n)/denom
        return t, p.pos + t*p.vel
    end
end

function collision(p::Particle{T}, w::Union{FiniteWall{T}, FiniteSplitterWall{T}}) where {T}
    n = normalvec(w, p.pos)
    denom = dot(p.vel, n)
    # case of velocity pointing away of wall:
    denom ≥ 0.0 && return nocollision(T)
    posdot = dot(w.sp-p.pos, n)
    # Case of particle starting behind finite wall:
    posdot ≥ 0.0 && return nocollision(T)
    colt = posdot/denom
    i = p.pos + colt * p.vel
    dfc = norm(i - w.center)
    if dfc > w.width/2
        return nocollision(T)
    else
        return colt, i
    end
end

function collision(p::Particle{T}, d::Circular{T}) where {T}

    dotp = dot(p.vel, normalvec(d, p.pos))
    dotp ≥ 0.0 && return nocollision(T)

    dc = p.pos - d.c
    B = dot(p.vel, dc)           #pointing towards circle center: B < 0
    C = dot(dc, dc) - d.r*d.r    #being outside of circle: C > 0
    Δ = B*B - C

    Δ ≤ 0.0 && return nocollision(T)
    sqrtD = sqrt(Δ)

    # Closest point:
    t = -B - sqrtD
    return t, p.pos + t*p.vel
end

function collision(p::Particle{T}, d::Antidot{T}) where {T}

    dotp = dot(p.vel, normalvec(d, p.pos))
    if d.pflag == true
        dotp ≥ 0 && return nocollision(T)
    end

    dc = p.pos - d.c
    B = dot(p.vel, dc)           #velocity towards circle center: B < 0
    C = dot(dc, dc) - d.r*d.r    #being outside of circle: C > 0
    Δ = B^2 - C
    Δ ≤ 0 && return nocollision(T)
    sqrtD = sqrt(Δ)

    # Closest point (may be in negative time):
    if dotp < 0
        t = -B - sqrtD
    else
        t = -B + sqrtD
    end
    t ≤ 0.0 ? nocollision(T) : (t, p.pos + t*p.vel)
end

function collision(p::Particle{T}, d::Semicircle{T}) where {T}

    dc = p.pos - d.c
    B = dot(p.vel, dc)         #velocity towards circle center: B > 0
    C = dot(dc, dc) - d.r*d.r    #being outside of circle: C > 0
    Δ = B^2 - C

    Δ ≤ 0 && return nocollision(T)
    sqrtD = sqrt(Δ)

    nn = dot(dc, d.facedir)
    if nn ≥ 0 # I am NOT inside semicircle
        # Return most positive time
        t = -B + sqrtD
    else # I am inside semicircle:
        t = -B - sqrtD
        # these lines make sure that the code works for ANY starting position:
        if t ≤ 0 || distance(p, d) ≤ accuracy(T)
            t = -B + sqrtD
        end
    end
    # This check is necessary to not collide with the non-existing side
    newpos = p.pos + p.vel * t
    if dot(newpos - d.c, d.facedir) ≥ 0 # collision point on BAD HALF;
        return nocollision(T)
    end
    # If collision time is negative, return Inf:
    t ≤ 0.0 ? nocollision(T) : (t, p.pos + t*p.vel)
end

function collision(p::Particle{T}, e::Ellipse{T}) where {T}
    # First check if particle is "looking at" eclipse if it is outside
    if e.pflag
        # These lines may be "not accurate enough" but so far all is good
        dotp = dot(p.vel, normalvec(e, p.pos))
        dotp ≥ 0.0 && return nocollision(T)
    end

    # http://www.ambrsoft.com/TrigoCalc/Circles2/Ellipse/EllipseLine.htm
    a = e.a; b = e.b
    # Translate particle with ellipse center (so that ellipse lies on [0, 0])
    pc = p.pos - e.c
    # Find μ, ψ for line equation y = μx + ψ describing particle
    μ = p.vel[2]/p.vel[1]
    ψ = pc[2] - μ*pc[1]

    # Determinant and intersection points follow from the link
    denomin = a*a*μ*μ + b*b
    Δ² = denomin - ψ*ψ
    Δ² ≤ 0 && return nocollision(T)
    Δ = sqrt(Δ²); f1 = -a*a*μ*ψ; f2 = b*b*ψ # just factors
    I1 = SV(f1 + a*b*Δ, f2 + a*b*μ*Δ)/denomin
    I2 = SV(f1 - a*b*Δ, f2 - a*b*μ*Δ)/denomin

    d1 = norm(pc - I1); d2 = norm(pc - I2)
    if e.pflag
        return d1 < d2 ? (d1, I1 + e.c) : (d2, I2 + e.c)
    else # inside the ellipse: one collision is _always_ valid
        if d1 < d2
            dmin, Imin = d1, I1
            dmax, Imax = d2, I2
        else
            dmin, Imin = d2, I2
            dmax, Imax = d1, I1
        end

        if dmin < accuracy(T) # special case for being very close to ellipse
            dotp = dot(p.vel, normalvec(e, Imin))
            dotp ≥ 0 && return (dmax, Imax + e.c)
        end
         # check which of the two points is ahead or behind the obstacle
        z1 = dot(pc - Imax, p.vel)
        return z1 < 0 ? (dmax, Imax + e.c) : (dmin, Imin + e.c)
    end
end

#######################################################################################
## Magnetic particle
#######################################################################################
function collision(p::MagneticParticle{T}, w::Wall{T}) where {T}
    ω = p.omega
    pc, pr = cyclotron(p)
    P0 = p.pos
    P2P1 = w.ep - w.sp
    P1P3 = w.sp - pc
    # Solve quadratic:
    a = dot(P2P1, P2P1)
    b = 2*dot(P2P1, P1P3)
    c = dot(P1P3, P1P3) - pr*pr
    Δ = b^2 -4*a*c
    # Check if line is completely outside (or tangent) of the circle:
    Δ ≤ 0.0 && return nocollision(T)
    # Intersection coefficients:
    u1 = (-b - sqrt(Δ))/2a
    u2 = (-b + sqrt(Δ))/2a
    cond1 = 0.0 ≤ u1 ≤ 1.0
    cond2 = 0.0 ≤ u2 ≤ 1.0
    # Check if the line (wall) is completely inside the circle:
    θ, I = nocollision(T)
    if cond1 || cond2
        dw = w.ep - w.sp
        for (u, cond) in ((u1, cond1), (u2, cond2))
            Y =  w.sp + u*dw
            if cond
                φ = realangle(p, w, Y)
                φ < θ && (θ = φ; I = Y)
            end
        end
    end
    # Collision time = arc-length until collision point
    return θ*pr, I
end

function collision(p::MagneticParticle{T}, o::Circular{T}) where {T}
    ω = p.omega
    pc, rc = cyclotron(p)
    p1 = o.c
    r1 = o.r
    d = norm(p1-pc)
    if (d >= rc + r1) || (d <= abs(rc-r1))
        return nocollision(T)
    end
    # Solve quadratic:
    a = (rc^2 - r1^2 + d^2)/2d
    h = sqrt(rc^2 - a^2)
    # Collision points (always 2):
    I1 = SVector{2, T}(
        pc[1] + a*(p1[1] - pc[1])/d + h*(p1[2] - pc[2])/d,
        pc[2] + a*(p1[2] - pc[2])/d - h*(p1[1] - pc[1])/d
    )
    I2 = SVector{2, T}(
        pc[1] + a*(p1[1] - pc[1])/d - h*(p1[2] - pc[2])/d,
        pc[2] + a*(p1[2] - pc[2])/d + h*(p1[1] - pc[1])/d
    )
    # Calculate real time until intersection:
    θ1 = realangle(p, o, I1)
    θ2 = realangle(p, o, I2)
    # Collision time, equiv. to arc-length until collision point:
    return θ1 < θ2 ? (θ1*rc, I1) : (θ2*rc, I2)
end

function collision(p::MagneticParticle{T}, o::Semicircle{T}) where {T}
    ω = p.omega
    pc, rc = cyclotron(p)
    p1 = o.c
    r1 = o.r
    d = norm(p1-pc)
    if (d >= rc + r1) || (d <= abs(rc-r1))
        return nocollision(T)
    end
    # Solve quadratic:
    a = (rc^2 - r1^2 + d^2)/2d
    h = sqrt(rc^2 - a^2)
    # Collision points (always 2):
    I1 = SVector{2, T}(
        pc[1] + a*(p1[1] - pc[1])/d + h*(p1[2] - pc[2])/d,
        pc[2] + a*(p1[2] - pc[2])/d - h*(p1[1] - pc[1])/d
    )
    I2 = SVector{2, T}(
        pc[1] + a*(p1[1] - pc[1])/d - h*(p1[2] - pc[2])/d,
        pc[2] + a*(p1[2] - pc[2])/d + h*(p1[1] - pc[1])/d
    )
    # Only consider intersections on the "correct" side of Semicircle:
    cond1 = dot(I1-o.c, o.facedir) < 0
    cond2 = dot(I2-o.c, o.facedir) < 0
    # Collision time, equiv. to arc-length until collision point:
    θ, I = nocollision(T)
    if cond1 || cond2
        for (Y, cond) in ((I1, cond1), (I2, cond2))
            if cond
                φ = realangle(p, o, Y)
                φ < θ && (θ = φ; I = Y)
            end
        end
    end
    # Collision time = arc-length until collision point
    return θ*rc, I
end

collision(p::MagneticParticle, e::Ellipse) = error(
"Magnetic propagation for Ellipse is not supported :( Consider contributing a "*
"method to `collision(p::MagneticParticle, e::Ellipse)`!")

"""
    realangle(p::MagneticParticle, o::Obstacle, I) -> θ
Given the intersection point `I` of the trajectory of a magnetic particle `p` with
some obstacle `o`, find the real angle that will be spanned until the particle
collides with the obstacle.

The function also takes care of problems that may arise when particles are very
close to the obstacle's boundaries, due to floating-point precision.
"""
function realangle(p::MagneticParticle{T}, o::Obstacle{T}, i::SV{T})::T where {T}

    pc = p.center; pr = p.r; ω = p.omega
    P0 = p.pos
    PC = pc - P0
    d2 = dot(i-P0,i-P0) #distance of particle from intersection point
    # Check dot product for close points:
    if d2 ≤ accuracy(T)*accuracy(T)
        dotp = dot(p.vel, normalvec(o,  p.pos))
        # Case where velocity points away from obstacle:
        dotp ≥ 0 && return T(Inf)
    end

    d2r = (d2/(2pr^2))
    d2r > 2 && (d2r = T(2.0))
    # Correct angle value for small arguments (between 0 and π/2):
    θprime = d2r < 1e-3 ? acos1mx(d2r) : acos(1.0 - d2r)

    # Get "side" of i:
    PI = i - P0
    side = cross2D(PI, PC)*ω
    # Get angle until i (positive number between 0 and 2π)
    side < 0 && (θprime = T(2π-θprime))
    return θprime
end

#######################################################################################
## next_collision
#######################################################################################
"""
    next_collision(p::AbstractParticle, bd::Billiard) -> i, tmin, cp
Compute the [`collision`](@ref) across all obstacles in `bd` and find the minimum
one. Return the index of colliding obstacle, the time and the collision point.
"""
function next_collision end

@generated function next_collision(p, bd::Billiard{T, L, BT}) where {T, L, BT}
    out = :(ind = 0; tmin = T(Inf); cp = SV{T}(0.0, 0.0))
    for j=1:L
        push!(out.args, quote
                            let x = bd[$j]
                                tcol, pcol = collision(p, x)
                                # Set minimum time:
                                if tcol < tmin
                                  tmin = tcol
                                  ind = $j
                                  cp = pcol
                                end
                            end
                        end
                        )
    end
    push!(out.args, :(return ind, tmin, cp))
    return out
end


text:
using Speechbox
using GMMbox
#using Distances
#using PyPlot
using Distributions
#using Clustering

ndata = convert(Array{Float32},rand(MvNormal(zeros(13),eye(13)),100));
g = gaussian(ndata);
gm = gmm(4,ndata);

data = rand(gm,5000000);

#=
fpath = "/home/vidhya/work/workspace/swb_mfcc";
mfc_flist = filter(x->endswith(x,".mfc"),readdir(fpath));
flist = mfc_flist;
for i=1:length(mfc_flist)
  flist[i] = joinpath(fpath,mfc_flist[i]);
end
=#

#=
ofile = "/home/vidhya/work/workspace/swb.mfc";
#combine_binary(ofile,flist);
data = read_binary(ofile);
=#

data = convert(Matrix{Float32},data);
g = gaussian(data);
gm = gmm(4,data);

p = zeros(eltype(data),size(data,2));

#=
@time p = logprob(gm,data);
@time p = logprob(gm,data);
@time logprob!(p,gm,data);
=#

@time @sync for i=1:10
  logprob!(p,gm,data);
end