text:
using BinaryBuilder, Pkg

name = "GLPK"
version = v"5.0"

# Collection of sources required to build GLPK
sources = [
    ArchiveSource("http://ftpmirror.gnu.org/gnu/glpk/glpk-$(version.major).$(version.minor).tar.gz",
                  "4a1013eebb50f728fc601bdd833b0b2870333c3b3e5a816eeba921d95bec6f15"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/glpk*
if [[ ${target} == *mingw* ]]; then
    export CPPFLAGS="-I${prefix}/include -D__WOE__=1"
else
    export CPPFLAGS="-I${prefix}/include"
fi
autoreconf -vi
./configure --prefix=${prefix} --host=${target} --build=${MACHTYPE} --with-gmp
make -j${nproc}
make install
"""

# Build for all platforms
platforms = supported_platforms()

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("GMP_jll", v"6.1.2"),
]

# Build the tarballs, and possibly a `build.jl` as well.
# Use the same preferred_gcc_version as GMP.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"6")

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
import REPL:ends_with_semicolon
import .Configuration
import .ExpressionExplorer: FunctionNameSignaturePair, is_joined_funcname, UsingsImports, external_package_names

Base.push!(x::Set{Cell}) = x

"Run given cells and all the cells that depend on them, based on the topology information before and after the changes."
function run_reactive!(session::ServerSession, notebook::Notebook, old_topology::NotebookTopology, new_topology::NotebookTopology, cells::Vector{Cell}; deletion_hook::Function=WorkspaceManager.delete_vars, persist_js_state::Bool=false)::TopologicalOrder
	# make sure that we're the only `run_reactive!` being executed - like a semaphor
	take!(notebook.executetoken)

	removed_cells = setdiff(keys(old_topology.nodes), keys(new_topology.nodes))
	cells = Cell[cells..., removed_cells...]

	# by setting the reactive node and expression caches of deleted cells to "empty", we are essentially pretending that those cells still exist, but now have empty code. this makes our algorithm simpler.
	new_topology = NotebookTopology(
		nodes=merge(
			new_topology.nodes,
			Dict(cell => ReactiveNode() for cell in removed_cells),
		),
		codes=merge(
			new_topology.codes,
			Dict(cell => ExprAnalysisCache() for cell in removed_cells)
		)
	)

	# save the old topological order - we'll delete variables assigned from it and re-evalutate its cells
	old_order = topological_order(notebook, old_topology, cells)

	old_runnable = old_order.runnable
	to_delete_vars = union!(Set{Symbol}(), defined_variables(old_topology, old_runnable)...)
	to_delete_funcs = union!(Set{Tuple{UUID,FunctionName}}(), defined_functions(old_topology, old_runnable)...)

	# get the new topological order
	new_order = topological_order(notebook, new_topology, union(cells, keys(old_order.errable)))
	to_run = setdiff(union(new_order.runnable, old_order.runnable), keys(new_order.errable))::Vector{Cell} # TODO: think if old error cell order matters

	# change the bar on the sides of cells to "queued"
	for cell in to_run
		cell.queued = true
	end
	for (cell, error) in new_order.errable
		cell.running = false
		cell.queued = false
		relay_reactivity_error!(cell, error)
	end

	# Send intermediate updates to the clients at most 20 times / second during a reactive run. (The effective speed of a slider is still unbounded, because the last update is not throttled.)
	# flush_send_notebook_changes_throttled,
	send_notebook_changes_throttled, flush_notebook_changes = throttled(1.0 / 20) do
		send_notebook_changes!(ClientRequest(session=session, notebook=notebook))
	end
	send_notebook_changes_throttled()

	# delete new variables that will be defined by a cell
	new_runnable = new_order.runnable
	to_delete_vars = union!(to_delete_vars, defined_variables(new_topology, new_runnable)...)
	to_delete_funcs = union!(to_delete_funcs, defined_functions(new_topology, new_runnable)...)

	# delete new variables in case a cell errors (then the later cells show an UndefVarError)
	new_errable = keys(new_order.errable)
	to_delete_vars = union!(to_delete_vars, defined_variables(new_topology, new_errable)...)
	to_delete_funcs = union!(to_delete_funcs, defined_functions(new_topology, new_errable)...)

	to_reimport = union(Set{Expr}(), map(c -> new_topology.codes[c].module_usings_imports.usings, setdiff(notebook.cells, to_run))...)

	deletion_hook((session, notebook), to_delete_vars, to_delete_funcs, to_reimport; to_run=to_run) # `deletion_hook` defaults to `WorkspaceManager.delete_vars`

	delete!.([notebook.bonds], to_delete_vars)

	local any_interrupted = false
	for (i, cell) in enumerate(to_run)

		cell.queued = false
		cell.running = true
		send_notebook_changes_throttled()

		if any_interrupted || notebook.wants_to_interrupt
			relay_reactivity_error!(cell, InterruptException())
		else
			run = run_single!(
				(session, notebook), cell,
				new_topology.nodes[cell], new_topology.codes[cell];
				persist_js_state=(persist_js_state || cell ∉ cells)
			)
			any_interrupted |= run.interrupted
		end

		cell.running = false
	end

	notebook.wants_to_interrupt = false
	flush_notebook_changes()
	# allow other `run_reactive!` calls to be executed
	put!(notebook.executetoken)
	return new_order
end

const lazymap = Base.Generator

function defined_variables(topology::NotebookTopology, cells)
	lazymap(cells) do cell
		topology.nodes[cell].definitions
	end
end

function defined_functions(topology::NotebookTopology, cells)
	lazymap(cells) do cell
		((cell.cell_id, namesig.name) for namesig in topology.nodes[cell].funcdefs_with_signatures)
	end
end

"Run a single cell non-reactively, set its output, return run information."
function run_single!(session_notebook::Union{Tuple{ServerSession,Notebook},WorkspaceManager.Workspace}, cell::Cell, reactive_node::ReactiveNode, expr_cache::ExprAnalysisCache; persist_js_state::Bool=false)
	run = WorkspaceManager.eval_format_fetch_in_workspace(
		session_notebook,
		expr_cache.parsedcode,
		cell.cell_id,
		ends_with_semicolon(cell.code),
		expr_cache.function_wrapped ? (filter(!is_joined_funcname, reactive_node.references), reactive_node.definitions) : nothing
	)
	set_output!(cell, run, expr_cache; persist_js_state=persist_js_state)
	if session_notebook isa Tuple && run.process_exited
		session_notebook[2].process_status = ProcessStatus.no_process
	end
	return run
end

function set_output!(cell::Cell, run, expr_cache::ExprAnalysisCache; persist_js_state::Bool=false)
	cell.output = CellOutput(
		body=run.output_formatted[1],
		mime=run.output_formatted[2],
		rootassignee=ends_with_semicolon(expr_cache.code) ? nothing : ExpressionExplorer.get_rootassignee(expr_cache.parsedcode),
		last_run_timestamp=time(),
		persist_js_state=persist_js_state,
	)
	cell.runtime = run.runtime
	cell.errored = run.errored
end

will_run_code(notebook::Notebook) = notebook.process_status != ProcessStatus.no_process && notebook.process_status != ProcessStatus.waiting_to_restart

###
# CONVENIENCE FUNCTIONS
###

"Do all the things!"
function update_save_run!(session::ServerSession, notebook::Notebook, cells::Array{Cell,1}; save::Bool=true, run_async::Bool=false, prerender_text::Bool=false, kwargs...)
	old = notebook.topology
	new = notebook.topology = updated_topology(old, notebook, cells)

	update_dependency_cache!(notebook)

	session.options.server.disable_writing_notebook_files || save_notebook(notebook)

	# _assume `prerender_text == false` if you want to skip some details_

	to_run_online = if !prerender_text
		cells
	else
		# this code block will run cells that only contain text offline, i.e. on the server process, before doing anything else
		# this makes the notebook load a lot faster - the front-end does not have to wait for each output, and perform costly reflows whenever one updates
		# "A Workspace on the main process, used to prerender markdown before starting a notebook process for speedy UI."
		original_pwd = pwd()
		offline_workspace = WorkspaceManager.make_workspace(
			(
				ServerSession(),
				notebook,
			),
			force_offline=true,
		)

		to_run_offline = filter(c -> !c.running && is_just_text(new, c) && is_just_text(old, c), cells)
		for cell in to_run_offline
			run_single!(offline_workspace, cell, new.nodes[cell], new.codes[cell])
		end

		cd(original_pwd)
		setdiff(cells, to_run_offline)
	end

	if will_run_code(notebook)
		run_task = @async begin
			run_reactive!(session, notebook, old, new, to_run_online; kwargs...)
		end
		if run_async
			run_task
		else
			fetch(run_task)
		end
	end
end

update_save_run!(session::ServerSession, notebook::Notebook, cell::Cell; kwargs...) = update_save_run!(session, notebook, [cell]; kwargs...)
update_run!(args...) = update_save_run!(args...; save=false)

"""
	throttled(f::Function, timeout::Real)

Return a function that when invoked, will only be triggered at most once
during `timeout` seconds.
The throttled function will run as much as it can, without ever
going more than once per `wait` duration.
Inspired by FluxML
See: https://github.com/FluxML/Flux.jl/blob/8afedcd6723112ff611555e350a8c84f4e1ad686/src/utils.jl#L662
"""
function throttled(f::Function, timeout::Real)
	tlock = ReentrantLock()
	iscoolnow = false
	run_later = false

	function flush()
		lock(tlock) do
			run_later = false
			f()
		end
	end

	function schedule()
		@async begin
			sleep(timeout)
			if run_later
				flush()
			end
			iscoolnow = true
		end
	end
	schedule()

	function throttled_f()
		if iscoolnow
			iscoolnow = false
			flush()
			schedule()
		else
			run_later = true
		end
	end

	return throttled_f, flush
end

ext: jl


text:
using Test
@time @testset "Section" begin
    include("test_section.jl")
end
@time @testset "Mesh" begin
    include("test_mesh.jl")
end
@time @testset "Modal" begin
    include("test_modal.jl")
end
@time @testset "Buckling" begin
    include("test_buckling.jl")
end
@time @testset "Transient dynamics" begin
    include("test_dyn.jl")
end


text:
#=
    The Objective Sense
=#

MOI.canget(m::LinQuadOptimizer, ::MOI.ObjectiveSense) = true
MOI.get(m::LinQuadOptimizer,::MOI.ObjectiveSense) = m.obj_sense

MOI.canset(::LinQuadOptimizer, ::MOI.ObjectiveSense) = true
function MOI.set!(m::LinQuadOptimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    if sense == MOI.MinSense
        change_objective_sense!(m, :min)
        m.obj_sense = MOI.MinSense
    elseif sense == MOI.MaxSense
        change_objective_sense!(m, :max)
        m.obj_sense = MOI.MaxSense
    elseif sense == MOI.FeasibilitySense
        # we set the objective sense to :min, and the objective to 0.0
        change_objective_sense!(m, :min)
        unsafe_set!(m, MOI.ObjectiveFunction{Linear}(), MOI.ScalarAffineFunction(VarInd[],Float64[],0.0))
        m.obj_type = AffineObjective
        m.obj_sense = MOI.FeasibilitySense
    else
        error("Sense $(sense) unknown.")
    end
end

#=
    The Objective Function
=#

function MOI.canset(m::LinQuadOptimizer, ::MOI.ObjectiveFunction{F}) where F<:MOI.AbstractFunction
    if MOI.get(m, MOI.ObjectiveSense()) == MOI.FeasibilitySense
        # it doesn't make sense to set an objective for a feasibility problem
        return false
    end
    return F in supported_objectives(m)
end

function MOI.set!(m::LinQuadOptimizer, ::MOI.ObjectiveFunction{MOI.SingleVariable}, var::MOI.SingleVariable)
    if m.obj_type == QuadraticObjective
        set_quadratic_objective!(m, Int[], Int[], Float64[])
    end
    m.obj_type = SingleVariableObjective
    m.single_obj_var = var.variable
    set_linear_objective!(m, [getcol(m, var.variable)], [1.0])
    m.objective_constant = 0
end

function MOI.set!(m::LinQuadOptimizer, ::MOI.ObjectiveFunction{F}, objf::Linear) where F
    cobjf = MOIU.canonical(objf)
    unsafe_set!(m, MOI.ObjectiveFunction{Linear}(), cobjf)
end

"""
    unsafe_set!(m, ::MOI.ObjectiveFunction{F}, objective::Linear) where F

Sets a linear objective function without cannonicalizing `objective`.
"""
function unsafe_set!(m::LinQuadOptimizer, ::MOI.ObjectiveFunction{F}, objf::Linear) where F
    if m.obj_type == QuadraticObjective
        # previous objective was quadratic...
        # zero quadratic part
        set_quadratic_objective!(m, Int[], Int[], Float64[])
    end
    m.obj_type = AffineObjective
    m.single_obj_var = nothing
    set_linear_objective!(m, getcol.(m, objf.variables), objf.coefficients)
    m.objective_constant = objf.constant
    nothing
end

function MOI.set!(m::LinQuadOptimizer, ::MOI.ObjectiveFunction, objf::Quad)
    m.obj_type = QuadraticObjective
    m.single_obj_var = nothing
    set_linear_objective!(m,
        getcol.(m, objf.affine_variables),
        objf.affine_coefficients
    )
    ri, ci, vi = reduceduplicates(
        getcol.(m, objf.quadratic_rowvariables),
        getcol.(m, objf.quadratic_colvariables),
        objf.quadratic_coefficients
    )
    set_quadratic_objective!(m,
        ri,
        ci,
        vi
    )
    m.objective_constant = objf.constant
    nothing
end

#=
    Get the objective function
=#

MOI.canget(m::LinQuadOptimizer, ::MOI.ObjectiveFunction{MOI.SingleVariable}) = m.obj_type == SingleVariableObjective
function MOI.get(m::LinQuadOptimizer, ::MOI.ObjectiveFunction{MOI.SingleVariable})
    SingleVariable(get(m.single_obj_var))
end

MOI.canget(m::LinQuadOptimizer, ::MOI.ObjectiveFunction{Linear}) = m.obj_type != QuadraticObjective
function MOI.get(m::LinQuadOptimizer, ::MOI.ObjectiveFunction{Linear})
    variable_coefficients = zeros(length(m.variable_references))
    get_linear_objective!(m, variable_coefficients)
    Linear(m.variable_references, variable_coefficients, m.objective_constant)
end

#=
    Modify objective function
=#

MOI.canmodifyobjective(m::LinQuadOptimizer, ::Type{MOI.ScalarCoefficientChange{Float64}}) = true
function MOI.modifyobjective!(m::LinQuadOptimizer, chg::MOI.ScalarCoefficientChange{Float64})
    if m.obj_type == SingleVariableObjective
        m.obj_type = AffineObjective
        m.single_obj_var = nothing
    end
    col = m.variable_mapping[chg.variable]
    # 0 row is the objective
    change_coefficient!(m, 0, col, chg.new_coefficient)
end


text:
struct OptimFitResults{T<:Optim.OptimizationResults, V<:Optim.AbstractOptimizerState}
    optres::T
    state::V
    # forwarddiff_hessian_callback::Function
end

# the optimizer interface <: Optimizer should implement
#  - minimize(objective, init, optimizer)
#  - minimizer(minimizationResult)
#  - measurement(minimizationResult)
#  - minimum(minimizationResult)
struct BFGSApproxHesse end

function minimize(func, init, optimizer::BFGSApproxHesse;
        algorithm = BFGS(), autodiff = :forwarddiff, show_trace::Bool = true, kws...)
    obj = OnceDifferentiable(func, init; autodiff)
    #
    options = Optim.Options(; show_trace, kws...)
    state = Optim.initial_state(algorithm, options, obj, init)
    #
    optres = optimize(obj, init, algorithm, options, state)
    #
    return OptimFitResults(optres, state)
end
#
#
minimizer(fr::OptimFitResults) = minimizer(fr.optres)
minimum(fr::OptimFitResults) = minimum(fr.optres)
#
#
invH(fr::OptimFitResults) = fr.state.invH
covmat(fr::OptimFitResults) = invH(fr)
errors(fr::OptimFitResults) = sqrt.(diag(fr.state.invH))
#
hessian_callback(p) = ForwardDiff.hessian(func, p)
invexacthessian(fr::OptimFitResults, p = minimizer(fr)) = inv(fr.forwarddiff_hessian_callback(p))
#
#
function measurements(fr, exacthessian::Bool=false)
    mv = minimizer(fr)
    δv = !(exacthessian) ? errors(fr) : sqrt.(diag(invexacthessian(fr)))
    return [±(m,δ) for (m,δ) in zip(mv, δv)]
end


text:
# SPDX-License-Identifier: X11
# 2020-12-07
# Day 3, Part 1

function main()
  input = readlines()
  # [cols][rows]

  cnt = 0
  x = 4
  bound = length(input[1])

  for i ∈ 2:length(input)
    input[i][x] == '#' && (cnt += 1)
    x += 3
    x > bound && (x = x - bound)
  end

  println(cnt)
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
module Resample

import LinearAlgebra: dot
import Statistics: mean

using Tables
using Dates

export resample
export Mean, Sum, First, None# TODO:, Last, Min, Max, None

abstract type AbstractSampleMethod end

abstract type AbstractWeightedSampleMethod <: AbstractSampleMethod end
abstract type AbstractSingleSampleMethod <: AbstractSampleMethod end

struct Mean <: AbstractWeightedSampleMethod end
struct Sum <: AbstractWeightedSampleMethod end

struct First <: AbstractSingleSampleMethod end
#= TODO:
struct Last <: AbstractSingleSampleMethod end
struct Min <: AbstractSingleSampleMethod end
struct Max <: AbstractSingleSampleMethod end
=#

struct None <: AbstractSampleMethod end

const WEIGHTED_SAMPLE_METHODS = (Mean(), Sum())
const SINGLE_SAMPLE_METHODS = (First(),)

include("tables.jl")

function resample(data::AbstractVector, org_inds, new_inds, method::AbstractWeightedSampleMethod = Mean())
    data, org_inds = sort_input(data, org_inds)
    new_inds = get_new_indices(org_inds, new_inds)
    inds, weights = resample_indices_and_weights(org_inds, new_inds, method)
    new_data = zeros(size(new_inds))
    for i in eachindex(new_data)
        new_data[i] = dot(data[inds[i]], weights[i])
    end

    return new_data
end

function resample(data::AbstractVector, org_inds, new_inds, method::AbstractSingleSampleMethod)
    data, org_inds = sort_input(data, org_inds)
    new_inds = get_new_indices(org_inds, new_inds)
    inds = resample_indices(org_inds, new_inds, method)
    new_data = similar(data, size(new_inds))
    for i in eachindex(new_data)
        new_data[i] = data[inds[i]]
    end

    return new_data
end

function resample(data::AbstractVector, org_inds, new_inds, ::None)
    @warn "The resample method `None()` is a placeholder to ignore DataFrame columns in resampling."
end

function sort_input(data, org_inds)
    sort_inds = issorted(org_inds) ? eachindex(org_inds) : sortperm(org_inds)
    return data[sort_inds], org_inds[sort_inds]
end

get_step(inds::StepRange) = step(inds)
get_step(inds) = inds[end] - inds[end-1]

function resample_indices_and_weights(org_inds, new_inds, method)
    Δ = get_step(new_inds)

    inds = [1:0 for i in eachindex(new_inds)]
    weights = [Float64[] for i in eachindex(new_inds)]

    for i in eachindex(new_inds)
        start = new_inds[i]
        stop = i == length(new_inds) ? new_inds[i] + Δ : new_inds[i + 1]
        inds[i], weights[i] = get_sample(org_inds, start, stop, method)
    end

    return inds, weights
end

function resample_indices(org_inds, new_inds, method)
    Δ = get_step(new_inds)
    inds = initialize_indices(length(new_inds), method)

    for i in eachindex(inds)
        start, stop = get_limits(new_inds, i, Δ)
        inds[i] = sample_indices(org_inds, start, stop, method)
    end

    return inds
end

function resample_indices(org_inds, new_inds, method::First)
    inds = initialize_indices(length(new_inds), method)

    for i in eachindex(inds)
        ind = findlast(j -> j <= new_inds[i], org_inds)
        if ind === nothing
            ind = 1
        end
        inds[i] = ind
    end

    return inds
end

function resample_weights(org_inds, new_inds, inds, method)
    Δ = get_step(new_inds)
    weights = [Float64[] for i in eachindex(new_inds)]

    for i in eachindex(new_inds)
        start, stop = get_limits(new_inds, i, Δ)
        weights[i] = sample_weights(org_inds, extrema(inds[i])..., start, stop, method)
    end

    return weights
end

initialize_indices(n, ::AbstractWeightedSampleMethod) = [1:0 for i in 1:n]
initialize_indices(n, ::AbstractSingleSampleMethod) = zeros(Int, n)

function get_limits(indices, i, Δ)
    start = indices[i]
    stop = i == length(indices) ? indices[i] + Δ : indices[i + 1]
    return start, stop
end

function sample_indices(org_inds, start, stop, ::AbstractWeightedSampleMethod)
    istart = findlast(i -> i <= start, org_inds)
    if istart === nothing
        istart = 1
    end
    istop = findlast(i -> i < stop, org_inds)
    if istop === nothing
        istop = 1
    end
    return istart:istop
end

function sample_weights(org_inds, istart, istop, start, stop, ::Mean)
    weights = diff([start; org_inds[istart + 1:istop]; stop]) ./ (stop - start)
    return weights
end

function sample_weights(org_inds, istart, istop, start, stop, ::Sum)
    start_prev = org_inds[istart]
    stop_next = istop == length(org_inds) ? stop : org_inds[istop + 1]

    weights = if istart == istop
        [(stop_next - start_prev) / (stop - start)]
    else
        start_next = org_inds[istart + 1]
        wstart = (start_next - start) / (start_next - start_prev)

        stop_prev = org_inds[istop]
        wstop = (stop - stop_prev) / (stop_next - stop_prev)

        [wstart; ones(istop - istart - 1); wstop]
    end
    return weights
end

function get_sample(org_inds, start, stop, method::Mean)
    istart, istop = extrema(sample_indices(org_inds, start, stop, method))
    weights = diff([start; org_inds[istart + 1:istop]; stop]) ./ (stop - start)
    return istart:istop, weights
end

function get_sample(org_inds, start, stop, method::Sum)
    istart, istop = extrema(sample_indices(org_inds, start, stop, method))

    start_prev = org_inds[istart]
    stop_next = istop == length(org_inds) ? stop : org_inds[istop + 1]

    weights = if istart == istop
        [(stop_next - start_prev) / (stop - start)]
    else
        start_next = org_inds[istart + 1]
        wstart = (start_next - start) / (start_next - start_prev)

        stop_prev = org_inds[istop]
        wstop = (stop - stop_prev) / (stop_next - stop_prev)

        [wstart; ones(istop - istart - 1); wstop]
    end
    return istart:istop, weights
end

get_new_indices(org_inds, new_inds::AbstractVector) = new_inds
get_new_indices(org_inds, step) = first(org_inds):step:last(org_inds)

end # module


text:
function get_xlim(data::DataFrame)
    num_particles_column = convert(Vector{Int},data[COL_NUM_PARTICLES])

    # add some padding on the left and right (in logspace)
    #return [10/1.1, 100*1.1] # TODO
    return [minimum(num_particles_column)/1.1, maximum(num_particles_column)*1.1]
end

function add_filled_circle(x, y, color)
    #ax = plt[:gca]()
    #xlim = ax[:get_xlim]()
    #ylim = ax[:get_ylim]()
    plt[:scatter]([x], [y], color=color, s=50, edgecolor="black", zorder=100)
    #ax[:set_xlim](xlim)
    #ax[:set_ylim](ylim)
end

function draw_arrow(x, y, xlen, ylen)
    ax = plt[:gca]()
    ax[:arrow](x, y, xlen, ylen, fc="black", ec="black",
        transform=ax[:transAxes], length_includes_head=true,
        head_width=0.05, head_length=0.05)
end

function text(x, y, content, fontsize)
    ax = plt[:gca]()
    plt[:text](x, y, content, fontsize=fontsize, transform=ax[:transAxes])
end


text:
using PiCamera
using Test


text:
#This file contains the incomplete factorization routine

#This function computes the inner procuct of two sparse vectors that are
#stores using the same nzind and nzval, but different offsets iter1:u1 resp.
#iter2:u2
function _innerprod( iter1::Int64,
                     u1::Int64,
                     iter2::Int64,
                     u2::Int64,
                     nzind::Array{Int64,1},
                     nzval::Array{Float64,1} )
  @inbounds begin
    res = zero(Float64)
    while ( iter1 <= u1 ) && ( iter2 <= u2 )
      while (nzind[iter1] == nzind[iter2])  &&( iter1 <= u1 ) && ( iter2 <= u2 )
        res += nzval[iter1] * nzval[iter2]
        iter1 += 1
        iter2 += 1
      end
      if nzind[iter1] < nzind[iter2]
        iter1 += 1
      else
        iter2 += 1
      end
    end
    return res

  end
end

#This function computes the in place incomplete Cholesky factorisation using
#an inner product wrapped in a function.
function _el_icholRWHL!( s::Array{Float64},
                          ind::Array{Int64, 1},
                          jnd::Array{Int64, 1},
                          )
  @inbounds begin
    #piv_i goes from 1 to N and designates the number, which is being treated
    for piv_i = 1 : ( size( ind, 1 ) - 1 )
      #piv_j ranges over pointers and designates where to find the present column
      for piv_j = ind[ piv_i ] : ( ind[ piv_i + 1 ] - 1 )
        #iter designates the pointer, at which the i-iterator starts
        iter::Int64 = ind[ piv_i ]
        #jter designates the pointer, at which the j-iterator starts
        jter::Int64 = ind[ jnd[ piv_j ] ]
        #The condition makes sure, that only the columns up to th pivot
        #are being summed.

        s[ piv_j ] -= _innerprod( iter, ind[ piv_i + 1 ] - 2,
                                  jter, ind[ jnd[ piv_j ] + 1 ] - 2,
                                  jnd,
                                  s)

        if jnd[ piv_j ] < piv_i
          s[ piv_j ] /= ( s[ ind[ jnd[ piv_j ] + 1 ] - 1 ] )
        elseif jnd[ piv_j ] == piv_i
          if ( s[ ind[ piv_i + 1 ] - 1 ] ) <= 0.
            s[ ind[ piv_i + 1 ] - 1 ] = 0
            #Treating the case where the matrix is near-low rank and our
            #factorization returns a low-rank matrix.
            maxj = piv_i
            for newpiv_i = ( piv_i + 1 ) : ( size( ind, 1 ) - 1 )
              newpiv_j::Int64 = ind[ newpiv_i ]
              while jnd[ newpiv_j ] < maxj
                #iter designates the pointer, at which the i-iterator starts
                iter = ind[ newpiv_i ]
                #jter designates the pointer, at which the j-iterator starts
                jter = ind[ jnd[ newpiv_j ] ]
                #The condition makes sure, that only the columns up to th pivot
                #are being summed.

                s[ newpiv_j ] -= _innerprod( iter, ind[ newpiv_i + 1 ] - 2,
                                          jter, ind[ jnd[ newpiv_j ] + 1 ] - 2,
                                          jnd,
                                          s)

                s[ newpiv_j ] /= ( s[ ind[ jnd[ newpiv_j ] + 1 ] - 1 ] )
                newpiv_j += 1
              end
              while newpiv_j <= ( ind[ newpiv_i + 1 ] - 1 )
                s[ newpiv_j ] = 0.
                newpiv_j += 1
              end
            end
            println( "Pivot in i = ", piv_i, " is nonpositive",
                     " returning low rank approximation" )
            return 0
          else
            s[ piv_j ] = sqrt( s[ piv_j ] )
          end
        else
          println( "Warning: jnd[ piv_j ] > piv_i" )
        end
      end
    end
  end
end

#In-place Cholesky factorization of an uppper triangular matrix
function icholU_high_level!( U::SparseMatrixCSC{Float64, Int64} )
  @fastmath _el_icholRWHL!( U.nzval, U.colptr, U.rowval )
end

max_stars_repo_name: f-t-s/Spelliptic

max_issues_repo_name: f-t-s/Spelliptic

max_forks_repo_name: f-t-s/Spelliptic


text:
"""
Optimization function to calculate warping for elastic regression

    regression_warp(beta, timet, q, y, alpha)
    :param beta: regression function
    :param timet: vector describing time samples
    :param q: vector describing srsf
    :param y: response value
    :param alpha: intercept

    Returns
    :return gamma_new: new gamma
"""
function regression_warp(beta::Vector, timet::Vector, q::Vector, y::Float64,
                         alpha::Float64)
    gam_M = optimum_reparam(beta, timet, q, method="DP2");
    qM = warp_q_gamma(timet, q, gam_M);
    y_M = trapz(timet, qM.*beta);

    gam_m = optimum_reparam(-1.*beta, timet, q, method="DP2");
    qm = warp_q_gamma(timet, q, gam_m);
    y_m = trapz(timet, qm.*beta);

    if y > alpha+y_M
        gamma_new = gam_M;
    elseif y < alpha + y_m
        gamma_new = gam_m;
    else
        gamma_new = zero_crossing(y-alpha, q, beta, timet, y_M, y_m,
                                  gam_M, gam_m);
    end

    return gamma_new

end

"""
Calculate warping for logistic regression

    logistic_warp(beta, timet, q, y)
    :param beta: regression function
    :param timet: time samples
    :param q: srsf
    :param y: response

    Returns
    :return gamma: new gamma
"""
function logistic_warp(beta::Vector, timet::Vector, q::Array, y)
    if y == 1
        gamma = optimum_reparam(beta, timet, q, method="DP2");
    elseif y == -1
        gamma = optimum_reparam(-1.*beta, timet, q, method="DP2");
    end

    return gamma

end

"""
Calculate logistic optimization function

    logit_optm(x::Vector, grad::Vector, Phi, y)
    :param x: samples
    :param grad: gradient
    :param Phi: coefficient matrix
    :param y: response
"""
function logit_optm(x::Vector, grad::Vector, Phi, y)
    if length(grad) > 0
        logit_gradient!(x, grad, Phi, y);
    end
    out = logit_loss(x, Phi, y);
    return out
end

"""
Logistic function

    phi(t)

    Returns
    :return out: phi(t)
"""
function phi(t)
    idx = t .> 0;
    out = Array(Float64, size(t));
    if sum(idx) > 0
        out[idx] = 1./(1+exp(-1*t[idx]));
    end
    exp_t = exp(t[~idx]);
    out[~idx] = exp_t ./ (1 + exp_t);
    return out
end

"""
Calculate logistic loss function

    logit_loss(b, X, y)
    :param b: coefficients
    :param X: matrix
    :param y: response

    Returns
    :return out: loss function
"""
function logit_loss(b, X, y)
    z = X * b;
    yz = y .* z;
    idx = yz .> 0;
    out = zeros(yz);
    out[idx] = log(1+exp(-1.*yz[idx]));
    out[~idx] = (-yz[~idx] + log(1+exp(yz[~idx])));
    out = sum(out);
    return out
end

"""
Calculate gradient of logistic optimization in place

    logit_gradient!(b, grad, X, y)
    :param b: coefficients
    :param grad: gradient
    :param X: matrix
    :param y: response
"""
function logit_gradient!(b, grad, X, y)
    z = X * b;
    z = phi(y.*z);
    z0 = (z-1) .* y;
    grad[:] = X.' * z0;
end

"""
Calculate Hessian of logistic optimization

    logit_hessian(s, b, X, y)
    :param s:
    :param b: coefficients
    :param X: matrix
    :param y: response
"""
function logit_hessian(s, b, X, y)
    z = X*b;
    z = phi(y.*z);
    d = z*(1-z);
    wa = d.*(X * s);
    Hs = X.' * wa;
    return Hs
end

"""
Calculate m-logistic warping using gradient method

    mlogit_warp_grad(alpha, beta, timet, q, y; max_itr=8000, tol=1e-10,
                     delt=0.008, display=0)
    :param alpha: intercept
    :param beta: regression function
    :param timet: vector describing time samples
    :param q: srsf
    :param y: response
    :param max_itr: maximum number of iterations
    :param tol: stopping tolerance
    :param delt: gradient step size
    :param display: display optimization iterations
"""
function mlogit_warp_grad(alpha, beta, timet, q, y; max_itr=8000,
                          tol=1e-10, delt=0.008, display=0)
    m1 = length(timet);
    m = size(beta,2);
    q /= norm(q);
    alpha /= norm(alpha);
    for i in 1:m
        beta[:, i] = beta[:, i] / norm(beta[:, i]);
    end
    gam1 = collect(linspace(0, 1, m1));
    gamout = zeros(m1);
    beta1 = reshape(beta, m1*m, 1);

    ccall((:mlogit_warp_grad, libfdasrsf), Void,
          (Ptr{Int32}, Ptr{Int32}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
          Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, Ptr{Float64},
          Ptr{Float64},  Ptr{Int32}, Ptr{Float64}),
          &m1, &m, alpha, beta1, timet, gam1, q, y, &max_itr, &tol, &delt,
          &display, gamout)

    return gamout
end

"""
Calculate warping for m-logistic elastic regression

    mlogit_optim(x, grad, Phi, Y)
    :param x: sample
    :param grad: gradient
    :param Phi: matrix
    :param Y: response matrix
"""
function mlogit_optm(x::Vector, grad::Vector, Phi, Y)
    if length(grad) > 0
        mlogit_gradient!(x, grad, Phi, Y);
    end
    out = mlogit_loss(x, Phi, Y);
    return out
end

"""
Calculate loss for m-logistic regression

    mlogit_loss(b, X, Y)
    :param b: coefficients
    :param X: matrix
    :param Y: response matrix
"""
function mlogit_loss(b, X, Y)
    N, m = size(Y);  # n_samples, n_classes
    M = size(X,2);   # n_features
    B = reshape(b, M, m);
    Yhat = X * B;
    Yhat -= repmat(minimum(Yhat, 2),1,size(Yhat,2));
    Yhat = exp(-1.*Yhat);
    # l1-normalize
    Yhat = Yhat./repmat(sum(Yhat, 2),1,size(Yhat,2));

    Yhat = Yhat .* Y;
    nll = sum(log(sum(Yhat,2)));
    nll /= -N;

    return nll
end

"""
Calculate m-logistic elastic regression loss function gradient in place

    mlogit_gradient!(b, grad, X, Y)
    :param b: coefficients
    :param grad: gradient
    :param X: matrix
    :param Y: response matrix
"""
function mlogit_gradient!(b, grad, X, Y)
    N, m = size(Y);  # n_samples, n_classes
    M = size(X,2); # n_features
    B = reshape(b, M, m);
    Yhat = X * B;
    Yhat -= repmat(minimum(Yhat, 2),1,size(Yhat,2));
    Yhat = exp(-1.*Yhat);
    # l1-normalize
    Yhat = Yhat./repmat(sum(Yhat, 2),1,size(Yhat,2));

    _Yhat = Yhat .* Y;
    _Yhat = _Yhat./repmat(sum(_Yhat,2),1,size(_Yhat,2));
    Yhat -= _Yhat;
    grad1 = X.' * Yhat;
    grad1 /= -N;
    grad[:] = reshape(grad1, M*m, 1);
end


text:
#=
Surprisingly there are only three numbers that can be written as the sum of fourth powers of their digits:

1634 = 14 + 64 + 34 + 44
8208 = 84 + 24 + 04 + 84
9474 = 94 + 44 + 74 + 44
As 1 = 14 is not a sum it is not included.

The sum of these numbers is 1634 + 8208 + 9474 = 19316.

Find the sum of all the numbers that can be written as the sum of fifth powers of their digits.
=#
#=
Note: Calculating limit:
Checking the number 111111111.
That number has 9 digits, so the maximum sum would be 9*59049 = 531441, which doesn't even come close to 111111111.
So any 9 digit number or greater can't ever reach a big enough sum.
=#
function sumofdigitpowers(a)
  total, i = 0, 2
  while i < length(digits(i))*9^a
    current = sum(map((x)->x^a, digits(i)))
    if current == i
      total += i
    end
    i+=1
  end
  total
end
@time println(sumofdigitpowers(5))

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
module MATLABCluster

using MATLAB

export MATLABManager
import Base.readline

immutable MATLABManager <: ClusterManager
    profile::ASCIIString
end
MATLABManager() = MATLABManager("")

type ConnectionInfoIOHack <: IO
    conninfo::ByteString
    line_buffered::Bool # unused
end
ConnectionInfoIOHack(conninfo) = ConnectionInfoIOHack(conninfo, true)
readline(x::ConnectionInfoIOHack) = x.conninfo

# Read the file written by diary and check if Julia is running
function get_worker_info(diary)
    for l in split(diary, '\n')
        private_hostname, port = Base.parse_connection_info(l)
        if private_hostname != ""
            return ConnectionInfoIOHack(l)
        end
    end
end

function Base.launch(cman::MATLABManager, np::Integer, config::Dict, instances_arr::Array, c::Condition)
    try
        exe = Base.shell_escape("$(config[:dir])/$(config[:exename])")
        exeflags = Base.shell_escape(config[:exeflags].exec...)
        cmd = replace("LD_LIBRARY_PATH= OMP_NUM_THREADS=1 $exe $exeflags", "'", "''")

        # Start MATLAB jobs
        jobvar = "jobs_$(randstring(12))"
        eval_string("""
            clust = parcluster($(cman.profile != "" ? "'$(cman.profile)'" : ""));
            username = clust.Username;
            host = clust.Host;

            j = createJob(clust);
            $jobvar = cell($np, 1);
            for i = 1:$(np)
                $jobvar{i} = createTask(j, @() system('$cmd', '-echo'), 0, {}, 'CaptureDiary', true);
            end
            submit(j);
        """)

        @mget username host

        print("Waiting for jobs to start...")
        for i = 1:np
            # Wait until MATLAB says job is running
            eval_string("""
                j = $jobvar{$i};
                wait(j, 'running')
            """)

            # Hack to wait until Julia is running
            local worker_info
            while true
                eval_string("""
                    state = j.State;
                    warning('off', 'all');
                    diary = j.Diary;
                    if isempty(diary)
                        diary = ' ';
                    end
                    warning('on', 'all');
                """)
                @mget state diary

                if state == "finished"
                    error("task failed:\n$diary")
                else
                    worker_info = get_worker_info(diary)
                    worker_info == nothing || break
                end
                # sleep(0.25)
            end
            inst = (worker_info, "$username@$host", merge(config, Dict(:jobvar => jobvar, :jobindex => i)))
            push!(instances_arr, (inst,))
            notify(c)
        end
        print("done\n")
    catch e
        showerror(STDERR, e, catch_backtrace())
        println()
        rethrow(e)
    end
    notify(c)
end

function Base.manage(cman::MATLABManager, id::Integer, config::Dict, op::Symbol)
    if op == :finalize
        eval_string("delete($(config[:jobvar]){$(config[:jobindex])});")
    end
end

end # module


text:
# Estimation of text extents using pango.

using Fontconfig

const libpangocairo = Cairo._jl_libpangocairo
const libpango = Cairo._jl_libpango

# Cairo text backend
const CAIRO_FONT_TYPE_TOY = 0
const CAIRO_FONT_TYPE_FT = 1
const CAIRO_FONT_TYPE_WIN32 = 2
const CAIRO_FONT_TYPE_QUARTZ = 3
const CAIRO_FONT_TYPE_USER = 4

# Mirroring a #define in the pango header.
const PANGO_SCALE = 1024.0

function pango_fmt_float(x::Float64)
    return @sprintf("%0.4f", x)
end

# Use the freetype/fontconfig backend to find the best match to a font
# description.
#
# Args:
#   desc: A string giving the font description. This can
#         also provide a comma-seperated list of families. E.g.,
#         "Helvetica, Arial 10"
#
# Returns:
#   A pointer to a PangoFontDescription with the closest match.
#
let available_font_families = Set{String}()
    for font_pattern in Fontconfig.list()
        push!(available_font_families, lowercase(format(font_pattern, "%{family}")))
    end

    const meta_families = Set(["serif", "sans", "sans-serif", "monospace",
                               "cursive", "fantasy"])

    global match_font
    function match_font(families::String, size::Float64)
        matched_family = "sans-serif"
        for family in [lowercase(strip(family, [' ', '"', '\''])) for family in split(families, ',')]
            if family in available_font_families || family in meta_families
                matched_family = family
                break
            end
        end
        family = format(match(Fontconfig.Pattern(family=family)), "%{family}")
        desc = @sprintf("%s %fpx", family, size)
        fd = ccall((:pango_font_description_from_string, libpango),
                   Ptr{Void}, (Ptr{Uint8},), bytestring(desc))
        return fd
    end
end

# Thin wrapper for a pango_layout object.
type PangoLayout
    layout::Ptr{Void}

    function PangoLayout()
        layout = ccall((:pango_layout_new, libpango),
                       Ptr{Void}, (Ptr{Void},), pango_cairo_ctx)
        # TODO: finalizer?

        new(layout)
    end
end

# Set the layout's font.
function pango_set_font(pangolayout::PangoLayout, family::String, pts::Number)
    fd = match_font(family, pts)
    ccall((:pango_layout_set_font_description, libpango),
          Void, (Ptr{Void}, Ptr{Void}), pangolayout.layout, fd)
end

# Find the width and height of a string.
#
# Args:
#   pangolayout: a pango layout object, with font, etc, set.
#   text: a string we might like to draw.
#
# Returns:
#   A (width, height) tuple in absolute units.
#
function pango_text_extents(pangolayout::PangoLayout, text::String)
    textarray = convert(Vector{Uint8}, bytestring(text))
    ccall((:pango_layout_set_markup, libpango),
          Void, (Ptr{Void}, Ptr{Uint8}, Int32),
          pangolayout.layout, textarray, length(textarray))

    extents = Array(Int32, 4)
    ccall((:pango_layout_get_extents, libpango),
          Void, (Ptr{Void}, Ptr{Int32}, Ptr{Int32}),
          pangolayout.layout, extents, C_NULL)

    width, height = (extents[3] / PANGO_SCALE)pt, (extents[4] / PANGO_SCALE)pt
end

# Find the minimum width and height needed to fit any of the given strings.
#
# (A "user-friendly" wrapper for pango_text_extents.)
#
# Args:
#   font_family: Something like a font name.
#   pts: Font size in points.
#   texts: One or more strings.
#
# Returns:
#   A (width, height) tuple in absolute units.
#
function max_text_extents(font_family::String, pts::Float64, texts::String...)
    pango_set_font(pangolayout::PangoLayout, font_family, pts)
    max_width  = 0mm
    max_height = 0mm
    for text in texts
        (width, height) = pango_text_extents(pangolayout::PangoLayout, text)
        max_width  = max_width.abs  < width.abs  ? width  : max_width
        max_height = max_height.abs < height.abs ? height : max_height
    end
    return (max_width, max_height)
end

# Same as max_text_extents but with font_size in arbitrary absolute units.
function max_text_extents(font_family::String, size::Measure,
                      texts::String...)
    if !isabsolute(size)
        error("text_extents requries font size be in absolute units")
    end
    return max_text_extents(font_family, size/pt, texts...)
end

# Return an array with the extents of each element
function text_extents(font_family::String, pts::Float64, texts::String...)
    pango_set_font(pangolayout::PangoLayout, font_family, pts)
    return [pango_text_extents(pangolayout::PangoLayout, text)
            for text in texts]
end

function text_extents(font_family::String, size::Measure, texts::String...)
    return text_extents(font_family, size/pt, texts...)
end

const pango_attrs = [
    (:PANGO_ATTR_LANGUAGE,        :PangoAttrLanguage),
    (:PANGO_ATTR_FAMILY,          :PangoAttrString),
    (:PANGO_ATTR_STYLE,           :PangoAttrInt),
    (:PANGO_ATTR_WEIGHT,          :PangoAttrInt),
    (:PANGO_ATTR_VARIANT,         :PangoAttrInt),
    (:PANGO_ATTR_STRETCH,         :PangoAttrInt),
    (:PANGO_ATTR_SIZE,            :PangoAttrSize),
    (:PANGO_ATTR_FONT_DESC,       :PangoAttrFontDesc),
    (:PANGO_ATTR_FOREGROUND,      :PangoAttrColor),
    (:PANGO_ATTR_BACKGROUND,      :PangoAttrColor),
    (:PANGO_ATTR_UNDERLINE,       :PangoAttrInt),
    (:PANGO_ATTR_STRIKETHROUGH,   :PangoAttrInt),
    (:PANGO_ATTR_RISE,            :PangoAttrInt),
    (:PANGO_ATTR_SHAPE,           :PangoAttrShape),
    (:PANGO_ATTR_SCALE,           :PangoAttrFloat),
    (:PANGO_ATTR_FALLBACK,        :PangoAttrFallback),
    (:PANGO_ATTR_LETTER_SPACING,  :PangoAttrInt),
    (:PANGO_ATTR_UNDERLINE_COLOR, :PangoAttrColor),
    (:PANGO_ATTR_ABSOLUTE_SIZE,   :PangoAttrSize),
    (:PANGO_ATTR_GRAVITY,         :PangoAttrInt),
    (:PANGO_ATTR_GRAVITY_HINT,    :PangoAttrInt)]

for (i, (attr, t)) in enumerate(pango_attrs)
    @eval begin
        const $attr = $i
    end
end

const PANGO_STYLE_NORMAL  = 0
const PANGO_STYLE_OBLIQUE = 1
const PANGO_STYLE_ITALIC  = 2

const PANGO_WEIGHT_THIN = 100
const PANGO_WEIGHT_ULTRALIGHT = 200
const PANGO_WEIGHT_LIGHT = 300
const PANGO_WEIGHT_BOOK = 380
const PANGO_WEIGHT_NORMAL = 400
const PANGO_WEIGHT_MEDIUM = 500
const PANGO_WEIGHT_SEMIBOLD = 600
const PANGO_WEIGHT_BOLD = 700
const PANGO_WEIGHT_ULTRABOLD = 800
const PANGO_WEIGHT_HEAVY = 900
const PANGO_WEIGHT_ULTRAHEAVY = 1000

# A Julia manifestation of a set of pango attributes
type PangoAttr
    rise::Maybe(Int)
    scale::Maybe(Float64)
    style::Maybe(Int)
    weight::Maybe(Int)

    function PangoAttr()
        new(nothing, nothing, nothing, nothing)
    end
end

function isempty(attr::PangoAttr)
    all([getfield(attr, name) === nothing for name in PangoAttr.names])
end

# Set an attribute in a PangoAttr
#
# Args:
#   attr: A PangoAttr to update.
#   attr_name: A pango attribute name (e.g., :PANGO_ATTR_RISE)
#   value: The value with which to update the attribute.
#
# Returns:
#   The attr.
function update_pango_attr(attr::PangoAttr, attr_name::Symbol, value)
    if attr_name == :PANGO_ATTR_RISE
        attr.rise = int64(value)
    elseif attr_name == :PANGO_ATTR_SCALE
        attr.scale = value
    elseif attr_name == :PANGO_ATTR_STYLE
        attr.style = int64(value)
    elseif attr_name == :PANGO_ATTR_WEIGHT
        attr.weight = int64(value)
    end
    attr
end

# Unpack the first part of a pango attribute
#
# Args:
#   ptr: A pointer to a PangoAttribute
#   t: The type of the attribute (e.g. PangoAttrInt)
#
# Returns:
#   A tuple of the form (start_idx, end_idx, value)
#
function unpack_pango_attr(ptr::Ptr{Void}, t::Symbol)
    ptr += sizeof(Ptr{Void}) # skip `klass` pointer
    ptr = convert(Ptr{Uint32}, ptr)
    idx = pointer_to_array(ptr, (2,))
    ptr += 2 * sizeof(Uint32)
    ptr = convert(Ptr{Void}, ptr)

    if t == :PangoAttrInt
        value = unpack_pango_int(ptr)
    elseif t == :PangoAttrFloat
        value = unpack_pango_float(ptr)
    else
        value = nothing
    end

    (idx[1], idx[2], value)
end

# Unpack a pango int attribute.
#
# Args:
#   ptr: A point to a PangoAttrInt plus sizeof(PangoAttribute)
#
# Returns:
#   And int value.
function unpack_pango_int(ptr::Ptr{Void})
    ptr = convert(Ptr{Int32}, ptr)
    pointer_to_array(ptr, (1,))[1]
end

function unpack_pango_float(ptr::Ptr{Void})
    ptr = convert(Ptr{Float64}, ptr)
    pointer_to_array(ptr, (1,))[1]
end

#function unpack_pango_size(ptr::Ptr{Void})
    #ptr = convert(Ptr{Int32}, ptr)
    #size = point_to_array(ptr, (1,))[1]
    #ptr = convert(Ptr{Uint32}, ptr)
    #absolute = point_to_array(ptr, (1,))[1] & 0x1

    #println(size, absolute)
    #nothing
#end

# TODO: unpacking other attributes

# Unpack a list of pango attributes
#
# Args:
#   ptr: A pointer to a PangoAttrList
#
# Returns:
#   A list of the form [(start_idx, attribute), ...] in which the start_idx
#   values are increasing and the attribute is a set of attributes that
#   should be applied starting at that position.
#
function unpack_pango_attr_list(ptr::Ptr{Void})
    attr_it = ccall((:pango_attr_list_get_iterator, libpango),
                    Ptr{Void}, (Ptr{Void},), ptr)

    # Alias some ugly C calls.
    attr_it_next = () -> ccall((:pango_attr_iterator_next, libpango),
                               Int32, (Ptr{Void},), attr_it)

    attr_it_get = attr_name -> ccall((:pango_attr_iterator_get, libpango),
                                     Ptr{Void}, (Ptr{Void}, Int32),
                                     attr_it, eval(attr_name))

    attr_it_range = () -> begin
        start_idx = Array(Int32, 1)
        end_idx = Array(Int32, 1)
        ccall((:pango_attr_iterator_range, libpango),
              Void, (Ptr{Void}, Ptr{Int32}, Ptr{Int32}),
              attr_it, start_idx, end_idx)
        (start_idx[1], end_idx[1])
    end

    attrs = Array((Int, PangoAttr), 0)

    while attr_it_next() != 0

        attr = PangoAttr()
        local start_idx

        for (attr_name, attr_type) in pango_attrs
            c_attr = attr_it_get(attr_name)
            (start_idx, end_idx) = attr_it_range()

            if c_attr != C_NULL
                (_, _, value) = unpack_pango_attr(c_attr, attr_type)
                update_pango_attr(attr, attr_name, value)
            end
        end

        push!(attrs, (start_idx, attr))
    end

    ccall((:pango_attr_iterator_destroy, libpango),
          Void, (Ptr{Void},), attr_it)

  attrs
end

function pango_to_svg(text::String)
    c_stripped_text = Array(Ptr{Uint8}, 1)
    c_attr_list = Array(Ptr{Void}, 1)

    ret = ccall((:pango_parse_markup, libpango),
                Int32, (Ptr{Uint8}, Int32, Uint32, Ptr{Ptr{Void}},
                        Ptr{Ptr{Uint8}}, Ptr{Uint32}, Ptr{Void}),
                bytestring(text), -1, 0, c_attr_list, c_stripped_text,
                C_NULL, C_NULL)

    if ret == 0
        error("Could not parse pango markup.")
    end

    # TODO: do c_stripped_text and c_attr_list need to be freed?

    bytearray =  str -> convert(Array{Uint8, 1}, str)

    text = bytearray(bytestring(c_stripped_text[1]))

    last_idx = 1
    open_tag = false
    baseline_shift = 0.0

    tagged_text = sprint() do io
        for (idx, attr) in unpack_pango_attr_list(c_attr_list[1])
            write(io, text[last_idx:idx])
            last_idx = idx + 1

            if open_tag
                write(io, "</tspan>")
            end

            if isempty(attr) && baseline_shift == 0.0
                open_tag = false
                continue
            end

            open_tag = true

            write(io, "<tspan style=\"dominant-baseline:inherit\"")

            # "baseline-shift" is not currently supported Firefox or IE.
            # if !(attr.rise === nothing)
            #     @printf(io, " baseline-shift=\"%s\"",
            #             fmt_float(((attr.rise / PANGO_SCALE)pt).abs))
            # end

            if !(attr.rise === nothing)
                bs = -((attr.rise / PANGO_SCALE)pt).abs
                @printf(io, " dy=\"%s\"", pango_fmt_float(bs))
                baseline_shift = bs
            elseif baseline_shift != 0.0
                @printf(io, " dy=\"%s\"", pango_fmt_float(-baseline_shift))
                baseline_shift = 0.0
            end

            if !(attr.scale === nothing)
                @printf(io, " font-size=\"%s%%\"",
                        pango_fmt_float(100.0 * attr.scale))
                baseline_shift *= attr.scale
            end

            if !(attr.style === nothing)
                if attr.style == PANGO_STYLE_NORMAL
                    @printf(io, " font-style=\"%s\"", "normal")
                elseif attr.style == PANGO_STYLE_OBLIQUE
                    @printf(io, " font-style=\"%s\"", "oblique")
                elseif attr.style == PANGO_STYLE_ITALIC
                    @printf(io, " font-style=\"%s\"", "italic")
                end
            end

            if !(attr.weight === nothing)
                @printf(io, " font-weight=\"%d\"", attr.weight)
            end

            write(io, ">")
        end

        write(io, text[last_idx:end])

        if open_tag
            write(io, "</tspan>")
        end
    end

    tagged_text
end


text:
"""
    timeseries(x::Node{{Union{Number, Point2}}})

Plots a sampled signal.
Usage:
```julia
signal = Node(1.0)
scene = timeseries(signal)
display(scene)
# @async is optional, but helps to continue evaluating more code
@async while isopen(scene)
    # aquire data from e.g. a sensor:
    data = rand()
    # update the signal
    signal[] = data
    # sleep/ wait for new data/ whatever...
    # It's important to yield here though, otherwise nothing will be rendered
    sleep(1/30)
end

```
"""
@recipe(TimeSeries, signal) do scene
    Attributes(
        history = 100;
        default_theme(scene, Lines)...
    )
end

signal2point(signal::Number, start) = Point2f0(time() - start, signal)
signal2point(signal::Point2, start) = signal
signal2point(signal, start) = error(""" Signal needs to be of type Number or Point.
Found: $(typeof(signal))
""")

function AbstractPlotting.plot!(plot::TimeSeries)
    # normal plotting code, building on any previously defined recipes
    # or atomic plotting operations, and adding to the combined `plot`:
    points = Node(fill(Point2f0(NaN), plot.history[]))
    buffer = copy(points[])
    lines!(plot, points)
    start = time()
    on(plot.signal) do x
        points[][end] = signal2point(x, start)
        circshift!(buffer, points[], 1)
        buff_ref = buffer
        buffer = points[]
        points[] = buff_ref
        update!(parent(plot))
    end
    plot
end


text:
"""
	OrderedPair

OrderedPair representation.
Ex:
```
julia> o = OrderedPair(1,3.2)
(1,3.2)
```
"""
mutable struct OrderedPair
	left::Real
	right::Real
	OrderedPair(left,right) = left > right ? error("OrderedPair must be ordered") : new(left,right)
end

"""
	string(p)

Returns a string representation of an `OrderedPair`.
"""
function string(p::OrderedPair)
	"($(p.left),$(p.right))"
end

"""
	show(io,p)

Prints the string representation of an `OrderedPair`.
"""
function show(io::IO, p::OrderedPair)
	println(string(p))
end

"""
	o1 == o2

Equality defined for `OrderedPair`.
"""
function ==(o1::OrderedPair, o2::OrderedPair)
	(o1.left == o2.left) && (o1.right == o2.right)
end

"""
	o1 < o2

Lexicographical order for `OrderedPair`.
"""
function isless(o1::OrderedPair, o2::OrderedPair)
	o1.left == o2.left ? o1.right < o2.right : o1.left < o2.left
end

"""
	Interval

`Interval` representation.
Attributes:
	- limits::OrderedPair : The left and right limits of the interval.
	- open_left::Bool : Set to true if the interval is open on the left.
	- open_right::Bool : Set to true if the interval is open on the right.
Ex:
```
julia> i = Interval(0,1)
[0,1]
julia> j = Interval(1,true,3.4)
]1,3.4]
julia> k = Interval(-1,true,3.12,true)
]-1,3.12[
```
"""
mutable struct Interval
    limits::OrderedPair
    open_left::Bool
    open_right::Bool
end

"""
	i1 == i2

Equality defined for `Interval`.
"""
function ==(i1::Interval, i2::Interval)
	(i1.limits == i2.limits) && (i1.open_left == i2.open_left) && (i1.open_right == i2.open_right)
end

"""
	i1 < i2

Lexicographical order for `Interval`.
Ex:
```
julia> Interval(0,2) < Interval(1,2,true)
true
julia> Interval(1,3) > Interval(-1,10)
true
julia> Interval(0,1) < Interval(0,1,true)
true
```
"""
function isless(i1::Interval, i2::Interval)
	if left(i1) < left(i2)
        return true
    elseif left(i1) == left(i2)
        if !i1.open_left & i2.open_left
        	return true
        elseif i1.open_left & !i2.open_left
        	return false
        else
        	if right(i1) < right(i2)
        		return true
        	elseif right(i1) == right(i2)
        		if !i1.open_right & i2.open_right
        			return true
        		else
        			return false
        		end
        	else
        		return false
        	end
        end
    else
        return false
    end
end

"""
	Interval(left, open_left, right, open_right)

Full constructor for `Interval`.
Requires:
	- left::Real Left limit of the interval.
	- open_left::Bool Set to true if the interval is open on the left.
	- right::Real Right limit of the interval.
	- open_right::Bool Set to true if the interval is open on the left.
Ex:
```
julia> i = Interval(0,false,1,true)
[0,1[
julia> j = Interval(1,false,2,false)
[1,2]
```
"""
function Interval(left::Real, open_left::Bool, right::Real, open_right::Bool)
	p = OrderedPair(left,right)
	Interval(p,open_left,open_right)
end

"""
	Interval(left, right)

Simplified constructor for closed `Interval`.
Requires:
	- left::Real Left limit of the interval.
	- right::Real Right limit of the interval.
Ex:
```
julia> i = Interval(0,2)
[0,2]
```
"""
function Interval(left::Real, right::Real)
	p = OrderedPair(left,right)
	Interval(p,false,false)
end

"""
	Interval(left, open_left, right)

Simplified constructor for `Interval` closed on the right.
Requires:
	- left::Real Left limit of the interval.
	- open_left::Bool Set to true if the interval is open on the left.
	- right::Real Right limit of the interval.
Ex:
```
julia> i = Interval(2,true,3)
]2,3]
```
"""
function Interval(left::Real, open_left::Bool, right::Real)
	p = OrderedPair(left,right)
	Interval(p,open_left,false)
end

"""
	Interval(left, right, open_right)

Simplified constructor for `Interval` closed on the left.
Requires:
	- left::Real Left limit of the interval.
	- right::Real Right limit of the interval.
	- open_right::Bool Set to true if the interval is open on the right.
Ex:
```
julia> i = Interval(-1,2,true)
[-1,2[
```
"""
function Interval(left::Real, right::Real, open_right::Bool)
	p = OrderedPair(left,right)
	Interval(p,false,open_right)
end

"""
	Interval()

Empty `Interval`.
Ex:
```
julia> i = Interval()
∅
```
"""
function Interval()
	Interval(0,true,0,true)
end

"""
	left(i)

Returns the left limit of the given `Interval`.
Ex:
```
julia> i = Interval(0,2)
[0,2]
julia> left(i)
0
```
"""
function left(i::Interval)
	i.limits.left
end

"""
	right(i)

Returns the right limit of the given `Interval`.
Ex:
```
julia> i = Interval(0,2)
[0,2]
julia> right(i)
2
```
"""
function right(i::Interval)
	i.limits.right
end

"""
	empty(i)

Returns true if the given `Interval` is empty.
Ex:
```
julia> i = Interval(0,1)
[0,1]
julia> empty(i)
false
julia> j = Interval()
∅
julia> empty(j)
true
```
"""
function empty(i::Interval)
	(left(i) == right(i)) && i.open_left && i.open_right
end

"""
	string(i)

Returns a string representation of an `Interval`.
Ex:
```
julia> i = Interval(0,1)
[0,1]
julia> string(i)
"[0,1]"
```
"""
function string(i::Interval)
	empty(i) && return "∅"
	if i.open_left
        i.open_right ? "]$(left(i)),$(right(i))[" : "]$(left(i)),$(right(i))]"
    else
        i.open_right ? "[$(left(i)),$(right(i))[" : "[$(left(i)),$(right(i))]"
    end
end

"""
	show(io,i)

Prints the string representation of an `Interval`.
"""
function show(io::IO, i::Interval)
	println(string(i))
end

"""
	x ∈ i

Returns true if x ∈ i and false if x ∉ i.
Ex:
```
julia> i = Interval(-1,true,2)
]-1,2]
julia> 0 ∈ i
true
julia> 2 ∈ i
true
julia> -1 ∈ i
false
julia> -0.9999999 ∈ i
true
```
"""
function ∈(x::Real, i::Interval)
    if i.open_left
        i.open_right ? left(i) < x < right(i) : left(i) < x <= right(i)
    else
        i.open_right ? left(i) <= x < right(i) : left(i) <= x <= right(i)
    end
end

"""
	i1 ⊆ i2

Returns true if the interval i1 is included in the interval i2.
Ex:
```
julia> i = Interval(0,3)
[0,3]
julia> j = Interval(1,2)
[1,2]
julia> j ⊆ i
true
julia> k = Interval(0,true,3,true)
]0,3[
julia> k ⊆ i
true
julia> i ⊈ k
true
```
"""
function ⊆(i1::Interval, i2::Interval)
    if i2.open_left
        if i2.open_right
        	if i1.open_left
        		c1 = left(i2) <= left(i1)
        	else
        		c1 = left(i2) < left(i1)
        	end
        	if i1.open_right
        		c2 = right(i1) <= right(i2)
        	else
        		c2 = right(i1) < right(i2)
        	end
        	return c1 && c2
        else
        	if i1.open_left
        		c1 = left(i2) <= left(i1)
        	else
        		c1 = left(i2) < left(i1)
        	end
        	return c1 && (right(i1) <= right(i2))
        end
    else
        if i2.open_right
        	if i1.open_right
        		c2 = right(i1) <= right(i2)
        	else
        		c2 = right(i1) < right(i2)
        	end
        	return (left(i2) <= left(i1)) && c2
        else
        	return (left(i2) <= left(i1)) && (right(i1) <= right(i2))
        end
    end
end

"""
	cardinal(i)

Returns the cardinal of the given `Interval` defined as:
```math
\\left| [a,b] \\right| = b - a
```
**Warning:** An `Interval` can be non empty and have a cardinal of zero (ex: [1,1]).
Ex:
```
julia> i = Interval(1,true,3)
]1,3]
julia> cardinal(i)
2
julia> cardinal(Interval(1,1))
0
```
"""
function cardinal(i::Interval)
	right(i) - left(i)
end

"""
	disjoint(i1,i2)

Returns true if interval i1 and interval i2 are disjoint, false otherwise.
Ex:
```
julia> i = Interval(0,1,true)
[0,1[
julia> j = Interval(1,true,2)
]1,2]
julia> disjoint(i,j)
true
julia> k = Interval(-1,0)
]-1,0]
julia> disjoint(i,k)
false
```
"""
function disjoint(i1::Interval, i2::Interval)
	i1 == i2 && return false
	(empty(i1) || empty(i2)) && return false
	if i1 < i2
		if right(i1) < left(i2)
			return true
		elseif right(i1) == left(i2)
			(i1.open_right || i2.open_left) ? (return true) : (return false)
		else
			return false
		end
	else
		disjoint(i2,i1)
	end
end

"""
	i1 ∩ i2

Returns the `Interval` corresponding to the intersection between `Interval` i1 and `Interval` i2.
If i1 and i2 are disjoint, this function returns an empty `Interval`.
Ex:
```
julia> i = Interval(0,1,true)
[0,1[
julia> j = Interval(0.8,true,1.2)
]0.8,1.2]
julia> i ∩ j
]0.8,1[
```
"""
function ∩(i1::Interval, i2::Interval)
    disjoint(i1,i2) && return Interval()
    (empty(i1) || empty(i2)) && return Interval()

    l::Real = left(i2)
    open_left::Bool = i2.open_left
    r::Real = right(i2)
    open_right::Bool = i2.open_right

    if left(i1) > left(i2)
        l = left(i1)
        open_left = i1.open_left
    elseif left(i1) == left(i2)
        open_left = i1.open_left | i2.open_left
    end

    if right(i1) < right(i2)
        r = right(i1)
        open_right = i1.open_right
    elseif right(i1) == right(i2)
        open_right = i1.open_right | i2.open_right
    end

    l == r ? Interval(OrderedPair(l,r),false,false) : Interval(OrderedPair(l,r),open_left,open_right)
end

"""
	i1 ∪ i2

Returns the `Interval` corresponding to the union between `Interval` i1 and `Interval` i2.
If i1 and i2 are disjoint, this function raises an error since the results wouldn't be an `Interval`.
Please use the `IntervalUnion` type to work with unions of disjoint intervals.
Ex:
```
julia> i = Interval(0,1,true)
[0,1[
julia> j = Interval(0.8,true,1.2)
]0.8,1.2]
julia> i ∪ j
[0,1.2]
```
"""
function ∪(i1::Interval,i2::Interval)
    if disjoint(i1,i2)
    	if right(i1)==left(i2) && (!i1.open_right || !i2.open_left)
    		return Interval(left(i1),i1.open_left,right(i2),i2.open_right)
    	elseif right(i2)==left(i1) && (!i2.open_right || !i1.open_left)
    		return Interval(left(i2),i2.open_left,right(i1),i1.open_right)
    	else
	        throw("Interval $(string(i1)) and $(string(i2)) are disjoints.")
	    end
    end
    l::Real = left(i2)
    open_left::Bool = i2.open_left
    r::Real = right(i2)
    open_right::Bool = i2.open_right

    if left(i1) < left(i2)
        l = left(i1)
        open_left = i1.open_left
    elseif left(i1) == left(i2)
        open_left = i1.open_left & i2.open_left
    end

    if right(i1) > right(i2)
        r = right(i1)
        open_right = i1.open_right
    elseif right(i1) == right(i2)
        open_right = i1.open_right & i2.open_right
    end

    l == r ? Interval(OrderedPair(l,r),false,false) : Interval(OrderedPair(l,r),open_left,open_right)
end

"""
	setdiff(i1,i2)

Setdiff is not implemented for `Interval`s.
Please use `IntervalUnion`s instead.
"""
function setdiff(i1::Interval, i2::Interval)
	throw("Not implemented.")
end

"""
	symdiff(i1,i2)

Symdiff is not implemented for `Interval`s.
Please use `IntervalUnion`s instead.
"""
function symdiff(i1::Interval, i2::Interval)
	throw("Not implemented.")
end

"""
	compact(i)

Returns the compact of the given `Interval`.
Ex:
```
julia> i = Interval(0,true,1)
]0,1]
julia> compact(i)
[0,1]
```
"""
function compact(i::Interval)
	empty(i) && return Interval()
	Interval(left(i),right(i))
end

"""
	compactness(i)

Retuns the compactness of a simple `Interval`, which is always 1 unless it is empty.
Ex:
```
julia> compactness(Interval(0,true,1))
1.0
julia> compactness(Interval())
0.0
```
"""
function compactness(i::Interval)
	empty(i) ? 0.0 : 1.0
end

"""
	superset(i)

Returns the superset of the given `Interval`, which is the `Interval` itself.
Ex:
```
julia> superset(Interval(0,true,1))
]0,1]
```
"""
function superset(i::Interval)
	return i
end

"""
	sample(i)

Returns a number drawn uniformly at random in the given `Interval`.
Ex:
```
julia> i = Interval(2,4,true)
[2,4[
julia> sample(i)
3.935796996224805
```
"""
function sample(i::Interval)
	if cardinal(i) == 0 && (i.open_left && i.open_right)
		throw("Cannot sample from empty Interval.")
	end
	s = rand() * (right(i)-left(i)) + left(i)
	if i.open_left && s == left(i)
		return sample(i)
	elseif i.open_right && s == right(i)
		return sample(i)
	else
		return s
	end
end

"""
	sample(i,n)

Returns an array of n random numbers drawn uniformly in the given `Interval`.
Ex:
```
julia> i = Interval(3,true,6)
]3,6]
julia> sample(i,4)
4-element Array{Float64,1}:
 5.543827369817867
 4.678054798740224
 3.1740822420010355
 3.6870186624440504
```
"""
function sample(i::Interval, n::Int64)
	[sample(i) for x in 1:n]
end

"""
	jaccard(i1,i2)

Returns the jaccard coefficient between two `Interval`s defined as:
```math
\\frac{\\left| A \\cap B \\right|}{\\left| A \\cup B \\right|}
```
Note: If |i1 ∪ i2|=0, this function returns 0.
Ex:
```
julia> i = Interval(3,true,6)
]3,6]
julia> j = Interval(4,7)
[4,7]
julia> jaccard(i,j)
0.5
```
"""
function jaccard(i1::Interval, i2::Interval)
	i1 == i2 && return 1.0
	u = i1 ∪ i2
	cardu = cardinal(u)
	cardu != 0 ? cardinal(i1 ∩ i2) / cardu : 0.0
end

"""
	overlap_coefficient(i1,i2)

Returns the overlap coefficient between two `Interval`s defined as:
```math
\\frac{\\left| A \\cap B \\right|}{min \\left( \\left|A \\right|, \\left|B \\right| \\right)}
```
This means that if setA is a subset of B or the converse then the overlap coefficient is equal to 1.
Note: If `|A|=0` and `|B|=0`, this function returns 0.
Ex:
```
julia> i = Interval(3,true,6)
]3,6]
julia> j = Interval(4,7)
[4,7]
julia> overlap_coefficient(i,j)
0.6666666666666666
```
"""
function overlap_coefficient(i1::Interval, i2::Interval)
	i1 == i2 && return 1.0
	denom = min(cardinal(i1),cardinal(i2))
	denom != 0 ? cardinal(i1 ∩ i2) / denom : 0
end

"""
	dice_coefficient(i1,i2)

Returns the Sørensen–Dice coefficient between two `Interval`s defined as:
```math
\\frac{2 \\left| A \\cap B \\right|}{\\left|A \\right| + \\left|B \\right|}
```
Note: If `|A|+|B|=0`, this function returns 0.
Ex:
```
julia> i = Interval(3,true,6)
]3,6]
julia> j = Interval(4,7)
[4,7]
julia> dice_coefficient(i,j)
0.6666666666666666
```
"""
function dice_coefficient(i1::Interval, i2::Interval)
	i1 == i2 && return 1.0
	denom = cardinal(i1) + cardinal(i2)
	denom != 0 ? 2 * cardinal(i1 ∩ i2) / denom : 0
end

"""
	IntervalUnion

`IntervalUnion` representation.
"""
mutable struct IntervalUnion
	components::Array{Interval,1}
	function IntervalUnion(intervals::Array{Interval,1})
		sorted_intervals = sort(intervals)
		cleaned_intervals = Interval[]
		for interval in sorted_intervals
			if !empty(interval)
				if length(cleaned_intervals) == 0
					push!(cleaned_intervals,interval)
				else
					if !disjoint(cleaned_intervals[end],interval) || ((right(cleaned_intervals[end]) == left(interval)) && (!cleaned_intervals[end].open_right || !interval.open_left))
						cleaned_intervals[end] = cleaned_intervals[end] ∪ interval
					else
						push!(cleaned_intervals,interval)
					end
				end
			end
		end
		new(cleaned_intervals)
	end
end

"""
	IntervalUnion()

Empty union of `Interval`s.
Ex:
```
julia> i = IntervalUnion()
∅
```
"""
function IntervalUnion()
	IntervalUnion(Interval[])
end

"""
	IntervalUnion(left,right)

Simplified constructor for `IntervalUnions` with only one `Interval` closed on both side.
Ex:
```
julia> i = IntervalUnion(2,3)
[2,3]
"""
function IntervalUnion(left::Real, right::Real)
	IntervalUnion([Interval(left,right)])
end

"""
	IntervalUnion(left,open_left,right)

Simplified constructor for `IntervalUnions` with only one `Interval` closed on the right side.
Ex:
```
julia> i = IntervalUnion(2,true,3)
]2,3]
julia> i = IntervalUnion(2,false,3)
[2,3]
"""
function IntervalUnion(left::Real, open_left::Bool, right::Real)
	IntervalUnion([Interval(left,open_left,right)])
end

"""
	IntervalUnion(left,right,open_right)

