text:
using DataFrames

struct DFBindInfo
    df::DataFrame
    input_columns::Vector
    scan_types::Vector{Type}
    result_types::Vector{Type}
    scan_functions::Vector{Function}

    function DFBindInfo(
        df::DataFrame,
        input_columns::Vector,
        scan_types::Vector{Type},
        result_types::Vector{Type},
        scan_functions::Vector{Function}
    )
        return new(df, input_columns, scan_types, result_types, scan_functions)
    end
end

function df_input_type(df, entry)
    return eltype(df[!, entry])
end

function df_result_type(df, entry)
    column_type = df_input_type(df, entry)
    if typeof(column_type) == Union
        # remove Missing type from the union
        column_type = Core.Compiler.typesubtract(column_type, Missing, 1)
    end
    return column_type
end

function df_julia_type(column_type)
    if column_type == Date
        column_type = Int32
    elseif column_type == Time
        column_type = Int64
    elseif column_type == DateTime
        column_type = Int64
    end
    return column_type
end

value_to_duckdb(val::T) where {T <: Date} =
    convert(Int32, Dates.date2epochdays(val) - ROUNDING_EPOCH_TO_UNIX_EPOCH_DAYS)
value_to_duckdb(val::T) where {T <: Time} = convert(Int64, Dates.value(val) / 1000)
value_to_duckdb(val::T) where {T <: DateTime} =
    convert(Int64, (Dates.datetime2epochms(val) - ROUNDING_EPOCH_TO_UNIX_EPOCH_MS) * 1000)
function value_to_duckdb(val::T) where {T <: AbstractString}
    throw(
        NotImplementedException(
            "Cannot use value_to_duckdb to convert string values - use DuckDB.assign_string_element on a vector instead"
        )
    )
end
function value_to_duckdb(val::T) where {T}
    return val
end

function df_scan_column(
    input_column::AbstractVector{DF_TYPE},
    df_offset::Int64,
    col_idx::Int64,
    result_idx::Int64,
    scan_count::Int64,
    output::DuckDB.DataChunk,
    ::Type{DUCK_TYPE},
    ::Type{DF_TYPE}
) where {DUCK_TYPE, DF_TYPE}
    vector::Vec = DuckDB.get_vector(output, result_idx)
    result_array::Vector{DUCK_TYPE} = DuckDB.get_array(vector, DUCK_TYPE)
    validity::ValidityMask = DuckDB.get_validity(vector)
    for i::Int64 in 1:scan_count
        val = getindex(input_column, df_offset + i)
        if val === missing
            DuckDB.setinvalid(validity, i)
        else
            result_array[i] = value_to_duckdb(val)
        end
    end
end

function df_scan_string_column(
    input_column::AbstractVector{DF_TYPE},
    df_offset::Int64,
    col_idx::Int64,
    result_idx::Int64,
    scan_count::Int64,
    output::DuckDB.DataChunk,
    ::Type{DUCK_TYPE},
    ::Type{DF_TYPE}
) where {DUCK_TYPE, DF_TYPE}
    vector::Vec = DuckDB.get_vector(output, result_idx)
    validity::ValidityMask = DuckDB.get_validity(vector)
    for i::Int64 in 1:scan_count
        val = getindex(input_column, df_offset + i)
        if val === missing
            DuckDB.setinvalid(validity, i)
        else
            DuckDB.assign_string_element(vector, i, val)
        end
    end
end

function df_scan_function(df, entry)
    result_type = df_result_type(df, entry)
    if result_type <: AbstractString
        return df_scan_string_column
    end
    return df_scan_column
end

function df_bind_function(info::DuckDB.BindInfo)
    # fetch the df name from the function parameters
    parameter = DuckDB.get_parameter(info, 0)
    name = DuckDB.getvalue(parameter, String)
    # fetch the actual df using the function name
    extra_data = DuckDB.get_extra_data(info)
    df = extra_data[name]

    # register the result columns
    input_columns = Vector()
    scan_types::Vector{Type} = Vector()
    result_types::Vector{Type} = Vector()
    scan_functions::Vector{Function} = Vector()
    for entry in names(df)
        result_type = df_result_type(df, entry)
        scan_function = df_scan_function(df, entry)
        push!(input_columns, df[!, entry])
        push!(scan_types, df_input_type(df, entry))
        push!(result_types, df_julia_type(result_type))
        push!(scan_functions, scan_function)

        DuckDB.add_result_column(info, entry, result_type)
    end
    return DFBindInfo(df, input_columns, scan_types, result_types, scan_functions)
end

mutable struct DFGlobalInfo
    pos::Int64
    global_lock::ReentrantLock

    function DFGlobalInfo()
        return new(0, ReentrantLock())
    end
end

mutable struct DFLocalInfo
    columns::Vector{Int64}
    current_pos::Int64
    end_pos::Int64

    function DFLocalInfo(columns)
        return new(columns, 0, 0)
    end
end

function df_global_init_function(info::DuckDB.InitInfo)
    bind_info = DuckDB.get_bind_info(info, DFBindInfo)
    # figure out the maximum number of threads to launch from the DF size
    row_count::Int64 = size(bind_info.df, 1)
    max_threads::Int64 = ceil(row_count / DuckDB.ROW_GROUP_SIZE)
    DuckDB.set_max_threads(info, max_threads)
    return DFGlobalInfo()
end

function df_local_init_function(info::DuckDB.InitInfo)
    columns = DuckDB.get_projected_columns(info)
    return DFLocalInfo(columns)
end

function df_scan_function(info::DuckDB.FunctionInfo, output::DuckDB.DataChunk)
    bind_info = DuckDB.get_bind_info(info, DFBindInfo)
    global_info = DuckDB.get_init_info(info, DFGlobalInfo)
    local_info = DuckDB.get_local_info(info, DFLocalInfo)

    if local_info.current_pos >= local_info.end_pos
        # ran out of data to scan in the local info: fetch new rows from the global state (if any)
        # we can in increments of 100 vectors
        lock(global_info.global_lock)
        row_count::Int64 = size(bind_info.df, 1)
        local_info.current_pos = global_info.pos
        total_scan_amount::Int64 = DuckDB.ROW_GROUP_SIZE
        if local_info.current_pos + total_scan_amount >= row_count
            total_scan_amount = row_count - local_info.current_pos
        end
        local_info.end_pos = local_info.current_pos + total_scan_amount
        global_info.pos += total_scan_amount
        unlock(global_info.global_lock)
    end
    scan_count::Int64 = DuckDB.VECTOR_SIZE
    current_row::Int64 = local_info.current_pos
    if current_row + scan_count >= local_info.end_pos
        scan_count = local_info.end_pos - current_row
    end
    local_info.current_pos += scan_count

    result_idx::Int64 = 1
    for col_idx::Int64 in local_info.columns
        if col_idx == 0
            result_idx += 1
            continue
        end
        bind_info.scan_functions[col_idx](
            bind_info.input_columns[col_idx],
            current_row,
            col_idx,
            result_idx,
            scan_count,
            output,
            bind_info.result_types[col_idx],
            bind_info.scan_types[col_idx]
        )
        result_idx += 1
    end
    DuckDB.set_size(output, scan_count)
    return
end

function register_data_frame(con::Connection, df::DataFrame, name::AbstractString)
    con.db.registered_objects[name] = df
    DBInterface.execute(
        con,
        string("CREATE OR REPLACE VIEW \"", name, "\" AS SELECT * FROM julia_df_scan('", name, "')")
    )
    return
end
register_data_frame(db::DB, df::DataFrame, name::AbstractString) = register_data_frame(db.main_connection, df, name)

function unregister_data_frame(con::Connection, name::AbstractString)
    pop!(con.db.registered_objects, name)
    DBInterface.execute(con, string("DROP VIEW IF EXISTS \"", name, "\""))
    return
end
unregister_data_frame(db::DB, name::AbstractString) = unregister_data_frame(db.main_connection, name)

function _add_data_frame_scan(db::DB)
    # add the data frame scan function
    DuckDB.create_table_function(
        db.main_connection,
        "julia_df_scan",
        [String],
        df_bind_function,
        df_global_init_function,
        df_scan_function,
        db.handle.registered_objects,
        true,
        df_local_init_function
    )
    return
end

max_stars_repo_name: taniabogatsch/duckdb

max_issues_repo_name: taniabogatsch/duckdb

max_forks_repo_name: taniabogatsch/duckdb


text:
struct IndexableGenerator{I,F}
    g::Base.Generator{I,F}
end

IndexableGenerator(args...) = IndexableGenerator(Base.Generator(args...))
Base.iterate(ig::IndexableGenerator, s...) = iterate(ig.g, s...)

Base.length(ig::IndexableGenerator) = length(ig.g)
Base.size(ig::IndexableGenerator) = size(ig.g)
Base.axes(ig::IndexableGenerator) = axes(ig.g)
Base.ndims(ig::IndexableGenerator) = ndims(ig.g)

Base.getindex(ig::IndexableGenerator, i::Integer) = ig.g.f(ig.g.iter[i])
Base.lastindex(ig::IndexableGenerator) = lastindex(ig.g.iter)

# Required for compatibility with: https://github.com/JuliaLang/julia/pull/42991
Base.last(ig::IndexableGenerator) = ig[end]

ext: jl


text:
using Reexport
@reexport module Embeddings

import Simplices.Delaunay.delaunay

using StaticArrays
using Distributions
using Statistics
using NearestNeighbors

# Abstract embedding type
include("AbstractEmbedding.jl")

# Composite embedding types holding different types of data depending on
# the application.
include("composite_types/Embedding.jl")
include("composite_types/SimpleEmbedding.jl")
include("composite_types/LinearlyInvariantEmbedding.jl")

# Invariantize embeddings (ensure last point lies within the convex hull
# of preceding points).
include("invariantizing/invariantize.jl")

# Interface with NearestNeighbors, allowing to create trees from embeddings.
include("NearestNeighbors/nearestneighbors.jl")

# Functions to perform delay embeddings
include("embed_functions/customembed.jl")

end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using QHull

function GLLar(points::Array{Float64,2}, cells::Array{Array{Int64, 1}, 1})

	points = convert(Array{Float64,2}, points')
	vertices=Vector{Float32}()
	normals =Vector{Float32}()
	for cell in cells

		ch = QHull.chull(points[cell,:])
		verts = ch.vertices
		vdict = Dict(zip(verts, 1:length(verts)))
		hfaces = [[vdict[u],vdict[v],vdict[w]] for (u,v,w) in ch.simplices]
		hdict = Dict(zip(1:length(cell), cell))
		trias = [[hdict[t1],hdict[t2],hdict[t3]] for (t1,t2,t3) in hfaces]
		pts = points[verts,:]

		for tria in trias
			p3,p2,p1 = points[tria[1],:],points[tria[2],:],points[tria[3],:]
			p1 = convert(Point3d, p1)
			p2 = convert(Point3d, p2)
			p3 = convert(Point3d, p3)
			n=0.5*computeNormal(p1::Point3d,p2::Point3d,p3::Point3d)

			append!(vertices,p1); append!(normals,n)
			append!(vertices,p2); append!(normals,n)
			append!(vertices,p3); append!(normals,n)
		end
	end

	ret=GLMesh(GL_TRIANGLES)
	ret.vertices = GLVertexBuffer(vertices)
	ret.normals  = GLVertexBuffer(normals)
	return ret
end

V,CV = Lar.cuboidGrid([20,10,1])
VIEW([
      #GLCuboid(Box3d(Point3d(0,0,0),Point3d(1,1,1)))
      GLLar(V,CV)
      GLAxis(Point3d(0,0,0),Point3d(1,1,1))
])


text:
using SparseGrids
using Test

@testset "combvec" begin
    @testset "Two 2D vectors" begin
        a1 = [[1; 3], [2; 4]]
        a = combvec(a1)

        @test a == [[1; 2], [3; 2], [1; 4], [3; 4]]
    end

    @testset "A 2D vector and two 1D vectors" begin
        a1 = [[1; 2], [3], [4]]
        a = combvec(a1)

        @test a == [[1; 3; 4], [2; 3; 4]]
    end
end


text:
# automatically generated by the FlatBuffers compiler, do not modify

# module: MyGame

module MyGame
    include("InParentNamespace.jl")
    include("Example2/Example2.jl")
    include("Example/Example.jl")
end


text:
export BrickletRemoteSwitchV2RemoteConfiguration
struct BrickletRemoteSwitchV2RemoteConfiguration
    remote_type::Integer
    minimum_repeats::Integer
    callback_enabled::Bool
end

export BrickletRemoteSwitchV2RemoteStatusA
struct BrickletRemoteSwitchV2RemoteStatusA
    house_code::Integer
    receiver_code::Integer
    switch_to::Integer
    repeats::Integer
end

export BrickletRemoteSwitchV2RemoteStatusB
struct BrickletRemoteSwitchV2RemoteStatusB
    address::Integer
    unit::Integer
    switch_to::Integer
    dim_value::Integer
    repeats::Integer
end

export BrickletRemoteSwitchV2RemoteStatusC
struct BrickletRemoteSwitchV2RemoteStatusC
    system_code::Char
    device_code::Integer
    switch_to::Integer
    repeats::Integer
end

export BrickletRemoteSwitchV2SPITFPErrorCount
struct BrickletRemoteSwitchV2SPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletRemoteSwitchV2Identity
struct BrickletRemoteSwitchV2Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletRemoteSwitchV2
"""
Controls remote mains switches and receives signals from remotes
"""
mutable struct BrickletRemoteSwitchV2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletRemoteSwitchV2(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_remote_switch_v2")
        deviceInternal = package.BrickletRemoteSwitchV2(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_switching_state
"""
    $(SIGNATURES)

Returns the current switching state. If the current state is busy, the
Bricklet is currently sending a code to switch a socket. It will not
accept any calls of switch socket functions until the state changes to ready.

How long the switching takes is dependent on the number of repeats, see
:func:`Set Repeats`.
"""
function get_switching_state(device::BrickletRemoteSwitchV2)
    return device.deviceInternal.get_switching_state()
end

export set_repeats
"""
    $(SIGNATURES)

Sets the number of times the code is sent when one of the Switch Socket
functions is called. The repeats basically correspond to the amount of time
that a button of the remote is pressed.

Some dimmers are controlled by the length of a button pressed,
this can be simulated by increasing the repeats.
"""
function set_repeats(device::BrickletRemoteSwitchV2, repeats)
    device.deviceInternal.set_repeats(repeats)
end

export get_repeats
"""
    $(SIGNATURES)

Returns the number of repeats as set by :func:`Set Repeats`.
"""
function get_repeats(device::BrickletRemoteSwitchV2)
    return device.deviceInternal.get_repeats()
end

export switch_socket_a
"""
    $(SIGNATURES)

To switch a type A socket you have to give the house code, receiver code and the
state (on or off) you want to switch to.

A detailed description on how you can figure out the house and receiver code
can be found :ref:`here <remote_switch_bricklet_type_a_house_and_receiver_code>`.
"""
function switch_socket_a(device::BrickletRemoteSwitchV2, house_code, receiver_code, switch_to)
    device.deviceInternal.switch_socket_a(house_code, receiver_code, switch_to)
end

export switch_socket_b
"""
    $(SIGNATURES)

To switch a type B socket you have to give the address, unit and the state
(on or off) you want to switch to.

To switch all devices with the same address use 255 for the unit.

A detailed description on how you can teach a socket the address and unit can
be found :ref:`here <remote_switch_bricklet_type_b_address_and_unit>`.
"""
function switch_socket_b(device::BrickletRemoteSwitchV2, address, unit, switch_to)
    device.deviceInternal.switch_socket_b(address, unit, switch_to)
end

export dim_socket_b
"""
    $(SIGNATURES)

To control a type B dimmer you have to give the address, unit and the
dim value you want to set the dimmer to.

A detailed description on how you can teach a dimmer the address and unit can
be found :ref:`here <remote_switch_bricklet_type_b_address_and_unit>`.
"""
function dim_socket_b(device::BrickletRemoteSwitchV2, address, unit, dim_value)
    device.deviceInternal.dim_socket_b(address, unit, dim_value)
end

export switch_socket_c
"""
    $(SIGNATURES)

To switch a type C socket you have to give the system code, device code and the
state (on or off) you want to switch to.

A detailed description on how you can figure out the system and device code
can be found :ref:`here <remote_switch_bricklet_type_c_system_and_device_code>`.
"""
function switch_socket_c(device::BrickletRemoteSwitchV2, system_code, device_code, switch_to)
    device.deviceInternal.switch_socket_c(system_code, device_code, switch_to)
end

export set_remote_configuration
"""
    $(SIGNATURES)

Sets the configuration for **receiving** data from a remote of type A, B or C.

* Remote Type: A, B or C depending on the type of remote you want to receive.
* Minimum Repeats: The minimum number of repeated data packets until the callback
  is triggered (if enabled).
* Callback Enabled: Enable or disable callback (see :cb:`Remote Status A` callback,
  :cb:`Remote Status B` callback and :cb:`Remote Status C` callback).
"""
function set_remote_configuration(device::BrickletRemoteSwitchV2, remote_type, minimum_repeats, callback_enabled)
    device.deviceInternal.set_remote_configuration(remote_type, minimum_repeats, callback_enabled)
end

export get_remote_configuration
"""
    $(SIGNATURES)

Returns the remote configuration as set by :func:`Set Remote Configuration`
"""
function get_remote_configuration(device::BrickletRemoteSwitchV2)
    return device.deviceInternal.get_remote_configuration()
end

export get_remote_status_a
"""
    $(SIGNATURES)

Returns the house code, receiver code, switch state (on/off) and number of
repeats for remote type A.

Repeats == 0 means there was no button press. Repeats >= 1 means there
was a button press with the specified house/receiver code. The repeats are the
number of received identical data packets. The longer the button is pressed,
the higher the repeat number.

Use the callback to get this data automatically when a button is pressed,
see :func:`Set Remote Configuration` and :cb:`Remote Status A` callback.
"""
function get_remote_status_a(device::BrickletRemoteSwitchV2)
    return device.deviceInternal.get_remote_status_a()
end

export get_remote_status_b
"""
    $(SIGNATURES)

Returns the address (unique per remote), unit (button number), switch state
(on/off) and number of repeats for remote type B.

If the remote supports dimming the dim value is used instead of the switch state.

If repeats=0 there was no button press. If repeats >= 1 there
was a button press with the specified address/unit. The repeats are the number of received
identical data packets. The longer the button is pressed, the higher the repeat number.

Use the callback to get this data automatically when a button is pressed,
see :func:`Set Remote Configuration` and :cb:`Remote Status B` callback.
"""
function get_remote_status_b(device::BrickletRemoteSwitchV2)
    return device.deviceInternal.get_remote_status_b()
end

export get_remote_status_c
"""
    $(SIGNATURES)

Returns the system code, device code, switch state (on/off) and number of repeats for
remote type C.

If repeats=0 there was no button press. If repeats >= 1 there
was a button press with the specified system/device code. The repeats are the number of received
identical data packets. The longer the button is pressed, the higher the repeat number.

Use the callback to get this data automatically when a button is pressed,
see :func:`Set Remote Configuration` and :cb:`Remote Status C` callback.
"""
function get_remote_status_c(device::BrickletRemoteSwitchV2)
    return device.deviceInternal.get_remote_status_c()
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletRemoteSwitchV2)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletRemoteSwitchV2, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletRemoteSwitchV2)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletRemoteSwitchV2, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletRemoteSwitchV2, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletRemoteSwitchV2, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletRemoteSwitchV2)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletRemoteSwitchV2)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletRemoteSwitchV2)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletRemoteSwitchV2, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletRemoteSwitchV2)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletRemoteSwitchV2)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletRemoteSwitchV2, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
function gputest(f, xs...; checkgrad=true, atol=1e-10)
    cpu_in = xs
    gpu_in = CuArray.(xs)

    cpu_out = f(cpu_in...)
    gpu_out = f(gpu_in...)
    @test collect(cpu_out) ≈ collect(gpu_out)

    if checkgrad
        cpu_grad = gradient((x...) -> sum(f(x...)), cpu_in...)
        gpu_grad = gradient((x...) -> sum(f(x...)), gpu_in...)
        for (cpu_g, gpu_g) in zip(cpu_grad, gpu_grad)
            if cpu_g === nothing
                @test gpu_g === nothing
            else
                @test collect(cpu_g) ≈ collect(gpu_g)  atol=atol
            end
        end
    end
end


text:
using Random: randn!

struct RandomFriction <: ElectronicFrictionProvider
    ndofs::Int
end

function friction!(::RandomFriction, F::AbstractMatrix, ::AbstractMatrix)
    randn!(F)
    F .= F'F
    F .= (F + F')/2
end

NQCModels.ndofs(model::RandomFriction) = model.ndofs


text:
function split_into_cols(df::AbstractDataFrame, col)
    new_cols = map(anymod_split, df[!, col])
    cols_len = map(length, new_cols)
    return new_cols, cols_len
end

function convert_splitted_cols_to_tuple_of_arrays(col, lens)
    return Tuple(
        map(x -> i <= lens[x] ? col[x][i] : col[x][lens[x]], 1:length(col)) for
        i in 1:maximum(lens)
    )
end

function expand_col!(df::T, col) where {T<:AbstractDataFrame}
    cols, lens = split_into_cols(df, col)
    tup_of_arr = convert_splitted_cols_to_tuple_of_arrays(cols, lens)
    for (i, c) in enumerate(tup_of_arr)
        df[!, col * "_$i"] = c
    end
    return df
end

anymod_split(x::AbstractString) = split(x, " < ")
anymod_split(x::AbstractString, i::Int) = split(x, " < ")[i]
function anymod_split(x)
    try
        y = string(x)
        return anymod_split(y)
    catch
        @error "Could not convert $x into an string."
    end
end


text:
# The Computer Language Benchmarks Game
# https://salsa.debian.org/benchmarksgame-team/benchmarksgame/

# contributed by Jarret Revels
# based on Mario Pernici Python's program
# modified by hanabi1224

using Printf

function pidigits(N::Int)
    i = k = 0
    k1 = 1
    n,a,d,t,u = map(BigInt,(1,0,1,0,0))

    while true
        k += 1
        t = n << 1
        n *= k
        a += t
        k1 += 2
        a *= k1
        d *= k1

        if a >= n
            t,u = divrem(n*3 +a, d)
            u += n
            if d > u
                i += 1
                print(t)
                mod_i = mod(i,10)
                if mod_i == 0
                    @printf("\t:%d\n", i)
                end
                if i >= N
                    if mod_i > 0
                        for n in 1:(10 - mod_i)
                            print(' ')
                        end
                        @printf("\t:%d\n", i)
                    end
                    return
                end
                a -= d * t
                a *= 10
                n *= 10
            end
        end
    end
end

function real_main()
    n = parse(Int,ARGS[1])
    pidigits(n)
end

if abspath(PROGRAM_FILE) == @__FILE__
    real_main()
end


text:
"""Utility functions"""

# """
# Compute the squared L2 norm.
# """
# norm²(x) = sum(abs2, x)
"""
Spectral flattening.
"""
function spectrumflatten(x::AbstractArray{T,1}, Nnorm::Int) where T <:Real
    if Nnorm >= length(x)
        xfilt = x.-median(x)
        xfilt[xfilt.<0] .= 0
        return xfilt
    end
    M = map(median, Subsequence(x, Nnorm, Nnorm-1))
    xfilt = x.-M
    xfilt[xfilt.<0] .= 0
    xfilt
end
function spectrumflatten(x::AbstractArray{T,2}, Nnorm::Int; dims::Int=2) where T <: Real
    mapslices(v->spectrumflatten(v, Nnorm), x, dims=dims)
end

"""
Get myriad constant given α and scale.
"""
function myriadconstant(α, scale)
    (α/(2-α+eps()))*(scale^2)
end

"""
Get myriad constant given `x`.
"""
function myriadconstant(x::AbstractVector{T}) where T<:Real
    d = fit(AlphaStable, x)
    myriadconstant(d.α, d.scale)
end

# """
# Get vector myriad constants given α and Rₘ.
# """
# vmyriadconstant(α::T, Rₘ::AbstractMatrix{T}) where T = (α/(2-α+eps())), Rₘ

# """
# Get vector myriad constants given `x`.
# """
# function vmyriadconstant(x::AbstractVector{T}, m::Integer=4) where T<:Real
#     d = fit(AlphaSubGaussian, x, m)
#     vmyriadconstant(d.α, d.R)
# end

"""
Convert a real signal `x` to an acoustic pressure signal in micropascal.
"""
function pressure(x::AbstractVector{T}, sensitivity::T, gain::T; voltparams::Union{Nothing, Tuple{Int, T}}=nothing) where T<:Real
    ν = exp10(sensitivity/20)
    G = exp10(gain/20)
    if voltparams !== nothing
        nbits, vref = voltparams
        x .*= vref/(2^(nbits-1))
    end
    x./(ν*G)
end

"""
Generate a Hilbert envelope of a real signal `x`.
"""
function envelope(x::AbstractVector{T}) where T<:Real
    abs.(hilbert(x))
end

"""
Get the normalized envelope of of a real signal `x`.
"""
function normalize_envelope(x::AbstractVector{T}) where T<:Real
    env = envelope(x)
    env/sum(env)
end

"""
Get the normalized spectrum of a real signal `x`.
"""
function normalize_spectrum(s::AbstractMatrix{T}) where T<:Real
    sf = sum(s, dims=2)
    sf/sum(sf)
end

"""
Compute ordinal patterns of a real signal `x`.
"""
function ordinalpatterns(x::AbstractVector{T}, m::Integer, τ::Integer=1) where T<:Real
    n = length(x) - τ*m + 1
    ps = Vector{Int}[]
    counts = Int[]
    @inbounds for t in 1:n
        s = @view x[t:τ:t+τ*(m-1)]
        p = sortperm(s)
#        v = [p == ptmp for ptmp in ps]
        cntindex = nothing
        for (i, ptmp) in enumerate(ps)
            if p == ptmp
                cntindex = i
                break
            end
        end
        if cntindex !== nothing #any(v)
            counts[cntindex] += 1
        else
            push!(ps, p)
            push!(counts, 1)
        end
    end
    counts ./ sum(counts)
end


text:
fatalerrors = length(ARGS) > 0 && ARGS[1] == "-f"
quiet = length(ARGS) > 0 && ARGS[1] == "-q"
anyerrors = false

my_tests = (
    "query/node.jl",
    # "query/query.jl",
    "collect.jl",
    "lift.jl"
)

println("Running tests:")

using Base.Test

@testset "All tests" begin
    for my_test in my_tests
        try
            include(my_test)
            println("\t\033[1m\033[32mPASSED\033[0m: $(my_test)")
        catch e
            anyerrors = true
            println("\t\033[1m\033[31mFAILED\033[0m: $(my_test)")
            if fatalerrors
                rethrow(e)
            elseif !quiet
                showerror(STDOUT, e, backtrace())
                println()
            end
        end
    end
end

if anyerrors
    throw("Tests failed")
end


text:
using Test
import Random
import Distributions
import StatsBase

import MatrixCompletion.Utilities.Indexing:_check_bernoulli, _check_gaussian,_check_poisson

bernoulliMatrixTest1 = Random.rand(Distributions.Bernoulli(0.5),100);
nonBernoulliMatrixTest1 = Random.rand(Distributions.Uniform(1,10),100);
poissonMatrixTest1 = Random.rand(Distributions.Poisson(3),100);
nonPoissonMatrixTest1 = Random.rand(Distributions.Uniform(1,10),100);
poissonBernoulliDifferenceTest1 = Random.rand(Distributions.Poisson(3),100);
poissonBernoulliDifferenceTest2 = Random.rand(Distributions.Bernoulli(0.4),100);

@test _check_bernoulli(bernoulliMatrixTest1) == true
@test _check_bernoulli(nonBernoulliMatrixTest1) == false
@test _check_poisson(poissonMatrixTest1) == true
@test _check_poisson(nonBernoulliMatrixTest1) == false
@test _check_bernoulli(poissonBernoulliDifferenceTest1) == false
@test _check_poisson(poissonBernoulliDifferenceTest1) == true
@test _check_poisson(poissonBernoulliDifferenceTest2) == false
@test _check_bernoulli(poissonBernoulliDifferenceTest2) == true

import MatrixCompletion.Utilities.Indexing:construct_type_matrix,construct_index_tracker,DIST_FLAGS
import MatrixCompletion.Utilities.Indexing:Bernoulli,
                                           Poisson,
                                           Gaussian,
                                           Gamma,
                                           NegativeBinomial

# test 1: bernoulli only
type_matrix_test_input1 = Random.rand(Distributions.Bernoulli(0.5),5,5)
type_matrix_test_expected1 = Array{Union{DIST_FLAGS,Missing}}(undef,5,5)
fill!(type_matrix_test_expected1,Bernoulli)
@test construct_type_matrix(type_matrix_test_input1) ==  type_matrix_test_expected1

construct_type_matrix(type_matrix_test_input1)
# test 2: poisson only
type_matrix_test_input2 = Random.rand(Distributions.Poisson(5),5,5)
type_matrix_test_expected2 = Array{Union{DIST_FLAGS,Missing}}(undef,5,5)
fill!(type_matrix_test_expected2,Poisson)
@test construct_type_matrix(type_matrix_test_input2) == type_matrix_test_expected2

# test 3: gaussian only
type_matrix_test_input3 = Random.rand(Distributions.Gaussian(0,1),5,5)
type_matrix_test_expected3 = Array{Union{DIST_FLAGS,Missing}}(undef,5,5)
fill!(type_matrix_test_expected3,Gaussian)
@test construct_type_matrix(type_matrix_test_input3) == type_matrix_test_expected3

import MatrixCompletion.Utilities.RandomMatrices.rand
# test4
type_matrix_test_input4 = rand([(Distributions.Bernoulli(0.5),50=>25,10),(Distributions.Gaussian(5,10),50=>25,10)])
type_matrix_test_expected4_part1 = Array{Union{DIST_FLAGS,Missing}}(undef,50,25)
fill!(type_matrix_test_expected4_part1,Bernoulli)
type_matrix_test_expected4_part2 = Array{Union{DIST_FLAGS,Missing}}(undef,50,25)
fill!(type_matrix_test_expected4_part2,Gaussian)
type_matrix_test_expected4 = hcat(type_matrix_test_expected4_part1,type_matrix_test_expected4_part2)
@test construct_type_matrix(type_matrix_test_input4) == type_matrix_test_expected4

import MatrixCompletion.Utilities.Indexing:IndexTracker
import MatrixCompletion.Utilities.Indexing:construct_index_tracker
# INDEX TRACKER TEST 1
id_tracker_test_input1 = construct_type_matrix(Random.rand(Distributions.Gaussian(0,1),10,10))
id_tracker_test_output1 = construct_index_tracker(input_type_matrix = id_tracker_test_input1)
id_tracker_test_expect1_gaussian = [CartesianIndex(i,j) for i in 1:10 for j in 1:10]
@test sort(id_tracker_test_expect1_gaussian) == sort(id_tracker_test_output1.Gaussian)
@test isempty(id_tracker_test_output1.Gamma)
@test isempty(id_tracker_test_output1.Bernoulli)
@test isempty(id_tracker_test_output1.Poisson)
@test isempty(id_tracker_test_output1.NegativeBinomial)
@test isempty(id_tracker_test_output1.Missing)

# INDEX TRACKER TEST 2
id_tracker_test_input2 = construct_type_matrix(rand([(Distributions.Bernoulli(0.5), 100=>25, 10),
                                                     (Distributions.Gaussian(5,10), 100=>25, 10),
                                                     (Distributions.Poisson(10),    100=>25, 10),
                                                     (Distributions.Gaussian(0,1),  100=>25, 10)]))
id_tracker_test_output2 = construct_index_tracker(input_type_matrix = id_tracker_test_input2)
id_tracker_test_expect2_gaussian = [CartesianIndex(i,j) for i in 1:100 for j in [26:50;76:100]]
id_tracker_test_expect2_bernoulli = [CartesianIndex(i,j) for i in 1:100 for j in 1:25]
id_tracker_test_expect2_poisson = [CartesianIndex(i,j) for i in 1:100 for j in 51:75]
@test sort(id_tracker_test_output2.Gaussian) == sort(id_tracker_test_expect2_gaussian)
@test sort(id_tracker_test_output2.Bernoulli) == sort(id_tracker_test_expect2_bernoulli)
@test sort(id_tracker_test_output2.Poisson) == sort(id_tracker_test_expect2_poisson)
@test isempty(id_tracker_test_output2.Missing)
@test isempty(id_tracker_test_output2.Gamma)
@test isempty(id_tracker_test_output2.NegativeBinomial)


text:
using Plots
include("experiment_helpers.jl")

# Create Roadway
roadway_opts = MergingRoadwayOptions(merging_point=1.4, lane_width=0.28)
roadway = build_roadway(roadway_opts)

# Create players
T = Float64
p = 3
model = BicycleGame(p=p)
n = model.n
m = model.m

# Define the horizon of the problem
N = 15 # N time steps
dt = 0.10 # each step lasts 0.1 second
probsize = ProblemSize(N,model) # Structure holding the relevant sizes of the problem

# Define the objective of each player
# We use a LQR cost
Q = [Diagonal(3*SVector{model.ni[i],T}([1., 1.2, 1., 1.])) for i=1:p] # Quadratic state cost
R = [Diagonal(0.3*ones(SVector{model.mi[i],T})) for i=1:p] # Quadratic control cost
# Desrired state
xf = [SVector{model.ni[1],T}([2.0, 0.05, 0.0, 0]),
      SVector{model.ni[2],T}([2.0, 0.05, 0.0, 0]),
      SVector{model.ni[2],T}([2.0, 0.05, 0.0, 0]),
      ]
# Desired control
uf = [zeros(SVector{model.mi[i],T}) for i=1:p]
# Objectives of the game
game_obj = Algames.GameObjective(Q,R,xf,uf,N,model)
# radius = 0.22*ones(p)
# μ = 1.0*ones(p)
# add_collision_cost!(game_obj, radius, μ)

# Define the constraints that each player must respect
game_con = Algames.GameConstraintValues(probsize)
radius = 0.08
add_collision_avoidance!(game_con, radius)

# Define the initial state of the system
x0 = SVector{model.n,T}([
    0.2,     0.2,    0.2,
   -0.45,   -0.13,   0.13,
    0.45,    0.60,   0.60,
    0.15,    0.0,    0.0,
    ])

# Define the Options of the solver
opts = Options()
opts.ls_iter = 15
opts.outer_iter = 20
opts.inner_iter = 20
opts.ρ_0 = 1e0
opts.reg_0 = 1e-5
opts.α_dual = 1.0
opts.λ_max = 1.0*1e7
opts.ϵ_dyn = 1e-6
opts.ϵ_sta = 1e-6
opts.ϵ_con = 1e-6
opts.ϵ_opt = 1e-6
opts.regularize = true
# Define the game problem
prob = Algames.GameProblem(N,dt,x0,model,opts,game_obj,game_con)

# Solve the problem
newton_solve!(prob)

plot!(prob.model, prob.pdtraj.pr)
plot!(prob.stats)

players = Vector{Player{T}}(undef, p)
players[1] = Player(model, roadway.lane[3])
players[2] = Player(model, roadway.lane[3])
players[3] = Player(model, roadway.lane[3])

# Create Scenario
sce = Scenario(model, roadway, players)

# Initialize visualizers
vis = Visualizer()
open(vis)

# Visualize trajectories
set_scenario!(vis, sce)
set_env!(vis, VehicleState(0.9, -0.2, 0.0, 0.0))
set_camera_birdseye!(vis, height=3.5)

build_waypoint!(vis, sce.player, N, key=0)
set_waypoint_traj!(vis, model, sce, prob.pdtraj.pr, key=0)
set_traj!(vis, model, sce, prob.pdtraj.pr)
set_line_traj!(vis, model, sce.player, prob.pdtraj.pr, vis_opts=LineVisualizationOptions(line_width=8.0))

get_num_active_constraint(prob)
prob_copy = deepcopy(prob)
ascore = ActiveSetCore(probsize)
active_vertical_mask!(ascore, prob_copy.game_con)
vmask = deepcopy(ascore.vmask)

function get_figure(prob::Algames.GameProblem)
	probize = prob.probsize
	ascore = ActiveSetCore(probsize)
	update_nullspace!(ascore, prob, prob.pdtraj)
	ns = length(ascore.null.vec)
	N = prob.probsize.N
	n = prob.probsize.n
	m = prob.probsize.m
	p = prob.probsize.p
	S = prob.probsize.S
	Sh = S + p*(p-1)*(N-1)
	Sv = S + Int(p*(p-1)*(N-1)/2)
	x0 = get_primal_dual(ascore, prob)
	v0 = deepcopy(ascore.null.vec[2])

 	pdtrajs = []

	plt = plot()
	for l = 1:100
		ascore = ActiveSetCore(probsize)
		active_vertical_mask!(ascore, prob.game_con)
		vmask = deepcopy(ascore.vmask)
		na = get_num_active_constraint(prob)

		β = 1.0
		α = 1e-3
		λ = 1e-10*rand(Sv)
		x = get_primal_dual(ascore, prob)

		xp, λ = simple_projection(α, λ, x, v0, β, ascore, prob, vmask)
		set_primal_dual!(ascore, prob, xp)
		push!(pdtrajs, deepcopy(prob.pdtraj))

		plot!(model, prob.pdtraj.pr, plt=plt)
		prob_c = deepcopy(prob)
		set_primal_dual!(ascore, prob_c, x+α*v0)
	end
	display(plt)

	plt = plot()
	set_primal_dual!(ascore, prob, x0)
	plot!(model, prob.pdtraj.pr, plt=plt)
	set_primal_dual!(ascore, prob, x0+100*1e-3*v0)
	plot!(model, prob.pdtraj.pr, plt=plt)
	display(plt)

	return prob, pdtrajs
end

prob_copy = deepcopy(prob)
prob_copy, pdtrajs = get_figure(prob_copy)

for (k,pdtraj) in enumerate(pdtrajs[[1,20,40,60,80,100]])
	vis_opts = LineVisualizationOptions(line_width=2.0, α=1.0)
	set_line_traj!(vis, model, sce.player, pdtraj.pr, vis_opts=vis_opts, key=k)
end

settransform!(vis["env/roadway"], Translation(0.0, -0.05, 0.0))

# for (k,pdtraj) in enumerate(pdtrajs[[(1:1:80); (80:-1:1)]])
for (k,pdtraj) in enumerate(pdtrajs[[(1:1:80); (80:-1:1)]])
	set_traj!(vis, model, sce, pdtraj.pr[end:end])
	set_line_traj!(vis, model, sce.player, pdtraj.pr, vis_opts=LineVisualizationOptions(line_width=8.0))
	sleep(0.3)
	save_image(vis)
end

prob.pdtraj.pr

[(1:5:100);]
[(1:5:100);]
plot([(1:5:100); (100:-5:1)])

using FFMPEG
using MeshCat
filename = "GNE_continuum"
MeshCat.convert_frames_to_video(
# convert_frames_to_video(
    "/home/simon/Downloads/$filename.tar",
    "/home/simon/Documents/$filename.mp4", overwrite=true)

convert_video_to_gif(
    "/home/simon/Documents/$filename.mp4",
    "/home/simon/Documents/$filename.gif", overwrite=true)

unpack_cmd = MeshCat.unpack_cmd
using FFMPEG

function convert_frames_to_video(tar_file_path::AbstractString, output_path::AbstractString="output.mp4"; framerate=60, overwrite=false)
    output_path = abspath(output_path)
    if !isfile(tar_file_path)
        error("Could not find the input file $tar_file_path")
    end
    if isfile(output_path) && !overwrite
        error("The output path $output_path already exists. To overwrite that file, you can pass `overwrite=true` to this function")
    end

    mktempdir() do tmpdir
        run(unpack_cmd(tar_file_path, tmpdir, ".zip", nothing))
        cmd = ["-r", string(framerate), "-i", "*.png", "-vcodec", "libx264", "-preset", "slow", "-crf", "18"]
        if overwrite
            push!(cmd, "-y")
        end
        push!(cmd, output_path)

        cd(tmpdir) do
            FFMPEG.exe(cmd...)
        end
    end
    @info("Saved output as $output_path")
    return output_path
end

function convert_video_to_gif(video_file_path::AbstractString, output_path::AbstractString="output.gif";
    framerate::Int=30, start_time=0., duration=1e3, overwrite=false, width::Int=1080, height::Int=-2, hq_colors::Bool=false)
    output_path = abspath(output_path)

    if !isfile(video_file_path)
        error("Could not find the input file $video_file_path")
    end
    if isfile(output_path) && !overwrite
        error("The output path $output_path already exists. To overwrite that file, you can pass `overwrite=true` to this function")
    end

    mktempdir() do tmpdir
        # run(MeshCat.unpack_cmd(video_file_path, tmpdir, ".mp4", nothing)) # unpack the .tar file
        # cmd = ["-r", string(framerate), "-i", "%07d.png", "-vcodec", "libx264", "-preset", "slow", "-crf", "18"]
        color_map = hq_colors ?
            "[0:v] fps=$framerate, scale=$width:$height,split [a][b];[a] palettegen=stats_mode=single [p];[b][p] paletteuse=new=1" :
            "[0:v] fps=$framerate, scale=$width:$height,split [a][b];[a] palettegen [p];[b][p] paletteuse"
        cmd = ["-ss", string(start_time), "-t", string(duration), "-i", video_file_path, "-filter_complex", color_map]
        if overwrite
            push!(cmd, "-y")
        end
        push!(cmd, output_path)

        cd(tmpdir) do
            FFMPEG.exe(cmd...)
        end
    end
    @info("Saved output as $output_path")
    return output_path
end


text:
using TensorBoardLogger #import the TensorBoardLogger package
using Logging #import Logging package
using TestImages
using Flux, Flux.Data.MNIST
using LinearAlgebra

logger = TBLogger("embeddinglogs", tb_append) #create tensorboard logger

################log embeddings example: MNIST data################
n_samples = 100
features = 100
data = MNIST.images()[1:n_samples]
n = size(data, 1)
data = hcat(data...)
data = reshape(data, (28, 28, n))
imgs = data
data = permutedims(data, (3, 1, 2))
data = reshape(data, n, 28*28)
data = convert(Array{Float64,2}, data)
svd_data = svd(data)
u = svd_data.U
s = Array(Diagonal(svd_data.S))[:, 1:features]
vt = svd_data.Vt[1:features, :]
data = u*s
metadata = MNIST.labels()[1:n]

#using explicit function interface
log_embeddings(logger, "embedding/mnist/explicitinterface", data, metadata = metadata, img_labels = TBImages(imgs, HWN), step = 1)


text:
# After reviewing the available paths, you realize you might have time to visit
# a single small cave twice. Specifically, big caves can be visited any number
# of times, a single small cave can be visited at most twice, and the remaining
# small caves can be visited at most once. However, the caves named start and
# end can only be visited exactly once each: once you leave the start cave, you
# may not return to it, and once you reach the end cave, the path must end
# immediately.
#
# Now, the 36 possible paths through the first example above are:
#
# start,A,b,A,b,A,c,A,end
# start,A,b,A,b,A,end
# start,A,b,A,b,end
# start,A,b,A,c,A,b,A,end
# start,A,b,A,c,A,b,end
# start,A,b,A,c,A,c,A,end
# start,A,b,A,c,A,end
# start,A,b,A,end
# start,A,b,d,b,A,c,A,end
# start,A,b,d,b,A,end
# start,A,b,d,b,end
# start,A,b,end
# start,A,c,A,b,A,b,A,end
# start,A,c,A,b,A,b,end
# start,A,c,A,b,A,c,A,end
# start,A,c,A,b,A,end
# start,A,c,A,b,d,b,A,end
# start,A,c,A,b,d,b,end
# start,A,c,A,b,end
# start,A,c,A,c,A,b,A,end
# start,A,c,A,c,A,b,end
# start,A,c,A,c,A,end
# start,A,c,A,end
# start,A,end
# start,b,A,b,A,c,A,end
# start,b,A,b,A,end
# start,b,A,b,end
# start,b,A,c,A,b,A,end
# start,b,A,c,A,b,end
# start,b,A,c,A,c,A,end
# start,b,A,c,A,end
# start,b,A,end
# start,b,d,b,A,c,A,end
# start,b,d,b,A,end
# start,b,d,b,end
# start,b,end
#
# The slightly larger example above now has 103 paths through it, and the even
# larger example now has 3509 paths through it.
#
# Given these new rules, how many paths through this cave system are there?

# load data into a dictionary
caves = Dict{String, Vector{String}}()
foreach(eachline("input.txt")) do l
  (k, v) = split(l, '-')
  if haskey(caves, k)
    push!(caves[k], v)
  else
    caves[k] = [v]
  end
  if haskey(caves, v)
    push!(caves[v], k)
  else
    caves[v] = [k]
  end
end

# A dictionary of small caves names to an arbitrary index. Used to map small
# caves to a BitVector of visited caves.
small_cave_idxs = Dict([(k,i) for (i,k) in enumerate(filter(c -> all(islowercase, c), keys(caves)))])

# Implements a BFS of caves - returns the number of ways to reach end
function visit_cave(cave::String, visited::BitVector, has_double_visited::Bool)
  # have we reached end?
  if cave == "end"
    return 1
  end

  # if cave is a small cave, mark it visited
  if haskey(small_cave_idxs, cave)
    visited = copy(visited)
    visited[small_cave_idxs[cave]] = true
  end

  # recurse
  foldl(caves[cave]; init=0) do acc, adjacent_cave
    if haskey(small_cave_idxs, adjacent_cave)
      if (has_double_visited && visited[small_cave_idxs[adjacent_cave]]) || adjacent_cave == "start"
        # we've already visited this small cave, and already double visited a
        # small cave, so we can't go this way
        return acc
      else
        # if either condition is false, it's either our first time visiting
        # this cave, or we can use it as our one revisit
        return acc + visit_cave(adjacent_cave, visited, has_double_visited || visited[small_cave_idxs[adjacent_cave]])
      end
    else
      # can always go into a large cave
      return acc + visit_cave(adjacent_cave, visited, has_double_visited)
    end
  end
end

# Do eet!
cnt = visit_cave("start", falses(length(small_cave_idxs)), false)
println("Result: ", cnt)


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

name = "FFMPEG"
version = v"4.1.0"

# Collection of sources required to build FFMPEG
sources = [
    "https://ffmpeg.org/releases/ffmpeg-4.1.tar.bz2" =>
    "b684fb43244a5c4caae652af9022ed5d85ce15210835bce054a33fb26033a1a5",
]

# Bash recipe for building across all platforms
# TODO: Theora and Opus once their releases are available
script = raw"""
cd $WORKSPACE/srcdir
cd ffmpeg-4.1/
sed -i 's/-lflite"/-lflite -lasound"/' configure
apk add coreutils
./configure            \
  --prefix=$prefix     \
  --enable-gpl         \
  --enable-version3    \
  --enable-nonfree     \
  --disable-static     \
  --enable-shared      \
  --disable-debug      \
  --enable-avresample  \
  --enable-libass      \
  --enable-libfdk-aac  \
  --enable-libfreetype \
  --enable-libmp3lame  \
  --enable-libvorbis   \
  --enable-libvpx      \
  --enable-libx264     \
  --enable-encoders    \
  --enable-decoders    \
  --enable-muxers      \
  --enable-demuxers    \
  --enable-parsers     \
  --extra-cflags="-I${prefix}/include" \
  --extra-ldflags="-L${prefix}/lib"
#--enable-libx265
make -j${nproc}
make install

"""
# For libx264rgb see discussion at https://stackoverflow.com/a/40409031/1364192

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libavformat", Symbol("")),
    LibraryProduct(prefix, "libavcodec", Symbol("")),
    LibraryProduct(prefix, "libavutil", Symbol("")),
    LibraryProduct(prefix, "libpostproc", Symbol("")),
    LibraryProduct(prefix, "libswresample", Symbol("")),
    LibraryProduct(prefix, "libavdevice", Symbol("")),
    LibraryProduct(prefix, "libavresample", Symbol("")),
    LibraryProduct(prefix, "libavfilter", Symbol("")),
    LibraryProduct(prefix, "libswscale", Symbol(""))
]

# Dependencies that must be installed before this package can be built
# TODO: Theora and Opus once their releases are available
dependencies = [
    "https://github.com/JuliaIO/LibassBuilder/releases/download/v0.14.0/build_libass.v0.14.0.jl",
    "https://github.com/SimonDanisch/FDKBuilder/releases/download/0.1.6/build_libfdk.v0.1.6.jl",
    "https://github.com/SimonDanisch/NASMBuilder/releases/download/2.13.3/build_nasm.v2.13.3.jl",
    "https://github.com/SimonDanisch/FribidiBuilder/releases/download/0.14.0/build_fribidi.v0.14.0.jl",
    "https://github.com/JuliaGraphics/FreeTypeBuilder/releases/download/v2.9.0-0/build.jl",
    "https://github.com/SimonDanisch/LAMEBuilder/releases/download/3.100.0/build_liblame.v3.100.0.jl",
    "https://github.com/JuliaIO/LibVorbisBuilder/releases/download/v1.3.6/build_libvorbis.v1.3.6.jl",
    "https://github.com/staticfloat/OggBuilder/releases/download/v1.3.3-7/build_Ogg.v1.3.3.jl",
    "https://github.com/jpsamaroo/LibVPXBuilder/releases/download/v5.0.0/build_LibVPX.v5.0.0.jl",
    "https://github.com/jpsamaroo/YasmBuilder/releases/download/v1.3.0-pre/build_YasmBuilder.v1.3.0.jl",
    "https://github.com/ianshmean/x264Builder/releases/download/v2019.5.25-noyasm/build_x264Builder.v2019.5.25-noyasm.jl",
    #"https://github.com/jpsamaroo/x265Builder/releases/download/v2.8-noyasm/build_x265Builder.v2.8.0.jl",
]

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

max_stars_repo_name: ianshmean/FFMPEGBuilder

max_issues_repo_name: ianshmean/FFMPEGBuilder

max_forks_repo_name: ianshmean/FFMPEGBuilder


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

module ImageInpainting

using ImageMorphology
using ImageFiltering
using FFTW: set_num_threads
using CpuId: cpucores

include("utils.jl")
include("pointgradients.jl")
include("inpaint.jl")

export
  # main function
  inpaint,

  # algorithms
  Criminisi

end


text:
function rosenbrock_init!(x0::Array{T,1}) where {T<:Real}
  x0[1:2:end] .= -1.2
  x0[2:2:end] .=  1.0
  return nothing
end

function rosenbrock_fg!(x::Array{T,1}, gx::Array{T,1}) where {T<:Real}
  local c1::T = 1
  local c2::T = 2
  local c10::T = 10
  local c200::T = 200
  x1 = x[1:2:end]
  x2 = x[2:2:end]
  t1 = c1 .- x1
  t2 = c10*(x2 - x1.*x1)
  g2 = c200*(x2 - x1.*x1)
  gx[1:2:end] = -c2*(x1 .* g2 + t1)
  gx[2:2:end] = g2
  return sum(t1.*t1) + sum(t2.*t2)
end

function rosenbrock_test(n::Integer=20, m::Integer=3; single::Bool=false)
  T = (single ? Float32 : Float64)
  x0 = Array{T}(undef, n)
  rosenbrock_init!(x0)
  lbfgs(rosenbrock_fg!, x0, m, verb=VERBOSE)
end

function run_tests(::Type{T}, n::Integer=20) where {T}
    prec = (T <: Float64 ? "double" : "single")
    x0 = Array{T}(undef, n)
    xsol = ones(T,n)
    atol = 1e-3
    rosenbrock_init!(x0)

    @printf("\nTesting NLCG in %s precision\n", prec)
    x1 = nlcg(rosenbrock_fg!, x0, verb=VERBOSE)
    err = maximum(abs.(x1 .- xsol))
    @printf("Maximum absolute error: %.3e\n", err)
    @test err < atol

    @printf("\nTesting VMLMB in %s precision with Oren & Spedicato scaling\n", prec)
    x2 = vmlmb(rosenbrock_fg!, x0, verb=VERBOSE)
    #scaling=OptimPack.SCALING_OREN_SPEDICATO)
    err = maximum(abs.(x2 .- xsol))
    @printf("Maximum absolute error: %.3e\n", err)
    @test err < atol

    @printf("\nTesting VMLMB in %s precision with Oren & Spedicato scaling\n", prec)
    x3 = vmlmb(rosenbrock_fg!, x0, verb=VERBOSE, mem=15)
    #scaling=OptimPack.SCALING_OREN_SPEDICATO)
    err = maximum(abs.(x3 .- xsol))
    @printf("Maximum absolute error: %.3e\n", err)
    @test err < atol

    @printf("\nTesting VMLMB in %s precision with nonnegativity\n", prec)
    x4 = vmlmb(rosenbrock_fg!, x0, verb=VERBOSE, lower=0)
    err = maximum(abs.(x4 .- xsol))
    @printf("Maximum absolute error: %.3e\n", err)
    @test err < atol

    #@printf("\nTesting VMLM in %s precision with Barzilai & Borwein scaling\n", prec)
    #x5 = vmlmb(rosenbrock_fg!, x0, verb=VERBOSE,
    #           scaling=OptimPack.SCALING_BARZILAI_BORWEIN)
    #err = maximum(abs.(x5 .- xsol))
    #@printf("Maximum absolute error: %.3e\n", err)
    #@test err < atol
end

# Run tests in double and single precisions.
run_tests(Float32)
run_tests(Float64)


text:
using Distributions
import Distributions.logpdf
using Base.Cartesian
using Base.Threads
using FillArrays
using CuArrays

include("gpulogpdfs.jl")

export logpdf
export rand

export For
struct For{D,P,X,N}
    params :: Array{P,N}
end

#########################################################
# T <: NTuple{N,J} where {J <: Integer}
#########################################################

For(f, dims::J...) where {J <: Integer} = For(f, dims)

function For(f, dims::T) where {N, J<:Integer, T<:NTuple{N,J}}
    θ = one(CartesianIndex{N})
    d = f(Tuple(θ)...)
    D = typeof(d)
    X = eltype(d)
    p = params(d)
    P = typeof(p)
    ps = Array{P,N}(undef, dims)
    ps[θ] = p
    @inbounds @simd for θ in CartesianIndices(dims)[2:end]
        d = f(Tuple(θ)...)
        ps[θ] = params(d)
    end
    For{D,P,X,N}(ps)
end

@inline function logpdf(d::For{D}, xs::AbstractArray) where D
    s = 0.0
    for (p, x) in zip(d.params, xs)
        s += logpdf(D(p...), x)
    end
    s
end

@inline function logpdf(d  :: For{Normal{X}, Tuple{X,X}, X, N},
                        xs :: CuArray{<:Real, N}) where {X<:Real, N}
    sum(cunormlogpdf.(d.params |> cu, xs))
end

function Base.rand(dist::For{D}) where D
    map(dist.params) do p
        (rand ∘ D)(p...)
    end
end

#########################################################
# T <: NTuple{N,J} where {J <: AbstractUnitRange}
#########################################################

# For(f, θ::J...) where {J <: AbstractUnitRange} = For(f,θ)
#
# function For(f::F, θ::T) where {F, N, J <: AbstractRange, T <: NTuple{N,J}}
#     d = f.(ones(Int, N)...)
#     D = typeof(d)
#     X = eltype(d)
#     For{F, NTuple{N,J}, D, X}(f,θ)
# end
#
#
# @inline function logpdf(d::For{F,T,D,X1},xs::AbstractArray{X2,N}) where {F, N, J <: AbstractRange,  T <: NTuple{N,J}, D, X1, X2}
#     s = 0.0
#     @inbounds @simd for θ in CartesianIndices(d.θ)
#         s += logpdf(d.f(Tuple(θ)...), xs[θ])
#     end
#     s
# end
#
#
# function Base.rand(dist::For{F,T}) where {F,  N, J <: AbstractRange, T <: NTuple{N,J}}
#     map(CartesianIndices(dist.θ)) do I
#         (rand ∘ dist.f)(Tuple(I)...)
#     end
# end

#########################################################
# T <: Base.Generator
#########################################################

# function For(f::F, θ::T) where {F, T <: Base.Generator}
#     d = f(θ.f(θ.iter[1]))
#     D = typeof(d)
#     X = eltype(d)
#     For{F, T, D, X}(f,θ)
# end
#
#
# @inline function logpdf(d :: For{F,T}, x) where {F,T <: Base.Generator}
#     s = 0.0
#     for (θj, xj) in zip(d.θ, x)
#         s += logpdf(d.f(θj), xj)
#     end
#     s
# end
#
# @inline function rand(d :: For{F,T,D,X}) where {F,T <: Base.Generator, D, X}
#     rand.(Base.Generator(d.θ.f, d.θ.iter))
# end

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

# export logpdf2
# @inline function logpdf2(d::For{F,N,X1},xs) where {F,N, X1, X2}
#     results = zeros(eltype(xs), nthreads())
#
#     θ = CartesianIndices(d.θ)
#
#     total = Threads.Atomic{Float64}(0.0)
#     @threads for tid in 1:nthreads()
#         # split work
#         start = 1 + ((tid - 1) * length(xs)) ÷ nthreads()
#         stop = (tid * length(xs)) ÷ nthreads()
#         domain = start:stop
#
#         s = 0.0
#         for j in domain
#             @inbounds θj = θ[j]
#             @inbounds s += logpdf(d.f(Tuple(θj)...), xs[θj])
#         end
#
#         Threads.atomic_add!(total, s)
#     end
#
#     total.value
# end

# @inline function importanceSample(p::For{F1,N,X}, q::For{F2,N,X}) where {F1,F2,N,X}
#     _ℓ = 0.0
#     x = Array{X, N}(undef, length.(q.θ))
#     @inbounds @simd for θ in CartesianIndices(q.θ)
#         I = Tuple(θ)
#         ℓx = importanceSample(p.f(I...), q.f(I...))
#         _ℓ += ℓx.ℓ
#         x[θ] = ℓx.val
#     end
#     Weighted(_ℓ,x)
# end

# using Transducers
# using Transducers: @next, complete
# function Transducers.__foldl__(rf, val, d::For)
#     for θ in d.θ
#         val = @next(rf, val, f(θ))
#     end
#     return complete(rf, val)
# end


text:
function t_interfaces()

N = 10	# variables
M = 25  # samples

vectordata = rand(M)
matrixdata = rand(N,M)
labels = round.(2*rand(M))
multilabels = round.(2*rand(2,M))

A01 = j4pr.datacell(vectordata)
A02 = j4pr.datacell(vectordata, labels)
A03 = j4pr.datacell(vectordata, multilabels)
A04 = j4pr.datacell(matrixdata)
A05 = j4pr.datacell(matrixdata, labels)
A06 = j4pr.datacell(matrixdata, multilabels)
AV = [A01, A02, A03, A04, A05, A06]

# Simple getindex type stability checks
#print("Checking getindex type stability for data cells... ")

# Index for ::Void and empty index
Test.@test nothing[1] isa Void
Test.@test nothing[:,1] isa Void
Test.@test nothing[:,[1,2]] isa Void
Test.@test nothing[:,1:2] isa Void

# Index with single integer (works for all data types e.g. vector and matrix)
for A in AV
	Test.@test A[5] isa typeof(A)
end

# Integer + Colon (works only for matrix data)
for A in AV[4:end]
	Test.@test A[5,:] isa typeof(A)
end

# Colon + Integer (works only for matrix data)
for A in AV[4:end]
	Test.@test A[:,5] isa typeof(A)
end

# Range/Vector (works only for vector data)
for A in AV[1:3]
	Test.@test A[[1,2,3]] isa typeof(A)
	Test.@test A[1:3] isa typeof(A)
end

# Range/Vector + Colon (works only for matrix data)
for A in AV[4:end]
	Test.@test A[[1,2,3],:] isa typeof(A)
	Test.@test A[1:3,:] isa typeof(A)
end

# Colon + Range/Vector (works only for matrix data)
for A in AV[4:end]
	Test.@test A[:,[1,2,3]] isa typeof(A)
	Test.@test A[:,1:3] isa typeof(A)
end
# Range/Vector + Range Vector (works only for matrix data)
for A in AV[4:end]
	Test.@test A[[1,2,3],[1,2,3]] isa typeof(A)
	Test.@test A[1:3,[1,2,3]] isa typeof(A)
	Test.@test A[1:3,[1,2,3]] isa typeof(A)
	Test.@test A[1:3,1:3] isa typeof(A)
end
#println("PASSED")

# Simple getindex type stability checks
#print("Checking setindex! for vector and matrix data... ")

Xvec = rand(3)
Xmat = rand(2,2)

for A in AV[1:3]
	A[1:length(Xvec)] = Xvec
	Test.@test A.x[1:length(Xvec)] == Xvec
end

for A in AV[4:end]
	A[1:size(Xmat,1),1:size(Xmat,2)] = Xmat
	Test.@test A.x[1:size(Xmat,1),1:size(Xmat,2)]  == Xmat
end
#println("PASSED")

# Simple pipe iteration checks
W1 = j4pr.ohenc("binary")
W2 = j4pr.scaler!(Dict(1=>"mean"))
W3 = j4pr.functioncell((x)->x)

PT = [W1;W2]
PP=[W1 W2]
PS = W1+W2

PC = [PT [PS;PP]]

#print("Checking pipe iteration... ")
for P in [PT, PP, PS]
	i = 1
	for W in P
		Test.@test W == P.x[i]
		i+=1
	end
end

for (i,P) in enumerate(PC)
	Test.@test P == PC.x[i]
end
#println("PASSED")

end


text:
@testset "645.set-mismatch.jl" begin
    @test find_error_nums([1, 2, 2, 4]) == (2, 3)
    @test find_error_nums([1, 1]) == (1, 2)
end


text:
"""
Bundle functions regarding (rooted) tree graphs.
"""
module Tree

include("Tree/root.jl")

"""
A rooted tree represanted by its parent vector, i.e. if `u` is the parent node of `v` then `parent[v] == u`.
We always require that root is its own parent (`parent[root] == root`) and `length(parent) == n`.

Often you might want to have a [`ChildrenIndex`](@ref).
For debugging, the [`hierarchy`](@ref) function might be useful.
"""
struct RootedTree
    root::Int
    parent::Vector{Int}
end

RootedTree(parent::Vector{Int}) =
    RootedTree(find_root(parent), parent)

Base.length(t::RootedTree) =  length(t.parent)

struct WeightedTree{W}
    tree::RootedTree
    lam::W
end

include("Tree/degrees.jl")
include("Tree/children.jl")
include("Tree/dfs.jl")
include("Tree/lca.jl")
include("Tree/inorder.jl")
include("Tree/hierarchy.jl")

end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
import Pkg
Pkg.activate("..")

using CrimsonSkyline
using Distributions: Normal, LogNormal, MvNormal, truncated
using Logging
using StatsBase: mean, std
using Random: seed!

seed!(2021)

const Maybe{T} = Union{T, Nothing}

function eye(d :: Int)
    mat = zeros(Float64, d, d)
    for i in 1:d
        mat[i, i] = 1.0
    end
    mat
end

function linear_model!(t :: Trace, out :: Vector{Float64}, X :: Matrix{Float64}, y :: Maybe{Vector{Float64}})
    size_X = size(X)
    D = size_X[1]
    N = size_X[2]
    β = sample(t, :β, MvNormal(zeros(D), eye(D)))
    scale = sample(t, :scale, LogNormal(0.0, 1.0))
    loc = X' * β
    if y === nothing
        for n in 1:N
            out[n] = sample(t, (:y, n), Normal(loc[n], scale))
        end
    else
        for n in 1:N
            out[n] = observe(t, (:y, n), Normal(loc[n], scale), y[n])
        end
    end
end

beta_proposal(old_t :: Trace, new_t :: Trace, _, _, _) = propose(new_t, :β, MvNormal(old_t[:β].value, 0.5 .* eye(length(old_t[:β].value))))
scale_proposal(old_t :: Trace, new_t :: Trace, _, _, _) = propose(new_t, :scale, truncated(Normal(old_t[:scale].value, 0.25), 0.0, Inf))

function main()
    @info "Bayesian linear regression"
    D = 10
    N = 200
    X = randn(D, N)
    out = Vector{Float64}(undef, N)
    @info "Using model to generate dataset"
    t = trace()
    linear_model!(t, out, X, nothing)
    y = [t[(:y, n)].value for n in 1:N]
    β_true = t[:β].value
    @info "True regression coefficients: $β_true"

    @info "Conducting inference using user-defined Metropolis kernels"
    @time results = mh(
        linear_model!,
        [beta_proposal, scale_proposal];
        params = (out, X, y),
        burn = 1000,
        thin = 50,
        num_iterations = 5000,
        inverse_verbosity = 500
    )
    β_mean = mean(results, :β)
    β_std = std(results, :β)
    @info "Posterior mean regression coefficients = $β_mean"
    @info "Posterior std regression coefficients = $β_std"

    # serve posterior predictive using empirical posterior
    new_X = randn(D, 10)
    post_pred_linear_model = update(linear_model!, results)
    @info "Serving predictions from posterior predictive"
    @time served_trace, _ = post_pred_linear_model(trace(), out, new_X, nothing)
end

main()


text:
##############################
# memory allocation
##############################

# # first start julia with:
# julia --track-allocation=user malloc.jl <args>
# where args is e.g. true
local_cluster = parse(Bool, ARGS[1])
@info "tracking memory allocation with local_cluster=$local_cluster"

using BenchmarkTools, Profile, Traceur, InteractiveUtils, Lasso

include("helpers.jl")

kwargs = Dict(:verbose=>false)

n = 300
p = 3
covars, counts = bdata(n, p, 10)
@assert sum(all(iszero.(counts), dims=1)) == 0 "some columns are always zero"

@info "First run is to ensure that everything is compiled (because compilation allocates memory)."
hdmr(covars, counts; local_cluster=local_cluster, kwargs...)

@info "clear stuff contaminated by compilation"
Profile.clear_malloc_data()

@info "Run your commands again"
hdmr(covars, counts; local_cluster=local_cluster, kwargs...)

@info "done tracking memory allocation"

@info "Quit julia!"
exit()

# Finally, navigate to the directory holding your source code.
# Start julia (without command-line flags), and analyze the results using

# using Coverage
# analyze_malloc(".")  # could be "." for the current directory, or "src", etc.

# This will return a vector of MallocInfo objects, specifying the number of bytes allocated
# , the file name, and the line number. These are sorted in increasing order of allocation size.


text:
using Knet, KnetModules
import KnetModules.forward

function vgg16(;trained=true, o...)
    model = VGG([64, 128, 256, 512, 512],
                [2,  2,   3,   3,   3];
                o...)
    # TODO: load pre-trained weights
    return model
end

function vgg19(;trained=true, o...)
    model = VGG([64, 128, 256, 512, 512],
                [2,  2,   4,   4,   4];
                o...)
    return model
end

type VGG <: KnetModule
    conv::Sequential
    fc::Union{Sequential, Void}
    top::Union{Linear, Void}
end

function VGG(channels::Array{Int,1}, repeats::Array{Int,1};
             fc=true, top=true, pdrop=0.5)
    @assert (fc || ~top) "Top cannot exist without fc layers"
    # Conv
    conv = Sequential()
    i = 3
    for (c, r) in zip(channels, repeats)
        block = Sequential()
        for _ = 1:r
            add!(block, Conv(3, 3, i, c; padding=1))
            add!(block, ReLU())
            i = c
        end
        add!(block, MaxPool())
        add!(conv, block)
    end
    # FC
    fcin = i * div(224, 2^length(channels))^2
    fc = fc ? Sequential(
        Linear(4096, fcin), ReLU(), Dropout(pdrop),
        Linear(4096, 4096), ReLU(), Dropout(pdrop)
    ) : nothing
    # output
    top = top ? Linear(1000, 4096) : nothing
    return VGG(conv, fc, top)
end

function forward(ctx, vgg::VGG, x)
    o = @mc vgg.conv(x)
    if vgg.fc !== nothing
        o = @mc vgg.fc(o)
        if vgg.top !== nothing
            o = @mc vgg.top(o)
        end
    end
    return o
end


text:
using Documenter, Literate
ENV["PLOTS_DEFAULT_BACKEND"] = "GR"
ENV["GKSwstype"]="nul"
using Plots
using Plots.PlotMeasures: mm

default(
    fontfamily = "Computer Modern", linewidth=1.5,
    titlefontsize=(16+8), guidefontsize=(11+5),
    tickfontsize=(8+4), legendfontsize=(8+4),
    left_margin=5mm, right_margin=5mm)

src = joinpath(@__DIR__, "src")
lit = joinpath(@__DIR__, "lit")

config = Dict(
    "credit" => false,  # credit is configured to render in Documenter instead
    "repo_root_url"=> "https://github.com/simonsobs/PSpipe/tree/planckcov/project/Planck_cov",
)

nonexecute_config = copy(config)
nonexecute_config["codefence"] = "```julia" => "```"
execution_exclusion = []

for (root, _, files) ∈ walkdir(lit), file ∈ files
    splitext(file)[2] == ".jl" || continue
    ipath = joinpath(root, file)
    opath = splitdir(replace(ipath, lit=>src))[1]
    if file ∉ execution_exclusion
        Literate.markdown(ipath, opath; config=config)
    else
        @warn("not executing unless forced for ", file)
        Literate.markdown(ipath, opath; config=nonexecute_config)
    end
end

