text:
function valid_super_cell(size::Integer=3)
    cell = rand(-3:3, (size, size))
    while det(cell) ≤ 1
        cell = rand(-3:3, (3, 3))
    end
    cell
end

cells = Any[
  [1.5 -0.5 -0.5;  0  1 -0.5; 0.5 -0.5  1]u"nm",
  [1.5 -0.5 -0.5;  0  1 -0.5; 0.5 -0.5  1]
]
@testset "> Potential equivalents regression $(unit(eltype(cell)))" for cell in cells
    cell = [1.5 -0.5 -0.5;  0  1 -0.5; 0.5 -0.5  1]u"nm"
    avecs, bvecs, cvecs = Crystals.SpaceGroup.potential_equivalents(cell)
    norms = mapreducedim(x -> x*x, +, cell, 1)
    @test size(avecs) == (3, 4)
    @test all(mapreducedim(x -> x * x, +, avecs, 1) .== norms[1])

    @test size(bvecs) == (3, 8)
    @test all(mapreducedim(x -> x * x, +, bvecs, 1) .== norms[2])

    @test size(cvecs) == (3, 8)
    @test all(mapreducedim(x -> x * x, +, cvecs, 1) .== norms[3])
end

@testset "> Point group operations" begin
    parameters = [
        ([0 0.5 0.5; 0.5 0 0.5; 0.5 0.5 0]u"nm", 48),
        ([-0.5 0.5 0.5; 0.5 -0.5 0.5; 0.5 0.5 -0.5], 48),
        ([-0.6 0.5 0.5; 0.6 -0.5 0.5; 0.6 0.5 -0.5], 4),
        ([-0.7 0.7 0.7; 0.6 -0.5 0.5; 0.6 0.5 -0.5], 8),
        ([-0.765 0.7 0.7; 0.665 -0.5 0.5; 0.6 0.5 -0.5], 2)
    ]
    allops = point_group([0 0.5 0.5; 0.5 0 0.5; 0.5 0.5 0])

    @testset "Cell: $cell" for (cell, numops) in parameters
        ops = point_group(cell)
        @test length(ops) == numops
        for op in ops
            @test size(op) == (3, 3)
            transformation = inv(cell) * op * cell
            @test all(isinteger, round.(transformation, 8))
            @test volume(transformation) ≈ 1.0

            if numops != 48
                found = 0
                for op in allops
                    transformation = inv(cell) * op * cell
                    all(isinteger, round.(transformation, 8)) || continue
                    isapprox(volume(transformation), 1, rtol=1e-8) || continue
                    found += 1
                end
                @test found == numops
            end
        end
    end
end

@testset "> Inner translations" begin
    @testset ">> Implementation details" begin
        a = DataFrame(species=["Al", "Al", "O", "O", "O"], momo=[1, 2, 3, 3, 1],
                      mama=[:+, :+, :+, :-, :-])

        species = Crystals.SpaceGroup.species_ids(a, [:species])
        @test species[1] == species[2]
        @test species[3] == species[4] == species[5]
        @test species[1] ≠ species[3]

        species = Crystals.SpaceGroup.species_ids(a, [:species, :momo])
        @test length(Set(species[[1, 2, 3, 5]])) == 4
        @test species[3] == species[4]
    end

    diamond = Lattices.diamond()
    diamond[:species] = ["A", "A"]
    @testset ">> No translations" begin
        translations = Crystals.SpaceGroup.inner_translations(diamond)
        @test size(translations, 2) == 0
        @test size(translations, 1) == size(diamond[:position], 1)
    end
    @testset ">> Supercell" begin
        cell = valid_super_cell(3)
        large = supercell(diamond, diamond.cell * cell)
        @test nrow(large) == round(Integer, det(cell)) * nrow(diamond)
        translations = Crystals.SpaceGroup.inner_translations(large)
        @test size(translations, 2) == round(Integer, det(cell) - 1)
        @test size(translations, 1) == size(diamond[:position], 1)
    end
end

@testset "> Make primitive" begin
    zinc_blende = Lattices.zinc_blende()
    @test is_primitive(zinc_blende)
    @test primitive(zinc_blende).cell ≈ zinc_blende.cell
    @test nrow(primitive(zinc_blende)) == nrow(zinc_blende)

    cell = valid_super_cell(3)
    large = supercell(zinc_blende, zinc_blende.cell * cell)
    @test !is_primitive(large)
    prim = primitive(large)
    @test niggly(prim.cell) ≈ niggly(zinc_blende.cell)
    @test nrow(prim) == nrow(zinc_blende)
    for i in 1:2
        prim_pos = prim[i, :position]
        d_pos = zinc_blende[i, :position]
        @test is_periodic(prim_pos, d_pos, prim.cell)
        @test is_periodic(prim_pos, d_pos, zinc_blende.cell)
    end
end

@testset "> Space group operations" begin
    @testset ">> fcc" begin
      fcc = Lattices.fcc()
      ops = space_group(fcc)
      @test length(ops) == 48
      for op in ops
          # No translations
          @test op([0, 0, 0]u"nm") ≈ [0, 0, 0]u"nm"
      end
    end

    @testset ">> b5" begin
        b5 = Lattices.b5()
        ops = space_group(b5)
        @test length(ops) == 48
        @test count(o -> all(abs.(ustrip.(o([0, 0, 0]u"nm"))) .< 1e-12), ops) == 12
    end
end


text:
using CUDA
using Adapt
using DataStructures
using Random

struct CollectPoints{A,B}
    maze_columns::Int64
    maze_rows::Int64
    maze_cell_size::Int32
    agent_radius::Int64
    point_radius::Int32
    agent_movement_range::Float32
    reward_per_collected_positive_point::Float32
    reward_per_collected_negative_point::Float32
    number_time_steps::Int32
    number_inputs::Int64
    number_outputs::Int64
    maze_walls_north::A
    maze_walls_east::A
    maze_walls_south::A
    maze_walls_west::A
    maze_randoms::B
    agents_positions::B
    positive_points_positions::B
    negative_points_positions::B

end

function CollectPoints(configuration::OrderedDict, number_individuals::Int)

    CollectPoints(
        configuration["maze_columns"],
        configuration["maze_rows"],
        convert(Int32, configuration["maze_cell_size"]),
        configuration["agent_radius"],
        convert(Int32, configuration["point_radius"]),
        convert(Float32, configuration["agent_movement_range"]),
        convert(Float32, configuration["reward_per_collected_positive_point"]),
        convert(Float32, configuration["reward_per_collected_negative_point"]),
        convert(Int32, configuration["number_time_steps"]),
        10,
        2,
        CUDA.fill(true, (configuration["maze_columns"], configuration["maze_rows"], number_individuals)),
        CUDA.fill(true, (configuration["maze_columns"], configuration["maze_rows"], number_individuals)),
        CUDA.fill(true, (configuration["maze_columns"], configuration["maze_rows"], number_individuals)),
        CUDA.fill(true, (configuration["maze_columns"], configuration["maze_rows"], number_individuals)),
        CUDA.fill(0, (configuration["maze_columns"] * configuration["maze_rows"], number_individuals)),
        CUDA.fill(0, (2, number_individuals)),
        CUDA.fill(0, (2, number_individuals)),
        CUDA.fill(0, (2, number_individuals)),
    )
end

Adapt.@adapt_structure CollectPoints

function initialize(threadID, blockID, input, environments::CollectPoints, offset, env_seed)

    if threadID == 1
        Random.seed!(Random.default_rng(), env_seed)
        create_maze(threadID, blockID, environments, offset)
    end

    return

end

function get_memory_requirements(environments::CollectPoints)
    total_amount_of_cells = environments.maze_columns * environments.maze_rows

    return sizeof(Int32) * (2 * total_amount_of_cells) + sizeof(Int32) * 8
end

function get_number_inputs(environments::CollectPoints)
    return environments.number_inputs
end

function get_number_outputs(environments::CollectPoints)
    return environments.number_outputs
end

function render(environments::CollectPoints)

    # Draw white background
    background("white")

    # Draw agent
    sethue("green")
    agent_position = Array(environments.agents_positions[:, 1])
    circle(agent_position[1], agent_position[2], environments.agent_radius, :fill)

    # Draw positive point
    sethue("blue")
    positive_points_position = Array(environments.positive_points_positions[:, 1])
    circle(positive_points_position[1], positive_points_position[2], environments.point_radius, :fill)

    # Draw negative point
    sethue("red")
    negative_points_position = Array(environments.negative_points_positions[:, 1])
    circle(negative_points_position[1], negative_points_position[2], environments.point_radius, :fill)

    # Draw maze
    sethue("black")
    setline(2)
    maze_walls_north = Array(environments.maze_walls_north[:, :, 1])
    maze_walls_east = Array(environments.maze_walls_east[:, :, 1])
    maze_walls_south = Array(environments.maze_walls_south[:, :, 1])
    maze_walls_west = Array(environments.maze_walls_west[:, :, 1])

    for cell_x = 1:environments.maze_columns
        for cell_y = 1:environments.maze_rows

            x_left, x_right, y_top, y_bottom = get_coordinates_maze_cell(cell_x, cell_y, environments.maze_cell_size)

            # Draw walls
            if maze_walls_north[cell_x, cell_y] == true
                line(Point(x_left, y_top), Point(x_right, y_top), :stroke)
            end

            if maze_walls_south[cell_x, cell_y] == true
                line(Point(x_left, y_bottom), Point(x_right, y_bottom), :stroke)
            end

            if maze_walls_east[cell_x, cell_y] == true
                line(Point(x_right, y_top), Point(x_right, y_bottom), :stroke)
            end

            if maze_walls_west[cell_x, cell_y] == true
                line(Point(x_left, y_top), Point(x_left, y_bottom), :stroke)
            end

        end
    end
end

function get_coordinates_maze_cell(cell_x::Int, cell_y::Int, maze_cell_size::Int32)

    x_left = maze_cell_size * (cell_x - 1)
    x_right = maze_cell_size * cell_x
    y_top = maze_cell_size * (cell_y - 1)
    y_bottom = maze_cell_size * cell_y

    return x_left, x_right, y_top, y_bottom

end

function create_maze(threadID, blockID, environments::CollectPoints, offset)

    # Total number of cells.
    total_amount_of_cells = environments.maze_columns * environments.maze_rows

    index_cell_stack = 0
    cell_stack = @cuDynamicSharedMem(Int32, (total_amount_of_cells, 2), offset)
    offset += sizeof(cell_stack)

    neighbours = @cuDynamicSharedMem(Int32, (4, 2), offset)

    # Total number of visited cells during maze construction.
    nv = 1

    current_cell_x = 1
    current_cell_y = 1

    while nv < total_amount_of_cells

        number_neighbours = find_valid_neighbors(blockID, current_cell_x, current_cell_y, neighbours, environments)

        if number_neighbours == 0
            # We've reached a dead end: backtrack, pop cell from cell stack
            current_cell_x = cell_stack[index_cell_stack, 1]
            current_cell_y = cell_stack[index_cell_stack, 2]
            index_cell_stack -= 1
            continue
        end

        # Choose a random neighbouring cell and move to it.
        k = rand(1:number_neighbours)
        next_cell_x = neighbours[k, 1]
        next_cell_y = neighbours[k, 2]

        # Store chosen random value to enable unit test for the create_maze function
        environments.maze_randoms[nv, blockID] = k

        # Knock down the wall between current_cell and next_cell
        knock_down_wall(blockID, current_cell_x, current_cell_y, next_cell_x, next_cell_y, environments)

        # Push current cell to cell stack
        index_cell_stack += 1
        cell_stack[index_cell_stack, 1] = current_cell_x
        cell_stack[index_cell_stack, 2] = current_cell_y

        current_cell_x = next_cell_x
        current_cell_y = next_cell_y

        nv += 1
    end

end

# Return a list of unvisited neighbours to cell
function find_valid_neighbors(blockID, cell_x, cell_y, neighbours, environments::CollectPoints)

    number_neighbours = 0

    neighbour_x = 0
    neighbour_y = 0

    for i = 1:4

        # Check West
        if i == 1
            neighbour_x = cell_x - 1
            neighbour_y = cell_y
        end

        # Check East
        if i == 2
            neighbour_x = cell_x + 1
            neighbour_y = cell_y
        end

        # Check South
        if i == 3
            neighbour_x = cell_x
            neighbour_y = cell_y + 1
        end

        # Check North
        if i == 4
            neighbour_x = cell_x
            neighbour_y = cell_y - 1
        end

        if is_valid_neighbor(blockID, neighbour_x, neighbour_y, environments)
            number_neighbours += 1
            neighbours[number_neighbours, 1] = neighbour_x
            neighbours[number_neighbours, 2] = neighbour_y
        end
    end

    return number_neighbours
end

function is_valid_neighbor(blockID, cell_x, cell_y, environments::CollectPoints)

    if (1 ≤ cell_x ≤ environments.maze_columns) && (1 ≤ cell_y ≤ environments.maze_rows)

        if has_all_walls(blockID, cell_x, cell_y, environments)
            return true
        end

    end

    return false
end

function has_all_walls(blockID, cell_x, cell_y, environments::CollectPoints)

    if environments.maze_walls_north[cell_x, cell_y, blockID] == false
        return false
    end

    if environments.maze_walls_south[cell_x, cell_y, blockID] == false
        return false
    end

    if environments.maze_walls_east[cell_x, cell_y, blockID] == false
        return false
    end

    if environments.maze_walls_west[cell_x, cell_y, blockID] == false
        return false
    end

    return true
end

# Knock down the wall between two cells
function knock_down_wall(blockID, cell1_x, cell1_y, cell2_x, cell2_y, environments::CollectPoints)

    # Both cells are in the same row, knock down north and south walls
    if cell1_x == cell2_x

        if cell1_y < cell2_y
            environments.maze_walls_south[cell1_x, cell1_y, blockID] = false
            environments.maze_walls_north[cell2_x, cell2_y, blockID] = false
        else
            environments.maze_walls_north[cell1_x, cell1_y, blockID] = false
            environments.maze_walls_south[cell2_x, cell2_y, blockID] = false
        end
    end

    # Both cells are in the same column, knock down east and west walls
    if cell1_y == cell2_y

        if cell1_x < cell2_x
            environments.maze_walls_east[cell1_x, cell1_y, blockID] = false
            environments.maze_walls_west[cell2_x, cell2_y, blockID] = false
        else
            environments.maze_walls_west[cell1_x, cell1_y, blockID] = false
            environments.maze_walls_east[cell2_x, cell2_y, blockID] = false
        end

    end
end

ext: jl


text:
"""
A `DiscreteDictionary{I,T}` is a dictionary where `I` is a discrete type, e.g.,
the integers. Expansions in this discrete set typically constitute a vector space.

The elements of the set are the Euclidean unit vectors, with a finite index set.
"""
abstract type DiscreteDictionary{I,T} <: Dictionary{I,T}
end

coefficienttype(::Type{DiscreteDictionary{I,T}}) where {I,T} = T

# The point x is in the support of d exactly when it is within the bounds of
# the index set, so we can do a checkbounds with Bool argument (which does not
# throw an error but returns true or false).
dict_in_support(d::DiscreteDictionary{I}, idx, x::I) where I = checkbounds(Bool, d, x)

function dict_in_support(d::DiscreteDictionary{I}, idx, x) where I
    try
        # Attempt to convert x to a native index. If that fails, the result is false.
        idxn = native_index(d, x)
        checkbounds(Bool, d, idxn)
    catch e
        false
    end
end

# Evaluation of discrete sets works as follows:
# -> eval_element: does bounds check on idx
#   -> unsafe_eval_element1: in_support corresponds to bounds check on x
#     -> unsafe_eval_element: return zero(T), or one(T) if idx == x
unsafe_eval_element(set::DiscreteDictionary{I,T}, idx::I, x::I) where {I,T} =
    idx == x ? one(T) : zero(T)

# In the routine above, we made sure that idx and x have the same type, so that
# they can be compared. If they do not have the same type, we can try to
# promote the indices.
unsafe_eval_element(set::DiscreteDictionary, idx, x) =
    unsafe_eval_element(set, native_index(set, idx), native_index(set, x))

isdiscrete(dict::Dictionary) = false
isdiscrete(dict::DiscreteDictionary) = true

#########################
# Concrete discrete sets
#########################

# We have to make these concrete sets because they store the size, and this
# depends on the type of the data. This is so in spite of the fact that type
# parameter I (the index type) already contains information.

"""
A `DiscreteVectorDictionary{T}` describes the linear space of vectors of finite length
`n` with element type `T`.
"""
struct DiscreteVectorDictionary{T} <: DiscreteDictionary{Int,T}
    n   ::  Int
end

# We set a default codomain type Float64
DiscreteVectorDictionary(n::Int) = DiscreteVectorDictionary{Float64}(n)

size(d::DiscreteVectorDictionary) = (d.n,)

similar(d::DiscreteVectorDictionary{T}, ::Type{Int}, n::Int) where {T} = DiscreteVectorDictionary{T}(n)

support(d::DiscreteVectorDictionary) = ClosedInterval{Int}(1, length(d))

"""
A `DiscreteArrayDictionary{N,T}` describes the linear space of arrays of finite size
`size(d)` with element type `T`.
"""
struct DiscreteArrayDictionary{T,N} <: DiscreteDictionary{ProductIndex{N},T}
    size    ::  NTuple{N,Int}
end

DiscreteArrayDictionary(size...) = DiscreteArrayDictionary{Float64}(size...)

DiscreteArrayDictionary{T}(size::Int...) where {T} = DicreteArrayDictionary{T}(size)

DiscreteArrayDictionary{T}(size::NTuple{N,Int}) where {T,N} =
    DiscreteArrayDictionary{T,N}(size)

DiscreteArrayDictionary(array::AbstractArray{T,N}) where {T,N} =
    DiscreteArrayDictionary{T,N}(size(array))

size(d::DiscreteArrayDictionary) = d.size

similar(d::DiscreteArrayDictionary{T}, ::Type{ProductIndex{N}}, dims::Vararg{Int,N}) where {T,N} = DiscreteArrayDictionary{T}(dims)

ordering(d::DiscreteArrayDictionary{N}) where {N} = ProductIndexList{N}(size(d))

native_index(d::DiscreteArrayDictionary, idx) = product_native_index(size(d), idx)

"""
A `DiscreteMultiArrayDictionary{A,T}` describes the linear space of multi-arrays
with element type `T`.
"""
struct DiscreteMultiArrayDictionary{T} <: DiscreteDictionary{MultilinearIndex,T}
    offsets ::  Vector{Int}
end

size(d::DiscreteMultiArrayDictionary) = (d.offsets[end],)

unsafe_offsets(d::DiscreteMultiArrayDictionary) = d.offsets

function similar(d::DiscreteMultiArrayDictionary, T::Type{MultilinearIndex}, size::Int...)
    @assert length(d) == prod(size)
    similar(d, T, d.offsets)
end

similar(d::DiscreteMultiArrayDictionary{T}, ::Type{MultilinearIndex}, offsets::Vector{Int}) where {T} = DiscreteMultiArrayDictionary{T}(offsets)

ordering(d::DiscreteMultiArrayDictionary) = MultilinearIndexList(unsafe_offsets(d))

support(d::DiscreteMultiArrayDictionary) = ClosedInterval{Int}(1, length(d))


text:
# Tests on Univariate LogitNormal distributions

using Distributions
using Random, Test
using StatsFuns

####### Core testing procedure
function test_logitnormal(g::LogitNormal, n_tsamples::Int=10^6,
                          rng::Union{AbstractRNG, Missing} = missing)
    d = length(g)
    #mn = mean(g)
    md = median(g)
    #mo = mode(g)
    #s = var(g)
    #e = entropy(g)
    @test partype(g) == Float64
    #@test isa(mn, Float64)
    @test isa(md, Float64)
    #@test isa(mo, Float64)
    #@test isa(s, Float64)
    @test md ≈ logistic(g.μ)
    #@test entropy(g) ≈ d*(1 + Distributions.log2π)/2 + logdetcov(g.normal)/2 + sum(mean(g.normal))
     @test insupport(g,1e-8)
     # corner cases of 0 and 1 handled as in support
     @test insupport(g,1.0)
     @test pdf(g,0.0) == 0.0
     @test insupport(g,0.0)
     @test pdf(g,1.0) == 0.0
     @test !insupport(g,-1e-8)
     @test pdf(g,-1e-8) == 0.0
     @test !insupport(g,1+1e-8)
     @test pdf(g,1+1e-8) == 0.0

    # sampling
    if ismissing(rng)
        X = rand(g, n_tsamples)
    else
        X = rand(rng, g, n_tsamples)
    end
    @test isa(X, Array{Float64,1})

    # evaluation of logpdf and pdf
    for i = 1:min(100, n_tsamples)
        @test logpdf(g, X[i]) ≈ log(pdf(g, X[i]))
    end
    @test logpdf.(g, X) ≈ log.(pdf.(g, X))
    @test isequal(logpdf(g, 0),-Inf)
    @test isequal(logpdf(g, 1),-Inf)
    @test isequal(logpdf(g, -eps()),-Inf)

    # test the location and scale functions
    @test location(g) == g.μ
    @test scale(g) == g.σ
    @test params(g) == (g.μ, g.σ)
end

###### General Testing
@testset "Logitnormal tests" begin
    test_logitnormal( LogitNormal() )
    test_logitnormal( LogitNormal(2,0.5) )
    d = LogitNormal(Float32(2))
    typeof(rand(d, 5)) # still Float64
    @test typeof(convert(LogitNormal{Float64}, d)) == typeof(LogitNormal(2,1))
end


text:
# Quadratic programming in MATLAB-like style
#
#   minimize x^2 + xy + y^2 + yz + z^2
#
#   s.t.    x + 2 y + 3 z >= 4
#           x +   y       >= 1
#

using Xpress, Compat.Test
using Compat.SparseArrays
using Compat.LinearAlgebra

@testset "QP 2" begin
    model = xpress_model(
    name = "qp_02",
    f = [0., 0., 0.],
    H = [2. 1. 0.; 1. 2. 1.; 0. 1. 2.],
    A = -[1. 2. 3.; 1. 1. 0.],
    b = -[4., 1.])

    @test Xpress.getq_upper(model) == triu(sparse([2. 1. 0.; 1. 2. 1.; 0. 1. 2.]))
    @test Xpress.getq(model) == sparse([2. .5 0.; .5 2. .5; 0. .5 2.])
    @test Xpress.get_obj(model) == [0,0,0]
    @test Xpress.get_constrmatrix(model) == sparse([-1 -2 -3; -1 -1 0])
    @test Xpress.get_rhs(model) == [-4,-1]
    @test Xpress.get_lb(model) == [-1e20, -1e20, -1e20]
    @test Xpress.get_ub(model) == [1e20, 1e20, 1e20]
    @test Xpress.get_sense(model) == [:(<=),:(<=)]

    optimize(model)

    ans = get_optiminfo(model)
    @test ans.status_lp == :optimal

    @test isapprox(get_solution(model), [0.5714, .4285, .8571]; atol = 1e-3)

    @test isapprox(get_objval(model), 1.857; atol = 1e-3)
end


text:
abstract type ConeNormalization end

struct TrivialNormalization <: ConeNormalization end
struct StandardNormalization <: ConeNormalization end

"""
    normalize_cglp!(::ConeNormalization, cglp, variables)

Add a constraint or set of constraints to the CGLP using the
`ConeNormalization` technique. Return the `ConstraintRef`
"""
function normalize_cglp! end

function normalize_cglp!(::TrivialNormalization, cglp::JuMP.Model, variables::NamedTuple)
    @constraint(cglp, variables.u0 + variables.v0 == 1)
end

function normalize_cglp!(::StandardNormalization, cglp::JuMP.Model, variables::NamedTuple)
    @constraint(cglp, sum(variables.u + variables.v) + variables.u0 + variables.v0 == 1)
end


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

function LRNDescriptor(n::Float64, alpha::Float64, beta::Float64, k::Float64)
  p = Ptr{Void}[0]
  cudnnCreateLRNDescriptor(p)
  cudnnSetLRNDescriptor(p[1], n, alpha, beta, k)
  desc = new(p[1])
  finalizer(desc, cudnnDestroyLRNDescriptor)
  desc
end

Base.unsafe_convert(::Type{Ptr{Void}}, desc::LRNDescriptor) = desc.ptr

function crosschanel()
  cudnnLRNCrossChannelForward()
end


text:
using Poisson

poisson(3)


text:
include("src/MyUtils.jl")
using MyUtils

@testset "project Tests" begin
    @test project(5.0,4.0,6.0,3.0,8.0) == 5.5
    @test project(0,0,1,0,1) == 0.0
    @test project(1.0,0.0,1.0,0.0,1.0) == 1.0
    @test project(1.0,0,1.0,0,2.0) == 2.0
    @test project(2,0.0,1.0,0.0,1.0) == 2.0
    @test project(0.25,1.0,0.0,0.0,1.0) == 0.75
end;


text:
include("runBatchSim.jl")
include("linspecer.jl")
using Printf

# dt = .1
# A = [1 0 0 dt;
#      0 1 dt -dt;
#      0 0 1 0;
#      0 0 0 1]
# C = [0. 1. 0 0;
#      0 0 0 1]
# muPrior = [0. 0; 10. 0; 0. 1; 1. 0]
# endT = 100.
# sigmas = [1., 100., 0., 1.]
#
# static = StaticWorld(A, C, muPrior, endT, dt)
#
# flexSigValues = [10, 1, .1]
# flexSigInd = 4
# flexVarName = "sigmaVest"

function compareSims(A::Array{Float64,2}, C::Array{Float64,2}, muPrior::Array{Float64,2}, endT::Number, dt::Number,
                    sigmas::Array{Float64,1}, flexSigValues::Array{Float64,1}, flexSigInd::Int64, flexVarName::String)
    static = StaticWorld(A, C, muPrior, endT, dt)
    compareSims(static, sigmas, flexSigValues, flexSigInd, flexVarName)
    return
end

function compareSims(A::Array{Float64,2}, C::Array{Float64,2}, muPrior::Array{Float64,2}, sigmas::Array{Float64,1},
                    flexSigValues::Array{Float64,1}, flexSigInd::Int64, flexVarName::String)
    static = StaticWorld(A, C, muPrior)
    compareSims(static, sigmas, flexSigValues, flexSigInd, flexVarName)
    return
end

function compareSims(static::StaticWorld, sigmas::Array{Float64,1}, flexSigValues::Array{Float64,1}, flexSigInd::Int64, flexVarName::String)
    set_zero_subnormals(true)
    variationNum = length(flexSigValues)
    colors = linspecer(variationNum)
    # figure()
    for v = 1:variationNum
        println(string("Variation #", v))
        plotOpts = PlotOpts(@sprintf("%s = %s", flexVarName, string.(flexSigValues[v])), colors[v,:])
        setindex!(sigmas, flexSigValues[v], flexSigInd)
        simOpts = SimOpts(sigmas, 1000)
        @time runBatchSim(plotOpts, static, simOpts)
    end
    legend()
    return
end

max_stars_repo_name: rhoadesScholar/Drugowitsch-Rotation

max_issues_repo_name: rhoadesScholar/Drugowitsch-Rotation

max_forks_repo_name: rhoadesScholar/Drugowitsch-Rotation


text:
using DataDeps

register(DataDep(
  "arrhythmia",
  "http://archive.ics.uci.edu/ml/datasets/Arrhythmia",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/arrhythmia/arrhythmia.data",
  "a7f0f4ca289a4c58b5ed85a9adb793189acd38425828ce3dfbb70adb45f30169",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=280,
        feature_indices=1:279,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
mutable struct FullTrackedState <: State
    queues::Vector{Vector{FullJob}} #A vector of queues holding the jobs waiting in buffer
    in_transit::BinaryMinHeap{FullJob} #Jobs in transit between queues
    left_system::Vector{FullJob} #Jobs that have left the system
    params::NetworkParameters #The parameters of the queueing system
    arrival_distr::Gamma{Float64}
    travel_distr::Gamma{Float64}
    service_distrs::Vector{Gamma{Float64}}
end

"""
Constructs an initilized tracked network state from the given parameters
"""
FullTrackedState(params::NetworkParameters) = FullTrackedState(
    [Vector{FullJob}[] for _ in 1:params.L], # Initial queues
    BinaryMinHeap{FullJob}(),                # Intial transit
    Vector{FullJob}[],                       # Initial left system
    params,
    rate_scv_gamma(params.λ, params.gamma_scv),
    rate_scv_gamma(params.η, params.gamma_scv),
    map((μ) -> rate_scv_gamma(μ, params.gamma_scv), params.μ_vector)
    )

queued_count(state::FullTrackedState)::Int = sum([length(state.queues[i]) for i in 1:length(state.queues)])

transit_count(state::FullTrackedState)::Int = length(state.in_transit)

num_in_queue(q::Int, state::FullTrackedState)::Int = length(state.queues[q])

function clear_left(state::FullTrackedState)
    empty!(state.left_system)
end

""" Prints a full history of all jobs that the system processes, sorted by entry time"""
function show(io::IO, state::FullTrackedState)
    all_jobs = FullJob[]
    for queue in state.queues
        append!(all_jobs, deepcopy(queue))
    end

    in_transit = deepcopy(state.in_transit)
    append!(all_jobs, extract_all!(in_transit))

    append!(all_jobs, deepcopy(state.left_system))

    sort!(all_jobs, by = (job) -> entry_time(job))
    println(io, "Tracked State: ")
    for job in all_jobs
        print(io, job.history)
        print(io, "\n")
    end
end

function new_job(q::Int, time::Float64, state::FullTrackedState)
    push!(state.queues[q], FullJob(InQueue(time, q)))
end

function new_transit(time::Float64, transit_time::Float64, state::FullTrackedState)
    push!(state.in_transit, FullJob(InTransit(time, transit_time)))
end

function pop_queue(q::Int, time::Float64, state::FullTrackedState)
    job = popfirst!(state.queues[q])
    push!(job, LeaveSystem(time))
    push!(state.left_system, job)
end

function queue_to_transit(q::Int, time::Float64, transit_time::Float64, state::FullTrackedState)
    job = popfirst!(state.queues[q])
    push!(job, InTransit(time, transit_time))
    push!(state.in_transit, job)
end

function pop_transit(time::Float64, state::FullTrackedState)
    job = pop!(state.in_transit)
    push!(job, LeaveSystem(time))
    push!(state.left_system, job)
end

function transit_to_queue(q::Int, time::Float64, state::FullTrackedState)
    job = pop!(state.in_transit)
    push!(job, InQueue(time, q))
    push!(state.queues[q], job)
end

function update_transit(time::Float64, transit_time::Float64, state::FullTrackedState)
    job = pop!(state.in_transit)
    push!(job, InTransit(time, transit_time))
    push!(state.in_transit, job)
end

function failed_arrival(time::Float64, state::FullTrackedState)
    push!(state.left_system, FullJob(LeaveSystem(time)))
end


text:
using AtomsBase
using Test
using Unitful
using UnitfulAtomic
using PeriodicTable

@testset "Interface" begin
    box = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]u"m"
    bcs = [Periodic(), Periodic(), DirichletZero()]
    positions = [[0.25, 0.25, 0.25], [0.75, 0.75, 0.75]]u"m"
    elements = [:C, :C]
    atoms = [Atom(elements[i], positions[i]) for i in 1:2]

    @testset "Atoms" begin
        @test position(atoms[1]) == [0.25, 0.25, 0.25]u"m"
        @test velocity(atoms[1]) == [0.0, 0.0, 0.0]u"bohr/s"
        @test element(atoms[1]) == PeriodicTable.elements[:C]
        @test atomic_symbol(atoms[1]) == :C
        @test atomic_number(atoms[1]) == 6
        @test atomic_mass(atoms[1])   == 12.011u"u"
    end

    @testset "System" begin
        flexible = FlexibleSystem(atoms, box, bcs)
        fast     = FastSystem(flexible)
        @test length(flexible) == 2
        @test size(flexible)   == (2, )

        @test bounding_box(flexible) == [[1, 0, 0], [0, 1, 0], [0, 0, 1]]u"m"
        @test boundary_conditions(flexible) == [Periodic(), Periodic(), DirichletZero()]
        @test periodicity(flexible) == [1, 1, 0]
        @test n_dimensions(flexible) == 3
        @test position(flexible) == [[0.25, 0.25, 0.25], [0.75, 0.75, 0.75]]u"m"
        @test position(flexible, 1) == [0.25, 0.25, 0.25]u"m"
        @test velocity(flexible)[1] == [0.0, 0.0, 0.0]u"bohr/s"
        @test velocity(flexible)[2] == [0.0, 0.0, 0.0]u"bohr/s"
        @test atomic_mass(flexible)   == [12.011, 12.011]u"u"
        @test atomic_number(fast) == [6, 6]
        @test atomic_number(fast, 1) == 6
        @test ismissing(velocity(fast, 2))
        @test atomic_symbol(flexible, 2) == :C
        @test atomic_number(flexible, 2) == 6
        @test atomic_mass(flexible, 1)   == 12.011u"u"

        @test ismissing(velocity(fast))
        @test all(position(fast)      .== position(flexible))
        @test all(atomic_symbol(fast) .== atomic_symbol(flexible))
    end
end


text:
# This file was generated, do not modify it. # hide
datac = coerce(datac, autotype(datac, rules=(:discrete_to_continuous,)));


text:
module ImageRegression
import ..AbstractOperator.get_misfit
using ..AbstractOperator
using PyPlot, LinearAlgebra, Statistics
using Random, Images, ImageIO, FileIO, ..CommonToAll

import ..Model, ..Options, ..OptionsStat, ..OptionsNonstat

export Img, get_image_data, calc_image_RMS, get_image_prediction_points,
        plot_image_data, plot_image_posterior, slice_image_posterior

mutable struct Img<:Operator2D

    filename          :: String
    x                 :: StepRangeLen
    y                 :: StepRangeLen
    fractrain         :: Float64
    dec               :: Int
    gausskernelwidth  :: Int
    f                 :: Array{Float64, 2}
    d                 :: Array{Float64, 2}
    useML             :: Bool
    σ                 :: Float64
    select            :: Array{Bool}
    N                 :: Int
    r                 :: Array{Float64, 1}
    xtrain            :: Array{Float64, 2}
end

function Img(;
             filename         = "",
             dx               = 0.01,
             fractrain        = 0.02,
             dec::Int         = 2,
             gausskernelwidth = 7,
             useML            = false,
             σ                = 1.0,
             select           = [true],
             N                = 1,
             r                = [-1.],
             xtrain           = [1 1]
             )

    @assert fractrain > 0 && fractrain < 1
    @assert !(filename == "")
    @assert dec > 1
    @assert σ > 0.0
    f = get_image(filename, gausskernelwidth, dec)
    x = 0:dx:dx*size(f,2)-1
    y = 0:dx:dx*size(f,1)-1
    Img(filename, x, y, fractrain, dec, gausskernelwidth, f, f, useML, σ, select, N, r, xtrain)
end

function get_image(filename::String, gausskernelwidth::Int, dec)
    f = Gray.(load(filename))
    f = convert(Array{Float64, 2}, f)[1:dec:end,:1:dec:end]
    # convert image into something that looks like log resistivity
    f = -1 .+ 3*f
    # smooth image
    imfilter(f,Kernel.gaussian(gausskernelwidth))
end

function get_image_data(img::Img;
                           sdmaxfrac = 0.05,
                           rseed     = 12,
                           ybreak    = 12321.0,
                           takeevery = 4
                          )
    @assert ybreak != 12321.0
    @assert sdmaxfrac > 0 && sdmaxfrac < 1
    Random.seed!(rseed)
    f, x, y = img.f, img.x, img.y
    noisyd = NaN .+ zeros(Float64, size(f))
    ntrain = round(Int, img.fractrain*length(f))
    δtry = sdmaxfrac*max(f...)
    Xtrain = zeros(2,0)
    linidx = randperm(length(f))[1:ntrain]
    for (i,l) in enumerate(linidx)
        row, col = Tuple(CartesianIndices(f)[l])
        if img.y[row] > ybreak
            noisyd[row, col] = f[row, col] + δtry*randn()
            Xtrain = hcat(Xtrain, [x[col]; y[row]])
        else
            if rem(row, takeevery) == 0 && rem(col, takeevery) == 0
                noisyd[row, col] = f[row, col] + δtry*randn()
                Xtrain = hcat(Xtrain, [x[col]; y[row]])
            end
        end
    end
    img.d = noisyd
    img.σ = δtry
    img.select = .!isnan.(noisyd[:])
    img.r = noisyd[img.select]
    img.xtrain = Xtrain
    img.N = sum(img.select)
    nothing
end

function get_image_prediction_points(img::Img)
    f, x, y = img.f, img.x, img.y
    Xall = zeros(2, length(f))
    for i in 1:length(f)
        yid, xid = Tuple(CartesianIndices(f)[i])
        Xall[:,i] = [x[xid]; y[yid]]
    end
    Xall
end

function plot_image_data(img::Img;  s=3, fsize=12)
    f, x, y = img.f, img.x, img.y
    ftrain, Xtrain = img.d[img.select], img.xtrain
    f1, ax1 = plt.subplots(1,2,figsize=(6.7,2.7), sharex=true, sharey=true)
    im1 = ax1[1].imshow(f, extent=[x[1],x[end],y[end],y[1]])
    cb1 = colorbar(im1, ax=ax1[1])
    ax1[2].imshow(f, extent=[x[1],x[end],y[end],y[1]], alpha=0.0)
    im2 = ax1[2].scatter(Xtrain[1,:], Xtrain[2,:], c=ftrain, s=s)
    cb2 = colorbar(im2)
    ax1[2].axis([x[1],x[end],y[end],y[1]])
    nicenup(gcf(), fsize=fsize)

end

function calc_image_RMS(img::Img)
    r = (img.d[img.select] - img.f[img.select])/img.σ
    @info "χ^2 error is $(r'*r) for $(img.N) points RMS: $(sqrt(norm(r)^2/img.N))"
    nothing
end

function get_misfit(m::Model, opt::Options, img::Img)
    chi2by2 = 0.0
    if !opt.debug
        d, r = img.d, img.r
        select = img.select
        r .= m.fstar[select] - d[select]
        if img.useML
            chi2by2 = 0.5*img.N*log(norm(r)^2)
        else
            chi2by2 = norm(r)^2/(2img.σ^2)
        end
    end
    return chi2by2
end

function slice_image_posterior( M::AbstractArray,
                                opt::Options,
                                img::Img, roworcol::Symbol;
                                rowcolnum = 1,
                                nbins = 50,
                                rhomin=Inf,
                                rhomax=-Inf,
                                qp1=0.05,
                                qp2=0.95,
                                islscale = false,
                                pdfnormalize=false,
                                temperaturenum=1)
    @assert temperaturenum == 1
    Mslices = Array{Array{Float64}, 1}(undef, length(M))
    for i in 1:length(M)
        m = reshape(M[i],length(img.y), length(img.x))
        if roworcol == :row
            Mslices[i] = m[rowcolnum,:]
        else
            Mslices[i] = m[:,rowcolnum]
        end
    end
    himage, edges, CI, = gethimage(img, Mslices, opt, burninfrac=0.0, temperaturenum=temperaturenum,
                nbins=nbins, rhomin=rhomin, rhomax=rhomax, qp1=qp1, qp2=qp2, islscale=islscale,
                pdfnormalize=pdfnormalize)
end

function get_image_quantile(M::AbstractArray, q=0.5)
    mall = hcat(M...)
    quantM = zeros(length(M[1]))
    for i = 1:length(M[1])
        quantM[i] = quantile(mall[i,:], q)
    end
    quantM
end

function plot_image_posterior(optns::OptionsNonstat, opts::OptionsStat,
                        img::Img;
                        rownum = 10,
                        colnum = 10,
                        temperaturenum = 1,
                        nbins = 50,
                        burninfrac=0.5,
                        qp1=0.05,
                        qp2=0.95,
                        cmappdf = "magma",
                        cmapmean = "bone",
                        figsizerows=[11,6],
                        figsizecols=[7.5,9.5],
                        pdfnormalize=false,
                        fsize=10,
                        aspect=1.0,
                        domean = false,
                        getquant = 0.5)
    @assert temperaturenum == 1
    @assert 0.0 < getquant < 1.0
    M = assembleTat1(optns, :fstar, burninfrac=burninfrac, temperaturenum=temperaturenum)
    if domean
        mns = mean(M)
    else
        mns = get_image_quantile(M, getquant)
    end
    mns = reshape(mns, length(img.y), length(img.x))
    himage_r_ns, edges_r_ns, CI_r_ns = slice_image_posterior(M, optns, img, :row, rowcolnum=rownum, nbins = nbins, qp1=qp1, qp2=qp2,
                                pdfnormalize=pdfnormalize, temperaturenum=temperaturenum)
    himage_c_ns, edges_c_ns, CI_c_ns = slice_image_posterior(M, optns, img, :col, rowcolnum=colnum, nbins = nbins, qp1=qp1, qp2=qp2,
                                pdfnormalize=pdfnormalize, temperaturenum=temperaturenum)
    M = assembleTat1(opts, :fstar, burninfrac=burninfrac, temperaturenum=temperaturenum)
    if domean
        ms = 0.5*log10.(mean(M))
    else
        a, b = [m[1,:] for m in M], [m[2,:] for m in M]
        ms = 0.5*log10.([get_image_quantile(a, getquant) get_image_quantile(b, getquant)])'
    end
    msx = reshape(ms[1,:], length(img.y), length(img.x))
    msy = reshape(ms[2,:], length(img.y), length(img.x))
    Mx, My = Array{Array{Float64}, 1}(undef, length(M)), Array{Array{Float64}, 1}(undef, length(M))
    for i = 1:length(M)
        Mx[i] = M[i][1,:]
        My[i] = M[i][2,:]
    end
    himage_x_c, edges_x_c, CI_x_c = slice_image_posterior(Mx, opts, img, :col, rowcolnum=colnum, nbins = nbins, qp1=qp1, qp2=qp2,
                                pdfnormalize=pdfnormalize, temperaturenum=temperaturenum, islscale=true)
    himage_y_c, edges_y_c, CI_y_c = slice_image_posterior(My, opts, img, :col, rowcolnum=colnum, nbins = nbins, qp1=qp1, qp2=qp2,
                                pdfnormalize=pdfnormalize, temperaturenum=temperaturenum, islscale=true)
    himage_x_r, edges_x_r, CI_x_r = slice_image_posterior(Mx, opts, img, :row, rowcolnum=rownum, nbins = nbins, qp1=qp1, qp2=qp2,
                                pdfnormalize=pdfnormalize, temperaturenum=temperaturenum, islscale=true)
    himage_y_r, edges_y_r, CI_y_r = slice_image_posterior(My, opts, img, :row, rowcolnum=rownum, nbins = nbins, qp1=qp1, qp2=qp2,
                                pdfnormalize=pdfnormalize, temperaturenum=temperaturenum, islscale=true)
    f = figure(figsize=figsizerows)
    s1 = subplot(231)
    im1 = s1.pcolormesh(img.x, img.y, mns, cmap=cmapmean)
    s1.plot([img.x[1], img.x[end]], [img.y[rownum], img.y[rownum]], "--k", alpha=0.5)
    central_tendency = domean ? ("Mean") : ("Percentile "*"$(round(Int, getquant*100))")
    s1.set_title(central_tendency*" pixel value")
    s1.set_ylabel("distance y")
    # s1.set_aspect(aspect)
    s1.set_xlabel("distance x")
    cb1 = colorbar(im1, ax=s1)
    s2 = subplot(232, sharex=s1, sharey=s1)
    im2 = s2.pcolormesh(img.x, img.y, msx, cmap=cmapmean)
    s2.plot([img.x[1], img.x[end]], [img.y[rownum], img.y[rownum]], "--k", alpha=0.5)
    s2.set_title(central_tendency*" "*L"\log_{10}λ_x")
    s2.set_ylabel("distance y")
    s2.set_xlabel("distance x")
    # s2.set_aspect(aspect)
    s3 = subplot(233, sharex=s1, sharey=s1)
    cb2 = colorbar(im2, ax=s2)
    im3 = s3.pcolormesh(img.x, img.y, msy, cmap=cmapmean)
    s3.plot([img.x[1], img.x[end]], [img.y[rownum], img.y[rownum]], "--k", alpha=0.5)
    # s3.set_aspect(aspect)
    s3.invert_yaxis()
    s3.set_title(central_tendency*" "*L"\log_{10}λ_y")
    s3.set_ylabel("distance y")
    s3.set_xlabel("distance x")
    cb3 = colorbar(im3, ax=s3)
    s4 = subplot(234, sharex=s1)
    im4 = s4.pcolormesh(img.x, edges_r_ns[:], himage_r_ns', cmap=cmappdf)
    s4.plot(img.x, mns[rownum,:][:],"--w", alpha=0.5)
    s4.plot(img.x, CI_r_ns,"--c", alpha=0.5)
    s4.set_title("Pixel value PDF")
    s4.set_ylabel("pixel value")
    s4.set_xlabel("distance x")
    cb4 = colorbar(im4, ax=s4)
    cb4.ax.set_title("pdf")
    s5 = subplot(235, sharex=s1)
    im5 = s5.pcolormesh(img.x, edges_x_r[:], himage_x_r', cmap=cmappdf)
    cb5 = colorbar(im5, ax=s5)
    cb5.ax.set_title("pdf")
    s5.plot(img.x, msx[rownum,:][:],"--w", alpha=0.5)
    s5.plot(img.x, CI_x_r,"--c", alpha=0.5)
    s5.set_title(L"\log_{10}λ_x"*" PDF")
    s5.set_ylabel(L"\log_{10}λ_x")
    s5.set_xlabel("distance x")
    s6 = subplot(236, sharex=s1)
    im6 = s6.pcolormesh(img.x, edges_y_r[:], himage_y_r', cmap=cmappdf)
    cb6 = colorbar(im6, ax=s6)
    cb6.ax.set_title("pdf")
    s6.plot(img.x, msy[rownum,:][:],"--w", alpha=0.5)
    s6.plot(img.x, CI_y_r,"--c", alpha=0.5)
    s6.set_title(L"\log_{10}λ_y"*" PDF")
    s6.set_ylabel(L"\log_{10}λ_y")
    s6.set_xlabel("distance x")
    nicenup(f, fsize=fsize)

    f = figure(figsize=figsizecols)
    s1 = subplot(321)
    im1 = s1.pcolormesh(img.x, img.y, mns, cmap=cmapmean)
    s1.plot([img.y[colnum], img.y[colnum]], [img.y[1], img.y[end]], "--k", alpha=0.5)
    s1.set_title(central_tendency*" pixel value")
    s1.set_ylabel("distance y")
    s1.set_xlabel("distance x")
    # s1.set_aspect(aspect)
    cb1 = colorbar(im1, ax=s1)
    s2 = subplot(323, sharex=s1, sharey=s1)
    im2 = s2.pcolormesh(img.x, img.y, msx, cmap=cmapmean)
    s2.plot([img.y[colnum], img.y[colnum]], [img.y[1], img.y[end]], "--k", alpha=0.5)
    s2.set_title(central_tendency*" "*L"\log_{10}λ_x")
    s2.set_ylabel("distance y")
    s2.set_xlabel("distance x")
    # s2.set_aspect(aspect)
    s3 = subplot(325, sharex=s1, sharey=s1)
    cb2 = colorbar(im2, ax=s2)
    im3 = s3.pcolormesh(img.x, img.y, msy, cmap=cmapmean)
    s3.plot([img.y[colnum], img.y[colnum]], [img.y[1], img.y[end]], "--k", alpha=0.5)
    # s3.set_aspect(aspect)
    s3.invert_yaxis()
    s3.set_title(central_tendency*" "*L"\log_{10}λ_y")
    s3.set_ylabel("distance y")
    s3.set_xlabel("distance x")
    cb3 = colorbar(im3, ax=s3)
    s4 = subplot(322, sharey=s1)
    im4 = s4.pcolormesh(edges_c_ns[:], img.y, himage_c_ns, cmap=cmappdf)
    s4.plot(mns[:,colnum][:], img.y, "--w", alpha=0.5)
    s4.plot(CI_c_ns, img.y, "--c", alpha=0.5)
    s4.set_title("Pixel value PDF")
    s4.set_xlabel("pixel value")
    s4.set_ylabel("distance y")
    cb4 = colorbar(im4, ax=s4)
    cb4.ax.set_title("pdf")
    s5 = subplot(324, sharey=s1)
    im5 = s5.pcolormesh(edges_x_c[:], img.y, himage_x_c, cmap=cmappdf)
    cb5 = colorbar(im5, ax=s5)
    cb5.ax.set_title("pdf")
    s5.plot(msx[:,colnum][:], img.y, "--w", alpha=0.5)
    s5.plot(CI_x_c, img.y, "--c", alpha=0.5)
    s5.set_title(L"\log_{10}λ_x"*" PDF")
    s5.set_xlabel(L"\log_{10}λ_x")
    s5.set_ylabel("distance y")
    s6 = subplot(326, sharey=s1)
    im6 = s6.pcolormesh(edges_y_c[:], img.y, himage_y_c, cmap=cmappdf)
    cb6 = colorbar(im6, ax=s6)
    cb6.ax.set_title("pdf")
    s6.plot(msy[:,colnum][:], img.y, "--w", alpha=0.5)
    s6.plot(CI_y_c, img.y, "--c", alpha=0.5)
    s6.set_title(L"\log_{10}λ_y"*" PDF")
    s6.set_xlabel(L"\log_{10}λ_y")
    s6.set_ylabel("distance y")
    nicenup(f, fsize=fsize)
end

function plot_image_posterior(opt::Options,
                        img::Img;
                        rownum = 10,
                        colnum = 10,
                        temperaturenum = 1,
                        nbins = 50,
                        burninfrac=0.5,
                        qp1=0.05,
                        qp2=0.95,
                        cmappdf = "magma",
                        cmapmean = "bone",
                        figsize=[7.5,5.9],
                        pdfnormalize=false,
                        fsize=10,
                        aspect=1.0,
                        domean = false,
                        getquant = 0.5)
    @assert temperaturenum == 1
    @assert 0.0 < getquant < 1.0
    if typeof(opt) == OptionsStat
        @assert opt.updatenonstat == false
        @assert opt.needλ²fromlog == false
    end
    M = assembleTat1(opt, :fstar, burninfrac=burninfrac, temperaturenum=temperaturenum)
    if domean
        m = mean(M)
    else
        a = [m' for m in M]
        m = get_image_quantile(a, getquant)
    end
    m = reshape(m, length(img.y), length(img.x))
    himage_r, edges_r, CI_r = slice_image_posterior(M, opt, img, :row, rowcolnum=rownum, nbins = nbins, qp1=qp1, qp2=qp2,
                                pdfnormalize=pdfnormalize, temperaturenum=temperaturenum)
    himage_c, edges_c, CI_c = slice_image_posterior(M, opt, img, :col, rowcolnum=colnum, nbins = nbins, qp1=qp1, qp2=qp2,
                                pdfnormalize=pdfnormalize, temperaturenum=temperaturenum)
    f = figure(figsize=figsize)
    s1 = subplot(223)
    central_tendency = domean ? ("Mean") : ("Percentile "*"$(round(Int, getquant*100))")
    im1 = s1.pcolormesh(img.x, img.y, m, cmap=cmapmean)
    s1.plot([img.x[1], img.x[end]], [img.y[rownum], img.y[rownum]], "--k", alpha=0.5)
    s1.plot([img.y[colnum], img.y[colnum]], [img.y[1], img.y[end]], "--k", alpha=0.5)
    s1.set_title(central_tendency*" pixel value")
    s1.set_ylabel("distance y")
    s1.invert_yaxis()
    # s1.set_aspect(aspect)
    s1.set_xlabel("distance x")
    cb1 = colorbar(im1, ax=s1)
    s2 = subplot(221, sharex=s1)
    im2 = s2.pcolormesh(img.x, edges_r[:], himage_r', cmap=cmappdf)
    s2.plot(img.x, m[rownum,:][:],"--w", alpha=0.5)
    s2.plot(img.x, CI_r, "--c", alpha=0.5)
    s2.set_title("Pixel value PDF")
    s2.set_ylabel("pixel value")
    s2.set_xlabel("distance x")
    cb2 = colorbar(im2, ax=s2)
    cb2.ax.set_title("pdf")
    s3 = subplot(224, sharey=s1)
    im3 = s3.pcolormesh(edges_c[:], img.y, himage_c, cmap=cmappdf)
    s3.plot(m[:,colnum][:], img.y, "--w", alpha=0.5)
    s3.plot(CI_c, img.y, "--c", alpha=0.5)
    s3.set_title("Pixel value PDF")
    s3.set_xlabel("pixel value")
    s3.set_ylabel("distance y")
    cb3 = colorbar(im3, ax=s3)
    cb3.ax.set_title("pdf")
    nicenup(f, fsize=fsize)
end

end

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
@safetestset "Pendulum Period Tests" begin
    include("pendulum_test.jl")
end

@safetestset "Nutation Behavior Test" begin
    include("nutation_test.jl")
end

@safetestset "Linearization Test" begin
    include("linearization_test.jl")
end

@safetestset "Generic Joint Test" begin
    include("genericjoint_test.jl")
end


text:
# Deadfish implementation in Julia

# Run fully compliant version (infinite input/output loop):
#                   deadfish() or deadfish(0)

# Run version that exits when entering 'h':
#                   deadfish(1)

# Run version that simulates manual input and
# appropriate output, given the input string as a whole:
#                   deadfish("input string")

function deadfish()
#= Run fully compliant version of deadfish if
   no parameters are given =#
    deadfish(0)
end

function deadfish(compliance::Int)
    x=0
    while true
        print(">> ")

        input=readline(STDIN)[1]

        # Make sure x is not greater then 256 ;)
        x == 256 ? x=0 : nothing
        x == -1  ? x=0 : nothing
        input == 'i' ? x=x+1 :
        input == 'd' ? x=x-1 :
        input == 'o' ? println(x) :
        input == 's' ? x=x^2 :
        input == 'h' && compliance ==1 ? break : println()
    end
end

function deadfish(code::AbstractString)
    x=0
    for i = 1:length(code)
        input=code[i]
        println(">> $input")

        # Make sure x is not greater then 256 ;)
        x == 256 ? x=0 : nothing
        x == -1  ? x=0 : nothing
        input == 'i' ? x=x+1 :
        input == 'd' ? x=x-1 :
        input == 'o' ? println(x) :
        input == 's' ? x=x^2 : println()
    end
end


text:
@inline iszero{T<:Real}(x::T) = (x == zero(T))
@inline isone{T<:Real}(x::T) = (x == one(T))

@inline function isfinite{T<:Rvl}(ivl::T)
    isfinite(ivl.lo) & isfinite(ivl.hi)
end
@inline function issubnormal{T<:Rvl}(ivl::T)
    issubnormal(ivl.lo) | issubnormal(ivl.hi)
end
@inline function isnormal{T<:Rvl}(ivl::T)
    isfinite(ivl) & !issubnormal(ivl)
end
@inline function isinf{T<:Rvl}(ivl::T)
    isinf(ivl.lo) | isinf(ivl.hi)
end
@inline function isnan{T<:Rvl}(ivl::T)
    isnan(ivl.lo) | isnan(ivl.hi)
end

@inline function isPoint{T<:Rvl}(ivl::T)
    (ivl.lo == ivl.hi) && isfinite(ivl.lo)
end
@inline function isNotPoint{T<:Rvl}(ivl::T)
    (ivl.lo != ivl.hi) || !isfinite(ivl.lo)
end
@inline function isBoundry{T<:Rvl}(ivl::T)
    (ivl.lo != ivl.hi) && isnan(ivl)
end
@inline function isNotBoundry{T<:Rvl}(ivl::T)
    (ivl.lo == ivl.hi) || !isnan(ivl)
end
@inline function isExact{T<:Rvl}(ivl::T)
    (isExact(ivl.lo) & isExact(ivl.hi)) && (isfinite(ivl.lo) & isfinite(ivl.hi))
end
@inline function isNotExact{T<:Rvl}(ivl::T)
    (isNotExact(ivl.lo) | isNotExact(ivl.hi)) || (!isfinite(ivl.lo) | !isfinite(ivl.hi))
end

@inline function isExactPoint{T<:Rvl}(ivl::T)
    isPoint(ivl) && isExact(ivl)
end
@inline function isNotExactPoint{T<:Rvl}(ivl::T)
    isNotPoint(ivl) || isNotExact(ivl)
end
@inline function isInexactPoint{T<:Rvl}(ivl::T)
    isPoint(ivl) && isNotExact(ivl)
end
@inline function isNotInexactPoint{T<:Rvl}(ivl::T)
    isNotPoint(ivl) || isExact(ivl)
end

@inline function isExactBoundry{T<:Rvl}(ivl::T)
    isBoundry(ivl) && isExact(ivl)
end
@inline function isNotExactBoundry{T<:Rvl}(ivl::T)
    isNotBoundry(ivl) || isNotExact(ivl)
end
@inline function isInexactBoundry{T<:Rvl}(ivl::T)
    isBoundry(ivl) && isNotExact(ivl)
end
@inline function isNotInexactBoundry{T<:Rvl}(ivl::T)
    isNotBoundry(ivl) || isExact(ivl)
end

@inline function isClosedOpen{T<:Rvl}(ivl::T)
    isExact(ivl.lo) & isInexact(ivl.hi)
end
@inline function isOpenClosed{T<:Rvl}(ivl::T)
    isInexact(ivl.lo) & isExact(ivl.hi)
end
@inline function isOpenOpen{T<:Rvl}(ivl::T)
    isInexact(ivl.lo) & isInexact(ivl.hi)
end
@inline function isClosedClosed{T<:Rvl}(ivl::T)
    isExact(ivl.lo) & isExact(ivl.hi)
end


text:
using DrWatson
@quickactivate
using ArgParse
using GenerativeAD
using GenerativeAD.Models: anomaly_score, generalized_anomaly_score_gpu
using BSON
using PyCall
using StatsBase: fit!, predict, sample

using Flux
using MLDataPattern
using OrderedCollections

s = ArgParseSettings()
@add_arg_table! s begin
   "max_seed"
		default = 1
		arg_type = Int
		help = "maximum number of seeds to run through"
	"dataset"
		default = "iris"
		arg_type = String
		help = "dataset"
	"sampling"
		default = "random"
		arg_type = String
		help = "sampling of hyperparameters"
    "contamination"
    	arg_type = Float64
    	help = "contamination rate of training data"
    	default = 0.0
end
parsed_args = parse_args(ARGS, s)
@unpack dataset, max_seed, sampling, contamination = parsed_args

modelname ="DeepSVDD"

function sample_params()
	parameters_rng = (
		hdim   		= 2 .^(4:9),
		zdim   		= 2 .^(3:8),
		nlayers   	= 3:4,
		activation  = ["relu", "swish", "tanh"],
		objective   = ["soft-boundary", "one-class"],
		nu   		= [0.01f0, 0.1f0, 0.5f0, 0.99f0],
		lr_svdd  	= 10f0 .^ (-4:-3),
		lr_ae   	= 10f0 .^ (-4:-3),
		batch_size   = 2 .^ (5:7),
		init_seed   = 1:Int(1e8),
	)

	return (;zip(keys(parameters_rng), map(x->sample(x, 1)[1], parameters_rng))...)
end

function create_space()
	pyReal = pyimport("skopt.space")["Real"]
	pyInt = pyimport("skopt.space")["Integer"]
	pyCat = pyimport("skopt.space")["Categorical"]

	parameters_rng = (
		hdim   		= pyInt(4, 9,						                name="log2_hdim"),
		zdim   		= pyInt(3, 8,						                name="log2_zdim"),
		nlayers   	= pyInt(3, 4,						                name="nlayers"),
		activation  = pyCat(categories=["relu", "swish", "tanh"],       name="activation"),
		objective   = pyCat(categories=["soft-boundary", "one-class"],  name="objective"),
		nu     	    = pyReal(1f-2, 0.99f0,                    	        name="nu"),
		lr_svdd   	= pyReal(1f-5, 1f-3, prior="log-uniform",	        name="lr_svdd"),
		lr_ae   	= pyReal(1f-5, 1f-3, prior="log-uniform",	        name="lr_ae"),
		batch_size  = pyInt(5, 7, 								        name="log2_batch_size"),
	)
end
function fit(data, parameters)
	all_parameters = merge(
		parameters,
		(
			iters=5000,
			check_every=30,
			ae_check_every = 30,
			patience=10,
			ae_iters=5000
		)
	)

    svdd = GenerativeAD.Models.svdd_constructor(;all_parameters...) |> gpu

	# define optimiser
	try
		global info, fit_t, _, _, _ = @timed fit!(svdd, data, all_parameters)
	catch e
		println("Error caught => $(e).")
		return (fit_t = NaN, model = nothing, history = nothing, n_parameters = NaN), []
	end

    training_info = (
		fit_t = fit_t,
		model = info[2]|>cpu,
		history_svdd = info[1][1], # losses through time
		history_ae = info[1][2],
		npars = info[3], # number of parameters
		iters = info[4] # optim iterations of model
		)

    return training_info, [(x -> GenerativeAD.Models.anomaly_score_gpu(svdd, x), parameters)]
    # if there is no gpu on pc anomaly_score will automaticly run on cpu
end

#_________________________________________________________________________________________________

try_counter = 0
max_tries = 10*max_seed
cont_string = (contamination == 0.0) ? "" : "_contamination-$contamination"
sampling_string = sampling == "bayes" ? "_bayes" : ""
prefix = "experiments$(sampling_string)/tabular$(cont_string)"
dataset_folder = datadir("$(prefix)/$(modelname)/$(dataset)")
while try_counter < max_tries
	if sampling == "bayes"
		parameters = GenerativeAD.bayes_params(
								create_space(),
								dataset_folder,
								sample_params; add_model_seed=true)
	else
		parameters = sample_params()
	end

    for seed in 1:max_seed
		savepath = joinpath(dataset_folder, "seed=$(seed)")
		mkpath(savepath)

		# get data
		data = GenerativeAD.load_data(dataset, seed=seed, contamination=contamination)

        @info "Trying to fit $modelname on $dataset with parameters $(parameters)..."
		@info "Train/validation/test splits: $(size(data[1][1], 2)) | $(size(data[2][1], 2)) | $(size(data[3][1], 2))"
		@info "Number of features: $(size(data[1][1])[1])"

		# update parameter
		parameters = merge(parameters, (idim=size(data[1][1],1), ))
		# here, check if a model with the same parameters was already tested
		if GenerativeAD.check_params(savepath, parameters)
			#(X_train,_), (X_val, y_val), (X_test, y_test) = data
			training_info, results = fit(data, parameters)
			# saving model separately
			if training_info.model !== nothing
				tagsave(joinpath(savepath, savename("model", parameters, "bson", digits=5)), Dict("model"=>training_info.model), safe = true)
				training_info = merge(training_info, (model = nothing,))
			end
			save_entries = merge(training_info, (modelname = modelname, seed = seed, dataset = dataset, contamination = contamination))
			# now loop over all anomaly score funs
			all_scores = [GenerativeAD.experiment(result..., data, savepath; save_entries...) for result in results]
			if sampling == "bayes" && length(all_scores) > 0
				@info("Updating cache with $(length(all_scores)) results.")
				GenerativeAD.update_bayes_cache(dataset_folder,
						all_scores;
						ignore=Set([:init_seed, :idim, :iters, :check_every, :patience, :ae_check_every, :ae_iters]))
			end
			global try_counter = max_tries + 1
		else
			@info "Model already present, sampling new hyperparameters..."
			global try_counter += 1
		end
	end
end
(try_counter == max_tries) ? (@info "Reached $(max_tries) tries, giving up.") : nothing


text:
# function ElectricField(sim::Simulation{1},m::Integer=1)
    # val = zeros(ComplexF64,3length(sim),m)
    # return ElectricField(sim,val)
# end

function ScalarField(sim::Simulation{1},val::AbstractVecOrMat) where {N,M}
    bl1, bl2 = sim.boundary.bls
    start = sim.x[1]
    stop = sim.x[end]
    typeof(bl1)<:AbstractComplexBL{1} ? start = start + (bl1.depth+8sim.dx) : nothing
    typeof(bl2)<:AbstractComplexBL{2} ? stop  = stop  - (bl2.depth+8sim.dx) : nothing
    start_ind = ceil(Int,Lattices.latticeindex(sim.lattice,start))
    stop_ind = floor(Int,Lattices.latticeindex(sim.lattice,stop))
    return VectorField{1,1}(sim.x,val,sim.lattice[start_ind],sim.lattice[stop_ind],[start_ind+1],[stop_ind+1])
end

function ElectricField(sim::Simulation{1},val::AbstractVecOrMat) where {N,M}
    bl1, bl2 = sim.boundary.bls
    start = sim.x[1]
    stop = sim.x[end]
    typeof(bl1)<:AbstractComplexBL{1} ? start = start + (bl1.depth+8sim.dx) : nothing
    typeof(bl2)<:AbstractComplexBL{2} ? stop  = stop  - (bl2.depth+8sim.dx) : nothing
    start_ind = ceil(Int,Lattices.latticeindex(sim.lattice,start))
    stop_ind = floor(Int,Lattices.latticeindex(sim.lattice,stop))
    return VectorField{1,3}(sim.x,val,sim.lattice[start_ind],sim.lattice[stop_ind],[start_ind+1],[stop_ind+1])
end


text:
"""
    readpng(pathname)

Read a PNG file.

This returns a image object suitable for placing on the current drawing with `placeimage()`.
You can access its `width` and `height` fields:

    image = readpng("/tmp/test-image.png")
    w = image.width
    h = image.height
"""
function readpng(pathname)
    return Cairo.read_from_png(pathname)
end

function paint_with_alpha(ctx::Cairo.CairoContext, a = 0.5)
    Cairo.paint_with_alpha(get_current_cr(), a)
end

"""
    placeimage(img, xpos, ypos; centered=false)

Place a PNG image on the drawing at (`xpos`/`ypos`). The image `img` has been previously
loaded using `readpng()`.

Use keyword `centered=true` to place the center of the image at the position.
"""
function placeimage(img::Cairo.CairoSurface, xpos, ypos; centered=false)
    if centered == true
        w, h = img.width, img.height
        xpos, ypos = xpos - (w/2), ypos - (h/2)
    end
    Cairo.set_source_surface(get_current_cr(), img, xpos, ypos)
    # no alpha
    Cairo.paint(get_current_cr())
end

"""
    placeimage(img, pos; centered=false)

Place the top left corner of the PNG image on the drawing at `pos`.

Use keyword `centered=true` to place the center of the image at the position.
"""
placeimage(img::Cairo.CairoSurface, pt::Point; kwargs...) = placeimage(img, pt.x, pt.y; kwargs...)

"""
    placeimage(img, xpos, ypos, a; centered=false)

Place a PNG image on the drawing at (`xpos`/`ypos`) with transparency `a`.

Use keyword `centered=true` to place the center of the image at the position.
"""
function placeimage(img::Cairo.CairoSurface, xpos, ypos, alpha; centered=false)
    if centered == true
        w, h = img.width, img.height
        xpos, ypos = xpos - (w/2), ypos - (h/2)
    end
    Cairo.set_source_surface(get_current_cr(), img, xpos, ypos)
    paint_with_alpha(get_current_cr(), alpha)
end

"""
    placeimage(img, pos, a; centered=false)

Place a PNG image on the drawing at `pos` with transparency `a`.

Use keyword `centered=true` to place the center of the image at the position.
"""
placeimage(img::Cairo.CairoSurface, pt::Point, alpha; kwargs...) =
  placeimage(img::Cairo.CairoSurface, pt.x, pt.y, alpha; kwargs...)


text:
function _sortperm_unstable!(idx, x, ranges, last_valid_range, ord, a; threads = true)
    @_threadsfor threads for i in 1:last_valid_range
        rangestart = ranges[i]
        i == last_valid_range ? rangeend = length(x) : rangeend = ranges[i+1] - 1
        if (rangeend - rangestart) == 0
            continue
        end
        ds_sort!(x, idx, rangestart, rangeend, a, ord)
    end
end

function _sortperm_unstable_using_hpsort!(idx, x, ranges, last_valid_range, ord, a)
    for i in 1:last_valid_range
        rangestart = ranges[i]
        i == last_valid_range ? rangeend = length(x) : rangeend = ranges[i+1] - 1
        if (rangeend - rangestart) == 0
            continue
        end
        if rangeend - rangestart + 1 > Threads.nthreads()*2
            hp_ds_sort!(x, idx, a, ord; lo = rangestart, hi = rangeend)
        else
            ds_sort!(x, idx, rangestart, rangeend, a, ord)
        end
    end
end

# NOT OK
# we should find starts here
function fast_sortperm_int_threaded!(x, original_P, copy_P, ranges, rangelen, minval, misatleft, last_valid_range, ::Val{T}) where T
    starts = [T[] for i in 1:Threads.nthreads()]
    Threads.@threads for i in 1:last_valid_range
        rangestart = ranges[i]
        i == last_valid_range ? rangeend = length(x) : rangeend = ranges[i+1] - 1
        # if (rangeend - rangestart) == 0
        #     continue
        # end
        if misatleft
            _start_vals = _ds_sort_int_missatleft_nopermx!(x, original_P, copy_P, rangestart, rangeend, rangelen, minval, Val(T))
        else
            _start_vals = _ds_sort_int_missatright_nopermx!(x, original_P, copy_P, rangestart, rangeend, rangelen, minval, Val(T))
        end
        _cleanup_starts!(_start_vals, rangeend - rangestart + 1)
        append!(starts[Threads.threadid()], _start_vals .+ rangestart .- 1)
    end
    cnt = 1
    flag = false
    @inbounds for i in 1:length(starts)
        for j in 1:length(starts[i])
            ranges[cnt] = starts[i][j]
            if cnt > 1 && ranges[cnt] < ranges[cnt - 1]
                flag = true
            end
            cnt += 1
        end
    end
    flag && sort!(view(ranges, 1:cnt-1))
    return cnt - 1
end

function fast_sortperm_int_lm!(x, original_P, copy_P, ranges, rangelen, minval, misatleft, last_valid_range, ::Val{T}) where T
    starts = T[]
    for i in 1:last_valid_range
        rangestart = ranges[i]
        i == last_valid_range ? rangeend = length(x) : rangeend = ranges[i+1] - 1
        # if (rangeend - rangestart) == 0
        #     continue
        # end
        if misatleft
            _start_vals = _ds_sort_int_missatleft_nopermx_threaded_lm!(x, original_P, copy_P, rangestart, rangeend, rangelen, minval, Val(T))
        else
            _start_vals = _ds_sort_int_missatright_nopermx_threaded_lm!(x, original_P, copy_P, rangestart, rangeend, rangelen, minval, Val(T))
        end
        _cleanup_starts!(_start_vals, rangeend - rangestart + 1)
        append!(starts, _start_vals .+ rangestart .- 1)
    end
    @inbounds for i in 1:length(starts)
        ranges[i] = starts[i]
    end
    return length(starts)
end

function fast_sortperm_int!(x, original_P, copy_P, ranges, rangelen, minval, misatleft, last_valid_range, ::Val{T}) where T
    starts = T[]
    for i in 1:last_valid_range
        rangestart = ranges[i]
        i == last_valid_range ? rangeend = length(x) : rangeend = ranges[i+1] - 1
        # if (rangeend - rangestart) == 0
        #     continue
        # end
        if misatleft
            _start_vals = _ds_sort_int_missatleft_nopermx_threaded!(x, original_P, copy_P, rangestart, rangeend, rangelen, minval, Val(T))
        else
            _start_vals = _ds_sort_int_missatright_nopermx_threaded!(x, original_P, copy_P, rangestart, rangeend, rangelen, minval, Val(T))
        end
        _cleanup_starts!(_start_vals, rangeend - rangestart + 1)
        append!(starts, _start_vals .+ rangestart .- 1)
    end
    @inbounds for i in 1:length(starts)
        ranges[i] = starts[i]
    end
    return length(starts)
end

function _sortperm_int!(idx, idx_cpy, x, ranges, where, last_valid_range, missingatleft, ord, a; threads = true)
    @_threadsfor threads for i in 1:last_valid_range
        rangestart = ranges[i]
        i == last_valid_range ? rangeend = length(x) : rangeend = ranges[i+1] - 1
        if (rangeend - rangestart + 1) == 1
            continue
        end
        _minval = stat_minimum(x, lo = rangestart, hi = rangeend)
        if ismissing(_minval)
            continue
        else
            minval::Int = _minval
        end
        maxval::Int = stat_maximum(x, lo = rangestart, hi = rangeend)
        # the overflow is check before calling _sortperm_int!
        rangelen = maxval - minval + 1
        if rangelen < div(rangeend - rangestart + 1, 2)
            if missingatleft
                ds_sort_int_missatleft!(x, idx, idx_cpy, where[Threads.threadid()], rangestart, rangeend, rangelen, minval)
            else
                ds_sort_int_missatright!(x, idx, idx_cpy, where[Threads.threadid()], rangestart, rangeend, rangelen, minval)
            end
        else
            ds_sort!(x, idx, rangestart, rangeend, a, ord)
        end
    end
end

function _apply_by_f_barrier(x::AbstractVector{T}, by, rev, threads) where T
    needrev = rev
    missat = :right
    CT = Core.Compiler.return_type(_date_value∘by, Tuple{nonmissingtype(T)})
    if CT == Bool
        CT = Int8
    end
    CT = Union{Missing, CT}
    # _temp = Vector{CT}(undef, length(x))
    _temp = _our_vect_alloc(CT, length(x))
    # we should make sure changing sign doesn't overflow
    if rev && nonmissingtype(CT) <: Union{Bool, Int8, Int16, Int32, Int64} && isless(typemin(nonmissingtype(CT)), threads ? hp_minimum(_date_value∘by, x) : stat_minimum(_date_value∘by, x))
        _by = x-> -_date_value(by(x))
        needrev = false
        missat = :left
    else
        _by = x-> _date_value(by(x))
    end
    _temp, _by, needrev, missat
end

missatleftless(x, y) = isless(x, y)
missatleftless(::Missing, y) = true
missatleftless(x, ::Missing) = false
missatleftless(::Missing, ::Missing) = false

function _apply_by!(_temp, x::AbstractVector{T}, idx, _by, rev, needrev, missat, threads) where T
    @_threadsfor threads for j in 1:length(x)
        # if by(x) is Date or DateTime only grab its value
        @inbounds _temp[j] = _by(x[idx[j]])
    end
    # it is beneficial to check if the value of floats can be converted to Int
    # TODO rev=true should also be considered
    intable = false
    if !rev && eltype(_temp) <: Union{Missing, Float64}
        intable = true
        @_threadsfor threads for j in 1:length(x)
            @inbounds _is_intable(_temp[j]) && !isequal(_temp[j], -0.0) ? true : (intable = false && break)
        end
    end
    if !rev
        return (_temp, ord(isless, identity, rev, Forward), :right, needrev, intable)
    else
        if missat == :left
            return (_temp, ord(missatleftless, identity, needrev, Forward), missat, needrev, intable)
        else
            return (_temp, ord(isless, identity, needrev, Forward), missat, needrev, intable)
        end
    end
end

function _apply_by(x::AbstractVector, idx, by, rev; threads = true)
    _temp, _by, needrev, missat = _apply_by_f_barrier(x, by, rev, threads)
    _apply_by!(_temp, x, idx, _by, rev, needrev, missat, threads)
end

function _modify_poolarray_to_integer!(refs, trans; threads = true)
    @inbounds @_threadsfor threads for i in 1:length(refs)
        refs[i] = trans[refs[i]]
    end
end

function _fast_path_modify_to_integer!(xrefs, perm; threads = true)
    @inbounds @_threadsfor threads for i in 1:length(xrefs)
        xrefs[i] = perm[xrefs[i]]
    end
end

function _fill_idx_for_sort!(idx; threads = true)
    @inbounds @_threadsfor threads for i in 1:length(idx)
        idx[i] = i
    end
end

function _check_memory_availability_for_hp_sort(x)
    return 4*sizeof(x) < Base.Sys.free_memory()
end

function _fill_ranges_for_fast_int_sort!(ranges, _starts_vals)
    for kk in 1:length(_starts_vals)
        ranges[kk] = _starts_vals[kk] # TODO we don't need ranges here
    end
end
# T is either Int32 or Int64 based on how many rows ds has
function ds_sort_perm(ds, colsidx, by::Vector{<:Function}, rev::Vector{Bool}, a::Base.Sort.Algorithm,  ::Val{T}; skipcol = 0, skipcol_mkcopy = true, notsortpaforjoin = false, givenrange = nothing, threads = true) where T
    @assert length(colsidx) == length(by) == length(rev) "each col should have all information about lt, by, and rev"
    stable = false
    # arrary to keep the permutation of rows
    if givenrange === nothing
        if skipcol>0 && !skipcol_mkcopy
            # make use of ds.perm space
            idx = index(ds).perm
        else
            idx = Vector{T}(undef, nrow(ds))
        end
        _fill_idx_for_sort!(idx; threads = threads)
    else
        givenrange.idx isa UnitRange ? idx = collect(givenrange.idx) : idx = givenrange.idx
    end

    if givenrange === nothing
        if skipcol>0
            last_valid_range = index(ds).ngroups[]
        else
            last_valid_range::T = 1
        end
    else
        last_valid_range = givenrange.lastvalid
    end
    # ranges keep the starts of blocks of sorted rows
    # rangescpy is a copy which will be help to rearrange ranges in place

    if givenrange === nothing
        if skipcol>0
            if skipcol_mkcopy
                ranges = resize!(copy(index(ds).starts), nrow(ds))
            else # make use of ds.starts
                ranges = resize!(index(ds).starts, nrow(ds))
            end
        else
            ranges = Vector{T}(undef, nrow(ds))
            ranges[1] = 1
        end
    else
        ranges = resize!(givenrange.starts, nrow(ds))
    end

    inbits = zeros(Bool, nrow(ds))

    # rangescpy[1] = 1
    # in case we have integer columns with few distinct values
    int_permcpy = T[]
    int_where = [T[] for _ in 1:Threads.nthreads()]
    int_count = [T[] for _ in 1:Threads.nthreads()]
    for i in 1:length(colsidx)
        if i <= skipcol
            continue
        end
        x = _columns(ds)[colsidx[i]]
        # pooledarray are treating differently (or general case of poolable data)
        if DataAPI.refpool(x) !== nothing
            if by[i] == identity
                _tmp = copy(x)
            else
                _tmp = map(by[i], x)
            end
            if givenrange === nothing
                if i > 1
                    _tmp.refs .= _threaded_permute(_tmp.refs, idx; threads = threads)
                end
            else
                _tmp.refs .= _threaded_permute(_tmp.refs, idx; threads = threads)
            end
            # collect here is to handle Categorical array.
            if !notsortpaforjoin
                if _tmp isa PooledArray
                    _fast_path_modify_to_integer!(_tmp.refs, invperm(sortperm(Dataset(x = DataAPI.refpool(_tmp), copycols = false), :, rev = rev[i], stable = false, alg = a)); threads = threads)
                else
                    # This path is for Categorical array and must be optimised
                    levs = allowmissing(DataAPI.levels(_tmp))
                    push!(levs, missing)
                    rev[i] ? reverse!(levs) : nothing
                    aaa = map(x->get(DataAPI.invrefpool(_tmp), x, missing), levs)
                    trans = Dict{eltype(aaa), Int32}(aaa .=> 1:length(aaa))
                    _modify_poolarray_to_integer!(_tmp.refs, trans; threads = threads)
                end
                _tmp = _tmp.refs
            else # this path is used for sorting (not correct sorting) refs which will be used in joining functions
                _tmp = _tmp.refs
            end
            _ordr = ord(isless, identity, false, Forward)
            _needrev = false
            intable = false
            _missat = :right
        else
            _tmp, _ordr, _missat, _needrev, intable=  _apply_by(x, idx, by[i], rev[i]; threads = threads)

        end
        if !_needrev && (eltype(_tmp) <: Union{Missing, INTEGERS} || intable) && eltype(_tmp) !== Missing
            # further check for fast integer sort
            n = length(_tmp)
            if n > 1
                if threads
                    _minval = hp_minimum(_tmp)
                else
                    _minval = stat_minimum(_tmp)
                end
                if ismissing(_minval)
                    continue
                else
                    minval::Integer = _minval
                end
                if threads
                    maxval::Integer = hp_maximum(_tmp)
                else
                    maxval = stat_maximum(_tmp)
                end
                rnglen = BigInt(maxval) - BigInt(minval) + 1
                o1 = false
                if rnglen < typemax(Int)
                    o1 = true
                    rangelen = Int(rnglen)
                end
                if o1 && maxval < typemax(Int) && (rangelen <= div(n,2))
                    if last_valid_range == 1
                      #for many levels but not too many we use parallel version
                      if threads && Threads.nthreads() > 1 && nrow(ds) > Threads.nthreads() && rangelen > 100_000 && rangelen*Threads.nthreads() < n

                          if _missat == :left
                              _starts_vals = _ds_sort_int_missatleft_nopermx_threaded!(_tmp, idx, rangelen, minval, Val(T))
                          else
                              _starts_vals = _ds_sort_int_missatright_nopermx_threaded!(_tmp, idx, rangelen, minval, Val(T))
                          end
                          _cleanup_starts!(_starts_vals, n)
                          last_valid_range = length(_starts_vals)
                          _fill_ranges_for_fast_int_sort!(ranges, _starts_vals)
                          if length(colsidx) == 1
                              stable = true
                          end
                          continue
                      elseif threads && Threads.nthreads() > 1 && nrow(ds) > Threads.nthreads() && rangelen > 100_000
                          if _missat == :left
                              _starts_vals = _ds_sort_int_missatleft_nopermx_threaded_lm!(_tmp, idx, rangelen, minval, Val(T))
                          else
                              _starts_vals = _ds_sort_int_missatright_nopermx_threaded_lm!(_tmp, idx, rangelen, minval, Val(T))
                          end
                          _cleanup_starts!(_starts_vals, n)
                          last_valid_range = length(_starts_vals)
                          _fill_ranges_for_fast_int_sort!(ranges, _starts_vals)
                          if length(colsidx) == 1
                              stable = true
                          end
                          continue
                      else

                          if _missat == :left
                              _starts_vals = _ds_sort_int_missatleft_nopermx!(_tmp, idx, rangelen, minval, Val(T))
                          else
                              _starts_vals = _ds_sort_int_missatright_nopermx!(_tmp, idx, rangelen, minval, Val(T))
                          end
                          _cleanup_starts!(_starts_vals, n)
                          last_valid_range = length(_starts_vals)
                          _fill_ranges_for_fast_int_sort!(ranges, _starts_vals)
                          if length(colsidx) == 1
                              stable = true
                          end
                          continue
                      end
                  else
                      resize!(int_permcpy, length(idx))
                      copy!(int_permcpy, idx)
                      if threads && last_valid_range < Threads.nthreads() && rangelen > 100_000

                          last_valid_range = fast_sortperm_int_lm!(_tmp, idx, int_permcpy, ranges, rangelen, minval, _missat == :left, last_valid_range, Val(T))
                          continue
                      end
                      # we should check if there are many observations in each group/ this version of int sort allocate where for every iteration
                      if threads && n/last_valid_range > max(2.0*rangelen,10^6) && rangelen*Threads.nthreads() < n/2 # if rangelen is not that much

                          last_valid_range = fast_sortperm_int_threaded!(_tmp, idx, int_permcpy, ranges, rangelen, minval, _missat == :left, last_valid_range, Val(T))
                          continue
                      else
                          for nt in 1:(threads ? Threads.nthreads() : 1)
                              resize!(int_where[nt], rangelen + 2)
                          end

                          _sortperm_int!(idx, int_permcpy, _tmp, ranges, int_where, last_valid_range, _missat == :left, _ordr, a; threads = threads)
                      end
                  end
                else
                    if threads && last_valid_range == 1 && Threads.nthreads() > 1 && nrow(ds) > Threads.nthreads()

                        hp_ds_sort!(_tmp, idx, a, _ordr)
                    else
                        if threads && last_valid_range < floor(Int, log2(Threads.nthreads()))/2
                            _sortperm_unstable_using_hpsort!(idx, _tmp, ranges, last_valid_range, _ordr, a)
                        else
                            _sortperm_unstable!(idx, _tmp, ranges, last_valid_range, _ordr, a; threads = threads)
                        end
                    end
                end
            end
        else
            if threads && last_valid_range == 1 && Threads.nthreads() > 1 && nrow(ds) > Threads.nthreads()
                hp_ds_sort!(_tmp, idx, a, _ordr)
            else
                # log2(thcnt) is used because of the current parallel sorting algorithm
                if threads && last_valid_range < floor(Int, log2(Threads.nthreads()))/2
                    _sortperm_unstable_using_hpsort!(idx, _tmp, ranges, last_valid_range, _ordr, a)
                else
                    _sortperm_unstable!(idx, _tmp, ranges, last_valid_range, _ordr, a; threads = threads)
                end
            end
        end
        # last_valid_range = _fill_starts!(ranges, _tmp, rangescpy, last_valid_range, _ordr, Val(T))
        last_valid_range = _fill_starts_v2!(ranges, inbits, _tmp, last_valid_range, _ordr, Val(T); threads = threads)
        last_valid_range == nrow(ds) && return (ranges, idx, last_valid_range, stable)
    end
    return (ranges, idx, last_valid_range, stable)
end

function _stablise_sort!(ranges, idx, last_valid_range, a; threads = true)
    @_threadsfor threads for i in 1:last_valid_range
        rangestart = ranges[i]
        i == last_valid_range ? rangeend = length(idx) : rangeend = ranges[i+1] - 1
        if (rangeend - rangestart) == 0
            continue
        end
        # if QuickSort is selected we make sure the worst case scenario is not that much bad
        if rangeend - rangestart + 1 > 1000
            sort!(idx, rangestart, rangeend, MergeSort, Forward)
        else
            sort!(idx, rangestart, rangeend, a, Forward)
        end
    end
end

function _sortperm(ds, cols::MultiColumnIndex, rev::Vector{Bool}; a = HeapSortAlg(), mapformats = true, stable = true, skipcol = 0, skipcol_mkcopy = true, notsortpaforjoin = false, givenrange = nothing, threads = true)
    colsidx = index(parent(ds))[cols]
    @assert length(colsidx) == length(rev) "`rev` argument must be the same as length of selected columns"
    if _check_for_fast_sort(ds, colsidx, rev, mapformats; notsortpaforjoin = notsortpaforjoin, givenrange = givenrange) == 0
        return copy(_group_starts(ds)), 1:nrow(ds), _ngroups(ds)
    end
    if _check_for_fast_sort(ds, colsidx, rev, mapformats; notsortpaforjoin = notsortpaforjoin, givenrange = givenrange) == 1
        colsidx, ranges, last_valid_index = _find_starts_of_groups(ds, colsidx, nrow(ds) < typemax(Int32) ? Val(Int32) : Val(Int64), mapformats = mapformats, threads = threads)
        if !skipcol_mkcopy
            _fill_idx_for_sort!(_get_sort_perms(ds); threads = threads)
            # make use of perm in original data (simply replace it with new one)
            return ranges, _get_sort_perms(ds), last_valid_index
        else
            return ranges, 1:nrow(ds), last_valid_index
        end
    end
    if skipcol == 0 && _check_for_fast_sort(ds, colsidx, rev, mapformats; notsortpaforjoin = notsortpaforjoin, givenrange = givenrange) == 2
        skipcol = length(_groupcols(ds))
    end
    if skipcol == -1 && _check_for_fast_sort(ds, colsidx, rev, mapformats; notsortpaforjoin = notsortpaforjoin) == 2
        skipcol = length(_groupcols(ds))
    end
    by = Function[]

    # TODO is it safe to refer to parents formats?
    if mapformats
        for j in 1:length(colsidx)
            push!(by, getformat(parent(ds), colsidx[j]))
        end
    else
        for j in 1:length(colsidx)
            push!(by, identity)
        end
    end
    ranges, idx, last_valid_range, isstable = ds_sort_perm(parent(ds), colsidx, by, rev, a, nrow(ds) < typemax(Int32) ? Val(Int32) : Val(Int64); skipcol = skipcol, skipcol_mkcopy = skipcol_mkcopy, notsortpaforjoin = notsortpaforjoin, givenrange = givenrange, threads = threads)
    if stable && !isstable
        if length(idx) == last_valid_range
            return ranges, idx, last_valid_range
        else
            _stablise_sort!(ranges, idx, last_valid_range, a, threads = threads)
            return ranges, idx, last_valid_range
        end
    else
        return ranges, idx, last_valid_range
    end
end

function _sortperm(ds, cols::MultiColumnIndex, rev::Bool = false; a = HeapSortAlg(), mapformats = true, stable = true, skipcol = 0, skipcol_mkcopy = true, notsortpaforjoin = false, givenrange = nothing, threads = true)
    colsidx = index(parent(ds))[cols]
    revs = repeat([rev], length(colsidx))
    if _check_for_fast_sort(ds, colsidx, revs, mapformats; notsortpaforjoin = notsortpaforjoin, givenrange = givenrange) == 0
        return copy(_group_starts(ds)), 1:nrow(ds), _ngroups(ds)
    end
    if _check_for_fast_sort(ds, colsidx, revs, mapformats; notsortpaforjoin = notsortpaforjoin, givenrange = givenrange) == 1
        colsidx, ranges, last_valid_index = _find_starts_of_groups(ds, colsidx, nrow(ds) < typemax(Int32) ? Val(Int32) : Val(Int64), mapformats = mapformats, threads = threads)

        if !skipcol_mkcopy
            _fill_idx_for_sort!(_get_sort_perms(ds); threads = threads)
            # make use of perm in original data (simply replace it with new one)
            return ranges, _get_sort_perms(ds), last_valid_index
        else
            return ranges, 1:nrow(ds), last_valid_index
        end
    end
    if skipcol == 0 && _check_for_fast_sort(ds, colsidx, revs, mapformats; notsortpaforjoin = notsortpaforjoin, givenrange = givenrange) == 2
        skipcol = length(_groupcols(ds))
    end
    if skipcol == -1 && _check_for_fast_sort(ds, colsidx, revs, mapformats; notsortpaforjoin = notsortpaforjoin) == 2
        skipcol = length(_groupcols(ds))
    end
    _sortperm(parent(ds), cols, revs; a = a , mapformats = mapformats, stable = stable, skipcol = skipcol, skipcol_mkcopy = skipcol_mkcopy, notsortpaforjoin = notsortpaforjoin, givenrange = givenrange, threads = threads)
end

_sortperm(ds, col::ColumnIndex, rev::Bool = false; a = HeapSortAlg(), mapformats = true, stable = true, skipcol = 0, skipcol_mkcopy = true, notsortpaforjoin = false, givenrange = nothing, threads = true) = _sortperm(ds, [col], rev; a = a,  mapformats = mapformats, stable = stable, skipcol = skipcol, skipcol_mkcopy = skipcol_mkcopy, notsortpaforjoin = notsortpaforjoin, givenrange = givenrange, threads = threads)

function _check_for_fast_sort(ds, colsidx, rev, mapformats; notsortpaforjoin = false, givenrange = nothing)
    scols = _sortedcols(ds)
    revs = _get_rev(ds)
    if revs isa Bool
        revs = repeat([revs], length(scols))
    end
    fmt = _get_fmt(ds)
    notsortpaforjoin && return -1
    givenrange !== nothing && return -1
    if colsidx == scols && rev == revs && mapformats == fmt
        return 0
    elseif length(colsidx) < length(scols)
        if colsidx == view(scols, 1:length(colsidx)) && rev == view(revs, 1:length(colsidx)) && mapformats == fmt
            return 1
        end
    elseif length(colsidx) > length(scols)
        if view(colsidx, 1:length(scols)) == scols && view(rev, 1:length(scols)) == revs && mapformats == fmt
            return 2
        end
    else
        return -1
    end
end

function _sortperm_v(ds::SubDataset, cols::MultiColumnIndex, rev = false; a = HeapSortAlg(), mapformats = true, stable = true, skipcol = 0, skipcol_mkcopy = true, notsortpaforjoin = false, givenrange = nothing, threads = true)
    colsidx = index(ds)[cols]
    if rev isa AbstractVector
        @assert length(rev) == length(colsidx) "length of rev and the number of selected columns must match"
        revs = rev
    else
        revs = repeat([rev], length(colsidx))
    end
    by = Function[]

    if mapformats
        for j in 1:length(colsidx)
            push!(by, getformat(ds, colsidx[j]))
        end
    else
        for j in 1:length(colsidx)
            push!(by, identity)
        end
    end
    ranges, idx, last_valid_range, isstable = ds_sort_perm(ds, colsidx, by, revs, a, nrow(ds) < typemax(Int32) ? Val(Int32) : Val(Int64); skipcol = skipcol, skipcol_mkcopy = skipcol_mkcopy, notsortpaforjoin = notsortpaforjoin, givenrange = givenrange, threads = threads)
    if stable && !isstable
        if length(idx) == last_valid_range
            return ranges, idx, last_valid_range
        else
            _stablise_sort!(ranges, idx, last_valid_range, a; threads = threads)
            return ranges, idx, last_valid_range
        end
    else
        return ranges, idx, last_valid_range
    end
end
_sortperm_v(ds::SubDataset, col::ColumnIndex, rev = false; a = HeapSortAlg(), mapformats = true, stable = true, skipcol = 0, skipcol_mkcopy = true, notsortpaforjoin = false, givenrange = nothing, threads = true) = _sortperm_v(ds, [col], rev, a = a, mapformats = mapformats, stable = stable, skipcol = skipcol, skipcol_mkcopy = skipcol_mkcopy, notsortpaforjoin = notsortpaforjoin, givenrange = givenrange, threads = threads)


text:
using Documenter
using StateSpace

makedocs(
    sitename = "StateSpace",
    format = Documenter.HTML(),
    modules = [StateSpace],
    pages = ["Documentation" => "index.md",
             "Interface" => "interface.md",
             "Examples" => "examples.md" ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo   = "github.com/npsm/StateSpace.jl.git",
)


text:
#
# Do work here ...
#

# This is a 1x3 Matrix{Float64}
M = [1.1 2.3 4.5]

# This is a 3-element Vector{Float64}
v = [1.1; 2.3; 4.5]
println("Julia side: I will return a Vector{Float64}: ", v)

return v

max_stars_repo_name: yundantianchang/dynload-julia

max_issues_repo_name: yundantianchang/dynload-julia

max_forks_repo_name: yundantianchang/dynload-julia


text:
using MathML, EzXML, Symbolics, SpecialFunctions, IfElse, AbstractTrees, Test

# these tests mimic the README examples
fn = "data/math.xml"
doc = readxml(fn)
num = MathML.parse_node(doc.root)
@variables compartment k1 S1
true_num = prod([compartment, k1, S1])
@test isequal(num, [true_num])

fn = "data/eq.xml"
doc = readxml(fn)
eqs = MathML.parse_node(doc.root)
@variables T I Par_94 Par_90
true_eqs = [
    T ~ Par_94,
    I ~ Par_90,
]
@test isequal(eqs, true_eqs)

@variables w x y z a b t

str = "<apply><compose/><ci>x</ci><ci>y</ci><ci>z</ci></apply>"
@test isequal(MathML.parse_str(str), x ∘ y ∘ z)

str = """<cn type="rational">22<sep/>7</cn>"""
@test isequal(MathML.parse_str(str), 22 // 7)

str = """<cn type="e-notation">5<sep/>2</cn>"""
@test isequal(MathML.parse_str(str), 500)

str = """<cn type="complex-polar"> 2 <sep/> 3.1415 </cn>"""
@test isapprox(MathML.parse_str(str), Complex(-2, 0), atol=1e-3)

str = """<cn type="complex-cartesian"> 12.3 <sep/> 5 </cn>"""
@test isequal(MathML.parse_str(str), Complex(12.3, 5))

# heaviside and nesting test
str = """
<piecewise>
  <piece>
    <apply><minus/><ci>x</ci></apply>
    <apply><lt/><ci>x</ci><cn>0</cn></apply>
  </piece>
  <piece>
    <cn>0</cn>
    <apply><eq/><ci>x</ci><cn>0</cn></apply>
  </piece>
  <piece>
    <ci>x</ci>
    <apply><gt/><ci>x</ci><cn>0</cn></apply>
  </piece>
</piecewise>
"""
# MathML.parse_str(str) # no <otherwise>

str = """
<piecewise>
<piece>
  <apply>
      <divide/>
      <apply>
        <times/>
        <ci>x</ci>
        <apply>
            <plus/>
            <ci>y</ci>
            <apply>
              <times/>
              <ci>z</ci>
              <ci>a</ci>
            </apply>
        </apply>
      </apply>
      <apply>
        <minus/>
        <cn>1</cn>
        <apply>
            <times/>
            <ci>z</ci>
            <ci>b</ci>
        </apply>
      </apply>
  </apply>
  <apply>
      <leq/>
      <ci>t</ci>
      <cn>1</cn>
  </apply>
</piece>
<otherwise>
  <apply>
      <times/>
      <ci>x</ci>
      <ci>y</ci>
  </apply>
</otherwise>
</piecewise>
"""
@test isequal(MathML.parse_str(str),
  IfElse.ifelse(
    IfElse.ifelse(1. - t >= 0, 1, 0) > 0.5,
     x * (y + a * z) * ((1.0 - (b * z))^-1),
      x * y))

# factorial
str = "<apply><factorial/><ci>x</ci></apply>"
@test isequal(MathML.parse_str(str), factorial(x))

str = "<apply><factorial/><cn>5</cn></apply>"
@test MathML.parse_str(str) == 120

# root
str = """
<apply><root/>
  <degree><ci type="integer">y</ci></degree>
  <ci>x</ci>
</apply>
"""
@test isequal(MathML.parse_str(str), x^y)

str = "<apply><root/><ci>x</ci></apply>"
@test isequal(MathML.parse_str(str), sqrt(x))

#
str = """
<apply><divide/>
<ci>x</ci>
<ci>y</ci>
</apply>
"""
@test isequal(MathML.parse_str(str), x / y)

#
str = "<apply><max/><cn>2</cn><cn>3</cn><cn>5</cn></apply>"
@test isequal(MathML.parse_str(str), 5)

#
str = "<apply><min/><ci>x</ci><ci>y</ci></apply>"
@test isequal(MathML.parse_str(str), min(x, y))

#
str = "<apply><minus/><cn>3</cn></apply>"
@test isequal(MathML.parse_str(str), -3)

#
str = "<apply><minus/><ci>x</ci><ci>y</ci></apply>"
@test isequal(MathML.parse_str(str), x - y)

#
str = "<apply><plus/><ci>x</ci><ci>y</ci><ci>z</ci></apply>"
@test isequal(MathML.parse_str(str), x + y + z)

#
str = "<apply><rem/><ci>x</ci><ci>y</ci></apply>"
@test isequal(MathML.parse_str(str), rem(x, y))

#
# str = "<apply><gcd/><ci>a</ci><ci>b</ci><ci>c</ci></apply>" # gcd is a ways off, gröbner bases
# MathML.parse_str(str)
# @test isequal(MathML.parse_str(str), )

#
str = "<apply><abs/><ci>x</ci></apply>"
@test isequal(MathML.parse_str(str), abs(x))

#
# str = """
# <apply><conjugate/>
#   <apply><plus/>
#     <ci>x</ci>
#     <apply><times/><cn>ⅈ<!--DOUBLE-STRUCK ITALIC SMALL I--></cn><ci>y</ci></apply>
#   </apply>
# </apply>
# """
# MathML.parse_str(str)
# @test isequal(MathML.parse_str(str), )

#  lcm
# str = "<apply><lcm/><ci>a</ci><ci>b</ci><ci>c</ci></apply>"
# MathML.parse_str(str)
# @test isequal(MathML.parse_str(str), )

# floor and ceil - RoundingMode
# str = "<apply><floor/><ci>a</ci></apply>"
# MathML.parse_str(str)
# @test isequal(MathML.parse_str(str), )

#  do we actually want it to evaluate, or maybe just return expr
str = "<apply><gt/><cn>3</cn><cn>2</cn></apply>"
@test isequal(MathML.parse_str(str), true)

#
str = "<apply><lt/><cn>2</cn><cn>3</cn><cn>4</cn></apply>"
@test isequal(MathML.parse_str(str), true)

#
str = "<apply><exp/><ci>x</ci></apply>"
@test isequal(MathML.parse_str(str), exp(x))

#
str = "<apply><ln/><ci>x</ci></apply>"
@test isequal(MathML.parse_str(str), log(x))

#  dopee
str = """
<vector>
  <apply><plus/><ci>x</ci><ci>y</ci></apply>
  <cn>3</cn>
  <cn>7</cn>
</vector>
"""
@test isequal(MathML.parse_str(str), [x + y, 3, 7])

str = "<bvar><ci>x</ci></bvar>"
@test isequal(MathML.parse_str(str), (Num(Variable(:x)), 1))

str = "<bvar><ci>x</ci><degree><cn>2</cn></degree></bvar>"
@test isequal(MathML.parse_str(str), (Num(Variable(:x)), 2))

str = """
<apply><diff/>
  <bvar><ci>x</ci></bvar>
  <apply><sin/><ci>x</ci></apply>
</apply>
"""
@test isequal(expand_derivatives(MathML.parse_str(str)), cos(Num(Variable(:x))))

str = """
<apply><diff/>
  <bvar><ci>x</ci><degree><cn>2</cn></degree></bvar>
  <apply><power/><ci>x</ci><cn>4</cn></apply>
</apply>
"""
@test isequal(expand_derivatives(MathML.parse_str(str)), 12 * Num(Variable(:x))^2)

# macro test
ml = MathML"""
<apply><diff/>
  <bvar><ci>x</ci><degree><cn>2</cn></degree></bvar>
  <apply><power/><ci>x</ci><cn>4</cn></apply>
</apply>
"""
@test isequal(expand_derivatives(ml), 12 * Num(Variable(:x))^2)

str = """
  <lambda>
    <bvar>
      <ci> x </ci>
    </bvar>
    <bvar>
      <ci> y </ci>
    </bvar>
    <apply>
      <power/>
      <ci> x </ci>
      <ci> y </ci>
    </apply>
  </lambda>
"""
f = parse_str(str)
xml = parsexml(str).root
g = parse_lambda(xml)
@test g(3,5) == f(3,5) == [243]
@test isequal(f(x, y), [x^y])

# parse_apply with <ci> as firstelement problem
@test_throws KeyError parse_file("data/err.xml")


text:
using Test
using GeometryTypes
using SurfaceTopology

@info "Topology function tests"

faces =   Face{3,Int}[
    [6,   8,  11],
    [8,   6,   7],
    [5,   1,   4],
    [1,   5,   7],
    [5,   8,   7],
    [5,  10,  11],
    [8,   5,  11],
    [1,   3,   2],
    [3,   1,   7],
    [3,   6,   2],
    [6,   3,   7],
    [9,   5,   4],
    [5,  12,  10],
    [9,  12,   5],
    [10,  12,   4],
    [12,   9,   4]
]

triangles = []
for i in FaceRing(5,faces)
    push!(triangles,i)
end

@test sort(triangles)==[3,4,5,6,7,12,13,14]

triverticies = []
for i in EdgeRing(5,faces)
    #println("i")
    push!(triverticies,i)
end

@test (1,4) in triverticies

verticies = []
for i in VertexRing(5,faces)
    push!(verticies,i)
end

@test sort(verticies)==[1,4,7,8,9,10,11,12] #[1,2,6,7]

@info "Testing FaceDS"

# points = zero(faces)
fb = FaceDS(faces)

triangles = []
for i in FaceRing(5,fb)
    push!(triangles,i)
end
@test sort(triangles)==[3,4,5,6,7,12,13,14]

triverticies = []
for i in EdgeRing(5,fb)
    push!(triverticies,i)
end

@test Face(1,4) in triverticies

verticies = []
for i in VertexRing(5,fb)
     push!(verticies,i)
end

@test sort(verticies)==[1,4,7,8,9,10,11,12]

@info "Testing EdgeDS"

eb = EdgeDS(faces)

verticies = []
for i in VertexRing(5,eb)
     push!(verticies,i)
end

@test sort(verticies)==[1,4,7,8,9,10,11,12]

triverticies = []
for i in EdgeRing(5,eb)
    push!(triverticies,i)
end

@test Face(1,4) in triverticies

@info "Topology tests for Cached DS"

# At the moment limited to a closed surfaces

faces = Face{3,Int64}[
    [1, 12, 6], [1, 6, 2], [1, 2, 8], [1, 8, 11], [1, 11, 12], [2, 6, 10], [6, 12, 5],
    [12, 11, 3], [11, 8, 7], [8, 2, 9], [4, 10, 5], [4, 5, 3], [4, 3, 7], [4, 7, 9],
    [4, 9, 10], [5, 10, 6], [3, 5, 12], [7, 3, 11], [9, 7, 8], [10, 9, 2]
]

cds = CachedDS(faces)

@test sort(collect(VertexRing(3,faces)))==sort(collect(VertexRing(3,cds)))

@test sort(collect(EdgeRing(3,faces)))==sort(collect(EdgeRing(3,cds)))


text:
# U.S. government securities/Treasury constant maturities/Nominal
function readUSTreasuries()

    ## set url
    url = "https://www.federalreserve.gov/datadownload/Output.aspx?rel=H15&series=bf17364827e38702b42a58cf8eaa3f78&lastobs=&from=&to=&filetype=csv&label=include&layout=seriescolumn&type=package"

    ## download and format data
    filename = download(url)
    data = readtable(filename, skipstart = 5,
                     nastrings = ["", "ND"])
    dats = Date[Date(data[ii, 1]) for ii=1:size(data, 1)]

    nObs, nCols = size(data)
    vals = zeros(Float64, nObs, nCols - 1)
    for ii=1:(nCols-1)
        xx = data[:, ii+1]
        thisvals = xx.data
        thisvals[xx.na] = NaN
        vals[:, ii] = thisvals
    end

    nams = names(data[:, 2:end])
    nams = [String(thisNam) for thisNam in nams]
    intRates = TimeSeries.TimeArray(dats, vals, nams)

    ## get metadata
    metaData = readtable(filename, nrows = 6, header = false,
                         nastrings = ["", "ND"])
    metaDf = EconDatasets.formatMetaData(metaData)

    return (intRates, metaDf)
end

function convertFedDate(dat::Date)
    yyyy = String(Dates.year(dat))
    m = String(Dates.month(dat))
    if length(m) == 1
        m = String("0", m)
    end
    d = String(Dates.day(dat))
    if length(d) == 1
        d = String("0", d)
    end
    return String(m, "/", d, "/", yyyy)
end

function formatMetaData(metaData::DataFrame)
    nVars = size(metaData, 2)-1
    nMeta = size(metaData, 1)
    varNames = Symbol[Symbol(metaData[ii, 1]) for ii=1:nMeta]
    df = DataFrame()
    for ii=nMeta:-1:1
        vals = DataArray(typeof(metaData[ii, 2]), nVars)
        for jj=1:nVars
            vals[jj] = metaData[ii, jj+1]
        end
        df[varNames[ii]] = vals
    end
    return df
end