Simplified constructor for `IntervalUnions` with only one `Interval` closed on the left side.
Ex:
```
julia> i = IntervalUnion(2,3,true)
[2,3[
julia> i = IntervalUnion(2,3,false)
[2,3]
"""
function IntervalUnion(left::Real, right::Real, open_right::Bool)
	IntervalUnion([Interval(left,right,open_right)])
end

"""
	IntervalUnion(left,open_left,right,open_right)

Simplified constructor for `IntervalUnions` with only one `Interval`.
Ex:
```
julia> i = IntervalUnion(2,true,3,true)
]2,3[
julia> i = IntervalUnion(2,false,3,false)
[2,3]
"""
function IntervalUnion(left::Real, open_left::Bool, right::Real, open_right::Bool)
	IntervalUnion([Interval(left,open_left,right,open_right)])
end

"""
	empty(iu)

Returns true if the given `IntervalUnion` is empty.
Ex:
```
julia> empty(IntervalUnion())
true
julia> empty(IntervalUnion([Interval(0,1)]))
false
```
"""
function empty(iu::IntervalUnion)
	(number_of_components(iu) == 0) || all([empty(c) for c in iu.components])
end

"""
	string(iu)

Returns a string representation of an `IntervalUnion`.
Ex:
```
julia> string(IntervalUnion([Interval(0,1,true),Interval(1,true,2)]))
"[0,1[ ∪ ]1,2]"
```
"""
function string(iu::IntervalUnion)
	if empty(iu)
		return "∅"
	elseif number_of_components(iu) == 1
		return string(iu.components[1])
	else
		return reduce((x,y)->"$(string(x)) ∪ $(string(y))", iu.components)
	end
end

"""
	show(io,iu)

Prints the string representation of an `IntervalUnion`.
"""
function show(io::IO, iu::IntervalUnion)
	println(string(iu))
end

"""
	iu1 == iu2

Equality defined for `IntervalUnion`s.
"""
function ==(iu1::IntervalUnion, iu2::IntervalUnion)
	if number_of_components(iu1) != number_of_components(iu2)
		return false
	end
	for (i,j) in zip(iu1.components,iu2.components)
		if i != j
			return false
		end
	end
	true
end

"""
	left(iu)

Returns a vector of left limits of the given `IntervalUnion`.
Ex:
```
julia> i = IntervalUnion([Interval(0,1,true),Interval(1,true,2)])
[0,1[ ∪ ]1,2]
julia> left(i)
2-element Array{Int64,1}:
 0
 1
```
"""
function left(iu::IntervalUnion)
	[left(i) for i in iu.components]
end

"""
	right(iu)

Returns a vector of right limits of the given `IntervalUnion`.
Ex:
```
julia> i = IntervalUnion([Interval(0,1,true),Interval(1,true,2)])
[0,1[ ∪ ]1,2]
julia> left(i)
2-element Array{Int64,1}:
 1
 2
```
"""
function right(iu::IntervalUnion)
	[right(i) for i in iu.components]
end

"""
	limits(iu)

Returns a sorted array of limits of the given `IntervalUnion`.
Ex:
```
julia> i = IntervalUnion([Interval(0,1,true),Interval(1,true,2)])
[0,1[ ∪ ]1,2]
julia> limits(i)
3-element Array{Real,1}:
 0
 1
 2
```
"""
function limits(iu::IntervalUnion)
	limit_set = Set{Real}()
	for i in iu.components
		push!(limit_set,left(i))
		push!(limit_set,right(i))
	end
	limit_array = Real[]
	for l in limit_set
		push!(limit_array,l)
	end
	sort(limit_array)
end

"""
	compact(iu)

Returns the compact of the given `IntervalUnion`.
Ex:
```
julia> i = IntervalUnion([Interval(0,true,1),Interval(2,3)])
]0,1] ∪ [2,3]
julia> compact(i)
[0,3]
```
"""
function compact(iu::IntervalUnion)
	empty(iu) && return IntervalUnion()
	l = minimum(left(iu))
	r = maximum(right(iu))
	if l == -Inf || r == Inf
		throw("IntervalUnion $(string(iu)) has no compact.")
	end
	IntervalUnion([Interval(l,r)])
end

"""
	compactness(iu)

Returns the compactness of the given `IntervalUnion` defined as:
```math
\\frac{ \\left| iu \\right|}{\\left| compact(iu) \\right|}
```
Note: If `|compact(iu)|=0`, this function returns 0.
Ex:
```
julia> i = IntervalUnion([Interval(0,true,1),Interval(3,5)])
]0,1] ∪ [3,5]
julia> compactness(i)
0.6
```
"""
function compactness(iu::IntervalUnion)
	cciu = cardinal(compact(iu))
	cciu != 0 ? cardinal(iu) / cciu : 0
end

"""
	superset(iu)

Returns the superset of the given `IntervalUnion`.
Ex:
```
julia> i = IntervalUnion([Interval(0,true,1),Interval(2,3)])
]0,1] ∪ [2,3]
julia> superset(i)
]0,3]
```
Note: `superset` and `compact` are strongly related but `compact` is always closed on both sides.
"""
function superset(iu::IntervalUnion)
	if empty(iu) || number_of_components(iu) < 2
		throw("IntervalUnion $(string(iu)) has not enough components.")
	end
	smallest = minimum(iu.components)
	largest  = maximum(iu.components)
	IntervalUnion([Interval(left(smallest), smallest.open_left, right(largest), largest.open_right)])
end

"""
	super_complement(iu)

Returns the complement of the superset of the given `IntervalUnion`.
Ex:
```
julia> i = IntervalUnion([Interval(0,true,1),Interval(2,3)])
]0,1] ∪ [2,3]
julia> super_complement(i)
]-Inf,0] ∪ ]3,Inf[
```
"""
function super_complement(iu::IntervalUnion)
	complement(superset(iu))
end

"""
	restricted_complement(iu)

Returns the complement of iu intersected with its superset.
Ex:
```
julia> i = IntervalUnion([Interval(0,1),Interval(4,5)])
[0,1] ∪ [4,5]
julia> restricted_complement(i)
]1,4[
julia> j = IntervalUnion([Interval(0,true,1,true),Interval(2.5,true,3),Interval(3.5,6)])
]0,1[ ∪ ]2.5,3] ∪ [3.5,6]
julia> restricted_complement(j)
[1,2.5] ∪ ]3,3.5[
```
"""
function restricted_complement(iu::IntervalUnion)
	superset(iu) ∩ complement(iu)
end

"""
	complement_cardinal(iu)

Returns the cardinal of the restricted complement of iu.
Ex:
```
julia> IntervalUnion([Interval(0,1),Interval(4,5)])
[0,1] ∪ [4,5]
julia> complement_cardinal(i)
3
julia> j = IntervalUnion([Interval(0,true,1,true),Interval(2.5,true,3),Interval(3.5,6)])
]0,1[ ∪ ]2.5,3] ∪ [3.5,6]
julia> complement_cardinal(j)
2
```
"""
function complement_cardinal(iu::IntervalUnion)
	cardinal(restricted_complement(iu))
end

"""
	x ∈ iu

Returns true if x ∈ iu and false if x ∉ iu.
Ex:
```
julia> i = IntervalUnion([Interval(0,true,1,true),Interval(2.5,true,3),Interval(3.5,6)])
]0,1[ ∪ ]2.5,3] ∪ [3.5,6]
julia> 0.5 ∈ i
true
julia> 1 ∈ i
false
julia> 0.9999999 ∈ i
true
```
"""
function ∈(x::Real, iu::IntervalUnion)
	for i in iu.components
		x ∈ i && return true
	end
	false
end

"""
	cardinal(iu)

Returns the cardinal of the given `IntervalUnion`.
```math
\\left| [a,b] ∪ ]c,d] ∪ [e,f[ ∪ ]g,h[ \\right| = (b - a) + (d-c) + (f-e) + (h-g)
```
Ex:
```
julia> i = IntervalUnion([Interval(0,true,1,true),Interval(2.5,true,3),Interval(3.5,6)])
]0,1[ ∪ ]2.5,3] ∪ [3.5,6]
julia> cardinal(i)
4.0
```
"""
function cardinal(iu::IntervalUnion)
	empty(iu) && return 0
	sum([cardinal(i) for i in iu.components])
end

"""
	number_of_components(iu)

Returns the number of components in the given `IntervalUnion`.
Ex:
```
julia> i = IntervalUnion([Interval(0,true,1,true),Interval(2.5,true,3),Interval(3.5,6)])
]0,1[ ∪ ]2.5,3] ∪ [3.5,6]
julia> number_of_components(i)
3
```
"""
function number_of_components(iu::IntervalUnion)
	length(iu.components)
end

"""
	iu1 ⊆ iu2

Returns true if the `IntervalUnion` iu1 is included in the `IntervalUnion` iu2.
Ex:
```
julia> i = IntervalUnion([Interval(0,true,1,true),Interval(2.5,true,3),Interval(3.5,6)])
]0,1[ ∪ ]2.5,3] ∪ [3.5,6]
julia> j = IntervalUnion([Interval(2.6,3),Interval(3.5,4,true)])
[2.6,3] ∪ [3.5,4[
julia> j ⊆ i
true
```
"""
function ⊆(iu1::IntervalUnion, iu2::IntervalUnion)
    for i in iu1.components
    	!any([i ⊆ j for j in iu2.components]) && return false
    end
    true
end

"""
	iu1 ∪ iu2

Returns the `IntervalUnion` corresponding to the union between
`IntervalUnion` iu1 and `IntervalUnion` iu2.
Ex:
```
julia> i = IntervalUnion([Interval(0,1,true),Interval(1,true,2.3)])
[0,1[ ∪ ]1,2.3]
julia> j = IntervalUnion([Interval(1,1),Interval(2.3,true,4,true)])
[1,1] ∪ ]2.3,4[
julia> i ∪ j
[0,4[
julia> k = IntervalUnion([Interval(-1,0,true),Interval(3,4)])
[-1,0[ ∪ [3,4]
julia> i ∪ k
[-1,1[ ∪ ]1,2.3] ∪ [3,4]
julia> i ∪ IntervalUnion()
[0,1[ ∪ ]1,2.3]
```
"""
function ∪(iu1::IntervalUnion, iu2::IntervalUnion)
	IntervalUnion(vcat(iu1.components,iu2.components))
end

"""
	iu1 ∩ iu2

Returns the `IntervalUnion` corresponding to the intersection between
`IntervalUnion` iu1 and `IntervalUnion` iu2.
Ex:
```
julia> i = IntervalUnion([Interval(0,1,true),Interval(1,true,2.3)])
[0,1[ ∪ ]1,2.3]
julia> j = IntervalUnion([Interval(1,1),Interval(2.3,true,4,true)])
[1,1] ∪ ]2.3,4[
julia> i ∩ j
∅
julia> k = IntervalUnion([Interval(-1,0,true),Interval(3,4)])
[-1,0[ ∪ [3,4]
julia> j ∩ k
[3,4[
julia> i ∩ IntervalUnion()
∅
```
"""
function ∩(iu1::IntervalUnion, iu2::IntervalUnion)
	components = Interval[]
	for c1 in iu1.components
		for c2 in iu2.components
			if !disjoint(c1,c2)
				push!(components,c1 ∩ c2)
			end
		end
	end
	IntervalUnion(components)
end

"""
	complement(iu)

Returns the complement of the given `IntervalUnion`.

Ex:
```
julia> i = IntervalUnion([Interval(0,1)])
[0,1]
julia> complement(i)
]-Inf,0[ ∪ ]1,Inf[
julia> j = IntervalUnion([Interval(0,0)])
[0,0]
julia> complement(j)
]-Inf,0[ ∪ ]0,Inf[
julia> k = IntervalUnion([Interval(0,true,0,true)])
]0,0[
julia> complement(k)
]-Inf,Inf[
julia> complement(IntervalUnion())
]-Inf,Inf[
julia> l = IntervalUnion([Interval(0,1,true),Interval(2,true,3)])
[0,1[ ∪ ]2,3]
julia> complement(l)
]-Inf,0[ ∪ [1,2] ∪ ]3,Inf[
```
"""
function complement(iu::IntervalUnion)
	components = Interval[]
	if empty(iu)
		return IntervalUnion([Interval(-Inf,true,Inf,true)])
	elseif number_of_components(iu) == 1
		if left(iu.components[1]) == -Inf && right(iu.components[1]) == Inf
			return IntervalUnion()
		elseif left(iu.components[1]) == -Inf && right(iu.components[1]) != Inf
			if iu.components[1].open_right
				return IntervalUnion([Interval(right(iu.components[1]),Inf,true)])
			else
				return IntervalUnion([Interval(right(iu.components[1]),true,Inf,true)])
			end
		elseif left(iu.components[1]) != -Inf && right(iu.components[1]) == Inf
			if iu.components[1].open_left
				return IntervalUnion([Interval(-Inf,true,left(iu.components[1]))])
			else
				return IntervalUnion([Interval(-Inf,true,left(iu.components[1]),true)])
			end
		else
			if iu.components[1].open_left
				if iu.components[1].open_right
					return IntervalUnion([Interval(-Inf,true,left(iu.components[1])),Interval(right(iu.components[1]),Inf,true)])
				else
					return IntervalUnion([Interval(-Inf,true,left(iu.components[1])),Interval(right(iu.components[1]),true,Inf,true)])
				end
			else
				if iu.components[1].open_right
					return IntervalUnion([Interval(-Inf,true,left(iu.components[1]),true),Interval(right(iu.components[1]),Inf,true)])
				else
					return IntervalUnion([Interval(-Inf,true,left(iu.components[1]),true),Interval(right(iu.components[1]),true,Inf,true)])
				end
			end
		end
	else
		if left(iu.components[1]) != -Inf
			if iu.components[1].open_left
				push!(components,Interval(-Inf,true,left(iu.components[1])))
			else
				push!(components,Interval(-Inf,true,left(iu.components[1]),true))
			end
		end
		for (c1,c2) in zip(iu.components[1:(end-1)],iu.components[2:end])
			if c1.open_right && c2.open_left
				push!(components,Interval(right(c1),left(c2)))
			elseif !c1.open_right && c2.open_left
				push!(components,Interval(right(c1),true,left(c2)))
			elseif c1.open_right && !c2.open_left
				push!(components,Interval(right(c1),left(c2),true))
			else
				push!(components,Interval(right(c1),true,left(c2),true))
			end
		end
		if right(iu.components[end]) != Inf
			if iu.components[end].open_right
				push!(components,Interval(right(iu.components[end]),Inf,true))
			else
				push!(components,Interval(right(iu.components[end]),true,Inf,true))
			end
		end
	end
	IntervalUnion(components)
end

"""
	setdiff(iu1,iu2)

Returns the `IntervalUnion` with elements in iu1 but not in iu2.
Ex:
```
julia> i = IntervalUnion([Interval(0,1,true),Interval(1,true,2.3)])
[0,1[ ∪ ]1,2.3]
julia> j = IntervalUnion([Interval(1,1),Interval(2.3,true,4,true)])
[1,1] ∪ ]2.3,4[
julia> k = IntervalUnion([Interval(-1,0,true),Interval(3,4)])
[-1,0[ ∪ [3,4]
julia> julia> setdiff(j,k)
[1,1] ∪ ]2.3,3[
```
"""
function setdiff(iu1::IntervalUnion, iu2::IntervalUnion)
	complement(iu2) ∩ iu1
end

"""
	symdiff(iu1,iu2)

Returns the `IntervalUnion` with elements in iu1 or iu2 but not in both.
Ex:
```
julia> i = IntervalUnion([Interval(0,1,true),Interval(1,true,2.3)])
[0,1[ ∪ ]1,2.3]
julia> j = IntervalUnion([Interval(1,1),Interval(2.3,true,4,true)])
[1,1] ∪ ]2.3,4[
julia> k = IntervalUnion([Interval(-1,0,true),Interval(3,4)])
[-1,0[ ∪ [3,4]
julia> symdiff(i,k)
[-1,1[ ∪ ]1,2.3] ∪ [3,4]
julia> symdiff(j,k)
[-1,0[ ∪ [1,1] ∪ ]2.3,3[ ∪ [4,4]
```
"""
function symdiff(iu1::IntervalUnion, iu2::IntervalUnion)
	setdiff(iu1 ∪ iu2, iu1 ∩ iu2)
end

"""
	sample(iu)

Returns a number drawn uniformly at random in the given `IntervalUnion`.
Ex:
```
julia> i = IntervalUnion([Interval(0,1,true),Interval(1,true,2.3)])
[0,1[ ∪ ]1,2.3]
julia> sample(i)
0.22557298488974253
```
"""
function sample(iu::IntervalUnion)
	if empty(iu)
		throw("Cannot sample from empty union of Intervals.")
	end
	size_components = [cardinal(i) for i in iu.components]
	full_size = sum(size_components)
	if full_size == 0
		return left(iu.components[rand(1:number_of_components(iu))])
	end
	size_components_normalized = size_components ./ full_size
	probabilities = cumsum(size_components_normalized)
	r = rand()
	idx_component = findfirst(x->x>r,probabilities)
	sample(iu.components[idx_component])
end

"""
	sample(iu,n)

Returns an array of n random numbers drawn uniformly in the given `IntervalUnion`.
Ex:
```
julia> i = IntervalUnion([Interval(0,1,true),Interval(1,true,2.3)])
[0,1[ ∪ ]1,2.3]
julia> sample(i,4)
4-element Array{Float64,1}:
 0.7974550091944592
 0.36823887828930935
 1.1796141777238074
 2.2731170537623138
```
"""
function sample(iu::IntervalUnion, n::Int64)
	[sample(iu) for x in 1:n]
end

"""
	jaccard(iu1,iu2)

Returns the jaccard similarity between two `IntervalUnion`s.
If |iu1 ∪ iu2|=0, this function returns 0.
Ex:
```
julia> i = IntervalUnion([Interval(0,1,true),Interval(1,true,2.3)])
[0,1[ ∪ ]1,2.3]
julia> j = IntervalUnion([Interval(1,1),Interval(2.3,true,4,true)])
[1,1] ∪ ]2.3,4[
julia> k = IntervalUnion([Interval(-1,0,true),Interval(3,4)])
[-1,0[ ∪ [3,4]
julia> jaccard(i,j)
0.0
julia> jaccard(j,k)
0.37037037037037035
```
"""
function jaccard(iu1::IntervalUnion, iu2::IntervalUnion)
	u = iu1 ∪ iu2
	cardu = cardinal(u)
	cardu != 0 ? cardinal(iu1 ∩ iu2) / cardu : 0.0
end

"""
	overlap_coefficient(iu1,iu2)

Returns the overlap coefficient between two `IntervalUnion`s defined as:
```math
\\frac{\\left| A \\cap B \\right|}{min \\left( \\left|A \\right|, \\left|B \\right| \\right)}
```
This means that if set A is a subset of B or the converse then the overlap coefficient is equal to 1.
Note: If `|A|=0` and `|B|=0`, this function returns 0.
Ex:
```
julia> i = IntervalUnion([Interval(0,1,true),Interval(1,true,2.3)])
[0,1[ ∪ ]1,2.3]
julia> j = IntervalUnion([Interval(1,1),Interval(2.3,true,4,true)])
[1,1] ∪ ]2.3,4[
julia> k = IntervalUnion([Interval(-1,0,true),Interval(3,4)])
[-1,0[ ∪ [3,4]
julia> overlap_coefficient(i,j)
0.0
julia> overlap_coefficient(j,k)
0.588235294117647
```
"""
function overlap_coefficient(iu1::IntervalUnion, iu2::IntervalUnion)
	denom = min(cardinal(iu1),cardinal(iu2))
	denom != 0 ? cardinal(iu1 ∩ iu2) / denom : 0
end

"""
	dice_coefficient(iu1,iu2)

Returns the Sørensen–Dice coefficient between two `IntervalUnion`s defined as:
```math
\\frac{2 \\left| A \\cap B \\right|}{\\left|A \\right| + \\left|B \\right|}
```
Note: If `|A|+|B|=0`, this function returns 0.
Ex:
```
julia> i = IntervalUnion([Interval(0,1,true),Interval(1,true,2.3)])
[0,1[ ∪ ]1,2.3]
julia> j = IntervalUnion([Interval(1,1),Interval(2.3,true,4,true)])
[1,1] ∪ ]2.3,4[
julia> k = IntervalUnion([Interval(-1,0,true),Interval(3,4)])
[-1,0[ ∪ [3,4]
julia> dice_coefficient(i,j)
0.0
julia> dice_coefficient(j,k)
0.5405405405405405
```
"""
function dice_coefficient(iu1::IntervalUnion, iu2::IntervalUnion)
	denom = cardinal(iu1) + cardinal(iu2)
	denom != 0 ? 2 * cardinal(iu1 ∩ iu2) / denom : 0
end


text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

#using Compat

const cublasCreate = cublasCreate_v2
const cublasDestroy = cublasDestroy_v2
const cublasGetVersion = cublasGetVersion_v2
const cublasSetStream = cublasSetStream_v2
const cublasGetStream = cublasGetStream_v2
const cublasGetPointerMode = cublasGetPointerMode_v2
const cublasSetPointerMode = cublasSetPointerMode_v2
const cublasSnrm2 = cublasSnrm2_v2
const cublasDnrm2 = cublasDnrm2_v2
const cublasScnrm2 = cublasScnrm2_v2
const cublasDznrm2 = cublasDznrm2_v2
const cublasSdot = cublasSdot_v2
const cublasDdot = cublasDdot_v2
const cublasCdotu = cublasCdotu_v2
const cublasCdotc = cublasCdotc_v2
const cublasZdotu = cublasZdotu_v2
const cublasZdotc = cublasZdotc_v2
const cublasSscal = cublasSscal_v2
const cublasDscal = cublasDscal_v2
const cublasCscal = cublasCscal_v2
const cublasCsscal = cublasCsscal_v2
const cublasZscal = cublasZscal_v2
const cublasZdscal = cublasZdscal_v2
const cublasSaxpy = cublasSaxpy_v2
const cublasDaxpy = cublasDaxpy_v2
const cublasCaxpy = cublasCaxpy_v2
const cublasZaxpy = cublasZaxpy_v2
const cublasScopy = cublasScopy_v2
const cublasDcopy = cublasDcopy_v2
const cublasCcopy = cublasCcopy_v2
const cublasZcopy = cublasZcopy_v2
const cublasSswap = cublasSswap_v2
const cublasDswap = cublasDswap_v2
const cublasCswap = cublasCswap_v2
const cublasZswap = cublasZswap_v2
const cublasIsamax = cublasIsamax_v2
const cublasIdamax = cublasIdamax_v2
const cublasIcamax = cublasIcamax_v2
const cublasIzamax = cublasIzamax_v2
const cublasIsamin = cublasIsamin_v2
const cublasIdamin = cublasIdamin_v2
const cublasIcamin = cublasIcamin_v2
const cublasIzamin = cublasIzamin_v2
const cublasSasum = cublasSasum_v2
const cublasDasum = cublasDasum_v2
const cublasScasum = cublasScasum_v2
const cublasDzasum = cublasDzasum_v2
const cublasSrot = cublasSrot_v2
const cublasDrot = cublasDrot_v2
const cublasCrot = cublasCrot_v2
const cublasCsrot = cublasCsrot_v2
const cublasZrot = cublasZrot_v2
const cublasZdrot = cublasZdrot_v2
const cublasSrotg = cublasSrotg_v2
const cublasDrotg = cublasDrotg_v2
const cublasCrotg = cublasCrotg_v2
const cublasZrotg = cublasZrotg_v2
const cublasSrotm = cublasSrotm_v2
const cublasDrotm = cublasDrotm_v2
const cublasSrotmg = cublasSrotmg_v2
const cublasDrotmg = cublasDrotmg_v2
const cublasSgemv = cublasSgemv_v2
const cublasDgemv = cublasDgemv_v2
const cublasCgemv = cublasCgemv_v2
const cublasZgemv = cublasZgemv_v2
const cublasSgbmv = cublasSgbmv_v2
const cublasDgbmv = cublasDgbmv_v2
const cublasCgbmv = cublasCgbmv_v2
const cublasZgbmv = cublasZgbmv_v2
const cublasStrmv = cublasStrmv_v2
const cublasDtrmv = cublasDtrmv_v2
const cublasCtrmv = cublasCtrmv_v2
const cublasZtrmv = cublasZtrmv_v2
const cublasStbmv = cublasStbmv_v2
const cublasDtbmv = cublasDtbmv_v2
const cublasCtbmv = cublasCtbmv_v2
const cublasZtbmv = cublasZtbmv_v2
const cublasStpmv = cublasStpmv_v2
const cublasDtpmv = cublasDtpmv_v2
const cublasCtpmv = cublasCtpmv_v2
const cublasZtpmv = cublasZtpmv_v2
const cublasStrsv = cublasStrsv_v2
const cublasDtrsv = cublasDtrsv_v2
const cublasCtrsv = cublasCtrsv_v2
const cublasZtrsv = cublasZtrsv_v2
const cublasStpsv = cublasStpsv_v2
const cublasDtpsv = cublasDtpsv_v2
const cublasCtpsv = cublasCtpsv_v2
const cublasZtpsv = cublasZtpsv_v2
const cublasStbsv = cublasStbsv_v2
const cublasDtbsv = cublasDtbsv_v2
const cublasCtbsv = cublasCtbsv_v2
const cublasZtbsv = cublasZtbsv_v2
const cublasSsymv = cublasSsymv_v2
const cublasDsymv = cublasDsymv_v2
const cublasCsymv = cublasCsymv_v2
const cublasZsymv = cublasZsymv_v2
const cublasChemv = cublasChemv_v2
const cublasZhemv = cublasZhemv_v2
const cublasSsbmv = cublasSsbmv_v2
const cublasDsbmv = cublasDsbmv_v2
const cublasChbmv = cublasChbmv_v2
const cublasZhbmv = cublasZhbmv_v2
const cublasSspmv = cublasSspmv_v2
const cublasDspmv = cublasDspmv_v2
const cublasChpmv = cublasChpmv_v2
const cublasZhpmv = cublasZhpmv_v2
const cublasSger = cublasSger_v2
const cublasDger = cublasDger_v2
const cublasCgeru = cublasCgeru_v2
const cublasCgerc = cublasCgerc_v2
const cublasZgeru = cublasZgeru_v2
const cublasZgerc = cublasZgerc_v2
const cublasSsyr = cublasSsyr_v2
const cublasDsyr = cublasDsyr_v2
const cublasCsyr = cublasCsyr_v2
const cublasZsyr = cublasZsyr_v2
const cublasCher = cublasCher_v2
const cublasZher = cublasZher_v2
const cublasSspr = cublasSspr_v2
const cublasDspr = cublasDspr_v2
const cublasChpr = cublasChpr_v2
const cublasZhpr = cublasZhpr_v2
const cublasSsyr2 = cublasSsyr2_v2
const cublasDsyr2 = cublasDsyr2_v2
const cublasCsyr2 = cublasCsyr2_v2
const cublasZsyr2 = cublasZsyr2_v2
const cublasCher2 = cublasCher2_v2
const cublasZher2 = cublasZher2_v2
const cublasSspr2 = cublasSspr2_v2
const cublasDspr2 = cublasDspr2_v2
const cublasChpr2 = cublasChpr2_v2
const cublasZhpr2 = cublasZhpr2_v2
const cublasSgemm = cublasSgemm_v2
const cublasDgemm = cublasDgemm_v2
const cublasCgemm = cublasCgemm_v2
const cublasZgemm = cublasZgemm_v2
const cublasSsyrk = cublasSsyrk_v2
const cublasDsyrk = cublasDsyrk_v2
const cublasCsyrk = cublasCsyrk_v2
const cublasZsyrk = cublasZsyrk_v2
const cublasCherk = cublasCherk_v2
const cublasZherk = cublasZherk_v2
const cublasSsyr2k = cublasSsyr2k_v2
const cublasDsyr2k = cublasDsyr2k_v2
const cublasCsyr2k = cublasCsyr2k_v2
const cublasZsyr2k = cublasZsyr2k_v2
const cublasCher2k = cublasCher2k_v2
const cublasZher2k = cublasZher2k_v2
const cublasSsymm = cublasSsymm_v2
const cublasDsymm = cublasDsymm_v2
const cublasCsymm = cublasCsymm_v2
const cublasZsymm = cublasZsymm_v2
const cublasChemm = cublasChemm_v2
const cublasZhemm = cublasZhemm_v2
const cublasStrsm = cublasStrsm_v2
const cublasDtrsm = cublasDtrsm_v2
const cublasCtrsm = cublasCtrsm_v2
const cublasZtrsm = cublasZtrsm_v2
const cublasStrmm = cublasStrmm_v2
const cublasDtrmm = cublasDtrmm_v2
const cublasCtrmm = cublasCtrmm_v2
const cublasZtrmm = cublasZtrmm_v2

# begin enum ANONYMOUS_1
typealias ANONYMOUS_1 UInt32
const CUBLAS_STATUS_SUCCESS = (UInt32)(0)
const CUBLAS_STATUS_NOT_INITIALIZED = (UInt32)(1)
const CUBLAS_STATUS_ALLOC_FAILED = (UInt32)(3)
const CUBLAS_STATUS_INVALID_VALUE = (UInt32)(7)
const CUBLAS_STATUS_ARCH_MISMATCH = (UInt32)(8)
const CUBLAS_STATUS_MAPPING_ERROR = (UInt32)(11)
const CUBLAS_STATUS_EXECUTION_FAILED = (UInt32)(13)
const CUBLAS_STATUS_INTERNAL_ERROR = (UInt32)(14)
const CUBLAS_STATUS_NOT_SUPPORTED = (UInt32)(15)
const CUBLAS_STATUS_LICENSE_ERROR = (UInt32)(16)
# end enum ANONYMOUS_1

# begin enum cublasStatus_t
typealias cublasStatus_t UInt32
const CUBLAS_STATUS_SUCCESS = (UInt32)(0)
const CUBLAS_STATUS_NOT_INITIALIZED = (UInt32)(1)
const CUBLAS_STATUS_ALLOC_FAILED = (UInt32)(3)
const CUBLAS_STATUS_INVALID_VALUE = (UInt32)(7)
const CUBLAS_STATUS_ARCH_MISMATCH = (UInt32)(8)
const CUBLAS_STATUS_MAPPING_ERROR = (UInt32)(11)
const CUBLAS_STATUS_EXECUTION_FAILED = (UInt32)(13)
const CUBLAS_STATUS_INTERNAL_ERROR = (UInt32)(14)
const CUBLAS_STATUS_NOT_SUPPORTED = (UInt32)(15)
const CUBLAS_STATUS_LICENSE_ERROR = (UInt32)(16)
# end enum cublasStatus_t

# begin enum ANONYMOUS_2
typealias ANONYMOUS_2 UInt32
const CUBLAS_FILL_MODE_LOWER = (UInt32)(0)
const CUBLAS_FILL_MODE_UPPER = (UInt32)(1)
# end enum ANONYMOUS_2

# begin enum cublasFillMode_t
typealias cublasFillMode_t UInt32
const CUBLAS_FILL_MODE_LOWER = (UInt32)(0)
const CUBLAS_FILL_MODE_UPPER = (UInt32)(1)
# end enum cublasFillMode_t

# begin enum ANONYMOUS_3
typealias ANONYMOUS_3 UInt32
const CUBLAS_DIAG_NON_UNIT = (UInt32)(0)
const CUBLAS_DIAG_UNIT = (UInt32)(1)
# end enum ANONYMOUS_3

# begin enum cublasDiagType_t
typealias cublasDiagType_t UInt32
const CUBLAS_DIAG_NON_UNIT = (UInt32)(0)
const CUBLAS_DIAG_UNIT = (UInt32)(1)
# end enum cublasDiagType_t

# begin enum ANONYMOUS_4
typealias ANONYMOUS_4 UInt32
const CUBLAS_SIDE_LEFT = (UInt32)(0)
const CUBLAS_SIDE_RIGHT = (UInt32)(1)
# end enum ANONYMOUS_4

# begin enum cublasSideMode_t
typealias cublasSideMode_t UInt32
const CUBLAS_SIDE_LEFT = (UInt32)(0)
const CUBLAS_SIDE_RIGHT = (UInt32)(1)
# end enum cublasSideMode_t

# begin enum ANONYMOUS_5
typealias ANONYMOUS_5 UInt32
const CUBLAS_OP_N = (UInt32)(0)
const CUBLAS_OP_T = (UInt32)(1)
const CUBLAS_OP_C = (UInt32)(2)
# end enum ANONYMOUS_5

# begin enum cublasOperation_t
typealias cublasOperation_t UInt32
const CUBLAS_OP_N = (UInt32)(0)
const CUBLAS_OP_T = (UInt32)(1)
const CUBLAS_OP_C = (UInt32)(2)
# end enum cublasOperation_t

# begin enum ANONYMOUS_6
typealias ANONYMOUS_6 UInt32
const CUBLAS_POINTER_MODE_HOST = (UInt32)(0)
const CUBLAS_POINTER_MODE_DEVICE = (UInt32)(1)
# end enum ANONYMOUS_6

# begin enum cublasPointerMode_t
typealias cublasPointerMode_t UInt32
const CUBLAS_POINTER_MODE_HOST = (UInt32)(0)
const CUBLAS_POINTER_MODE_DEVICE = (UInt32)(1)
# end enum cublasPointerMode_t

# begin enum ANONYMOUS_7
typealias ANONYMOUS_7 UInt32
const CUBLAS_ATOMICS_NOT_ALLOWED = (UInt32)(0)
const CUBLAS_ATOMICS_ALLOWED = (UInt32)(1)
# end enum ANONYMOUS_7

# begin enum cublasAtomicsMode_t
typealias cublasAtomicsMode_t UInt32
const CUBLAS_ATOMICS_NOT_ALLOWED = (UInt32)(0)
const CUBLAS_ATOMICS_ALLOWED = (UInt32)(1)
# end enum cublasAtomicsMode_t

# begin enum ANONYMOUS_8
typealias ANONYMOUS_8 UInt32
const CUBLAS_DATA_FLOAT = (UInt32)(0)
const CUBLAS_DATA_DOUBLE = (UInt32)(1)
const CUBLAS_DATA_HALF = (UInt32)(2)
const CUBLAS_DATA_INT8 = (UInt32)(3)
# end enum ANONYMOUS_8

# begin enum cublasDataType_t
typealias cublasDataType_t UInt32
const CUBLAS_DATA_FLOAT = (UInt32)(0)
const CUBLAS_DATA_DOUBLE = (UInt32)(1)
const CUBLAS_DATA_HALF = (UInt32)(2)
const CUBLAS_DATA_INT8 = (UInt32)(3)
# end enum cublasDataType_t

typealias cublasContext Void
typealias cublasHandle_t Ptr{cublasContext}


text:
# ---
# title: 307. Range Sum Query - Mutable
# id: problem307
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Binary Indexed Tree, Segment Tree
# link: <https://leetcode.com/problems/range-sum-query-mutable/description/>
# hidden: true
# ---
#
# Given an integer array _nums_ , find the sum of the elements between indices
# _i_ and _j_ ( _i_ ≤ _j_ ), inclusive.
#
# The _update(i, val)_ function modifies _nums_ by updating the element at index
# _i_ to _val_.
#
# **Example:**
#
#
#
#     Given nums = [1, 3, 5]
#
#     sumRange(0, 2) -> 9
#     update(1, 2)
#     sumRange(0, 2) -> 8
#
#
#
#
# **Constraints:**
#
#   * The array is only modifiable by the _update_ function.
#   * You may assume the number of calls to _update_ and _sumRange_ function is distributed evenly.
#   * `0 <= i <= j <= nums.length - 1`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
    Scene{E}

Container to store a list of entities.
The main difference from a regular array is that its size is defined at construction and is fixed.
(`push!` is O(1))

# Constructors

- `Scene(arr::AbstractVector; capacity::Int=length(arr))`
- `Scene(::Type{E}, capacity::Int=100) where {E}`

# Fields

To interact with `Scene` object it is preferable to use functions rather than accessing the fields directly.

- `entities::Vector{E}`
- `n::Int` current number of entities in the scene
"""
mutable struct Scene{E}
    entities::Vector{E} # NOTE: I tried StaticArrays; was not faster
    n::Int
end
function Scene(arr::AbstractVector{E}; capacity::Int=length(arr)) where {E}
    capacity ≥ length(arr) || error("capacity cannot be less than entitiy count! (N ≥ length(arr))")
    entities = Array{E}(undef, capacity)
    copyto!(entities, arr)
    return Scene{E}(entities, length(arr))
end
function Scene(::Type{E}, capacity::Int=100) where {E}
    entities = Array{E}(undef, capacity)
    return Scene{E}(entities, 0)
end

Base.show(io::IO, scene::Scene{E}) where {E}= @printf(io, "Scene{%s}(%d entities)", string(E), length(scene))

"""
    capacity(scene::Scene)
returns the maximum number of entities that can be put in the scene.
To get the current number of entities use `length` instead.
"""
capacity(scene::Scene) = length(scene.entities)

Base.length(scene::Scene) = scene.n
Base.getindex(scene::Scene, i::Int) = scene.entities[i]
Base.eltype(scene::Scene{E}) where {E} = E

Base.lastindex(scene::Scene) = scene.n
function Base.setindex!(scene::Scene{E}, entity::E, i::Int) where {E}
    scene.entities[i] = entity
    return scene
end
function Base.empty!(scene::Scene)
    scene.n = 0
    return scene
end
function Base.deleteat!(scene::Scene, entity_index::Int)
    for i in entity_index : scene.n - 1
        scene.entities[i] = scene.entities[i+1]
    end
    scene.n -= 1
    scene
end

function Base.iterate(scene::Scene{E}, i::Int=1) where {E}
    if i > length(scene)
        return nothing
    end
    return (scene.entities[i], i+1)
end

function Base.copyto!(dest::Scene{E}, src::Scene{E}) where {E}
    for i in 1 : src.n
        dest.entities[i] = src.entities[i]
    end
    dest.n = src.n
    return dest
end
Base.copy(scene::Scene{E}) where {E} = copyto!(Scene(E, capacity(scene)), scene)

function Base.push!(scene::Scene{E}, entity::E) where {E}
    scene.n += 1
    scene.entities[scene.n] = entity
    return scene
end

####
"""
    EntityScene{S,D,I} = Scene{Entity{S,D,I}}
Alias for `Scene` when the entities in the scene are of type `Entity`

# Constructors
- `EntityScene(::Type{S},::Type{D},::Type{I}) where {S,D,I}`
- `EntityScene(::Type{S},::Type{D},::Type{I},capacity::Int)`

"""
const EntityScene{S,D,I} = Scene{Entity{S,D,I}}
EntityScene(::Type{S},::Type{D},::Type{I}) where {S,D,I} = Scene(Entity{S,D,I})
EntityScene(::Type{S},::Type{D},::Type{I},capacity::Int) where {S,D,I} = Scene(Entity{S,D,I}, capacity)

Base.in(id::I, scene::EntityScene{S,D,I}) where {S,D,I} = findfirst(id, scene) !== nothing

function Base.findfirst(id::I, scene::EntityScene{S,D,I}) where {S,D,I}
    for entity_index in 1 : scene.n
        entity = scene.entities[entity_index]
        if entity.id == id
            return entity_index
        end
    end
    return nothing
end

function id2index(scene::EntityScene{S,D,I}, id::I) where {S,D,I}
    entity_index = findfirst(id, scene)
    if entity_index === nothing
        throw(BoundsError(scene, [id]))
    end
    return entity_index
end

"""
    get_by_id(scene::EntityScene{S,D,I}, id::I) where {S,D,I}
Retrieve the entity by its `id`. This function uses `findfirst` which is O(n).
"""
get_by_id(scene::EntityScene{S,D,I}, id::I) where {S,D,I} = scene[id2index(scene, id)]

function get_first_available_id(scene::EntityScene{S,D,I}) where {S,D,I}
    ids = Set{I}(entity.id for entity in scene)
    id_one = one(I)
    id = id_one
    while id ∈ ids
        id += id_one
    end
    return id
end
function Base.push!(scene::EntityScene{S,D,I}, s::S) where {S,D,I}
    id = get_first_available_id(scene)
    entity = Entity{S,D,I}(s, D(), id)
    push!(scene, entity)
end

Base.delete!(scene::EntityScene{S,D,I}, entity::Entity{S,D,I}) where {S,D,I} = deleteat!(scene, findfirst(entity.id, scene))
function Base.delete!(scene::EntityScene{S,D,I}, id::I) where {S,D,I}
    entity_index = findfirst(id, scene)
    if entity_index != nothing
        deleteat!(scene, entity_index)
    end
    return scene
end

###

function Base.write(io::IO, frames::Vector{EntityScene{S,D,I}}) where {S,D,I}
    println(io, length(frames))
    for scene in frames
        println(io, length(scene))
        for entity in scene
            write(io, entity.state)
            print(io, "\n")
            write(io, entity.def)
            print(io, "\n")
            write(io, string(entity.id))
            print(io, "\n")
        end
    end
end
function Base.read(io::IO, ::Type{Vector{EntityScene{S,D,I}}}) where {S,D,I}

    n = parse(Int, readline(io))
    frames = Array{EntityScene{S,D,I}}(undef, n)

    for i in 1 : n
        m = parse(Int, readline(io))
        scene = Scene(Entity{S,D,I}, m)
        for j in 1 : m
            state = read(io, S)
            def = read(io, D)
            strid = readline(io)
            id = try
                parse(I, strid)
            catch e
                if e isa MethodError
                    id = try
                        convert(I, strid)
                    catch e
                        error("Cannot parse $strid as $I")
                    end
                else
                    error("Cannot parse $strid as $I")
                end
            end
            push!(scene, Entity(state,def,id))
        end
        frames[i] = scene
    end

    return frames
end


text:
"""Options for the `modify_predictor_fn` keyword argument when defining a `Problem`."""

function append_tke(x, time_index, state_variables)
    vcat(x, state_variables.tke_avg[time_index])
end

function partial_temp_profile(z_set)
    function f(x, time_index, 𝒟)
        x[z_set]
    end
    f
end

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection