text:
using GLFW

struct Frame
    mWindow

    function Frame(parent)
        this = new(parent)
        return this
    end
end

@class Frame

mutable struct GlfwWindow <: IFeature
    mDevice::Device
    mChainFormat::vk.VkSurfaceFormatKHR

    mWidth::UInt32
    mHeight::UInt32
    mResizable::Bool
    mTitle::String
    mWindow::GLFW.Window
    mSurface::vk.VkSurfaceKHR
    mChain::vk.VkSwapchainKHR

    mImageReady::vk.VkSemaphore
    mRenderingComplete::vk.VkSemaphore

    mPresentIndex::UInt32
    mQueue::Queue

    mChainImages::Vector{Image}
    mChainViews::Vector{ImageView}
    mSwapchainHandler#::SwapchainBuildHandler

    mSwapchainInfo::vk.VkSwapchainCreateInfoKHR

    mFrameProxy::Frame

    function GlfwWindow(device::Device, format::vk.VkSurfaceFormatKHR,
                        width::UInt32, height::UInt32, resizable::Bool, title::String)
        this = new()
        this.mDevice = device
        this.mChain = C_NULL
        this.mChainFormat = format
        this.mWidth = width
        this.mHeight = height
        this.mResizable = resizable
        this.mTitle = title

        this.mChainImages = Vector{Image}()
        this.mChainViews = Vector{ImageView}()
        this.mSwapchainHandler = nothing

        GLFW.WindowHint(GLFW.CLIENT_API, GLFW.NO_API)
        GLFW.WindowHint(GLFW.VISIBLE, true)
        GLFW.WindowHint(GLFW.RESIZABLE, this.mResizable ? true : false)
        this.mWindow = GLFW.CreateWindow(width, height, title)
        GLFW.SetWindowSizeCallback(this.mWindow, (winHandle, w, h) -> begin
            # TODO: if check does not work
            # if this.mWindow == winHandle
                this.onResize(UInt32(w), UInt32(h))
            # end
        end)

        this.mQueue = namedQueue(this.mDevice, "present")
        this.mSurface = GLFW.CreateWindowSurface(getInstance(device), this.mWindow)
        vkDevice = getLogicalDevice(this.mDevice)
        this.mImageReady = VkExt.createSemaphore(vkDevice)
        this.mRenderingComplete = VkExt.createSemaphore(vkDevice)

        this.mFrameProxy = Frame(this)

        return this
    end
end

@class GlfwWindow

function LavaCore.:destroy(this::GlfwWindow)
    destroy!(this.mChainViews)

    vkDevice = getLogicalDevice(this.mDevice)
    if isdefined(this, :mImageReady)
        vk.vkDestroySemaphore(vkDevice, this.mImageReady, C_NULL)
    end
    if isdefined(this, :mRenderingComplete)
        vk.vkDestroySemaphore(vkDevice, this.mRenderingComplete, C_NULL)
    end

    if this.mChain != C_NULL
        vk.vkDestroySwapchainKHR(vkDevice, this.mChain, C_NULL)
    end

    vk.vkDestroySurfaceKHR(getInstance(this.mDevice), this.mSurface, C_NULL)

    if isdefined(this, :mWindow)
        GLFW.DestroyWindow(this.mWindow)
    end
end

function imageIndex(this::Frame)::UInt32
    return this.mWindow.mPresentIndex
end

function imageReady(this::Frame)::vk.VkSemaphore
    return this.mWindow.mImageReady
end

function renderingComplete(this::Frame)::vk.VkSemaphore
    return this.mWindow.mRenderingComplete
end

function setSize(this::GlfwWindow, width::UInt32, height::UInt32)
    this.mWidth = width
    this.mHeight = height
end

function getSurface(this::GlfwWindow)::vk.VkSurfaceKHR
    return this.mSurface
end

function getSurfaceFormat(this::GlfwWindow)::vk.VkSurfaceFormatKHR
    return this.mChainFormat
end

function getWidth(this::GlfwWindow)
    return this.mWidth
end

function getHeight(this::GlfwWindow)
    return this.mHeight
end

function buildSwapchainWith(this::GlfwWindow, handler)
    this.mSwapchainHandler = handler
    buildSwapchain(this)
end

function onResize(this::GlfwWindow, w::UInt32, h::UInt32)
    this.mWidth = w
    this.mHeight = h
    # TODO Executing build while dragging will get bad performance
    this.buildSwapchain()
end

function buildSwapchain(this::GlfwWindow)
    destroy!(this.mChainImages)
    destroy!(this.mChainViews)
    vkDevice = getLogicalDevice(this.mDevice)

    if this.mChain != C_NULL
        vk.vkDestroySwapchainKHR(vkDevice, this.mChain, C_NULL)
    end
    phyDevice = getPhysicalDevice(this.mDevice)
    capRef = Ref{vk.VkSurfaceCapabilitiesKHR}()
    vk.vkGetPhysicalDeviceSurfaceCapabilitiesKHR(phyDevice, this.mSurface, capRef)
    cap = capRef[]
    setSize(this, cap.currentExtent.width, cap.currentExtent.height)

    pres = bestMode(phyDevice, this.mSurface)

    supp = Ref{vk.VkBool32}(false)
    vk.vkGetPhysicalDeviceSurfaceSupportKHR(phyDevice, family(this.mQueue), this.mSurface, supp)
    if(supp[] != vk.VK_TRUE)
        error("The selected queue family can't present to this device.")
    end

    imageSharingMode = 0
    queueFamilyIndexCount = 0
    queueFamilyIndices = C_NULL
    families = 0
    if family(graphicsQueue(this.mDevice)) == family(this.mQueue)
        imageSharingMode = vk.VK_SHARING_MODE_EXCLUSIVE
    else
        imageSharingMode = vk.VK_SHARING_MODE_CONCURRENT
        families = [family(graphicsQueue(this.mDevice)), family(this.mQueue)]
        queueFamilyIndexCount = 2
        queueFamilyIndices = pointer(families)
    end

    createInfo = Ref(vk.VkSwapchainCreateInfoKHR(
        vk.VK_STRUCTURE_TYPE_SWAPCHAIN_CREATE_INFO_KHR, #sType
        C_NULL, #pNext
        vk.VkFlags(0), #flags
        this.mSurface, #surface
        cap.minImageCount, #minImageCount
        this.mChainFormat.format, #imageFormat
        this.mChainFormat.colorSpace, #imageColorSpace
        cap.currentExtent, #imageExtent
        1, #imageArrayLayers
        vk.VK_IMAGE_USAGE_COLOR_ATTACHMENT_BIT | vk.VK_IMAGE_USAGE_TRANSFER_DST_BIT, #imageUsage
        imageSharingMode, #imageSharingMode::VkSharingMode
        queueFamilyIndexCount, #queueFamilyIndexCount::UInt32
        queueFamilyIndices, #pQueueFamilyIndices::Ptr{UInt32}
        cap.currentTransform, #preTransform::VkSurfaceTransformFlagBitsKHR
        vk.VK_COMPOSITE_ALPHA_OPAQUE_BIT_KHR, #compositeAlpha::VkCompositeAlphaFlagBitsKHR
        pres, #presentMode::VkPresentModeKHR
        vk.VK_TRUE, #clipped::VkBool32
        vk.VK_NULL_HANDLE #oldSwapchain::VkSwapchainKHR
    ))

    vkDevice = getLogicalDevice(this.mDevice)
    swapChain = Ref{vk.VkSwapchainKHR}()
    GC.@preserve createInfo families begin
        if (vk.vkCreateSwapchainKHR(vkDevice, createInfo, C_NULL, swapChain) != vk.VK_SUCCESS)
            error("Failed to create swap chain!")
        end
    end
    this.mChain = swapChain[]
    chainImageHandles = VkExt.getSwapchainImagesKHR(vkDevice, this.mChain)
    imgCreateInfo = attachment2D(phyDevice, this.mWidth, this.mHeight, this.mChainFormat.format)
    for handle::vk.VkImage in chainImageHandles
        image = Image(this.mDevice, handleRef(imgCreateInfo)[], handle, vk.VK_IMAGE_VIEW_TYPE_2D)
        push!(this.mChainViews, createView(image))
        push!(this.mChainImages, image)
    end

    # TODO Why duplicate?
    # for img in this.mChainImages
    #     push!(this.mChainViews, createView(img))
    # end

    this.mSwapchainInfo = createInfo[]
    if this.mSwapchainHandler != nothing
        this.mSwapchainHandler(this.mChainViews)
    end
end

function bestMode(phyDev::vk.VkPhysicalDevice, surf::vk.VkSurfaceKHR)::vk.VkPresentModeKHR
    prios = [vk.VK_PRESENT_MODE_MAILBOX_KHR,
            # Fifo-Modes are bugged on nvidia currently, so rather accept tearing
            # than the whole system freezing
            vk.VK_PRESENT_MODE_IMMEDIATE_KHR,
            vk.VK_PRESENT_MODE_FIFO_KHR,
            vk.VK_PRESENT_MODE_FIFO_RELAXED_KHR]

    modes = VkExt.getSurfacePresentModesKHR(phyDev, surf)
    for target in prios
        for candidate in modes
            if (target == candidate)
                return target
            end
        end
    end
    return modes[1]
end

function startFrame(this::GlfwWindow)::Frame
    @assert this.mChain != C_NULL "You need to provide a handler for swapchain creation."

    while (true)
        vkDevice = getLogicalDevice(this.mDevice)
        imageIndex = Ref{UInt32}()
        res = vk.vkAcquireNextImageKHR(vkDevice, this.mChain, 1e9, this.mImageReady, C_NULL, imageIndex)
        this.mPresentIndex = imageIndex[]
        if (res == vk.VK_TIMEOUT)
            error("GlfwWindow::startFrame(): acquireNextImage timed out (>1s)")
            continue
        end

        if (res == vk.VK_ERROR_OUT_OF_DATE_KHR || res == vk.VK_SUBOPTIMAL_KHR)
            vk.vkDeviceWaitIdle(vkDevice)
            buildSwapchain(this)
            continue;
        end
        return this.mFrameProxy
    end
end

# When ~Frame is called
function endFrame(this::Frame)
    if (this.mWindow == nothing)
        return
    end

    chains = [this.mWindow.mChain]
    semaphores = [renderingComplete(this)]
    indices = [imageIndex(this)]
    preserves = [chains, semaphores, indices]

    info = Ref(vk.VkPresentInfoKHR(
        vk.VK_STRUCTURE_TYPE_PRESENT_INFO_KHR, #sType::VkStructureType
        C_NULL, #pNext::Ptr{Cvoid}
        1, #waitSemaphoreCount::UInt32
        pointer(semaphores), #pWaitSemaphores::Ptr{VkSemaphore}
        1, #swapchainCount::UInt32
        pointer(chains), #pSwapchains::Ptr{VkSwapchainKHR}
        pointer(indices), #pImageIndices::Ptr{UInt32}
        C_NULL #pResults::Ptr{VkResult}
    ))

    GC.@preserve info preserves semaphores indices chains begin
        err = vk.vkQueuePresentKHR(LavaCore.handle(this.mWindow.mQueue), info)
        vk.vkDeviceWaitIdle(getLogicalDevice(this.mWindow.mDevice))
        if err == vk.VK_ERROR_OUT_OF_DATE_KHR
            buildSwapchain(this.mWindow)
        end
    end
end


text:
module NetworkDynamicalSystems

using Graphs
using DynamicalSystems
using NetworkDynamics

end

ext: jl


text:
println("here")

using MetaGen
using JSON
import YAML
using Pipe: @pipe
using Random

include("useful_functions.jl")

path = "/Users/marleneberke/Documents/03_Yale/Projects/001_Mask_RCNN/scratch_work_07_16_21/09_18/"

dict = @pipe (path * "data_labelled.json") |> open |> read |> String |> JSON.parse

################################################################################
#Set up the observations
num_videos = 4 #total, including test set
num_frames = 20
threshold = 0.05

n_top = 5
params = Video_Params(n_possible_objects = 5)

receptive_fields = make_receptive_fields(params)
objects_observed, camera_trajectories = make_observations_office(dict, receptive_fields, num_videos, num_frames, threshold, n_top)

################################################################################
#Set up for MCMC
num_particles = 10
mcmc_steps_outer = 1
mcmc_steps_inner = 1

################################################################################
#Online MetaGen
num_videos_train = convert(Int64, num_videos/2)
@assert num_videos_train==50

if shuffle_type==0
	order = collect(1:num_videos_train)
elseif shuffle_type==1
	order = vcat(reverse(collect(26:50)), reverse(collect(1:25)))
elseif shuffle_type==2
	order = vcat(collect(26:50), collect(1:25))
elseif shuffle_type==3
	order = vcat(reverse(collect(1:25)), reverse(collect(26:50)))
end

training_objects_observed = objects_observed[order, :]
training_camera_trajectories = camera_trajectories[order, :]

#Set up the output files
online_V_file = open(path * "online_V.csv", "w")
file_header_V(online_V_file, params)
online_ws_file = open(path * "online_ws.csv", "w")
file_header_ws(online_ws_file, params, num_particles)

traces, inferred_world_states, avg_v = unfold_particle_filter(nothing,
	num_particles, mcmc_steps_outer, mcmc_steps_inner, training_objects_observed,
	training_camera_trajectories, params, online_V_file, online_ws_file, order)
close(online_V_file)
close(online_ws_file)
#
println("avg_v ", avg_v)
println("done with pf for online")

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

# #Retrospective MetaGen
#
#Set up the output file
retro_V_file = open(path * "retro_V.csv", "w")
file_header_V(retro_V_file, params)
retro_ws_file = open(path * "retro_ws.csv", "w")
file_header_ws(retro_ws_file, params, num_particles)

#training set and test set
order = vcat(order, 51:100)
input_objects_observed = vcat(objects_observed[order, :])
input_camera_trajectories = vcat(camera_trajectories[order, :])

traces, inferred_world_states, avg_v = unfold_particle_filter(avg_v, num_particles, mcmc_steps_outer, mcmc_steps_inner,
	input_objects_observed, input_camera_trajectories, params, retro_V_file, retro_ws_file, order)
close(retro_V_file)
close(retro_ws_file)

# # ################################################################################
# # #Lesioned MetaGen
# #
# #Set up the output file
# lesioned_V_file = open(path * "lesioned_V.csv", "w")
# file_header_V(lesioned_V_file, params)
# lesioned_ws_file = open(path * "lesioned_ws.csv", "w")
# file_header_ws(lesioned_ws_file, params, num_particles)
#
# v = zeros(length(params.possible_objects), 2)
# v[:,1] .= 1.0
# v[:,2] .= 0.5
# unfold_particle_filter(v, num_particles, mcmc_steps_outer, mcmc_steps_inner,
# 	objects_observed, camera_trajectories, params, lesioned_V_file, lesioned_ws_file)
# close(lesioned_V_file)
# close(lesioned_ws_file)
#
# println("done with pf for lesioned metagen")

################################################################################
#for writing an output file for a demo using Retro MetaGen. will only make sense for non-mixed up version

#undor re-ordering of inferred_world_states
inferred_world_states = inferred_world_states[order]

###### add to dictionary
out = write_to_dict(dict, camera_trajectories, inferred_world_states, num_videos, num_frames)

#open("../../scratch_work_07_16_21/output_tiny_set_detections.json","w") do f
open(path * "output.json","w") do f
	JSON.print(f,out)
end

println("finished writing json")

max_stars_repo_name: zhangir-azerbayev/MetaGen

max_issues_repo_name: zhangir-azerbayev/MetaGen

max_forks_repo_name: zhangir-azerbayev/MetaGen


text:
using HDF5;
using PyPlot;

#*******************************************************************************

primeira_camada = h5open("dados_split_operator.h5","r") do file
    names(file)
end

segunda_camada = h5open("dados_split_operator.h5","r") do file
    segunda_camada = [];
    for name in primeira_camada
        append!(segunda_camada,names(file[name]))
    end
    return segunda_camada
end

terceira_camada = h5open("dados_split_operator.h5","r") do file
    terceira_camada = [];
    for name in primeira_camada
        g = file[name];
        for name_grupo in segunda_camada
            append!(terceira_camada,names(g[name_grupo]))
        end
    end
    return terceira_camada
end

#*******************************************************************************

Atr_grupo = h5readattr("dados_split_operator.h5", "split_1D/potencial_degrau")

Atr_potencial = h5readattr("dados_split_operator.h5","split_1D/potencial_degrau/Potencial")

#*******************************************************************************

V_lido = h5open("dados_split_operator.h5","r") do file
    group = file[primeira_camada[1]];
    subgroup = group[segunda_camada[1]];
    read(subgroup["Potencial"])
end

x_lido = h5open("dados_split_operator.h5","r") do file
    group = file[primeira_camada[1]];
    subgroup = group[segunda_camada[1]];
    read(subgroup["grid"])
end

# plot(x_lido, V_lido)

#*******************************************************************************

terceira_camada[6]

nome_dados_real = terceira_camada[6]
nome_dados_imag = terceira_camada[5]

name_path_real = string("split_1D/potencial_degrau/",nome_dados_real)
name_path_imag = string("split_1D/potencial_degrau/",nome_dados_imag)

psi_lido_real = h5read("dados_split_operator.h5",name_path_real)
psi_lido_imag = h5read("dados_split_operator.h5",name_path_imag)

psi_lido_total = complex(psi_lido_real,psi_lido_imag)

PyPlot.matplotlib[:rc]("text", usetex=true) # allow tex rendering
PyPlot.matplotlib[:rc]("font", family="serif")

figure(2)
plot(x_lido, abs(psi_lido_total).^2)
plot(x_lido, 0.5*10.0^(-2)*V_lido, linewidth=2.0);
title("Gaussian distribution");
xlabel(L"$x$");
ylabel(L"$|\psi(x)|^2$");
# ylim(-0.05,1.05);
# xlim(-20,20);

max_stars_repo_name: mhlmedeiros/Split-Operator

max_issues_repo_name: mhlmedeiros/Split-Operator

max_forks_repo_name: mhlmedeiros/Split-Operator


text:
export i_mapfoldl, i_filter!, i_map!

"""
    i_mapfoldl(map, fold, out!, iter)

Reversible `mapfoldl` function, `map` can be irreversible, but `fold` should be reversible.
"""
@i function i_mapfoldl(map, fold, out!::T, iter) where T
    anc ← zero(T)
    for i=1:length(iter)
        anc += map(iter[i])
        fold(out!, anc)
        anc -= map(iter[i])
    end
end

"""
    i_filter!(f, out!, iter)

Reversible `filter` function, `out!` is an emptied vector.
"""
@i function i_filter!(f, out!::AbstractVector, x::AbstractVector{T}) where T
    @invcheckoff @inbounds for i = 1:length(x)
        if (f(x[i]), ~)
            COPYPUSH!(out!, x[i])
        end
    end
end


text:
function mkscale(scalevec::Vector{Float64})
    mult = eye(4)
    for i in 1:size(scalevec, 1)
        mult[i, i] = scalevec[i]
    end
    mult
end

function mktranslate(transvec::Vector{Float64})
    mult = eye(4)
    for i in 1:3
        mult[i, 4] = transvec[i]
    end
    mult
end

function mkrotate(angle::Float64, direction::Symbol)
    x = [cosd(angle), -sind(angle), sind(angle), cosd(angle)]
    indices = ROTATIONINDICES[direction]
    mult = eye(4)
    for i in 1:size(indices, 1)
        mult[indices[i, 1], indices[i, 2]] = x[i]
    end
    mult
end

max_stars_repo_name: AnAverageHuman/jlgraphics

max_issues_repo_name: AnAverageHuman/jlgraphics

max_forks_repo_name: AnAverageHuman/jlgraphics


text:
## LINEAR LEARNING NETWORKS (FOR INTERNAL USE ONLY)

# returns Model supertype - or `missing` if arguments are incompatible
function kind_(is_supervised, is_probabilistic)
    if is_supervised
        if ismissing(is_probabilistic) || !is_probabilistic
            return :DeterministicComposite
        else
            return :ProbabilisticComposite
        end
    else
        if ismissing(is_probabilistic) || !is_probabilistic
            return :UnsupervisedComposite
        else
            return missing
        end
    end
end

# Here `model` is a function or Unsupervised model (possibly
# `Static`).  This function wraps `model` in an appropriate `Node`,
# `node`, and returns `(node, mode.machine)`.
function node_(model, X)
    if model isa Unsupervised
        if model isa Static
            mach = machine(model)
        else
            mach = machine(model, X)
        end
        return transform(mach, X), mach
    else
        n = node(model, X)
        return n, n.machine
    end
end

# `models_and_functions` can include both functions (static
# operatations) and bona fide models (including `Static`
# transformers). If `ys == nothing` the learning network is assumed to
# be unsupervised. Otherwise: If `target === nothing`, then no target
# transform is applied; if `target !== nothing` and `inverse ===
# nothing`, then corresponding `transform` and `inverse_transform` are
# applied; if neither `target` nor `inverse` are `nothing`, then both
# `target` and `inverse` are assumed to be StaticTransformations, to
# be applied respectively to `ys` and to the output of the
# `models_and_functions` pipeline. Note that target inversion is
# applied to the output of the *last* nodal machine in the pipeline,
# corresponding to the last element of `models_and_functions`.

# No checks whatsoever are performed. Returns a learning network.
function linear_learning_network(Xs,
                                 ys,
                                 ws,
                                 target,
                                 inverse,
                                 models_and_functions...)

    n = length(models_and_functions)

    if ys !== nothing && target !== nothing
        yt, target_machine = node_(target, ys)
    else
        yt  = ys
    end

    if ws !== nothing
        tail_args = (yt, ws)
    else
        tail_args = (yt,)
    end

    nodes = Vector(undef, n + 1)
    nodes[1] = Xs

    for i = 2:(n + 1)
        m = models_and_functions[i-1]
        if m isa Supervised
            supervised_machine = machine(m, nodes[i-1], tail_args...)
            nodes[i] = predict(supervised_machine, nodes[i-1])
       else
            nodes[i] = node_(m, nodes[i-1]) |> first
        end
    end

    if target === nothing
        terminal_node=nodes[end]
    else
        if inverse === nothing
            terminal_node = inverse_transform(target_machine, nodes[end])
        else
            terminal_node = node_(inverse, nodes[end]) |> first
        end
    end

    return terminal_node

end

## PREPROCESSING

function eval_and_reassign_deprecated(modl, ex)
    s = gensym()
    evaluated = modl.eval(ex)
    modl.eval(:($s = $evaluated))
    return s, evaluated
end

# does expression processing, syntax and semantic
# checks. is_probabilistic is `true`, `false` or `missing`.
function pipeline_preprocess_deprecated(
    modl, ex, is_probabilistic::Union{Missing,Bool})

    ex isa Expr || pipe_alert(2)
    length(ex.args) > 1 || pipe_alert(3)
    ex.head == :call || pipe_alert(4)
    pipetype_ = ex.args[1]
    pipetype_ isa Symbol || pipe_alert(5)

    fieldnames_ = []               # for `from_network`
    models_ = []                   # for `from_network`
    models_and_functions_ = []     # for `linear_learning_network`
    models_and_functions  = []     # for `linear_learning_network`
    trait_value_given_name_ = Dict{Symbol,Any}()
    for ex in ex.args[2:end]
        if ex isa Expr
            if ex.head == :kw
                variable_ = ex.args[1]
                variable_ isa Symbol || pipe_alert(8)
                value_, value = eval_and_reassign_deprecated(modl, ex.args[2])
                if variable_ == :target
                    if value isa Function
                        value_, value =
                            eval_and_reassign_deprecated(modl,
                                              :(MLJBase.WrappedFunction($value)))
                    end
                    value isa Unsupervised ||
                        pipe_alert("Got $value where a function or "*
                                   "Unsupervised model instance was "*
                                   "expected. ")
                    target_ = value_
                    target = value
                    push!(fieldnames_, :target)
                elseif variable_ == :inverse
                    if value isa Function
                        value_, value =
                            eval_and_reassign_deprecated(modl,
                                              :(MLJBase.WrappedFunction($value)))
                    else
                        pipe_alert(10)
                    end
                    inverse_ = value_
                    inverse = value
                    push!(fieldnames_, :inverse)
                else
                    value isa Model ||
                        throw(ArgumentError("$value given where `Model` "*
                                            "instance expected. "))
                    push!(models_and_functions_, value_)
                    push!(models_and_functions, value)
                    push!(fieldnames_, variable_)
                end
                push!(models_, value_)
            else
                f = modl.eval(ex)
                f isa Function ||
                    pipe_alert("Perhaps a missing name, as `name=$f'?")
                push!(models_and_functions_, ex)
                push!(models_and_functions, f)
            end
        else
            f = modl.eval(ex)
            f isa Function || pipe_alert(7)
            push!(models_and_functions_, ex)
            push!(models_and_functions, f)
        end
    end

    (@isdefined target)  || (target = nothing; target_ = :nothing)
    (@isdefined inverse) || (inverse = nothing; inverse_ = :nothing)
    inverse !== nothing && target === nothing &&
        pipe_alert("You have specified `inverse=...` but no `target`.")

    supervised(m) = !(m isa Unsupervised) && !(m isa Function)

    supervised_components = filter(supervised, models_and_functions)

    length(supervised_components) < 2 ||
        pipe_alert("More than one component of the pipeline is a "*
                   "supervised model .")

    if length(supervised_components) == 1
        model = supervised_components[1]
        trait_value_given_name_[:supports_weights] =
            supports_weights(model)
    end

    is_supervised  =
        length(supervised_components) == 1

    # `kind_` is defined in /src/composition/models/from_network.jl
    kind = kind_(is_supervised, is_probabilistic)

    ismissing(kind) &&
        pipe_alert("Composite has no supervised components and so "*
                  "`prediction_type=:probablistic` (or "*
                   "`is_probabilistic=true`) "*
                  "declaration is not allowed. ")

    target isa MLJBase.WrappedFunction && inverse == nothing &&
        pipe_alert("It appears `target` is a function. "*
                   "You must therefore specify `inverse=...` .")

    target == nothing || is_supervised ||
        pipe_alert("`target=...` has been specified but no "*
                   "supervised components have been specified. ")

    return (pipetype_, fieldnames_, models_,
            models_and_functions_, target_, inverse_,
            kind, trait_value_given_name_)

end

function pipeline_preprocess_deprecated(modl, ex, kw_ex)
    kw_ex isa Expr || pipe_alert(10)
    kw_ex.head == :(=) || pipe_alert(11)
    if kw_ex.args[1] == :is_probabilistic
        value_ = kw_ex.args[2]
        if value_ == :missing
            value = missing
        else
            value = value_
            value isa Bool ||
                pipe_alert("`is_probabilistic` can only be `true` or `false`.")
        end
    elseif kw_ex.args[1] == :prediction_type
        value_ = kw_ex.args[2].value
        if value_ == :probabilistic
            value = true
        elseif value_ == :deterministic
            value =false
        else
            pipe_alert("`prediction_type` can only be `:probabilistic` "*
                      "or `:deterministic`.")
        end
    else
        pipe_alert("Unrecognized keywork `$(kw_ex.args[1])`.")
    end
    return pipeline_preprocess_deprecated(modl, ex, value)
end

function pipeline_deprecated_(modl, ex, kw_ex)

    Base.depwarn("Using deprecated @pipeline syntax. "*
                 "New syntax:\n   @pipeline model1 model2 ... \n "*
                 "Query the docstring or see manual for options and "*
                 "further details. ", Symbol("@pipeline"))

    (pipetype_, fieldnames_, models_, models_and_functions_,
     target_, inverse_, kind, trait_value_given_name_) =
         pipeline_preprocess_deprecated(modl, ex, kw_ex)

    if kind === :UnsupervisedComposite
        ys_ = :nothing
    else
        ys_ = :(source(kind=:target))
    end

    if haskey(trait_value_given_name_, :supports_weights) &&
        trait_value_given_name_[:supports_weights]
        ws_ = :(source(kind=:weights))
    else
        ws_ = nothing
    end

    if kind == :DeterministicComposite
        super_type = Deterministic
    elseif kind == :ProbabilisticComposite
        super_type = Probabilistic
    else
        super_type = Unsupervised
    end

    mach_ = quote
        MLJBase.linear_learning_network_machine($super_type,
                                                source(nothing), $ys_, $ws_,
                                                $target_, $inverse_,
                                                true,
                                                predict,
                                                $(models_and_functions_...))
    end

    field_declarations =
        [:($(fieldnames_[j])=$(models_[j])) for j in eachindex(models_)]

    supertype_ = _exported_type(modl.eval(mach_).model)

    struct_ = :(mutable struct $pipetype_ <: $supertype_
                        $(field_declarations...)
                    end)

    no_fields = length(models_) == 0

    from_network_(modl,
                  mach_,
                  pipetype_,
                  struct_,
                  no_fields,
                  trait_value_given_name_)

    # N_ex = quote

    #     MLJBase.linear_learning_network(source(nothing), $ys_, $ws_,
    #                                 $target_, $inverse_,
    #                                 $(models_and_functions_...))
    # end

    # from_network_(modl, pipetype_, fieldnames_, models_, N_ex,
    #               kind, trait_value_given_name_)

    return pipetype_

end

pipeline_deprecated_(modl, ex) =
    pipeline_deprecated_(modl, ex, :(is_probabilistic=missing))

macro pipeline_deprecated(exs...)

    pipetype_ = pipeline_deprecated_(__module__, exs...)

    esc(quote
        $pipetype_()
        end)

end


text:
"""
Configuration structure which allow the selection of specific parts of the Vulkan API.
"""
Base.@kwdef struct WrapperConfig
    "Include core API (with core extensions)."
    wrap_core::Bool = true
    "Include beta (provisional) exensions. Provisional extensions may break between patch releases."
    include_provisional_exts::Bool = false
    "Platform-specific families of extensions to include."
    include_platforms::Vector{PlatformType} = []
    "Path the wrapper will be written to."
    destfile::String
end

include_provisional_exts(config::WrapperConfig) = config.include_provisional_exts || PLATFORM_PROVISIONAL in config.include_platforms

function extensions(config::WrapperConfig)
    exts = filter(x -> x.is_provisional && include_provisional_exts(config) || x.platform in config.include_platforms || x.platform == PLATFORM_NONE && config.wrap_core, spec_extensions_supported)
end

function _filter_specs(specs, extensions, wrap_core)
    filter(specs) do spec
        ext = extension(spec)
        isnothing(ext) && wrap_core || ext in extensions
    end
end

filter_specs(config::WrapperConfig) = x -> _filter_specs(x, extensions(config), config.wrap_core)

abstract type Platform end

struct Linux <: Platform end
struct MacOS <: Platform end
struct BSD <: Platform end
struct Windows <: Platform end

WrapperConfig(p::Platform, destfile; kwargs...) = WrapperConfig(; include_platforms = platform_extensions(p), destfile, kwargs...)

platform_extensions(::Linux) = [PLATFORM_WAYLAND, PLATFORM_XCB, PLATFORM_XLIB, PLATFORM_XLIB_XRANDR]
platform_extensions(::MacOS) = [PLATFORM_MACOS, PLATFORM_METAL]
platform_extensions(::BSD) = [PLATFORM_WAYLAND, PLATFORM_XCB, PLATFORM_XLIB, PLATFORM_XLIB_XRANDR]
platform_extensions(::Windows) = [PLATFORM_WIN32]


text:
"""
Materialize a `GroupBy` into a `GroupedDataFrame`. Currently only eager
materialization is allowed, which requires two steps:

1. Materialize the source node into a `DataFrame`, then make a shallow copy of
    it that shares columns with `materialize(node.source)`.
2. Use `DataFrames.groupby` to generate the `GroupedDataFrame`.
"""
function materialize(node::GroupBy)
    # TODO: Performance optimization opportunity:
    #     Don't make a copy if all group keys are existing columns.
    tmp = copy(materialize(node.source), copycols=false)
    # TODO: Handle possible conflicts with existing columns here.
    for e in node.expressions
        tmp[e.alias] = generate_column(node.source, e)
    end
    DataFrames.groupby(tmp, [e.alias for e in node.expressions])
end


text:
module bases

import Base.==

export Basis, GenericBasis, CompositeBasis,
       tensor, ⊗, dualbasis, ptrace,
       equal_shape, equal_bases, multiplicable,
       IncompatibleBases,
       check_equal, check_multiplicable

"""
Abstract base class for all specialized bases.

The Basis class is meant to specify a basis of the Hilbert space of the
studied system. Besides basis specific information all subclasses must
implement a shape variable which indicates the dimension of the used
Hilbert space. For a spin 1/2 Hilbert space this would for example be the
vector Int[2]. A system composed of 2 spins could for example have a
shape vector Int[2 2].

Composite systems can be easily defined with help of the CompositeBasis class.
"""
abstract Basis