makedocs(
    sitename = "PSPipe Planck",
    modules = Module[],
    clean=false,
    format=Documenter.HTML(
        prettyurls=get(ENV, "CI", nothing) == "true",
        assets=["assets/so.css"],
    ),
    pages = [
        "Introduction" => "index.md",
        "setup.jl" => "setup.md",
        "util.jl" => "util.md",
        "rawspectra.jl" => "rawspectra.md",
        "fitnoisemodel.jl" => "fitnoisemodel.md",
        "slurmgen.jl" => "slurmgen.md",
        "signalsim.jl" => "signalsim.md",

        "Pipeline" => "pipeline.md"
        ]
    )

deploydocs(;
    repo="github.com/xzackli/PSPipePlanckRender.jl",
    devbranch="main"
)

max_stars_repo_name: xzackli/PSPipePlanckRender

max_issues_repo_name: xzackli/PSPipePlanckRender

max_forks_repo_name: xzackli/PSPipePlanckRender


text:
using RecipesBase
using Statistics

@recipe function f(spots::AbstractVector{Spot})
    nday = [s.nday for s in spots]
    lat = [s.lat for s in spots]
    Φmax = [s.Φmax for s in spots]

    seriestype --> :scatter
    xlabel --> "Time (days)"
    ylabel --> "Latitude (deg)"
    markersize --> 3Φmax ./ median(Φmax)(spots.Φmax)
    alpha --> 0.5
    color --> "#996699"
    markerstrokealpha --> 0

    nday, rad2deg.(lat)
end

@recipe function f(r::SpotDynamics)
    spots = r.spots
    nday = [s.nday for s in spots]
    lat = [s.lat for s in spots]
    Φmax = [s.Φmax for s in spots]

    seriestype --> :scatter
    xlabel --> "Time (days)"
    ylabel --> "Latitude (deg)"
    markersize --> 3Φmax ./ median(Φmax)
    alpha --> 0.5
    color --> "#996699"
    markerstrokealpha --> 0

    nday, rad2deg.(lat)
end

max_stars_repo_name: zclaytor/butterpy

max_issues_repo_name: zclaytor/butterpy

max_forks_repo_name: zclaytor/butterpy


text:
#=
def sample_actual_demand(price):
    avg_demand = 65 + (-0.8) * price
    theta = 0.1/4
    k = avg_demand / theta
    return np.random.gamma(k*theta, k*theta**2, 1)[0]
=#

function sample_actual_demand(price)
    avg_demand = 65 + (-0.8) * price
    theta = 0.1/4
    k = avg_demand / theta
    return rand(Gamma(k*theta, k*theta^2), 1)
end

#=
def emperical_mean(sampler, n):
    mean = 0
    for i in range(1, n):
        mean = mean + sampler()
    return mean/n
=#

function emperical_mean(sampler, n)
    mean = 0
    for i in 1:n-1
        mean += sampler()[1]
    end
    return mean / n
end

#=
def emperical_demand_curve(min_price, max_price, n):
    prices = np.linspace(min_price, max_price, n)
    sampling = 5000
    demands = map(lambda p: emperical_mean(functools.partial(sample_actual_demand, p), sampling), prices)
    return np.dstack((prices, list(demands)))[0]
=#

function emperical_demand_curve(min_price, max_price, n)
    prices = range(min_price, max_price, length=n)
    sampling = 5000
    demands = []
    for p in prices
        sampler() = sample_actual_demand(p)
        demands = append!(demands, emperical_mean(sampler, sampling))
    end
    return DataFrame(
        :prices => prices,
        :demands => demands
        )
end


text:
typealias HTMLTable Matrix{HTMLCell}

function renderformatspecific(io::IO, ::Type{HTMLCell}, cells::HTMLTable)
    n, m = size(cells)
    println(io, "
")
    for i=1:n
        println(io, "\t
")
        for j=1:m
            showcell(io, cells[i,j])
        end
        println(io, "\t
")
    end
    println(io,"
")
end


text:
"""
Package de fonctions de pénalité

liste des fonctions :
function Quadratic(mp::NLPModels.AbstractNLPModel,G::Function,H::Function,nb_comp::Int64,x::Vector,yg::Vector,yh::Vector,rho::Vector,usg::Vector,ush::Vector,uxl::Vector,uxu::Vector,ucl::Vector,ucu::Vector)
Lagrangian(mp::NLPModels.AbstractNLPModel,G::Function,H::Function,nb_comp::Int64,x::Vector,yg::Vector,yh::Vector,rho::Vector,usg::Vector,ush::Vector,uxl::Vector,uxu::Vector,ucl::Vector,ucu::Vector)

fonctions additionnelles :
RhoDetail(rho::Vector,mp::NLPModels.AbstractNLPModel,nb_comp::Int64)
"""
module Penalty

using NLPModels
using ForwardDiff

"""
Fonction de pénalité quadratique
"""
function Quadratic(mp::NLPModels.AbstractNLPModel,err::Vector,nb_comp::Int64,x::Vector,yg::Vector,
                   yh::Vector,rho::Vector,usg::Vector,ush::Vector,uxl::Vector,uxu::Vector,ucl::Vector,ucu::Vector)
 n=length(mp.meta.x0)

 rho_eqg,rho_eqh,rho_ineq_lvar,rho_ineq_uvar,rho_ineq_lcons,rho_ineq_ucons=RhoDetail(rho,mp,nb_comp)
 err_eq_g,err_eq_h,err_in_lv,err_in_uv,err_in_lc,err_in_uc=ErrDetail(err,n,nb_comp,mp.meta.ncon)

# Pen_eq2=rhof*norm((alas.mod.G(x)-yg))^2+rhof*norm((alas.mod.H(x)-yh))^2
if nb_comp>0
 Pen_eq=dot(rho_eqg.*err_eq_g,err_eq_g)+dot(rho_eqh.*err_eq_h,err_eq_h)
else
 Pen_eq=0.0
end

 Pen_in_lv=dot(rho_ineq_lvar.*err_in_lv,err_in_lv)
 Pen_in_uv=dot(rho_ineq_uvar.*err_in_uv,err_in_uv)

 if mp.meta.ncon!=0
  Pen_in_lc=dot(rho_ineq_lcons.*err_in_lc,err_in_lc)
  Pen_in_uc=dot(rho_ineq_ucons.*err_in_uc,err_in_uc)
 else
  Pen_in_uc=0.0
  Pen_in_lc=0.0
 end

 return 0.5*(Pen_eq+Pen_in_lv+Pen_in_uv+Pen_in_lc+Pen_in_uc)
end

function Quadratic(mp::NLPModels.AbstractNLPModel,err::Vector,nb_comp::Int64,x::Vector,yg::Vector,
                   yh::Vector,rho::Vector,usg::Vector,ush::Vector,uxl::Vector,uxu::Vector,ucl::Vector,ucu::Vector,f::Float64,rho_update::Float64)
 return rho_update*f
end

function Quadratic(mp::NLPModels.AbstractNLPModel,G::NLPModels.AbstractNLPModel,
                   H::NLPModels.AbstractNLPModel,err::Vector,nb_comp::Int64,
                   x::Vector,yg::Vector,yh::Vector,
                   rho::Vector,
                   usg::Vector,ush::Vector,
                   uxl::Vector,uxu::Vector,
                   ucl::Vector,ucu::Vector,
                   g::Vector)

 n=length(mp.meta.x0)

 rho_eqg,rho_eqh,rho_ineq_lvar,rho_ineq_uvar,rho_ineq_lcons,rho_ineq_ucons=RhoDetail(rho,mp,nb_comp)
 err_eq_g,err_eq_h,err_in_lv,err_in_uv,err_in_lc,err_in_uc=ErrDetail(err,n,nb_comp,mp.meta.ncon)

  G_in_lv=-rho_ineq_lvar.*err_in_lv
  G_in_uv=rho_ineq_uvar.*err_in_uv

  if mp.meta.ncon!=0
   Jtv(z)=NLPModels.jtprod(mp,x,z)
   G_in_lc=-Jtv(rho_ineq_lcons.*err_in_lc)
   G_in_uc=Jtv(rho_ineq_ucons.*err_in_uc)
  else
   G_in_uc=zeros(n)
   G_in_lc=zeros(n)
  end

  if nb_comp>0
   G_eq=NLPModels.jtprod(G,x,rho_eqg.*err_eq_g)+NLPModels.jtprod(H,x,rho_eqh.*err_eq_h)
   return vec([(G_eq+G_in_lv+G_in_uv+G_in_lc+G_in_uc)' (-rho_eqg.*err_eq_g)' (-rho_eqh.*err_eq_h)'])
  else
   return G_in_lv+G_in_uv+G_in_lc+G_in_uc
  end
end

function Quadratic(mp::NLPModels.AbstractNLPModel,G::NLPModels.AbstractNLPModel,
                   H::NLPModels.AbstractNLPModel,err::Vector,nb_comp::Int64,
                   x::Vector,yg::Vector,yh::Vector,
                   rho::Vector,
                   usg::Vector,ush::Vector,
                   uxl::Vector,uxu::Vector,
                   ucl::Vector,ucu::Vector,
                   g::Vector,rho_update::Float64)
 return rho_update*g
end

function Quadratic(mp::NLPModels.AbstractNLPModel,G::NLPModels.AbstractNLPModel,
                   H::NLPModels.AbstractNLPModel,err::Vector,nb_comp::Int64,
                   x::Vector,yg::Vector,yh::Vector,
                   rho::Vector,
                   usg::Vector,ush::Vector,
                   uxl::Vector,uxu::Vector,
                   ucl::Vector,ucu::Vector,
                   Hess::Array{Float64,2})

 n=length(mp.meta.x0)

 rho_eqg,rho_eqh,rho_ineq_lvar,rho_ineq_uvar,rho_ineq_lcons,rho_ineq_ucons=RhoDetail(rho,mp,nb_comp)
 err_eq_g,err_eq_h,err_in_lv,err_in_uv,err_in_lc,err_in_uc=ErrDetail(err,n,nb_comp,mp.meta.ncon)

  lv=zeros(n);ilv=find(x->x>0,err_in_lv);lv[ilv]=ones(length(ilv))
  uv=zeros(n);iuv=find(x->x>0,err_in_uv);uv[iuv]=ones(length(iuv))

  if nb_comp>0 && mp.meta.ncon!=0
   ilc=find(x->x>0,err_in_lc);rho_ineq_lcons[ilc]=zeros(length(ilc))
   iuc=find(x->x>0,err_in_uc);rho_ineq_ucons[iuc]=zeros(length(iuc))
   Hlc=hess(mp, x, obj_weight=1.0, y=rho_ineq_lcons.*err_in_lc)+(diagm(rho_ineq_lcons)*NLPModels.jac(mp,x))'*NLPModels.jac(mp,x)
   Huc=hess(mp, x, obj_weight=1.0, y=rho_ineq_ucons.*err_in_uc)+(diagm(rho_ineq_ucons)*NLPModels.jac(mp,x))'*NLPModels.jac(mp,x)

   HCG=hess(G, x, obj_weight=0.0, y=err_eq_g)+(diagm(rho_eqg)*NLPModels.jac(G,x))'*NLPModels.jac(G,x)
   HCH=hess(H, x, obj_weight=0.0, y=err_eq_h)+(diagm(rho_eqh)*NLPModels.jac(H,x))'*NLPModels.jac(H,x)

   return tril([diagm(rho_ineq_lvar.*lv)+diagm(rho_ineq_uvar.*uv)+HCG+HCH+Hlc+Huc zeros(n,2*nb_comp);zeros(2*nb_comp,n) diagm([rho_eqg;rho_eqh])])
  elseif nb_comp>0
   HCG=hess(G, x, obj_weight=0.0, y=err_eq_g)+(diagm(rho_eqg)*NLPModels.jac(G,x))'*NLPModels.jac(G,x)
   HCH=hess(H, x, obj_weight=0.0, y=err_eq_h)+(diagm(rho_eqh)*NLPModels.jac(H,x))'*NLPModels.jac(H,x)

   return tril([diagm(rho_ineq_lvar.*lv)+diagm(rho_ineq_uvar.*uv)+HCG+HCH zeros(n,2*nb_comp);zeros(2*nb_comp,n) diagm([rho_eqg;rho_eqh])])
  elseif mp.meta.ncon!=0
   ilc=find(x->x>0,err_in_lc);rho_ineq_lcons[ilc]=zeros(length(ilc))
   iuc=find(x->x>0,err_in_uc);rho_ineq_ucons[iuc]=zeros(length(iuc))

   Hlc=hess(mp, x, obj_weight=1.0, y=rho_ineq_lcons.*err_in_lc)+(diagm(rho_ineq_lcons)*NLPModels.jac(mp,x))'*NLPModels.jac(mp,x)
   Huc=hess(mp, x, obj_weight=1.0, y=rho_ineq_ucons.*err_in_uc)+(diagm(rho_ineq_ucons)*NLPModels.jac(mp,x))'*NLPModels.jac(mp,x)
   return tril(diagm(rho_ineq_lvar.*lv)+diagm(rho_ineq_uvar.*uv)+Hlc+Huc)
  else
   return tril(diagm(rho_ineq_lvar.*lv)+diagm(rho_ineq_uvar.*uv))
  end
end

function Quadratic(mp::NLPModels.AbstractNLPModel,G::NLPModels.AbstractNLPModel,
                   H::NLPModels.AbstractNLPModel,err::Vector,nb_comp::Int64,
                   x::Vector,yg::Vector,yh::Vector,
                   rho::Vector,
                   usg::Vector,ush::Vector,
                   uxl::Vector,uxu::Vector,
                   ucl::Vector,ucu::Vector,
                   Hess::Array{Float64,2},rho_update::Float64)
 return rho_update*Hess
end

"""
Fonction de pénalité lagrangienne
"""
function Lagrangian(mp::NLPModels.AbstractNLPModel,err::Vector,nb_comp::Int64,x::Vector,yg::Vector,
                   yh::Vector,rho::Vector,usg::Vector,ush::Vector,uxl::Vector,uxu::Vector,ucl::Vector,ucu::Vector)
 n=length(mp.meta.x0)
 #mod::MPCCmod.MPCC,r::Float64,s::Float64,t::Float64
 rho_eqg,rho_eqh,rho_ineq_lvar,rho_ineq_uvar,rho_ineq_lcons,rho_ineq_ucons=RhoDetail(rho,mp,nb_comp)
 err_eq_g,err_eq_h,err_in_lv,err_in_uv,err_in_lc,err_in_uc=ErrDetail(err,n,nb_comp,mp.meta.ncon)

 Lagrangian=dot(err_eq_g,usg)+dot(err_eq_h,ush)
# Pen_eq2=rhof*norm((G(x)-yg))^2+rhof*norm((H(x)-yh))^2
 Pen_eq=dot(rho_eqg.*err_eq_g,err_eq_g)+dot(rho_eqh.*err_eq_h,err_eq_h)
# Pen_in_lv=norm(sqrt(rho_ineq_lvar).*max(mp.meta.lvar-x+uxl./rho_ineq_lvar,0.0))^2
 Pen_in_lv=dot(rho_ineq_lvar.*err_in_lv,err_in_lv)
# Pen_in_uv=norm(sqrt(rho_ineq_uvar).*max(x-mp.meta.uvar+uxu./rho_ineq_uvar,0.0))^2
 Pen_in_uv=dot(rho_ineq_uvar.*err_in_uv,err_in_uv)
# Pen_in_lc=norm(sqrt(rho_ineq_lcons).*max(mp.meta.lcon-mp.c(x)+ucl./rho_ineq_lcons,0.0))^2
 Pen_in_lc=dot(rho_ineq_lcons.*err_in_lc,err_in_lc)
 Pen_in_uc=dot(rho_ineq_ucons.*err_in_uc,err_in_uc)

 return Lagrangian+0.5*(Pen_eq+Pen_in_lv+Pen_in_uv+Pen_in_lc+Pen_in_uc)
end

function Lagrangian(mp::NLPModels.AbstractNLPModel,G::NLPModels.AbstractNLPModel,
                   H::NLPModels.AbstractNLPModel,err::Vector,nb_comp::Int64,
                   x::Vector,yg::Vector,yh::Vector,
                   rho::Vector,
                   usg::Vector,ush::Vector,
                   uxl::Vector,uxu::Vector,
                   ucl::Vector,ucu::Vector,
                   g::Vector)
 n=length(mp.meta.x0)

 rho_eqg,rho_eqh,rho_ineq_lvar,rho_ineq_uvar,rho_ineq_lcons,rho_ineq_ucons=RhoDetail(rho,mp,nb_comp)
 err_eq_g,err_eq_h,err_in_lv,err_in_uv,err_in_lc,err_in_uc=ErrDetail(err,n,nb_comp,mp.meta.ncon)

  G_in_lv=-rho_ineq_lvar.*err_in_lv
  G_in_uv=rho_ineq_uvar.*err_in_uv

  if mp.meta.ncon!=0
   Jtv(z)=NLPModels.jtprod(mp,x,z)
   G_in_lc=-Jtv(rho_ineq_lcons.*err_in_lc)
   G_in_uc=Jtv(rho_ineq_ucons.*err_in_uc)
  else
   G_in_uc=zeros(n)
   G_in_lc=zeros(n)
  end

  if nb_comp>0
   G_eq=NLPModels.jtprod(G,x,rho_eqg.*err_eq_g+usg)+NLPModels.jtprod(H,x,rho_eqh.*err_eq_h+ush)
   return vec([(G_eq+G_in_lv+G_in_uv+G_in_lc+G_in_uc)' (-rho_eqg.*err_eq_g)' (-rho_eqh.*err_eq_h)'])
  else
   return G_in_lv+G_in_uv+G_in_lc+G_in_uc
  end
end

function Lagrangian(mp::NLPModels.AbstractNLPModel,G::NLPModels.AbstractNLPModel,
                   H::NLPModels.AbstractNLPModel,err::Vector,nb_comp::Int64,
                   x::Vector,yg::Vector,yh::Vector,
                   rho::Vector,
                   usg::Vector,ush::Vector,
                   uxl::Vector,uxu::Vector,
                   ucl::Vector,ucu::Vector,
                   Hess::Array{Float64,2})

 n=length(mp.meta.x0)

 rho_eqg,rho_eqh,rho_ineq_lvar,rho_ineq_uvar,rho_ineq_lcons,rho_ineq_ucons=RhoDetail(rho,mp,nb_comp)
 err_eq_g,err_eq_h,err_in_lv,err_in_uv,err_in_lc,err_in_uc=ErrDetail(err,n,nb_comp,mp.meta.ncon)

  lv=zeros(n);ilv=find(x->x>0,err_in_lv);lv[ilv]=ones(length(ilv))
  uv=zeros(n);iuv=find(x->x>0,err_in_uv);uv[iuv]=ones(length(iuv))

  if nb_comp>0 && mp.meta.ncon!=0
   ilc=find(x->x>0,err_in_lc);rho_ineq_lcons[ilc]=zeros(length(ilc))
   iuc=find(x->x>0,err_in_uc);rho_ineq_ucons[iuc]=zeros(length(iuc))
   Hlc=hess(mp, x, obj_weight=1.0, y=rho_ineq_lcons.*err_in_lc)+(diagm(rho_ineq_lcons)*NLPModels.jac(mp,x))'*NLPModels.jac(mp,x)
   Huc=hess(mp, x, obj_weight=1.0, y=rho_ineq_ucons.*err_in_uc)+(diagm(rho_ineq_ucons)*NLPModels.jac(mp,x))'*NLPModels.jac(mp,x)

   HCG=hess(G, x, obj_weight=1.0, y=err_eq_g+usg)+(diagm(rho_eqg)*NLPModels.jac(G,x))'*NLPModels.jac(G,x)
   HCH=hess(H, x, obj_weight=1.0, y=err_eq_h+ush)+(diagm(rho_eqh)*NLPModels.jac(H,x))'*NLPModels.jac(H,x)

   return tril([diagm(rho_ineq_lvar.*lv)+diagm(rho_ineq_uvar.*uv)+HCG+HCH+Hlc+Huc zeros(n,2*nb_comp);zeros(2*nb_comp,n) diagm([rho_eqg;rho_eqh])])
  elseif nb_comp>0
   HCG=hess(G, x, obj_weight=1.0, y=err_eq_g+usg)+(diagm(rho_eqg)*NLPModels.jac(G,x))'*NLPModels.jac(G,x)
   HCH=hess(H, x, obj_weight=1.0, y=err_eq_h+ush)+(diagm(rho_eqh)*NLPModels.jac(H,x))'*NLPModels.jac(H,x)
   return tril([diagm(rho_ineq_lvar.*lv)+diagm(rho_ineq_uvar.*uv)+HCG+HCH zeros(n,n);zeros(2*nb_comp,n) diagm([rho_eqg;rho_eqh])])
  elseif mp.meta.ncon!=0
   ilc=find(x->x>0,err_in_lc);rho_ineq_lcons[ilc]=zeros(length(ilc))
   iuc=find(x->x>0,err_in_uc);rho_ineq_ucons[iuc]=zeros(length(iuc))

   Hlc=hess(mp, x, obj_weight=1.0, y=rho_ineq_lcons.*err_in_lc)+(diagm(rho_ineq_lcons)*NLPModels.jac(mp,x))'*NLPModels.jac(mp,x)
   Huc=hess(mp, x, obj_weight=1.0, y=rho_ineq_ucons.*err_in_uc)+(diagm(rho_ineq_ucons)*NLPModels.jac(mp,x))'*NLPModels.jac(mp,x)
   return tril(diagm(rho_ineq_lvar.*lv)+diagm(rho_ineq_uvar.*uv)+Hlc+Huc)
  else
   return tril(diagm(rho_ineq_lvar.*lv)+diagm(rho_ineq_uvar.*uv))
  end
end

"""
Renvoie : rho_eq, rho_ineq_var, rho_ineq_cons
"""
function RhoDetail(rho::Vector,mp::NLPModels.AbstractNLPModel,nb_comp::Int64)
 nb_ineq_lvar=length(mp.meta.lvar)
 nb_ineq_uvar=length(mp.meta.uvar)
 nb_ineq_lcons=length(mp.meta.lcon)
 nb_ineq_ucons=length(mp.meta.ucon)

 return rho[1:nb_comp],rho[nb_comp+1:2*nb_comp],rho[2*nb_comp+1:2*nb_comp+nb_ineq_lvar],
        rho[2*nb_comp+nb_ineq_lvar+1:2*nb_comp+nb_ineq_lvar+nb_ineq_uvar],
        rho[2*nb_comp+nb_ineq_lvar+nb_ineq_uvar+1:2*nb_comp+nb_ineq_lvar+nb_ineq_uvar+nb_ineq_lcons],
        rho[2*nb_comp+nb_ineq_lvar+nb_ineq_uvar+nb_ineq_lcons+1:2*nb_comp+nb_ineq_lvar+nb_ineq_uvar+nb_ineq_lcons+nb_ineq_ucons]
end

function ErrDetail(err::Vector,n::Int64,nb_comp::Int64,ncon::Int64)

 if nb_comp>0
  err_eq_g=err[1:nb_comp]
  err_eq_h=err[nb_comp+1:2*nb_comp]
 else
  err_eq_g=Array{Float64}(0)
  err_eq_h=Array{Float64}(0)
 end

 err_in_lv=err[2*nb_comp+1:2*nb_comp+n]
 err_in_uv=err[2*nb_comp+n+1:2*nb_comp+2*n]

 if ncon>0
  err_in_lc=err[2*nb_comp+2*n+1:2*nb_comp+2*n+ncon]
  err_in_uc=err[2*nb_comp+2*n+ncon+1:2*nb_comp+2*n+2*ncon]
 else
  err_in_lc=Array{Float64}(0)
  err_in_uc=Array{Float64}(0)
 end

 return err_eq_g,err_eq_h,err_in_lv,err_in_uv,err_in_lc,err_in_uc
end

#end of module
end


text:
using PyCall
skl = pyimport("sklearn.manifold")
function skl_smacof(Δ; anchors = nothing, eps = ϵ / 2)
    Yhat = skl.smacof(Δ, eps = eps)[1]
    return align(Yhat, anchors)
end

# @testset "misc types" begin
#     X = 3 * rand(Float32, 500, 2)
#     Δ = dist(X)
#     # params = SmacofParams(Δ, anchors = X)
#     @time Y = smacof(Δ, anchors = X)
#     @show mse(Y, X)
#     @test mse(Y, X) < ϵ

#     X = 3 * rand(Float32, 500, 2)
#     Δ = dist(X)
#     # params = SmacofParams(Δ, anchors = X)
#     @time Y = smacof(Δ, anchors = X)
#     @test mse(Y, X) < ϵ
# end

# @testset "regular float 64" begin
#     for i in 1:2
#         X = 3 * rand(500, 2)
#         Δ = dist(X)

#         @time Y = smacof(Δ, anchors = X)
#         @show mse(Y, X)

#         @time Yhat = skl_smacof(Δ, anchors = X)
#         @show mse(Yhat, X)
#     end
# end

@testset "embedded square" begin
    for i in 1:2
        X = rand(500, 3)
        X[:, 3] .= 0.0
        U = X[:, 1:2]
        X = rotate(X)
        Δ = dist(X)
        @time Y = smacof(Δ, anchors = U, itmax = 1000)
        @test stress(Δ, dist(Y)) < ϵ

    end
end


text:
"Hermite expansion coefficients for Gaussian product"
function E(a::Float64, b::Float64, Rx::Float64, i::Int, j::Int, t::Int=0)
    if i < 0 || j < 0 || t < 0 || t > (i + j)
        # out of bounds for i, j, or t
        return 0.0
    end

    p = a + b

    if i == j == t == 0 # base case
        return exp(-a*b/p*Rx*Rx)
    end

    if j == 0 # decrement index i
        return (0.5/p)*E(a,b,Rx,i-1,0,t-1) -
               (b/p*Rx)*E(a,b,Rx,i-1,0,t)  +
               (t+1)*E(a,b,Rx,i-1,0,t+1)
    else # decrement index j
        return (0.5/p)*E(a,b,Rx,i,j-1,t-1) +
               (a/p*Rx)*E(a,b,Rx,i,j-1,t)  +
               (t+1)*E(a,b,Rx,i,j-1,t+1)
    end
end

function ∇E(a::Float64, b::Float64, Rx::Float64, i::Int, j::Int, t::Int=0)
    if Rx == 0.0
        return 0.0
    end

    return 2*a*E(a,b,Rx,i+1,j,t) - i*E(a,b,Rx,i-1,j,t)
end

function ∇∇E(a::Float64, b::Float64, Rx::Float64, i::Int, j::Int, t::Int=0)
    return -i*j*E(a,b,Rx,i-1,j-1,t)   +
            2*i*b*E(a,b,Rx,i-1,j+1,t) +
            2*a*j*E(a,b,Rx,i+1,j-1,t) -
            4*a*b*E(a,b,Rx,i+1,j+1,t)
end

function int1e_ovlp(
        a::Float64,
        LA::Vec3{Int},
        b::Float64,
        LB::Vec3{Int},
        RAB::Vec3{Float64},
    )
    p = a + b
    Sx = E(a, b, RAB.x, LA.x, LB.x, 0)
    Sy = E(a, b, RAB.y, LA.y, LB.y, 0)
    Sz = E(a, b, RAB.z, LA.z, LB.z, 0)
    return (π / p)^(1.5) * Sx * Sy * Sz
end

function int1e_kin(
        a::Float64,
        LA::Vec3{Int},
        b::Float64,
        LB::Vec3{Int},
        RAB::Vec3{Float64},
    )
    p = a + b
    Sx = E(a, b, RAB.x, LA.x, LB.x)
    Sy = E(a, b, RAB.y, LA.y, LB.y)
    Sz = E(a, b, RAB.z, LA.z, LB.z)
    ∇∇Sx = ∇∇E(a, b, RAB.x, LA.x, LB.x)
    ∇∇Sy = ∇∇E(a, b, RAB.y, LA.y, LB.y)
    ∇∇Sz = ∇∇E(a, b, RAB.z, LA.z, LB.z)
    return -0.5 * (π / p)^(1.5) *  (∇∇Sx * Sy * Sz + Sx * ∇∇Sy * Sz + Sx * Sy * ∇∇Sz)
end

"Auxiliary Hermite integrals"
function R(t::Int, u::Int, v::Int, n::Int, p::Float64, RPC::Vec3{Float64})
    if t == u == v == 0 # base case
        return (-2 * p)^n * boys(n, p * (RPC[1]*RPC[1] + RPC[2]*RPC[2] + RPC[3]*RPC[3]))
    elseif u == v == 0 # decrement index t
        if t == 1
            return RPC[1] * R(0, 0, 0, n+1, p, RPC)
        end

        return (t-1) * R(t-2, 0, 0, n+1, p, RPC) +
               RPC[1] * R(t-1, 0, 0, n+1, p, RPC)
    elseif v == 0 # decrement index u
        if u == 1
            return RPC[2] * R(t, 0, 0, n+1, p ,RPC)
        end

        return (u-1) * R(t, u-2, 0, n+1, p, RPC) +
               RPC[2] * R(t, u-1, 0, n+1, p, RPC)
    else # decrement index v
        if v == 1
            return RPC[3] * R(t, u, 0, n+1, p, RPC)
        end

        return (v-1) * R(t, u, v-2, n+1, p, RPC) +
               RPC[3] * R(t, u, v-1, n+1, p, RPC)
    end
end

function int1e_rinv(
        a::Float64,
        LA::Vec3{Int},
        b::Float64,
        LB::Vec3{Int},
        RAB::Vec3{Float64},
        RPC::Vec3{Float64},
    )
    p = a + b
    ret = 0.0
    for t = 0:LA.x + LB.x, u = 0:LA.y + LB.y, v = 0:LA.z + LB.z
        ret += E(a, b, RAB.x, LA.x, LB.x, t) *
               E(a, b, RAB.y, LA.y, LB.y, u) *
               E(a, b, RAB.z, LA.z, LB.z, v) *
               R(t, u, v, 0, p, RPC)
    end
    return 2 * π / p * ret
end

function int2e(
        a::Float64,
        LA::Vec3{Int},
        b::Float64,
        LB::Vec3{Int},
        c::Float64,
        LC::Vec3{Int},
        d::Float64,
        LD::Vec3{Int},
        RAB::Vec3{Float64},
        RCD::Vec3{Float64},
        RPQ::Vec3{Float64},
    )
    p = a + b
    q = c + d
    ω = p * q / (p + q)
    ret = 0.0
    for t = 0:LA.x + LB.x, u = 0:LA.y + LB.y, v = 0:LA.z + LB.z
        for τ = 0:LC.x + LD.x, μ = 0:LC.y + LD.y, υ = 0:LC.z + LD.z
            ret += E(a, b, RAB.x, LA.x, LB.x, t) *
                   E(a, b, RAB.y, LA.y, LB.y, u) *
                   E(a, b, RAB.z, LA.z, LB.z, v) *
                   E(c, d, RCD.x, LC.x, LD.x, τ) *
                   E(c, d, RCD.y, LC.y, LD.y, μ) *
                   E(c, d, RCD.z, LC.z, LD.z, υ) *
                   R(t+τ, u+μ, v+υ, 0, ω, RPQ) *
                   (-1)^(τ+μ+υ)
        end
    end
    ret *= 2 * π^2.5 / (p * q * sqrt(p + q))
    return ret
end


text:
module Salsa

# Document this Module via the README.md file.
@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), "```julia" => "```jldoctest")
end Salsa

# This is the entirety of the high-level Salsa API: Users create inputs and derived
# functions, and access them through a Runtime instance.
export @derived, @declare_input, Runtime, DerivedFunctionException

import MacroTools
# For registering DerivedFunctionException as a wrapped exception.
import ExceptionUnwrapping

using Base: @lock

include("Debug.jl")
using .Debug

# Computational results and demand are identified via dependency keys.
include("dependency_keys.jl")

# Things go wrong, it happens.
include("exceptions.jl")

# For debugging, Salsa can keep its own stack trace with additional information. This is
# also used to detect dependency cycles.
include("stack.jl")

# A trace data structure is used to keep track of dependencies as a derived function
# executes. Traces are somewhat heavy to allocate, so we keep a pre-allocated pool of them.
include("trace.jl")
include("trace_pool.jl")

# Salsa supports pluggable storage backends. These are responsible for storing the inputs,
# derived values, and dependency information discovered by Salsa.
include("storage_api.jl")

# The runtime API is Salsa's internal API. Different runtime types exist for different
# purposes.
include("runtime_api.jl")
include("runtime_generic.jl")
include("runtime_top_level.jl")
include("runtime_tracing.jl")

# --- Macro utils -----

# Return an array of the Symbol names of every argument in a function definition. To ensure
# every argument is named, for unnamed arguments (e.g. `_` or `::Int`), we _generate_ a name
# via `gensym`.
function _argnames(args)
    [
        name === nothing || name === :_ ? gensym("_$i") : name
        for (i, name) in enumerate(first.(map(MacroTools.splitarg, args)))
    ]
end

# Return an array of the Types of every argument in a function definition.
function _argtypes(args)
    getindex.(map(MacroTools.splitarg, args), Ref(2))
end

# ---------------

# Salsa's high-level API consists of two macros, one for declaring derived functions, the
# other for declaring inputs.
include("derived_macro.jl")
include("declare_input_macro.jl")

# Provide Storage-major view on Runtime types for ease of implementing Storage backends.
const RuntimeWithStorage{ST,CT} = Runtime{CT,ST}
const _TopLevelRuntimeWithStorage{ST,CT} = _TopLevelRuntime{CT,ST}
const _TracingRuntimeWithStorage{ST,CT} = _TracingRuntime{CT,ST}

# Default context, manual storage
RuntimeWithStorage{ST}(st = ST()) where {ST} = Runtime{EmptyContext,ST}(EmptyContext(), st)

include("inspect.jl")

# The Salsa package provides a default storage layer to back Salsa Runtimes, which
# implements a simple key-value store using in-memory dictionaries.
include("default_storage.jl")
using ._DefaultSalsaStorage: DefaultStorage, DefaultRuntime

# Simpler print function for default Runtime implementation.
function Base.show(io::IO, rt::_TopLevelRuntime{EmptyContext,DefaultStorage})
    print(io, "Salsa.Runtime($(rt.storage))")
end

function __init__()
    # Init the freelists at runtime based on the number of threads julia is configured with.
    _init_thread_local_pools_and_freelists()
end

end  # module


text:
using Colors, Base.Test

# showcompact
tomato_bisque =
           compose(context(),
                   (context(), circle(), fill(colorant"bisque")),
                   (context(), rectangle(), fill(colorant"tomato")))

io = IOBuffer()
showcompact(io, tomato_bisque)
str = String(take!(io))
@test str == "Context(Context(R,f),Context(C,f))"

# Tagging
function points(xa, ya)
    [(x, y) for (x, y) in zip(xa, ya)]
end

pnts = points(rand(5), rand(5))
p = polygon(pnts, :mypoints)
@test p.tag == :mypoints

r = rectangle(0, 1, 0.5, 0.8, :box)
@test r.tag == :box
r = rectangle(rand(5),rand(5),rand(5),rand(5),:manybox)
@test r.tag == :manybox

c = circle(0, 0.8, 1.2, :circle)
@test c.tag == :circle
c = circle(rand(5), rand(5), rand(5), :data)
@test c.tag == :data

elps = ellipse(0, 0.8, 1.2, 1.5, :ellipse)
@test elps.tag == :ellipse
elps = ellipse(rand(5),rand(5),rand(5),rand(5),:manyellipse)
@test elps.tag == :manyellipse

txt = text(1.5, 15, "hello", tag=:hello)
@test txt.tag == :hello
txt = text(rand(5),rand(5),map(x->randstring(5), 1:5), tag=:random)
@test txt.tag == :random

ln = line(pnts, :line)
@test ln.tag == :line

crv = curve((0,0), (1,0.5), (0.2,0.3), (0.7,-2.4), :curve)
@test crv.tag == :curve
crv = curve(pnts, pnts, pnts, pnts, :manycurve)
@test crv.tag == :manycurve

bm = bitmap("fake", rand(UInt8,10), 0, 1, 0.8, 0.7, :image)
@test bm.tag == :image

# type definitions & constructors (issue #149)
@test isa(Compose.polygon(), Compose.Polygon)
@test isa(Compose.polygon([(1,2),(3,5),(4,2)]), Compose.Polygon)

@test isa(Compose.rectangle(), Compose.Rectangle)
@test isa(Compose.rectangle(0,1,0.3,0.8), Compose.Rectangle)
@test isa(Compose.rectangle(rand(5),rand(5),rand(5),rand(5)), Compose.Rectangle)

@test isa(Compose.circle(), Compose.Circle)
@test isa(Compose.circle(3.2,1.4,0.8), Compose.Circle)
@test isa(Compose.circle(rand(5), rand(5), rand(5)), Compose.Circle)

@test isa(Compose.ellipse(), Compose.Ellipse)
@test isa(Compose.ellipse(0.5,0.5,0.3,0.2), Compose.Ellipse)
@test isa(Compose.ellipse(rand(5), rand(5), rand(5), rand(5)), Compose.Ellipse)

@test isa(Compose.text(0.5,0.4,"hello"), Compose.Text)
@test isa(Compose.text(rand(5),rand(5),["hello","there"]), Compose.Text)

@test isa(Compose.line(), Compose.Line)
@test isa(Compose.line([(1,2),(3,5),(4,2)]), Compose.Line)

# issue #172: default circle(xs, ys, rs) radius measure is context units
@test isequal(circle([0.5], [0.5], [0.1]).primitives[1].radius, 0.1cx)

# Gadfly issue 857 and 436
# make sure that newlines are respected by `text_extents`
font_family = "'PT Sans Caption','Helvetica Neue','Helvetica',sans-serif"
@test text_extents(font_family, 8pt, "test test")[1][2]*2 ==
                        text_extents(font_family, 8pt, "test\ntest")[1][2]


text:
Base.mean(X::NullableArray; skipnull::Bool = false) =
    sum(X; skipnull = skipnull) /
        Nullable(length(X.isnull) - (skipnull * countnz(X.isnull)))

function Base.mean{T, W, V}(X::NullableArray{T}, w::WeightVec{W, V};
                            skipnull::Bool=false)
    if skipnull
        _X = NullableArray(X.values .* w.values, X.isnull)
        _w = NullableArray(w.values, X.isnull)
        return sum(_X; skipnull=true) / sum(_w; skipnull=true)
    else
        anynull(X) ? Nullable{T}() : Nullable(mean(X.values, w))
    end
end

function Base.mean{T, W, V<:NullableArray}(X::NullableArray{T},
                                           w::WeightVec{W, V};
                                           skipnull::Bool=false)
    if skipnull
        _X = X .* w.values
        _w = NullableArray(w.values, _X.isnull)
        return sum(_X; skipnull=true) / sum(_w; skipnull=true)
    else
        anynull(X) || anynull(w) ? Nullable{T}() :
                                   Nullable(mean(X.values, w.values.values))
    end
end

function Base.varm{T}(X::NullableArray{T}, m::Number; corrected::Bool=true,
                      skipnull::Bool=false)
    if skipnull
        n = length(X)

        nnull = countnz(X.isnull)
        nnull == n && return Nullable(convert(Base.momenttype(T), NaN))
        nnull == n-1 && return Nullable(
            convert(Base.momenttype(T),
                    abs2(X.values[Base.findnextnot(X.isnull, 1)] - m)/(1 - Int(corrected))
            )
        )
        /(nnull == 0 ? Nullable(Base.centralize_sumabs2(X.values, m, 1, n)) :
                       mapreduce_impl_skipnull(Base.CentralizedAbs2Fun(m),
                                               Base.AddFun(), X),
          Nullable(n - nnull - Int(corrected))
        )
    else
        any(X.isnull) && return Nullable{T}()
        Nullable(Base.varm(X.values, m; corrected=corrected))
    end
end

function Base.varm{T, U<:Number}(X::NullableArray{T}, m::Nullable{U};
                                 corrected::Bool=true, skipnull::Bool=false)
    m.isnull && throw(NullException())
    return varm(X, m.value; corrected=corrected, skipnull=skipnull)
end

function Base.varzm{T}(X::NullableArray{T}; corrected::Bool=true,
                       skipnull::Bool=false)
    n = length(X)
    nnull = skipnull ? countnz(X.isnull) : 0
    (n == 0 || n == nnull) && return Nullable(convert(Base.momenttype(T), NaN))
    return sumabs2(X; skipnull=skipnull) /
           Nullable((n - nnull - Int(corrected)))
end

function Base.var(X::NullableArray; corrected::Bool=true, mean=nothing,
         skipnull::Bool=false)

    (anynull(X) & !skipnull) && return Nullable{eltype(X)}()

    if mean == 0 || isequal(mean, Nullable(0))
        return Base.varzm(X; corrected=corrected, skipnull=skipnull)
    elseif mean == nothing
        return varm(X, Base.mean(X; skipnull=skipnull); corrected=corrected,
             skipnull=skipnull)
    elseif isa(mean, Union{Number, Nullable})
        return varm(X, mean; corrected=corrected, skipnull=skipnull)
    else
        error()
    end
end

function Base.stdm(X::NullableArray, m::Number;
                   corrected::Bool=true, skipnull::Bool=false)
    return sqrt(varm(X, m; corrected=corrected, skipnull=skipnull))
end

function Base.stdm{T<:Number}(X::NullableArray, m::Nullable{T};
                              corrected::Bool=true, skipnull::Bool=false)
    return sqrt(varm(X, m; corrected=corrected, skipnull=skipnull))
end

function Base.std(X::NullableArray; corrected::Bool=true,
                              mean=nothing, skipnull::Bool=false)
    return sqrt(var(X; corrected=corrected, mean=mean, skipnull=skipnull))
end

max_stars_repo_name: davidagold/NullableStats

max_issues_repo_name: davidagold/NullableStats

max_forks_repo_name: davidagold/NullableStats