"""
An easy and fast to use general purpose basis.

Should only be used rarely since it defeats the purpose of checking that the
bases of state vectors and operators are correct for algebraic operations.
The cleaner way is to specify special bases for different systems, i.e.
there are FockBasis, MomentumBasis and PositionBasis.
"""
type GenericBasis <: Basis
    shape::Vector{Int}
end

GenericBasis(N::Int) = GenericBasis(Int[N])

"""
Basis for composite Hilbert spaces.

Simply stores the subbases in a vector and creates the shape vector directly
from the shape vectors of these subbases.
"""
type CompositeBasis <: Basis
    shape::Vector{Int}
    bases::Vector{Basis}
    CompositeBasis(bases::Basis...) = new([prod(b.shape) for b=bases], [bases...])
end

"""
Create composite bases.

Any given CompositeBasis is expanded so that the resulting CompositeBasis never
has another CompositeBasis as subbasis.
"""
tensor() = error("Tensor function needs at least one argument.")
tensor(b1::Basis, b2::Basis) = CompositeBasis(b1, b2)
tensor(b1::CompositeBasis, b2::CompositeBasis) = CompositeBasis(b1.bases..., b2.bases...)
tensor(b1::CompositeBasis, b2::Basis) = CompositeBasis(b1.bases..., b2)
tensor(b1::Basis, b2::CompositeBasis) = CompositeBasis(b1, b2.bases...)
tensor(bases::Basis...) = reduce(tensor, bases)
⊗(a,b) = tensor(a,b)

"""
Return the dual basis to the given basis.

If the basis is orthogonal the dual basis and the basis itself are identitcal.
"""
dualbasis(b::Basis) = b

"""
Total dimension of the Hilbert space.
"""
Base.length(b::Basis) = prod(b.shape)

"""
Check if two shape vectors are the same.

This implementation handles special cases in a very fast way so that their
overhead is very low in often performed operations.
"""
function equal_shape(a::Vector{Int}, b::Vector{Int})
    if a === b
        return true
    end
    if length(a) != length(b)
        return false
    end
    for i=1:length(a)
        if a[i]!=b[i]
            return false
        end
    end
    return true
end

"""
Check if two subbases vectors are identical.

This implementation handles special cases in a very fast way so that their
overhead is very low in often performed operations.
"""
function equal_bases(a::Vector{Basis}, b::Vector{Basis})
    if a===b
        return true
    end
    for i=1:length(a)
        if a[i]!=b[i]
            return false
        end
    end
    return true
end

==(b1::Basis, b2::Basis) = false
==(b1::GenericBasis, b2::GenericBasis) = equal_shape(b1.shape,b2.shape)
==(b1::CompositeBasis, b2::CompositeBasis) = equal_shape(b1.shape,b2.shape) && equal_bases(b1.bases,b2.bases)

"""
Check if two objects given in the specified bases are multiplicable.

In the case of orthogonal bases this mostly means that the two bases have to be
the same. For nonorthogonal bases this function can be overridden to account
for their behavior.
"""
multiplicable(b1::Basis, b2::Basis) = b1==dualbasis(b2)

function multiplicable(b1::CompositeBasis, b2::CompositeBasis)
    if !equal_shape(b1.shape,b2.shape)
        return false
    end
    for i=1:length(b1.shape)
        if !multiplicable(b1.bases[i], b2.bases[i])
            return false
        end
    end
    return true
end

"""
Exception that should be raised for an illegal algebraic operation.
"""
type IncompatibleBases <: Exception end

"""
Error throwing version of equality checking.

For a boolean version use the equality operator ==
"""
check_equal(b1::Basis, b2::Basis) = (b1==b2 ? true : throw(IncompatibleBases()))

"""
Error throwing version of multiplicativity checking.

For a boolean version use the function multiplicable.
"""
check_multiplicable(b1::Basis, b2::Basis) = (multiplicable(b1, b2) ? true : throw(IncompatibleBases()))

"""
Partial trace of a composite basis.
"""
function ptrace(b::CompositeBasis, indices::Vector{Int})
    reduced_basis = Basis[]
    for (i, subbasis) in enumerate(b.bases)
        if !(i in indices)
            push!(reduced_basis, subbasis)
        end
    end
    if length(reduced_basis)==0
        error("Nothing left.")
    elseif length(reduced_basis)==1
        return reduced_basis[1]
    else
        return CompositeBasis(reduced_basis...)
    end
end

end # module


text:
#=

  Fill-a-pix problem in Julia ConstraintSolver.jl

  From http://www.conceptispuzzles.com/index.aspx?uri=puzzle/fill-a-pix/basiclogic
  """
  Each puzzle consists of a grid containing clues in various places. The
  object is to reveal a hidden picture by painting the squares around each
  clue so that the number of painted squares, including the square with
  the clue, matches the value of the clue.
  """

  http://www.conceptispuzzles.com/index.aspx?uri=puzzle/fill-a-pix/rules
  """
  Fill-a-Pix is a Minesweeper-like puzzle based on a grid with a pixilated
  picture hidden inside. Using logic alone, the solver determines which
  squares are painted and which should remain empty until the hidden picture
  is completely exposed.
  """

  Fill-a-pix History:
  http://www.conceptispuzzles.com/index.aspx?uri=puzzle/fill-a-pix/history

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function fill_a_pix(problem,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))
    n,_ = size(problem)

    @variable(model, x[1:n,1:n], Bin)

    for i in 1:n, j in 1:n
        if problem[i,j] >= 0
            @constraint(model, problem[i,j] ==
                            sum([x[i+a,j+b] for a in -1:1, b in -1:1
                            if 1 <= i+a <= n &&
                               1 <= j+b <= n
                            ])
            )
        end
    end

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        # println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                # println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                print_matrix_str(x_val,".X")

            end
        end
    else
        println("status:$status")
    end

    return status
end

M = -1
fill_a_pix_problems = Dict(
    # Puzzle 1 from
    # http://www.conceptispuzzles.com/index.aspx?uri=puzzle/fill-a-pix/rules
    #
    :1 => resize_matrix([[M,M,M,M,M,M,M,M,0,M],
               [M,8,8,M,2,M,0,M,M,M],
               [5,M,8,M,M,M,M,M,M,M],
               [M,M,M,M,M,2,M,M,M,2],
               [1,M,M,M,4,5,6,M,M,M],
               [M,0,M,M,M,7,9,M,M,6],
               [M,M,M,6,M,M,9,M,M,6],
               [M,M,6,6,8,7,8,7,M,5],
               [M,4,M,6,6,6,M,6,M,4],
               [M,M,M,M,M,M,3,M,M,M]]
          ),

    # Puzzle 2 from
    # http://www.conceptispuzzles.com/index.aspx?uri=puzzle/fill-a-pix/rules
    #
    :2 => resize_matrix([[0,M,M,M,M,M,3,4,M,3],
                [M,M,M,4,M,M,M,7,M,M],
                [M,M,5,M,2,2,M,4,M,3],
                [4,M,6,6,M,2,M,M,M,M],
                [M,M,M,M,3,3,M,M,3,M],
                [M,M,8,M,M,4,M,M,M,M],
                [M,9,M,7,M,M,M,M,5,M],
                [M,M,M,7,5,M,M,3,3,0],
                [M,M,M,M,M,M,M,M,M,M],
                [4,4,M,M,2,3,3,4,3,M]]),

    # Puzzle from
    # http://www.conceptispuzzles.com/index.aspx?uri=puzzle/fill-a-pix/basiclogic
    #
    # Code: 030.15x15
    # ID: 03090000000
    #
    :3 => resize_matrix([[M,5,M,6,M,M,M,M,M,M,6,M,M,M,M],
                [M,M,7,6,M,4,M,M,4,M,M,8,9,M,5],
                [5,M,M,5,M,5,M,3,M,6,M,7,M,M,6],
                [4,M,2,M,4,M,4,M,3,M,2,M,M,9,M],
                [M,M,M,5,M,4,M,3,M,4,M,4,5,M,6],
                [M,4,3,3,4,M,M,M,4,M,2,M,M,M,M],
                [M,M,M,M,M,M,M,M,M,5,M,M,M,4,M],
                [3,M,3,M,M,3,M,M,M,5,M,4,4,M,M],
                [M,M,M,4,3,M,3,3,M,M,5,7,6,M,M],
                [4,M,M,M,2,M,3,3,2,M,8,9,M,5,M],
                [M,M,3,M,M,M,M,5,M,M,7,M,8,M,M],
                [4,M,M,3,2,M,M,M,M,M,7,M,M,6,M],
                [M,M,4,M,5,4,4,M,M,9,6,M,M,M,M],
                [M,3,5,7,M,6,M,M,M,M,M,M,7,M,M],
                [M,M,4,6,6,M,M,M,6,5,M,M,M,4,M]]),

)

for p in sort(collect(keys(fill_a_pix_problems)))
    println("\nproblem $p")
    @time fill_a_pix(fill_a_pix_problems[p])
end

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
using BenchmarkTools, SparseArrays, LinearAlgebra
using ExpmV

SUITE = BenchmarkGroup()

dimensions = 2 .^ (5:10)
density = 10. .^ (-4:-1)

for p in density
    SUITE[p] = BenchmarkGroup(["density"])
    for d in dimensions
        SUITE[p][d] = @benchmarkable ExpmV.expmv(rt,r,rv) setup=((rt,r,rv)=(rand(), sprandn(ComplexF64,$d,$d,$p), normalize(randn(ComplexF64,$d))))
    end
end


text:
function generate_action_defs(domain::Domain, state::State,
                              domain_type::Symbol, state_type::Symbol)
    available_def =
        generate_available(domain, state, domain_type, state_type)
    execute_def =
        generate_execute(domain, state, domain_type, state_type)
    all_action_defs = [available_def, execute_def]
    all_action_names = sortedkeys(get_actions(domain))
    all_action_types = Symbol[]
    for name in all_action_names
        action_defs =
            generate_action_defs(domain, state, domain_type, state_type, name)
        push!(all_action_types, action_defs[1])
        append!(all_action_defs, action_defs[2:end])
    end
    action_map = Expr(:tuple, (Expr(:(=), n, Expr(:call, act)) for (n, act)
                               in zip(all_action_names, all_action_types))...)
    get_actions_def = :(get_actions(::$domain_type) = $action_map)
    push!(all_action_defs, get_actions_def)
    return Expr(:block, all_action_defs...)
end

function generate_action_defs(domain::Domain, state::State,
                              domain_type::Symbol, state_type::Symbol,
                              action_name::Symbol)
    action_type = gensym("Compiled" * pddl_to_type_name(action_name) * "Action")
    action_typedef = :(struct $action_type <: CompiledAction end)
    groundargs_def =
        generate_groundargs(domain, state, domain_type, state_type,
                            action_name, action_type)
    available_def =
        generate_available(domain, state, domain_type, state_type,
                           action_name, action_type)
    execute_def =
        generate_execute(domain, state, domain_type, state_type,
                         action_name, action_type)
    method_defs =
        generate_action_methods(domain, state, domain_type, state_type,
                                action_name, action_type)
    return (action_type, action_typedef, method_defs,
            groundargs_def, available_def, execute_def)
end

function generate_action_methods(domain::Domain, state::State,
                              domain_type::Symbol, state_type::Symbol,
                              action_name::Symbol, action_type::Symbol)
    action = get_action(domain, action_name)
    get_name_def =
        :(get_name(::$action_type) = $(QuoteNode(get_name(action))))
    get_argvars_def =
        :(get_argvars(::$action_type) = $(QuoteNode(get_argvars(action))))
    get_argtypes_def =
        :(get_argtypes(::$action_type) = $(QuoteNode(get_argtypes(action))))
    get_precond_def =
        :(get_precond(::$action_type) = $(QuoteNode(get_precond(action))))
    get_effect_def =
        :(get_effect(::$action_type) = $(QuoteNode(get_effect(action))))
    method_defs = Expr(:block, get_name_def, get_argvars_def, get_argtypes_def,
                       get_precond_def, get_effect_def)
    return method_defs
end

function generate_groundargs(domain::Domain, state::State,
                             domain_type::Symbol, state_type::Symbol,
                             action_name::Symbol, action_type::Symbol)
    action = get_action(domain, action_name)
    argtypes = get_argtypes(action)
    if get_requirements(domain)[:typing]
        objs_exprs = [:(get_objects(state, $(QuoteNode(ty)))) for ty in argtypes]
    else
        objs_exprs = [:(get_objects(state)) for ty in argtypes]
    end
    iter_expr = :(Iterators.product($(objs_exprs...)))
    groundargs_def = quote
        function groundargs(domain::$domain_type, state::$state_type,
                            action::$action_type)
            return $iter_expr
        end
    end
    return groundargs_def
end


text:
@inline function add_dbfp_db(x::DoubleFloat{T}, y::T) where {T<:IEEEFloat}
    (isnan(LO(x)) | isnan(y)) && return add_dbfp_db_nonfinite(x,y)
    return DoubleFloat{T}(add_ddfp_dd(HILO(x), y))
end

@inline function sub_dbfp_db(x::DoubleFloat{T}, y::T) where {T<:IEEEFloat}
    (isnan(LO(x)) | isnan(y)) && return sub_dbfp_db_nonfinite(x,y)
    return DoubleFloat{T}(sub_ddfp_dd(HILO(x), y))
end

@inline function mul_dbfp_db(x::DoubleFloat{T}, y::T) where {T<:IEEEFloat}
    (isnan(LO(x)) | isnan(y)) && return mul_dbfp_db_nonfinite(x,y)
    return DoubleFloat{T}(mul_ddfp_dd(HILO(x), y))
end

@inline function dvi_dbfp_db(x::DoubleFloat{T}, y::T) where {T<:IEEEFloat}
    (isnan(LO(x)) | isnan(y)) && return dvi_dbfp_db_nonfinite(x,y)
    return DoubleFloat{T}(dvi_ddfp_dd(HILO(x), y))
end

@inline function add_dbfp_db_nonfinite(x::DoubleFloat{T}, y::T) where {T<:IEEEFloat}
    z = HI(x) + y
    return DoubleFloat{T}(z, T(NaN))
end

@inline function sub_dbfp_db_nonfinite(x::DoubleFloat{T}, y::T) where {T<:IEEEFloat}
    z = HI(x) - y
    return DoubleFloat{T}(z, T(NaN))
end

@inline function mul_dbfp_db_nonfinite(x::DoubleFloat{T}, y::T) where {T<:IEEEFloat}
    z = HI(x) * y
    return DoubleFloat{T}(z, T(NaN))
end

@inline function dvi_dbfp_db_nonfinite(x::DoubleFloat{T}, y::T) where {T<:IEEEFloat}
    z = div(HI(x), y)
    return DoubleFloat{T}(z, T(NaN))
end


text:
using PackageStream
using Base.Test

# Test functions
pack_unpack(a::Any) = unpack(pack(a)) == a
unpack_pack(b::Vector{UInt8}) = pack(unpack(b)) == b
function test_packing(a::Any, b::Vector{UInt8})
    return pack(a) == b && unpack(b) == a && pack_unpack(a) && unpack_pack(b)
end

# Test : Nil, Bool
@test test_packing(nothing, b"\xc0")
@test test_packing(false, b"\xc2")
@test test_packing(true, b"\xc3")

# Test : Integer
@test pack_unpack((r = rand(-2^4:typemax(Int8)); println("tiny = $r"); r))
@test pack_unpack((r = rand(typemin(Int8):-2^4-1); println("int8 = $r"); r))
@test pack_unpack((r = rand(typemax(Int8)+1:typemax(Int16)); println("int16+ = $r"); r))
@test pack_unpack((r = rand(typemin(Int16):typemin(Int8)-1); println("int16- = $r"); r))
@test pack_unpack((r = rand(typemax(Int16)+1:typemax(Int32)); println("int32+ = $r"); r))
@test pack_unpack((r = rand(typemin(Int32):typemin(Int16)-1); println("int32- = $r"); r))
@test pack_unpack((r = rand(typemax(Int32)+1:typemax(Int64)); println("int64+ = $r"); r))
@test pack_unpack((r = rand(typemin(Int64):typemin(Int32)-1); println("int64- = $r"); r))

# Test : String
strings = [
    str_tiny = "tiny",
    str_8 = "My string is not tiny but within 8 bytes",
    str_16 = "My string is long enough, not tiny, above 8 bytes and below 16. "^4,
    str_32 = "My string is long enough, not tiny, above 16 bytes and below 32. "^1024
]
@test mapreduce(pack_unpack, &, true, strings)
alphabet_bytes = b"\xd0\x1a\x41\x42\x43\x44\x45\x46\x47\x48\x49\x4a\x4b\x4c\x4d\x4e\x4f\x50\x51\x52\x53\x54\x55\x56\x57\x58\x59\x5a"
@test test_packing("ABCDEFGHIJKLMNOPQRSTUVWXYZ", alphabet_bytes)

# Test : List (Vector)
@test pack_unpack(Vector(1:10))
@test pack_unpack(Vector(1:100))
@test pack_unpack(Vector(1:1000))
@test pack_unpack(Vector(1:100000))
heterogeneous_list_bytes =
    b"\x93\x01\xC1\x40\x00\x00\x00\x00\x00\x00\x00\x85\x74\x68\x72\x65\x65"
@test test_packing([1, 2., "three"], heterogeneous_list_bytes)

# Test : Map (Dict)
@test pack_unpack(Dict([(i, 2i) for i in 1:10]))
@test pack_unpack(Dict([(i, 2i) for i in 1:100]))
@test pack_unpack(Dict([(i, 2i) for i in 1:1000]))
@test pack_unpack(Dict([(i, 2i) for i in 1:100000]))

# Test : Structure (Tuple)
@test pack_unpack(Tuple(1:10))
@test pack_unpack(Tuple(1:100))
@test pack_unpack(Tuple(1:1000))


text:
# This file was a part of Julia. License is MIT: http://julialang.org/license

import Base: AbstractSet, similar, copy, copy!, eltype, push!, pop!, delete!,
             empty!, isempty, union, union!, intersect, intersect!,
             setdiff, setdiff!, symdiff, symdiff!, in, iterate,
             last, length, show, display, hash, issubset, ==, <=, <, +, *, -, !, unsafe_getindex,
             unsafe_setindex!, findnextnot, first, getindex, rand, xor
using Random

if !isdefined(Base, :complement)
    export complement, complement!
else
    import Base: complement, complement!
end

struct CardSet32 <: AbstractSet{UInt32}
    cs::UInt32
end

#IDEA offset as type parameter -> CardSet32{:offset}, CardSet32{:5} -> get it with eltype

const Card = CardSet32 #a card is a card set with one element

CardSet32() = CardSet32(UInt32(0))

function CardSet32(cards::Array{Card, 1})
    result = CardSet32()
    for card in cards
        result = union(result, card)
    end
    return result
end

const emptyCardSet32 = CardSet32()

function <(cs1::Card, cs2::Card)
    cs1.cs < cs2.cs
end

eltype(cs1::CardSet32) = UInt32
eltype(::Type{CardSet32}) = UInt32

union(cs1::CardSet32, cs2::CardSet32) = CardSet32(cs1.cs | cs2.cs)
union(cs1::CardSet32, css...) = union(cs1, union(css...))
+(cs1::CardSet32, cs2::CardSet32) = union(cs1, cs2)

intersect(cs1::CardSet32, cs2::CardSet32) = CardSet32(cs1.cs & cs2.cs)
intersect(cs1::CardSet32, css...) = intersect(cs1, intersect(css...))
*(cs1::CardSet32, cs2::CardSet32) = intersect(cs1, cs2)

in(cs1::CardSet32, cs2::CardSet32) = isequal(cs1.cs, cs1.cs & cs2.cs)
hash(cs1::CardSet32) = cs1.cs

function in(cardIndex::Int64, cs2::CardSet32)
    if 0 <= cardIndex < 32
        return in(CardSet32(0x00000001 << cardIndex), cs2)
    end
    throw(BoundsError("needed: 0 <= card < 32"))
end

length(cs1::CardSet32) = count_ones(cs1.cs)

function toArray(cs1::CardSet32)
    result = Vector{Card}()
    for card in cs1 push!(result, card) end
    return result
end

complement(cs1::CardSet32) = CardSet32(~cs1.cs)
!(cs1::CardSet32) = complement(cs1)

setdiff(cs1::CardSet32, cs2::CardSet32) = intersect(cs1, complement(cs2))
-(cs1::CardSet32, cs2::CardSet32) = setdiff(cs1, cs2)

symdiff(cs1::CardSet32, cs2::CardSet32) = CardSet32(xor(cs1.cs, cs2.cs)) #xor / flip
xor(cs1::CardSet32, cs2::CardSet32) = symdiff(cs1, cs2)

#add card(s)
function add(cs1::CardSet32, items::CardSet32)
    if !UNSAFE && !isempty(intersect(items, cs1)) throw(ArgumentError("items already in cs1, cannot remove")) end
    union(cs1, items)
end

#remove card(s)
function remove(cs1::CardSet32, items::CardSet32)
    if !UNSAFE && !in(items, cs1) throw(ArgumentError("items not in cs1, cannot remove")) end
    setdiff(cs1, items)
end

isempty(cs1::CardSet32) = isequal(cs1.cs, UInt32(0))

const CardSet32_iterstate = Tuple{UInt32, UInt32}

function iterate(cs1::CardSet32)
    return iterate(cs1, (cs1.cs, 0x80000000)) #state is the set and an index
end

function iterate(cs1::CardSet32, state::CardSet32_iterstate)
    isequal(state[1], UInt32(0)) && return nothing
    leadzeros = leading_zeros(state[1])
    nextValue = state[2] >>> leadzeros
    return (CardSet32(nextValue), (state[1] << (leadzeros + 1), nextValue >>> 1))
end

function first(cs1::CardSet32)
    isempty(cs1) && throw(ArgumentError("collection must be non-empty"))
    CardSet32(0x80000000 >>> leading_zeros(cs1.cs))
end

function last(cs1::CardSet32)
    isempty(cs1) && throw(ArgumentError("collection must be non-empty"))
    CardSet32(0x00000001 << trailing_zeros(cs1.cs))
end

function getindex(cs1::CardSet32, index::Int)
    if  index <= 0 || index > length(cs1) throw(BoundsError("index out of bounds")) end
    state = start(cs1)
    val = CardSet32()
    for i in 1:index
        (val, state) = next(cs1, state)
    end
    return val
end

#pick a random card
function rand(rng::AbstractRNG, cs1::CardSet32)
    cs1[rand(rng, 1:length(cs1))]
    # rand(rng, toArray(cs1)) #TODO more efficiency: make set indexible by getindex
end

max_stars_repo_name: albapa/ulti

max_issues_repo_name: albapa/ulti

max_forks_repo_name: albapa/ulti


text:
# A tuple of (Pix,Spin,Basis) defines the generic behavior of our fields
abstract type Pix end
abstract type Spin end
abstract type Basis end

# All fields are a subtype of this.
abstract type Field{B<:Basis, S<:Spin, P<:Pix, T} <: AbstractVector{T} end

# Spin types, "S0" is spin-0, i.e. a scalar map. "S2" is spin-2 like QU, and S02
# is a tuple of S0 and S2 like TQU.
abstract type S0 <: Spin end
abstract type S2 <: Spin end
abstract type S02 <: Spin end

# the basis type for a FieldTuple can be a BasisTuple which just holds the bases
# of the sub-fields
abstract type BasisTuple{T} <: Basis end
promote_type(::Type{BasisTuple{BT1}}, ::Type{BasisTuple{BT2}}) where {BT1,BT2} = BasisTuple{Tuple{map_tupleargs(promote_type,BT1,BT2)...}}

# Basis types
abstract type Map       <: Basis end
abstract type Fourier   <: Basis end
abstract type QUMap     <: Basis end
abstract type EBMap     <: Basis end
abstract type QUFourier <: Basis end
abstract type EBFourier <: Basis end
const IQUFourier = BasisTuple{Tuple{Fourier,QUFourier}}
const IEBFourier = BasisTuple{Tuple{Fourier,EBFourier}}
const IQUMap     = BasisTuple{Tuple{Map,QUMap}}
const IEBMap     = BasisTuple{Tuple{Map,EBMap}}

basis_promotion_rules = Dict(
    # S0
    (Map,       Fourier)    => Map,
    # S2
    (QUMap,     QUFourier)  => QUMap,
    (EBMap,     EBFourier)  => EBFourier,
    (QUMap,     EBMap)      => QUMap,
    (QUFourier, EBFourier)  => QUFourier,
    (QUMap,     EBFourier)  => QUMap,
    (QUFourier, EBMap)      => QUFourier
)
for ((B1,B2),B′) in basis_promotion_rules
    @eval promote_rule(::Type{$B1}, ::Type{$B2}) = $B′
end

# A "basis-like" object, e.g. the lensing basis Ł or derivative basis Ð. For any
# particular types of fields, these might be different actual bases, e.g. the
# lensing basis is Map for S0 but QUMap for S2.
abstract type Basislike <: Basis end
(::Type{B})(f::F) where {F<:Field,B<:Basislike} = B(F)(f)
(::Type{B})(f′::Field, f::F) where {F<:Field,B<:Basislike} = B(F)(f′,f)
(::Type{B})(a::AbstractArray{<:Field}...) where {B<:Basis} = B.(a...)

# Basis in which derivatives are sparse for a given field
abstract type DerivBasis <: Basislike end
const Ð = DerivBasis
Ð!(args...) = Ð(args...)

# Basis in which lensing is a pixel remapping for a given field
abstract type LenseBasis <: Basislike end
const Ł = LenseBasis
Ł!(args...) = Ł(args...)

# B(f) where B is a basis converts f to that basis. This is the fallback if the
# field is already in the right basis.
(::Type{B})(f::Field{B}) where {B} = f

# The abstract `Basis` type means "any basis", hence this conversion rule:
Basis(f::Field) = f

# B(f′, f) converts f to basis B and stores the result inplace in f′. If f is
# already in basis B, we just return f (but note, we never actually set f′ in
# this case, which is more efficient, but necessitates some care when using this
# construct)
(::Type{B})(f′::Field{B}, f::Field{B}) where {B} = f

### ImplicitOp
# ImplicitOps and ImplicitFields are operators and fields which may not be
# explicitly stored, but can be multiplied, etc... and generally implicitly take
# on whatever size they need. The specification of Float32 as the the Array type
# is only for inference, such that its widened to whatever its needed to be the
# end.
abstract type ImplicitOp{B<:Basis, S<:Spin, P<:Pix} <: AbstractMatrix{Float32} end
abstract type ImplicitField{B<:Basis, S<:Spin, P<:Pix} <: Field{B,S,P,Float32} end

# no size and 0-length represents the fact that ImplicitOps aquire the size of
# the fields they're applied to. it also helps makes the generic printing
# methods for AbstractMatrix work well.
size(::Union{ImplicitOp,ImplicitField}) = ()
length(::Union{ImplicitOp,ImplicitField}) = 0

adapt_structure(to, x::Union{ImplicitOp,ImplicitField}) = x

# printing
show(io::IO, ::MIME"text/plain", L::ImplicitOp) = show(io,L)
show(io::IO, ::MIME"text/plain", L::Adjoint{<:Any,<:ImplicitOp}) = show(io,L)
show(io::IO, L::Adjoint{<:Any,<:ImplicitOp}) = (print(io,"Adjoint{"); show(io,parent(L)); print(io,"}"))
# this is the main function ImplicitOps should specialize if this default behavior isn't enough:
show(io::IO, L::ImplicitOp) = showarg(io, L, true)

# all CMBLensing operators are then either Diagonals or ImplicitOps
# ImplicitOrAdjOp are things for which algebra is done lazily, including Diagonal{<:ImplicitField}
const DiagOp{F<:Field,T} = Diagonal{T,F}
const LinOp{B,S,P} = Union{ImplicitOp{B,S,P},DiagOp{<:Field{B,S,P}}}
const ImplicitOrAdjOp{B,S,P} = Union{ImplicitOp{B,S,P}, Adjoint{<:Any,<:ImplicitOp{B,S,P}}, DiagOp{<:ImplicitField{B,S,P}}}
const LinOrAdjOp{B,S,P} = Union{ImplicitOrAdjOp{B,S,P},DiagOp{<:Field{B,S,P}}}

### Scalars
# scalars which are allowed in our expressions must be real because we
# implicitly assume our maps are real, and addition/multiplication by a complex
# number, even of the fourier transform, would break this.
const Scalar = Real
const FieldOrOp = Union{Field,LinOp}
const FieldOpScal = Union{Field,LinOp,Scalar}

"""
    logdet(L::LinOp, θ)

If L depends on θ, evaluates `logdet(L(θ))` offset by its fiducial value at
`L()`. Otherwise, returns 0.
"""
logdet(L::LinOp, θ) = depends_on(L,θ) ? logdet(L()\L(θ)) : 0

#
# ### Matrix conversion
#
# # We can build explicit matrix representations of linear operators by applying
# # them to a set of basis vectors which form a complete basis (in practice this
# # is prohibitive except for fairly small maps, e.g. 16x16 pixels, but is quite
# # useful)
#
# """
#     full(::Type{Fin}, ::Type{Fout}, L::LinOp; progress=true)
#     full(::Type{F}, L::LinOp; progress=true)
#
# Construct an explicit matrix representation of the linear operator `L` by
# applying it to a set of vectors which form a complete basis. The `Fin` and
# `Fout` types should be fields which specify the input and output bases for the
# representation (or just `F` if `L` is square and we want the same input/output
# bases)
#
# The name `full` is to be consistent with Julia's SparseArrays where `full`
# builds a full matrix from a sparse one.
# """
# full(::Type{Fin}, ::Type{Fout}, L::LinOp; progress=true) where {Fin<:Field, Fout<:Field} =
#     hcat(@showprogress(progress ? 1 : Inf, [(Fout(L*(x=zeros(length(Fin)); x[i]=1; x)[Tuple{Fin}]))[:] for i=1:length(Fin)])...)
# full(::Type{F}, L::LinOp; kwargs...) where {F<:Field} = full(F,F,L; kwargs...)
#
#

### Other generic stuff

# convenience "getter" functions for the Basis/Spin/Pix
# basis of UnionAlls like basis(Field) will return Basis (which means any Basis)
basis(::Type{<:Field}) = Basis
basis(::Type{<:Field{B,S,P}}) where {B<:Basis,S<:Spin,P<:Pix} = B
basis(::F) where {F<:Field} = basis(F)
spin(::Type{<:Field}) = Spin
spin(::Type{<:Field{B,S,P}}) where {B<:Basis,S<:Spin,P<:Pix} = S
spin(::F) where {F<:Field} = spin(F)
pix(::Type{<:Field}) = Pix
pix(::Type{<:Field{B,S,P}}) where {B<:Basis,S<:Spin,P<:Pix} = P
pix(::F) where {F<:Field} = pix(F)

# we use Field cat'ing mainly for plotting, e.g. plot([f f; f f]) plots a 2×2
# matrix of maps. the following definitions make it so that Fields aren't
# splatted into a giant matrix when doing [f f; f f] (which they would othewise
# be since they're Arrays)
hvcat(rows::Tuple{Vararg{Int}}, values::Field...) = hvcat(rows, ([x] for x in values)...)
hcat(values::Field...) = hcat(([x] for x in values)...)

### printing
print_array(io::IO, f::Field) = print_array(io, f[:])
show_vector(io::IO, f::Field) = show_vector(io, f[:])
show_vector(io::IO, f::ImplicitField) = print(io, "[…]")

# addition/subtraction works between any fields and scalars, promotion is done
# automatically if fields are in different bases
for op in (:+,:-), (T1,T2) in ((:Field,:Scalar),(:Scalar,:Field),(:Field,:Field))
    @eval ($op)(a::$T1, b::$T2) = broadcast($op,($T1==$T2 ? promote : tuple)(a,b)...)
end

# multiplication/division is not strictly defined for abstract vectors, but
# make it work anyway if the two fields are exactly the same type, in which case
# its clear we wanted broadcasted multiplication/division.
for op in (:*, :/)
    @eval ($op)(a::Field{B,S,P}, b::Field{B,S,P}) where {B,S,P} = broadcast($op, a, b)
end

(::Type{T})(f::Field{<:Any,<:Any,<:Any,<:Real}) where {T<:Real} = T.(f)
(::Type{T})(f::Field{<:Any,<:Any,<:Any,<:Complex}) where {T<:Real} = Complex{T}.(f)

one(f::Field) = fill!(similar(f), one(eltype(f)))

invalid_broadcast_error(B1,B2) =
    error("""Can't broadcast fields in different bases. ($B1, $B2)
    Try the same operation without broadcasting (which will do an automatic basis conversion).""")

@init @require PyCall="438e738f-606a-5dbb-bf0a-cddfbfd45ab0" begin
    # never try to auto-convert Fields or LinOps to Python arrays
    PyCall.PyObject(x::Union{LinOp,Field}) = PyCall.pyjlwrap_new(x)
end