text:
Vk.@bitmask_flag ResourceType::UInt32 begin
    RESOURCE_TYPE_VERTEX_BUFFER = 1
    RESOURCE_TYPE_INDEX_BUFFER = 2
    RESOURCE_TYPE_COLOR_ATTACHMENT = 4
    RESOURCE_TYPE_DEPTH_ATTACHMENT = 8
    RESOURCE_TYPE_STENCIL_ATTACHMENT = 16
    RESOURCE_TYPE_INPUT_ATTACHMENT = 32
    RESOURCE_TYPE_TEXTURE = 64
    RESOURCE_TYPE_BUFFER = 128
    RESOURCE_TYPE_IMAGE = 256
    RESOURCE_TYPE_DYNAMIC = 512
    RESOURCE_TYPE_STORAGE = 1024
    RESOURCE_TYPE_TEXEL = 2048
    RESOURCE_TYPE_UNIFORM = 4096
    RESOURCE_TYPE_SAMPLER = 8192
end

@enum ResourceClass::Int8 begin
    RESOURCE_CLASS_BUFFER = 1
    RESOURCE_CLASS_IMAGE = 2
    RESOURCE_CLASS_ATTACHMENT = 3
end

abstract type ResourceInfo end

struct BufferResourceInfo <: ResourceInfo
    size::Int
    usage::Vk.BufferUsageFlag
end

BufferResourceInfo(size) = BufferResourceInfo(size, Vk.BufferUsageFlag(0))

struct ImageResourceInfo <: ResourceInfo
    dims::Vector{Int}
    size_unit::SizeUnit
    format::Vk.Format
    usage::Vk.ImageUsageFlag
    aspect::Vk.ImageAspectFlag
end

function ImageResourceInfo(format; dims = [1., 1.], size_unit = SIZE_SWAPCHAIN_RELATIVE, usage = Vk.ImageUsageFlag(0), aspect = Vk.ImageAspectFlag(0))
    ImageResourceInfo(dims, size_unit, format, usage, aspect)
end

struct AttachmentResourceInfo <: ResourceInfo
    image_info::ImageResourceInfo
    AttachmentResourceInfo(image_info::ImageResourceInfo) = new(image_info)
end

AttachmentResourceInfo(args...; kwargs...) = AttachmentResourceInfo(ImageResourceInfo(args...; kwargs...))

struct Pass
    name::Symbol
    stages::Vk.PipelineStageFlag
end

Pass(stages) = Pass(gensym("Pass"), stages)

struct ResourceUsage
    type::ResourceType
    access::MemoryAccess
end

struct SynchronizationRequirements
    stages::Vk.PipelineStageFlag
    access::Vk.AccessFlag
    wait_semaphores::Vector{Vk.Semaphore}
end

SynchronizationRequirements(stages, access) = SynchronizationRequirements(stages, access, [])

"""
Frame graph implementation.

A frame graph has a list of virtual resources (buffers, images, attachments) that are
referenced by passes. They are turned into physical resources for the actual execution of those passes.

The frame graph uses two graph structures: a resource graph and an execution graph.

## Resource graph (bipartite, directed)

This bipartite graph has two types of vertices: passes and resources.
An edge from a resource to a pass describes a read dependency. An edge from a pass to a resource describes a write dependency.

Graph attributes:
- Resources:
    - `:name`: name of the resource
    - `:format` (if the resource describes either an image or an attachment)
    - `:usage`
    - `:size` (if the resource describes a buffer)
    - `:vresource`: description as a virtual resource
    - `:presource`: physical resource
    - `:current_layout` (if the resource describes an image)
    - `:last_write`: access bits and pipeline stages associated to the latest pass that wrote to the resource
    - `:synchronization_state`: access bits and pipeline stages of all passes that required to synchronize with this resource
        since the last write.
- Passes:
    - `:name`: name of the pass
- Edge between a resource and a pass (all directions)
    - `:image_layout` (if the resource describes an image)
    - `:usage`
    - `:aspect`
    - `:stage`: stages in which the resource is used
    - `:clear_value`: clear value, if applicable. If this property is absent and a load operation takes place on this resource in this pass,
        then the contents are preserved.

## Execution graph (directed, acyclic)

In this graph, vertices represent passes, and edges are resource dependencies between passes.
A topological sort of this graph represents a possible sequential execution order that respects execution dependencies.

Reusing the example above, the graph has three vertices: `gbuffer`, `lighting` and `adapt_luminance`.
`gbuffer` has five outgoing edges to `lighting`, each edge being labeled with a resource.
`lighting` has one outgoing edge to `adapt_luminance`.

This graph is generated just-in-time, to convert the resource graph into a linear sequence of passes.
"""
mutable struct FrameGraph
    device::Device
    resource_graph::MetaGraph{Int}
    resources::Dictionary{Symbol,Int}
    passes::Dictionary{Symbol,Int}
    frame::Frame
end

device(fg::FrameGraph) = fg.device

attribute(rg::AbstractMetaGraph, i, key) = get_prop(rg, i, key)
pass_attribute(rg::AbstractMetaGraph, i::Integer, key::Symbol) = attribute(rg, i, key)
pass_attribute(fg::FrameGraph, name::Symbol, key::Symbol) = attribute(fg.resource_graph, fg.passes[name], key)
resource_attribute(rg::AbstractMetaGraph, i::Integer, key::Symbol) = attribute(rg, i, key)
resource_attribute(fg::FrameGraph, name::Symbol, key::Symbol) = attribute(fg.resource_graph, fg.resources[name], key)
function set_attribute(fg::FrameGraph, pass::Symbol, resource::Symbol, key::Symbol, value)
    set_prop!(fg.resource_graph, fg.passes[pass], fg.resources[resource], key, value)
end

virtual_image(g, idx) = resource_attribute(g, idx, :virtual_resource)::ImageResourceInfo
virtual_buffer(g, idx) = resource_attribute(g, idx, :virtual_resource)::BufferResourceInfo
virtual_attachment(g, idx) = resource_attribute(g, idx, :virtual_resource)::AttachmentResourceInfo

image(g, idx) = resource_attribute(g, idx, :physical_resource)::Image
buffer(g::Union{FrameGraph, MetaGraph}, idx) = resource_attribute(g, idx, :physical_resource)::Buffer
attachment(g, idx) = resource_attribute(g, idx, :physical_resource)::Attachment

resource_type(g, i, j) = attribute(g, Edge(i, j), :type)::ResourceType
access(g, i, j) = attribute(g, Edge(i, j), :access)::MemoryAccess
access_bits(g, i, j) = attribute(g, Edge(i, j), :access_bits)::Vk.AccessFlag
image_layout(g, i, j) = attribute(g, Edge(i, j), :image_layout)::Vk.ImageLayout
clear_value(g, i, j) = attribute(g, Edge(i, j), :clear_value)::Vk.ClearValue
aspect(g, i, j) = attribute(g, Edge(i, j), :aspect)::Vk.ImageAspectFlag
# currently unused (not exposed to the user)
stages(g, i, j) = attribute(g, Edge(i, j), :stages)::Vk.PipelineStageFlag

resource_class(g, idx) = resource_attribute(g, idx, :class)::ResourceClass
current_layout(g, idx) = resource_attribute(g, idx, :current_layout)::Vk.ImageLayout
format(g, idx) = resource_attribute(g, idx, :format)::Vk.Format
buffer_usage(g, idx) = resource_attribute(g, idx, :usage)::Vk.BufferUsageFlag
image_usage(g, idx) = resource_attribute(g, idx, :usage)::Vk.ImageUsageFlag
_size(g, idx) = resource_attribute(g, idx, :size)::Int
_dims(g, idx) = resource_attribute(g, idx, :dims)::Vector{Int}
name(g, idx) = attribute(g, idx, :name)::Symbol
last_write(g, idx) = resource_attribute(g, idx, :last_write)::Pair{Vk.AccessFlag, Vk.PipelineStageFlag}
synchronization_state(g, idx) = resource_attribute(g, idx, :synchronization_state)::Dictionary{Vk.AccessFlag,Vk.PipelineStageFlag}
access(g, idx) = resource_attribute(g, idx, :access)::MemoryAccess
stages(g, idx) = pass_attribute(g, idx, :stages)::Vk.PipelineStageFlag
render_function(g, idx) = pass_attribute(g, idx, :render_function)
render_pass(g, idx) = pass_attribute(g, idx, :pass)::RenderPass

function FrameGraph(device::Device, frame::Frame = Frame(device))
    FrameGraph(device, MetaGraph(), Dictionary(), Dictionary(), frame)
end

function add_pass!(render_function, fg::FrameGraph, name::Symbol, pass::RenderPass, stages::Vk.PipelineStageFlag = Vk.PIPELINE_STAGE_ALL_GRAPHICS_BIT)
    !haskey(fg.passes, name) || error("Pass '$name' was already added to the frame graph. Passes can only be provided once.")
    g = fg.resource_graph
    add_vertex!(g)
    v = nv(g)
    set_prop!(g, v, :name, name)
    set_prop!(g, v, :stages, stages)
    set_prop!(g, v, :render_function, render_function)
    set_prop!(g, v, :pass, pass)
    insert!(fg.passes, name, v)
    nothing
end

add_pass!(fg::FrameGraph, pass::Pass) = add_pass!(fg, pass.name, pass.stages)

function add_resource!(fg::FrameGraph, name::Symbol, info::ResourceInfo)
    !haskey(fg.resources, name) || error("Resource '$name' was already added to the frame graph. Resources can only be provided once.")
    g = fg.resource_graph
    add_vertex!(g)
    v = nv(g)
    set_prop!(g, v, :name, name)
    set_prop!(g, v, :access, MemoryAccess(0))
    insert!(fg.resources, name, v)
    if info isa BufferResourceInfo
        set_prop!(g, v, :size, info.size)
        set_prop!(g, v, :usage, info.usage)
        set_prop!(g, v, :class, RESOURCE_CLASS_BUFFER)
    else
        image_info = if info isa AttachmentResourceInfo
            set_prop!(g, v, :class, RESOURCE_CLASS_ATTACHMENT)
            info.image_info
        else
            set_prop!(g, v, :class, RESOURCE_CLASS_IMAGE)
            info
        end
        set_prop!(g, v, :dims, image_info.dims)
        set_prop!(g, v, :usage, image_info.usage)
        set_prop!(g, v, :aspect, image_info.aspect)
        set_prop!(g, v, :format, image_info.format)
    end
    if haskey(fg.frame.resources, name)
        set_prop!(g, v, :physical_resource, fg.frame.resources[name].data)
        if resource_class(g, v) in (RESOURCE_CLASS_IMAGE, RESOURCE_CLASS_ATTACHMENT)
            set_prop!(g, v, :current_layout, image_layout(fg.frame.resources[name].data))
        end
    end
    nothing
end

function add_resource_usage!(fg::FrameGraph, iter)
    for (pass, usage) in pairs(iter)
        add_resource_usage!(fg, pass, usage)
    end
end

function add_resource_usage!(fg::FrameGraph, pass::Symbol, usage::Dictionary{Symbol,ResourceUsage})
    haskey(fg.passes, pass) || error("Unknown pass '$(pass)'")
    v = fg.passes[pass]
    g = fg.resource_graph
    for (name, usage) in pairs(usage)
        haskey(fg.resources, name) || error("Unknown resource '$name'")
        i = fg.resources[name]
        add_edge!(g, i, v) || error("A resource usage between pass '$pass' and resource '$name' was already added to the frame graph. Resource usages can only be provided once.")
        set_prop!(g, i, v, :type, usage.type)
        set_prop!(g, i, v, :access, usage.access)
    end
end

function resource_usages(ex::Expr)
    lines = @match ex begin
        Expr(:block, _...) => ex.args
        _ => [ex]
    end

    filter!(Base.Fix2(!isa, LineNumberNode), lines)

    usages = Dictionary{Symbol,Dictionary{Symbol,ResourceUsage}}()
    for line in lines
        (f, reads, writes) = @match line begin
            :($writes = $f($(reads...))) => @match writes begin
                Expr(:tuple, _...) => (f, reads, writes.args)
                ::Expr => (f, reads, [writes])
            end
            _ => error("Malformed expression, expected :(a, b = f(c, d)), got $line")
        end

        dict = Dictionary{Symbol,ResourceUsage}()
        for (name, type) in extract_resource_spec.(reads)
            !haskey(dict, name) || error("Resource $name for pass $f specified multiple times in read access")
            insert!(dict, name, ResourceUsage(type, READ))
        end
        for (name, type) in extract_resource_spec.(writes)
            usage = if haskey(dict, name)
                WRITE ∉ dict[name].access || error("Resource $name for pass $f specified multiple times in write access")
                ResourceUsage(type | dict[name].type, WRITE | READ)
            else
                ResourceUsage(type, WRITE)
            end
            set!(dict, name, usage)
        end
        !haskey(usages, f) || error("Pass $f is specified more than once")
        insert!(usages, f, dict)
    end
    usages
end

macro resource_usages(ex)
    :($(resource_usages(ex)))
end

function extract_resource_spec(ex::Expr)
    @match ex begin
        :($r::Buffer::Vertex) => (r => RESOURCE_TYPE_VERTEX_BUFFER)
        :($r::Buffer::Index) => (r => RESOURCE_TYPE_INDEX_BUFFER)
        :($r::Buffer::Storage) => (r => RESOURCE_TYPE_BUFFER | RESOURCE_TYPE_STORAGE)
        :($r::Buffer::Uniform) => (r => RESOURCE_TYPE_BUFFER | RESOURCE_TYPE_UNIFORM)
        :($r::Buffer) => (r => RESOURCE_TYPE_BUFFER)
        :($r::Color) => (r => RESOURCE_TYPE_COLOR_ATTACHMENT)
        :($r::Depth) => (r => RESOURCE_TYPE_DEPTH_ATTACHMENT)
        :($r::Stencil) => (r => RESOURCE_TYPE_STENCIL_ATTACHMENT)
        :($r::Depth::Stencil) || :($_::Stencil::Depth) => (r => RESOURCE_TYPE_DEPTH_ATTACHMENT | RESOURCE_TYPE_STENCIL_ATTACHMENT)
        :($r::Texture) => (r => RESOURCE_TYPE_TEXTURE)
        :($r::Image::Storage) => (r => RESOURCE_TYPE_IMAGE | RESOURCE_TYPE_STORAGE)
        :($r::Input) => (r => RESOURCE_TYPE_INPUT_ATTACHMENT)
        ::Symbol => error("Resource type annotation required for $ex")
        _ => error("Invalid or unsupported resource type annotation for $ex")
    end
end

function extract_resource_name(ex::Expr)
    @match ex begin
        :($r::$_::$_) => r
        :($r::$_) => r
        _ => error("Cannot extract resource name from $ex")
    end
end

function clear_attachments(fg::FrameGraph, pass::Symbol, color_clears, depth_clears = [], stencil_clears = [])
    clears = Dictionary{Symbol,Vk.ClearValue}()
    for (resource, color) in color_clears
        clear_value = Vk.ClearValue(Vk.ClearColorValue(convert(NTuple{4,Float32}, color)))
        insert!(clears, resource, clear_value)
    end
    for (resource, depth) in depth_clears
        clear_value = Vk.ClearValue(Vk.ClearDepthStencilValue(depth, 0))
        insert!(clears, resource, clear_value)
    end
    for (resource, stencil) in stencil_clears
        clear_value = Vk.ClearValue(Vk.ClearDepthStencilValue(0., stencil))
        insert!(clears, resource, clear_value)
    end
    for (resource, clear_value) in pairs(clears)
        set_attribute(fg, pass, resource, :clear_value, clear_value)
    end
end

function execution_graph(fg::FrameGraph)
    g = fg.resource_graph
    eg = SimpleDiGraph(length(fg.passes))
    for pass in fg.passes
        resources = neighbors(g, pass)
        for resource in resources
            WRITE in access(g, pass, resource) || continue
            for dependent_pass in neighbors(g, resource); dependent_pass == pass && continue
                add_edge!(eg, dependent_pass, pass)
            end
        end
    end
    eg
end

"""
Submit rendering commands to a device.

A command buffer is recorded, which may be split into multiple ones to take advantage of multithreading,
and is then submitted them to the provided device. A custom primary command buffer can be optionally passed as a keyword,
mostly intended for debugging purposes.

A semaphore to wait for can be provided to synchronize with other commands.
"""
function render(fg::FrameGraph; semaphore = nothing, command_buffer = nothing, submit = true)
    prepare!(fg)
    passes = sort_passes(fg)
    for pass in passes
        # create render pass & framebuffer objects
        prepare_render_pass(fg, pass)
    end
    records, pipeline_hashes = record_commands!(fg, passes)
    !isempty(pipeline_hashes) || error("No draw calls detected.")
    device = Lava.device(fg)
    create_pipelines!(device)

    # fill command buffer with synchronization commands & recorded commands
    if isnothing(command_buffer)
        command_buffer = request_command_buffer(device, Vk.QUEUE_GRAPHICS_BIT)
    end
    first_pipeline = device.pipeline_ht[first(values(pipeline_hashes))]
    initialize(command_buffer, device, fg.frame.gd, first_pipeline)
    flush(command_buffer, fg, passes, records, pipeline_hashes)
    submit || return command_buffer

    # submit rendering work
    wait_semaphores = device.transfer_ops
    !isnothing(semaphore) && push!(wait_semaphores, semaphore)
    submit_info = Vk.SubmitInfo2KHR(wait_semaphores, [Vk.CommandBufferSubmitInfoKHR(command_buffer)], [])
    Lava.submit(device, command_buffer.queue_family_index, [submit_info]; signal_fence = true, release_after_completion = [Ref(fg)])
end

function prepare!(fg::FrameGraph)
    resolve_attributes!(fg)
    create_physical_resources!(fg)
end

function record_commands!(fg::FrameGraph, passes)
    records = CompactRecord[]
    pipeline_hashes = Dictionary{ProgramInstance,UInt64}()
    g = fg.resource_graph

    # record commands and submit pipelines for creation
    for pass in passes
        record = CompactRecord(fg, pass)
        f = render_function(g, pass)
        f(record)
        push!(records, record)
        merge!(pipeline_hashes, submit_pipelines!(device(fg), render_pass(g, pass), record))
    end

    records, pipeline_hashes
end

function Base.flush(cb::CommandBuffer, fg::FrameGraph, passes, records, pipeline_hashes)
    binding_state = BindState()
    for (pass, record) in zip(passes, records)
        synchronize_before(cb, fg, pass)
        begin_render_pass(cb, fg, pass)
        binding_state = flush(cb, record, device(fg), binding_state, pipeline_hashes)
        Vk.cmd_end_render_pass_2(cb, Vk.SubpassEndInfo())
        synchronize_after(cb, fg, pass)
    end
end

function sort_passes(fg::FrameGraph)
    indices = topological_sort_by_dfs(execution_graph(fg))
    collect(values(fg.passes))[indices]
end

"""
Build barriers for all resources that require it.

Requires the extension `VK_KHR_synchronization2`.
"""
function synchronize_before(cb, fg::FrameGraph, pass::Integer)
    g = fg.resource_graph
    deps = Vk.DependencyInfoKHR([], [], [])
    for resource in neighbors(g, pass)
        type = resource_type(g, pass, resource)
        class = resource_class(g, resource)
        (req_access, req_stages) = access(g, pass, resource) => stages(g, pass)
        # if the resource was not written to recently, no synchronization is required
        if has_prop(g, resource, :last_write)
            (r_access, p_stages) = last_write(g, resource)
            req_access_bits = access_bits(type, req_access, req_stages)
            sync_state = synchronization_state(g, resource)
            synced_stages = Vk.PipelineStageFlag(0)
            barrier_needed = true
            for (access, stages) in pairs(sync_state)
                if covers(access, req_access_bits)
                    synced_stages |= (stages & req_stages)
                end
                if synced_stages == req_stages
                    barrier_needed = false
                    break
                else
                    # keep only stages that haven't been synced
                    req_stages &= ~synced_stages
                end
            end
            if barrier_needed
                @switch class begin
                    @case &RESOURCE_CLASS_BUFFER
                        buff = buffer(g, resource)
                        barrier = Vk.BufferMemoryBarrier2KHR(p_stages, req_stages, r_access, req_access_bits, 0, 0, buff, offset(buff), size(buff))
                        push!(deps.buffer_memory_barriers, barrier)
                    @case &RESOURCE_CLASS_IMAGE || &RESOURCE_CLASS_ATTACHMENT
                        view = if class == RESOURCE_CLASS_IMAGE
                            View(image(g, resource))
                        else
                            attachment(g, resource).view
                        end
                        new_layout = image_layout(g, pass, resource)
                        range = subresource_range(view)
                        barrier = Vk.ImageMemoryBarrier2KHR(C_NULL, p_stages, req_stages, r_access, req_access_bits, current_layout(g, resource), new_layout, 0, 0, view.image, range)
                        push!(deps.image_memory_barriers, barrier)
                        set_prop!(g, resource, :current_layout, new_layout)
                        view.image.layout[] = new_layout
                end
                set!(sync_state, req_access_bits, req_stages)
            end
        elseif class in (RESOURCE_CLASS_IMAGE, RESOURCE_CLASS_ATTACHMENT) && current_layout(g, resource) ≠ image_layout(g, pass, resource)
            # perform the required layout transition without further synchronization
            view = if class == RESOURCE_CLASS_IMAGE
                View(image(g, resource))
            else
                attachment(g, resource).view
            end
            new_layout = image_layout(g, pass, resource)
            barrier = Vk.ImageMemoryBarrier2KHR(C_NULL, Vk.PipelineStageFlag(0), Vk.PipelineStageFlag(0), Vk.AccessFlag(0), Vk.AccessFlag(0), current_layout(g, resource), new_layout, 0, 0, view.image, subresource_range(view))
            push!(deps.image_memory_barriers, barrier)
            set_prop!(g, resource, :current_layout, new_layout)
            view.image.layout[] = new_layout
        end
        if WRITE in req_access
            set_prop!(g, resource, :last_write, access_bits(resource_type(g, pass, resource), req_access, req_stages) => req_stages)
            set_prop!(g, resource, :synchronization_state, Dictionary{Vk.AccessFlag,Vk.PipelineStageFlag}())
        end
    end
    (!isempty(deps.memory_barriers) || !isempty(deps.image_memory_barriers) || !isempty(deps.buffer_memory_barriers)) && Vk.cmd_pipeline_barrier_2_khr(cb, deps)
end

function prepare_render_pass(fg::FrameGraph, pass::Integer)
    g = fg.resource_graph
    attachment_nodes = (resource for resource in neighbors(g, pass) if ResourceClass(resource_type(g, pass, resource)) == RESOURCE_CLASS_ATTACHMENT)
    attachms = Attachment[]
    attach_descs = Vk.AttachmentDescription2[]
    clear_values = Vk.ClearValue[]

    subpass_descriptions = Vk.SubpassDescription2[]
    subpass_dependencies = Vk.SubpassDependency2[]

    for (i, node) in enumerate(attachment_nodes)
        # render pass (global) attachment info
        attachm = attachment(g, i)
        push!(attachms, attachm)
        clear = has_prop(g, pass, node, :clear_value)
        push!(attach_descs, Vk.AttachmentDescription2(attachm, clear, image_layout(g, pass, node), image_layout(g, pass, node), aspect(g, pass, node)))
        clear_val = if clear
            clear_value(g, pass, node)
        else
            # need a filler clear value
            Vk.ClearValue(Vk.ClearColorValue(ntuple(Returns(0f0), 4)))
        end
        push!(clear_values, clear_val)

        # subpass (local) attachment info
        # only one subpass per pass is currently supported
        # we could later look up an internal subpass graph to support more
        for subpass in 1:1
            stage_flags = stages(g, pass)
            pipeline_bind_point = Vk.PIPELINE_STAGE_COMPUTE_SHADER_BIT in stage_flags ? Vk.PIPELINE_BIND_POINT_COMPUTE : Vk.PIPELINE_BIND_POINT_GRAPHICS
            color_attachments = Vk.AttachmentReference2[]
            depth_stencil_attachment = C_NULL
            input_attachments = Vk.AttachmentReference2[]
            type = resource_type(g, pass, node)
            ref = Vk.AttachmentReference2(i - 1, image_layout(g, pass, node), aspect(g, pass, node))
            if RESOURCE_TYPE_COLOR_ATTACHMENT in type
                push!(color_attachments, ref)
            end
            if RESOURCE_TYPE_DEPTH_ATTACHMENT in type || RESOURCE_TYPE_STENCIL_ATTACHMENT in type
                depth_stencil_attachment = ref
            end
            if RESOURCE_TYPE_INPUT_ATTACHMENT in type
                push!(input_attachments, ref)
            end
            push!(subpass_descriptions, Vk.SubpassDescription2(pipeline_bind_point, 0, input_attachments, color_attachments, []; depth_stencil_attachment))

            acc_bits = access_bits(g, pass, node)
            push!(subpass_dependencies, Vk.SubpassDependency2(Vk.SUBPASS_EXTERNAL, subpass - 1, 0; dst_stage_mask = stage_flags, dst_access_mask = acc_bits))
        end
    end
    rp = render_pass(g, pass)
    rp_handle = Vk.RenderPass(fg.device, attach_descs, subpass_descriptions, subpass_dependencies, [])
    (; width, height) = rp.area.extent
    fb = Vk.Framebuffer(fg.device, rp_handle, [att.view for att in attachms], width, height, 1)
    set_prop!(g, pass, :render_pass_handle, rp_handle)
    set_prop!(g, pass, :framebuffer, fb)
    set_prop!(g, pass, :clear_values, clear_values)
end

function begin_render_pass(cb, fg::FrameGraph, pass::Integer)
    g = fg.resource_graph
    Vk.cmd_begin_render_pass_2(cb,
        Vk.RenderPassBeginInfo(
            pass_attribute(g, pass, :render_pass_handle),
            pass_attribute(g, pass, :framebuffer),
            render_pass(g, pass).area,
            pass_attribute(g, pass, :clear_values)
        ),
        Vk.SubpassBeginInfo(Vk.SUBPASS_CONTENTS_INLINE)
    )
end

function synchronize_after(cb, fg, pass)
    nothing
end

"""
Deduce the Vulkan usage, layout and access flags form a resource given its type, stage and access.

The idea is to reconstruct information like `Vk.ACCESS_COLOR_ATTACHMENT_READ_BIT` and `Vk.IMAGE_LAYOUT_SHADER_READ_ONLY_OPTIMAL` from a more decoupled description.
"""
function image_layout(type::ResourceType, access::MemoryAccess, stage::Vk.PipelineStageFlag)
    @match (type, access) begin
        (&RESOURCE_TYPE_COLOR_ATTACHMENT, &READ) => Vk.IMAGE_LAYOUT_SHADER_READ_ONLY_OPTIMAL
        (&RESOURCE_TYPE_COLOR_ATTACHMENT, &WRITE) => Vk.IMAGE_LAYOUT_COLOR_ATTACHMENT_OPTIMAL
        (&RESOURCE_TYPE_COLOR_ATTACHMENT || &RESOURCE_TYPE_IMAGE || RESOURCE_TYPE_TEXTURE, &(READ | WRITE)) => Vk.IMAGE_LAYOUT_GENERAL
        (&RESOURCE_TYPE_DEPTH_ATTACHMENT, &READ) => Vk.IMAGE_LAYOUT_DEPTH_READ_ONLY_OPTIMAL
        (&RESOURCE_TYPE_DEPTH_ATTACHMENT, &WRITE) => Vk.IMAGE_LAYOUT_DEPTH_ATTACHMENT_OPTIMAL
        (&RESOURCE_TYPE_STENCIL_ATTACHMENT, &READ) => Vk.IMAGE_LAYOUT_STENCIL_READ_ONLY_OPTIMAL
        (&RESOURCE_TYPE_STENCIL_ATTACHMENT, &WRITE) => Vk.IMAGE_LAYOUT_STENCIL_ATTACHMENT_OPTIMAL
        (&(RESOURCE_TYPE_DEPTH_ATTACHMENT | RESOURCE_TYPE_STENCIL_ATTACHMENT), &READ) => Vk.IMAGE_LAYOUT_DEPTH_STENCIL_READ_ONLY_OPTIMAL
        (&(RESOURCE_TYPE_DEPTH_ATTACHMENT | RESOURCE_TYPE_STENCIL_ATTACHMENT), &WRITE) => Vk.IMAGE_LAYOUT_DEPTH_STENCIL_ATTACHMENT_OPTIMAL
        (&RESOURCE_TYPE_INPUT_ATTACHMENT || &RESOURCE_TYPE_TEXTURE || &RESOURCE_TYPE_IMAGE, &READ) => Vk.IMAGE_LAYOUT_SHADER_READ_ONLY_OPTIMAL
        _ => error("Unsupported combination of type $type and access $access")
    end
end

function buffer_usage_bits(type::ResourceType, access::MemoryAccess)
    bits = Vk.BufferUsageFlag(0)

    RESOURCE_TYPE_BUFFER | RESOURCE_TYPE_STORAGE in type && (bits |= Vk.BUFFER_USAGE_STORAGE_BUFFER_BIT)
    RESOURCE_TYPE_VERTEX_BUFFER in type && (bits |= Vk.BUFFER_USAGE_VERTEX_BUFFER_BIT)
    RESOURCE_TYPE_INDEX_BUFFER in type && (bits |= Vk.BUFFER_USAGE_INDEX_BUFFER_BIT)
    RESOURCE_TYPE_BUFFER in type && access == READ && (bits |= Vk.BUFFER_USAGE_UNIFORM_BUFFER_BIT)

    bits
end

function image_usage_bits(type::ResourceType, access::MemoryAccess)
    bits = Vk.ImageUsageFlag(0)

    RESOURCE_TYPE_COLOR_ATTACHMENT in type && (bits |= Vk.IMAGE_USAGE_COLOR_ATTACHMENT_BIT)
    (RESOURCE_TYPE_DEPTH_ATTACHMENT in type || RESOURCE_TYPE_STENCIL_ATTACHMENT in type) && (bits |= Vk.IMAGE_USAGE_DEPTH_STENCIL_ATTACHMENT_BIT)
    RESOURCE_TYPE_INPUT_ATTACHMENT in type && (bits |= Vk.IMAGE_USAGE_INPUT_ATTACHMENT_BIT)
    RESOURCE_TYPE_TEXTURE in type && (bits |= Vk.IMAGE_USAGE_SAMPLED_BIT)
    RESOURCE_TYPE_IMAGE in type && WRITE in access && (bits |= Vk.IMAGE_USAGE_STORAGE_BIT)

    bits
end

const SHADER_STAGES = |(
    Vk.PIPELINE_STAGE_VERTEX_SHADER_BIT,
    Vk.PIPELINE_STAGE_TESSELLATION_CONTROL_SHADER_BIT,
    Vk.PIPELINE_STAGE_TESSELLATION_EVALUATION_SHADER_BIT,
    Vk.PIPELINE_STAGE_GEOMETRY_SHADER_BIT,
    Vk.PIPELINE_STAGE_FRAGMENT_SHADER_BIT,
    Vk.PIPELINE_STAGE_COMPUTE_SHADER_BIT,
)

function access_bits(type::ResourceType, access::MemoryAccess, stage::Vk.PipelineStageFlag)
    bits = Vk.AccessFlag(0)
    RESOURCE_TYPE_VERTEX_BUFFER in type && (bits |= Vk.ACCESS_VERTEX_ATTRIBUTE_READ_BIT)
    RESOURCE_TYPE_INDEX_BUFFER in type && (bits |= Vk.ACCESS_INDEX_READ_BIT)
    if RESOURCE_TYPE_COLOR_ATTACHMENT in type
        READ in access && (bits |= Vk.ACCESS_COLOR_ATTACHMENT_READ_BIT)
        WRITE in access && (bits |= Vk.ACCESS_COLOR_ATTACHMENT_WRITE_BIT)
    end
    if (RESOURCE_TYPE_DEPTH_ATTACHMENT in type || RESOURCE_TYPE_STENCIL_ATTACHMENT in type)
        #TODO: support mixed access modes (depth write, stencil read and vice-versa)
        READ in access && (bits |= Vk.ACCESS_DEPTH_STENCIL_ATTACHMENT_READ_BIT)
        WRITE in access && (bits |= Vk.ACCESS_DEPTH_STENCIL_ATTACHMENT_WRITE_BIT)
    end
    RESOURCE_TYPE_INPUT_ATTACHMENT in type && (bits |= Vk.ACCESS_INPUT_ATTACHMENT_READ_BIT)
    if RESOURCE_TYPE_BUFFER in type && !iszero(stage & SHADER_STAGES)
        access == READ && (bits |= Vk.ACCESS_UNIFORM_READ_BIT)
        WRITE in access && (bits |= Vk.ACCESS_SHADER_WRITE_BIT)
    end
    RESOURCE_TYPE_TEXTURE in type && READ in access && (bits |= Vk.ACCESS_SHADER_READ_BIT)
    RESOURCE_TYPE_TEXTURE in type && WRITE in access && (bits |= Vk.ACCESS_SHADER_WRITE_BIT)
    bits
end

function aspect_bits(type::ResourceType)
    bits = Vk.ImageAspectFlag(0)
    RESOURCE_TYPE_COLOR_ATTACHMENT in type && (bits |= Vk.IMAGE_ASPECT_COLOR_BIT)
    RESOURCE_TYPE_DEPTH_ATTACHMENT in type && (bits |= Vk.IMAGE_ASPECT_DEPTH_BIT)
    RESOURCE_TYPE_STENCIL_ATTACHMENT in type && (bits |= Vk.IMAGE_ASPECT_STENCIL_BIT)
    bits
end

function descriptor_type(type::ResourceType, access::MemoryAccess)
    @match access, type begin
        &(RESOURCE_TYPE_IMAGE | RESOURCE_TYPE_STORAGE) => Vk.DESCRIPTOR_TYPE_STORAGE_IMAGE
        &(RESOURCE_TYPE_TEXEL | RESOURCE_TYPE_BUFFER) => Vk.DESCRIPTOR_TYPE_STORAGE_TEXEL_BUFFER
        &(RESOURCE_TYPE_UNIFORM | RESOURCE_TYPE_BUFFER) => Vk.DESCRIPTOR_TYPE_UNIFORM_BUFFER
        &(RESOURCE_TYPE_UNIFORM | RESOURCE_TYPE_BUFFER | RESOURCE_TYPE_DYNAMIC) => Vk.DESCRIPTOR_TYPE_UNIFORM_BUFFER_DYNAMIC
        &(RESOURCE_TYPE_UNIFORM | RESOURCE_TYPE_TEXEL | RESOURCE_TYPE_BUFFER) => Vk.DESCRIPTOR_TYPE_UNIFORM_TEXEL_BUFFER
        &(RESOURCE_TYPE_STORAGE | RESOURCE_TYPE_BUFFER) => Vk.DESCRIPTOR_TYPE_STORAGE_BUFFER
        &(RESOURCE_TYPE_STORAGE | RESOURCE_TYPE_BUFFER | RESOURCE_TYPE_DYNAMIC) => Vk.DESCRIPTOR_TYPE_STORAGE_BUFFER_DYNAMIC
        &RESOURCE_TYPE_TEXTURE => Vk.DESCRIPTOR_TYPE_SAMPLED_IMAGE
        &(RESOURCE_TYPE_TEXTURE | ESOURCE_TYPE_SAMPLER) => Vk.DESCRIPTOR_TYPE_COMBINED_IMAGE_SAMPLER
        &RESOURCE_TYPE_SAMPLER => Vk.DESCRIPTOR_TYPE_SAMPLER
        _ => error("Unsupported combination of type $type and access $access")
    end
end

function resolve_attributes!(fg::FrameGraph)
    g = fg.resource_graph
    for pass in fg.passes
        for resource in neighbors(g, pass)
            type = resource_type(g, pass, resource)
            resource_access = access(g, pass, resource)
            set_prop!(g, resource, :access, access(g, resource) | resource_access)
            pipeline_stages = if has_prop(g, pass, resource, :stages)
                stages(g, pass, resource)
            else
                stages(g, pass)
            end
            # class-independent attributes
            acc_bits = access_bits(type, resource_access, pipeline_stages)
            set_prop!(g, pass, resource, :access_bits, acc_bits)

            # class-dependent attributes
            class = resource_class(g, resource)
            if class == RESOURCE_CLASS_IMAGE || class == RESOURCE_CLASS_ATTACHMENT
                usage = image_usage_bits(type, resource_access)
                set_prop!(g, pass, resource, :usage, usage)
                layout = image_layout(type, resource_access, pipeline_stages)
                set_prop!(g, pass, resource, :image_layout, layout)
                aspect = aspect_bits(type)
                set_prop!(g, pass, resource, :aspect, aspect)
                set_prop!(g, resource, :usage, usage | image_usage(g, resource))
            else
                usage = buffer_usage_bits(type, resource_access)
                set_prop!(g, pass, resource, :usage, usage)
                set_prop!(g, resource, :usage, usage | buffer_usage(g, resource))
            end
        end
    end
end

function create_physical_resources!(fg::FrameGraph)
    g = fg.resource_graph
    for resource in fg.resources |> Filter(x -> !has_prop(g, x, :physical_resource))
        c = resource_class(g, resource)
        if c == RESOURCE_CLASS_IMAGE || c == RESOURCE_CLASS_ATTACHMENT
            image = ImageBlock(device(fg), Tuple(_dims(g, resource)), format(g, resource), image_usage(g, resource))
            allocate!(image, MEMORY_DOMAIN_DEVICE)
            if c == RESOURCE_CLASS_ATTACHMENT
                set_prop!(g, resource, :physical_resource, Attachment(View(image), access(g, resource)))
            else
                set_prop!(g, resource, :physical_resource, image)
            end
        else
            buffer = BufferBlock(device(fg), _size(g, resource); usage = buffer_usage(g, resource))
            allocate!(buffer, MEMORY_DOMAIN_DEVICE)
            set_prop!(g, resource, :physical_resource, buffer)
        end
    end
end

"""
Return whether `x` covers accesses of type `y`; that is, whether guarantees about memory access operations in `x` induce guarantees about memory access operations in `y`.

```jldoctest
julia> covers(Vk.ACCESS_MEMORY_WRITE_BIT, Vk.ACCESS_SHADER_WRITE_BIT)
true

julia> covers(Vk.ACCESS_SHADER_READ_BIT, Vk.ACCESS_UNIFORM_READ_BIT)
true

julia> covers(Vk.ACCESS_SHADER_WRITE_BIT, Vk.ACCESS_MEMORY_WRITE_BIT)
false
```

"""
function covers(x::Vk.AccessFlag, y::Vk.AccessFlag)
    y in x && return true
    if Vk.ACCESS_MEMORY_READ_BIT in x
        if Vk.ACCESS_MEMORY_READ_BIT ∉ y
            x &= ~Vk.ACCESS_MEMORY_READ_BIT
        end
        x |= |(
            Vk.ACCESS_ACCELERATION_STRUCTURE_READ_BIT_KHR,
            Vk.ACCESS_COLOR_ATTACHMENT_READ_BIT,
            Vk.ACCESS_COLOR_ATTACHMENT_READ_NONCOHERENT_BIT_EXT,
            Vk.ACCESS_COMMAND_PREPROCESS_READ_BIT_NV,
            Vk.ACCESS_CONDITIONAL_RENDERING_READ_BIT_EXT,
            Vk.ACCESS_DEPTH_STENCIL_ATTACHMENT_READ_BIT,
            Vk.ACCESS_FRAGMENT_DENSITY_MAP_READ_BIT_EXT,
            Vk.ACCESS_HOST_READ_BIT,
            Vk.ACCESS_INDEX_READ_BIT,
            Vk.ACCESS_INDIRECT_COMMAND_READ_BIT,
            Vk.ACCESS_INPUT_ATTACHMENT_READ_BIT,
            Vk.ACCESS_SHADER_READ_BIT,
            Vk.ACCESS_SHADING_RATE_IMAGE_READ_BIT_NV,
            Vk.ACCESS_TRANSFER_READ_BIT,
            Vk.ACCESS_TRANSFORM_FEEDBACK_COUNTER_READ_BIT_EXT,
            Vk.ACCESS_UNIFORM_READ_BIT,
            Vk.ACCESS_VERTEX_ATTRIBUTE_READ_BIT,
        )
        y in x && return true
    end
    if Vk.ACCESS_MEMORY_WRITE_BIT in x
        if Vk.ACCESS_MEMORY_WRITE_BIT ∉ y
            x &= ~Vk.ACCESS_MEMORY_WRITE_BIT
        end
        x |= |(
            Vk.ACCESS_ACCELERATION_STRUCTURE_WRITE_BIT_KHR,
            Vk.ACCESS_COLOR_ATTACHMENT_WRITE_BIT,
            Vk.ACCESS_COMMAND_PREPROCESS_WRITE_BIT_NV,
            Vk.ACCESS_DEPTH_STENCIL_ATTACHMENT_WRITE_BIT,
            Vk.ACCESS_HOST_WRITE_BIT,
            Vk.ACCESS_SHADER_WRITE_BIT,
            Vk.ACCESS_TRANSFER_WRITE_BIT,
            Vk.ACCESS_TRANSFORM_FEEDBACK_COUNTER_WRITE_BIT_EXT,
            Vk.ACCESS_TRANSFORM_FEEDBACK_WRITE_BIT_EXT,
        )
        y in x && return true
    end
    if Vk.ACCESS_SHADER_READ_BIT in x
        if Vk.ACCESS_SHADER_READ_BIT ∉ y
            x &= ~Vk.ACCESS_SHADER_READ_BIT
        end
        x |= Vk.ACCESS_UNIFORM_READ_BIT
        y in x && return true
    end
    false
end

function ResourceClass(type::ResourceType)
    RESOURCE_TYPE_TEXTURE in type && return RESOURCE_CLASS_IMAGE
    RESOURCE_TYPE_IMAGE in type && return RESOURCE_CLASS_IMAGE
    RESOURCE_TYPE_INPUT_ATTACHMENT in type && return RESOURCE_CLASS_ATTACHMENT
    RESOURCE_TYPE_COLOR_ATTACHMENT in type && return RESOURCE_CLASS_ATTACHMENT
    RESOURCE_TYPE_DEPTH_ATTACHMENT in type && return RESOURCE_CLASS_ATTACHMENT
    RESOURCE_TYPE_STENCIL_ATTACHMENT in type && return RESOURCE_CLASS_ATTACHMENT
    RESOURCE_TYPE_BUFFER in type && return RESOURCE_CLASS_BUFFER
    RESOURCE_TYPE_VERTEX_BUFFER in type && return RESOURCE_CLASS_BUFFER
    RESOURCE_TYPE_INDEX_BUFFER in type && return RESOURCE_CLASS_BUFFER
    error("Resource type '$type' does not belong to any resource class.")
end


text:
using StringParserPEG

function createRecord(rule,value,first,last,children)
  if length(children) == 0
    return String[]
  elseif length(children) == 1
    return children[1]
  else
    return append!([children[1]],children[2].children)
  end
end
csv = Grammar("""
  start => *((record & -(crlf)) {liftchild}) {childrenarray}
  record => ( field & *((-(',') & field){liftchild}) ) {Main.createRecord}
  field => ?(escaped_field | unescaped_field) {liftchild}

  escaped_field => (-('"') & escaped_field_value & -('"')) {liftchild}
  escaped_field_value => r([^"]+)r {nodevalue}
  unescaped_field => r([^\\n\\r,]+)r {nodevalue}

  crlf => r([\\n\\r]+)r
""")

data = """
1,2,3
4,5,6

this,is,a,"test and

only a test"
"""

(ast, pos, error) = parse(csv, data)
println("ast = $ast")

ext: jl


text:
# This method is called when a SemidiscretizationHyperbolic is constructed.
# It constructs the basic `cache` used throughout the simulation to compute
# the RHS etc.
function create_cache(mesh::P4estMesh, equations::AbstractEquations, dg::DG, ::Any, ::Type{uEltype}) where {uEltype<:Real}
  # Make sure to balance the p4est before creating any containers
  # in case someone has tampered with the p4est after creating the mesh
  p4est_balance(mesh.p4est, P4EST_CONNECT_FACE, C_NULL)
  # Due to a bug in p4est, the forest needs to be rebalanced twice sometimes
  # See https://github.com/cburstedde/p4est/issues/112
  p4est_balance(mesh.p4est, P4EST_CONNECT_FACE, C_NULL)

  elements   = init_elements(mesh, equations, dg.basis, uEltype)
  interfaces = init_interfaces(mesh, equations, dg.basis, elements)
  boundaries = init_boundaries(mesh, equations, dg.basis, elements)
  mortars    = init_mortars(mesh, equations, dg.basis, elements)

  cache = (; elements, interfaces, boundaries, mortars)

  # Add specialized parts of the cache required to compute the volume integral etc.
  cache = (;cache..., create_cache(mesh, equations, dg.mortar, uEltype)...)

  return cache
end

# Extract outward-pointing normal vector (contravariant vector ±Ja^i, i = index) as SVector
# Note that this vector is not normalized
@inline function get_normal_vector(direction, cache, indices...)
  @unpack contravariant_vectors, inverse_jacobian = cache.elements

  # If the mapping is orientation-reversing, the contravariant vectors' orientation
  # is reversed as well
  sign_jacobian = sign(inverse_jacobian[indices...])

  orientation = div(direction + 1, 2)
  normal = sign_jacobian * get_contravariant_vector(orientation, contravariant_vectors, indices...)

  # Contravariant vectors at interfaces in negative coordinate direction are pointing inwards
  # (after normalizing with sign(J) above)
  if direction in (1, 3, 5)
    normal *= -1
  end

  return normal
end

@inline ndofs(mesh::P4estMesh, dg::DG, cache) = nelements(cache.elements) * nnodes(dg)^ndims(mesh)

include("containers.jl")
include("dg_2d.jl")


text:
module PoreMatMod

using DataFrames, FIGlet, Graphs, LinearAlgebra, MetaGraphs, Pluto, Reexport, StatsBase
@reexport using Xtals

import Base.(∈), Base.show, Base.replace

function __init__()
    add_bonding_rules(tagged_bonding_rules())
end

export
    # search.jl
    Search, substructure_search, nb_isomorphisms, nb_locations, nb_ori_at_loc, isomorphic_substructures,

    # replace.jl
    substructure_replace,

    # moiety.jl
    moiety,

    # misc.jl
    PoreMatModGO

include("Ullmann.jl")
include("moiety.jl")
include("search.jl")
include("replace.jl")
include("ExampleHelpers.jl")
include("misc.jl")

end


text:
using Test, Metalhead
using Flux
using Flux: Zygote

function gradtest(model, input)
  y, pb = Zygote.pullback(() -> model(input), Flux.params(model))
  gs = pb(ones(Float32, size(y)))

  # if we make it to here with no error, success!
  return true
end

# CNN tests
@testset verbose = true "ConvNets" begin
  include("convnets.jl")
end

GC.gc()

# Other tests
@testset verbose = true "Other" begin
  include("other.jl")
end

GC.gc()

# ViT tests
@testset verbose = true "ViTs" begin
  include("vit-based.jl")
end


text:
push!(LOAD_PATH,"../src/")
using Documenter, GPLinearODEMaker

makedocs(
    modules=[GPLinearODEMaker],
    sitename="GPLinearODEMaker.jl",
    authors = "Christian Gilbertson")
    # pages = ["Home" => "index.md"])

deploydocs(repo = "github.com/christiangil/GPLinearODEMaker.jl.git",
           target = "build",
           push_preview = true)


text:
#=
Description:  Bubble Sort
Author: Tushar Gautam
Date: 8th Jan 2022
=#

using BenchmarkTools

# Length of array
LEN = 100
# Min/Max range for array elements
min = 0.
max = 100.

# Initialising Arrays
arr = rand(min:eps():max, LEN)

println("Array: ", round.(arr, digits=3))

# Bubble sort
for i = 0:LEN-2
    for j = 2:LEN-i
        if (arr[j-1] > arr[j])
            tmp = arr[j-1]
            arr[j-1] = arr[j]
            arr[j] = tmp
        end
    end
end

println("Sorted Array: ", round.(arr, digits=3))


text:
using Rootionals
using Test

R = Rootional

@testset "Init" begin
    @test R(1) == R([1], [1])
    @test R(0) == R([0], [1])

    @test R(3//5) == R([3, 5], [1, -1])

    @test one(R) == R(1)
    @test zero(R) == R(0)
end

@testset "Access" begin
    @test numerator(R([2], [1//2])) == R([2], [1//2])
    @test numerator(R([2], [-1//2])) == one(R)
    @test denominator(R([2], [1//2])) == one(R)
    @test denominator(R([2], [-1//2])) == R([2], [1//2])
end

@testset "Simplify" begin
    @test simplify!(R([3,3], [4, -1//2])) == R([3], [7//2])
    @test simplify!(R([2], [0])) == one(R)
    @test simplify!(R([7], [-0//2])) == one(R)
end

@testset "Conversions" begin
    @test one(R) == 1
    @test 5 == R(5)
    @test R(3//5) == 3//5
    @test 7//7 == one(R)

    @test float(R(2)) == 2.0
    @test float(R([2], [1//2])) ≈ sqrt(2)
    @test float(R([2], [-1//2])) ≈ 1/sqrt(2)
end

@testset "Arithmetic" begin
    @testset "*" begin
        @test R(2)*R(1//3) == R(2//3)
        @test R(2)*R(1//2) == one(R)
        @test R(2)*R([2], [1//2]) == R([2], [3//2])
        @test R(2)*R([2], [-1//2]) == R([2], [1//2])
        @test 2R([2], [-1//2]) == R([2], [1//2])
        @test R([2], [-1//2])*1//2 == R([2], [-3//2])
    end

    @testset "/" begin
        @test R(2)/R(2) == one(R)
        @test R([2], [1//2])/2 == R([2], [-1//2])
        @test 2/R([4, 5], [1, 1//2]) == R([2, 5], [-1, -1//2])
        @test 2/R(16) == R([2], [-3])
    end

    @testset "^" begin
        @test R([2], [1//2])^2 == 2
        @test (2R([2], [1//2]))^2 == 8
        @test R(4//3)^(1//2) == R([2, 3], [1, -1//2])
    end

    @testset "√" begin
        @test sqrt(R([2], [1//2])) == R([2], [1//4])
        @test cbrt(R([2], [1//2])) == R([2], [1//6])
    end
end

@testset "Strings" begin
    @test string(R([2], [1//2])) == "√2"
    @test string(R([2], [-1//2])) == "1/√2"
    @test string(R([3, 2, 6], [3//4, -2, -1//3])) == "(∜3)³/2²∛6"
    @test string(R([5], [17//8])) == "(⁸√5)¹⁷"
end


text:
using Test
using Carnival

@testset "My test" begin
    @test 1 + 1 == 2
    @test foo() == -1
end

max_stars_repo_name: charleskawczynski/Carnival

max_issues_repo_name: charleskawczynski/Carnival

max_forks_repo_name: charleskawczynski/Carnival


text:
@testset "Flux variability analysis" begin
    cp = test_simpleLP()
    optimizer = GLPK.Optimizer
    fluxes = flux_variability_analysis(cp, optimizer)

    @test size(fluxes) == (2, 2)
    @test fluxes ≈ [
        1.0 1.0
        2.0 2.0
    ]

    fluxes = flux_variability_analysis(cp, [2], optimizer)

    @test size(fluxes) == (1, 2)
    @test fluxes == Matrix{Float64}([2 2])

    # a special testcase for slightly sub-optimal FVA (gamma<1)
    cp = CoreModel(
        [-1.0 -1.0 -1.0],
        [0.0],
        [1.0, 0.0, 0.0],
        [0.0, 0.0, -1.0],
        1.0 * ones(3),
        ["r$x" for x = 1:3],
        ["m1"],
    )
    fluxes = flux_variability_analysis(cp, optimizer)
    @test fluxes ≈ [
        1.0 1.0
        0.0 0.0
        -1.0 -1.0
    ]
    fluxes = flux_variability_analysis(cp, optimizer; bounds = gamma_bounds(0.5))
    @test fluxes ≈ [
        0.5 1.0
        0.0 0.5
        -1.0 -0.5
    ]
    fluxes = flux_variability_analysis(cp, optimizer; bounds = _ -> (0, Inf))
    @test fluxes ≈ [
        0.0 1.0
        0.0 1.0
        -1.0 0.0
    ]

    @test isempty(flux_variability_analysis(cp, Vector{Int}(), GLPK.Optimizer))
    @test_throws DomainError flux_variability_analysis(cp, [-1], GLPK.Optimizer)
    @test_throws DomainError flux_variability_analysis(cp, [99999999], GLPK.Optimizer)
end

@testset "Parallel FVA" begin
    cp = test_simpleLP()
    pids = addprocs(2, topology = :master_worker)
    @everywhere using COBREXA, GLPK
    fluxes = flux_variability_analysis(cp, [1, 2], GLPK.Optimizer; workers = pids)
    @test fluxes ≈ [
        1.0 1.0
        2.0 2.0
    ]
    rmprocs(pids)
end

@testset "Flux variability analysis with StandardModel" begin
    model_path = download_data_file(
        "http://bigg.ucsd.edu/static/models/e_coli_core.json",
        joinpath("data", "e_coli_core.json"),
        "7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8",
    )

    model = load_model(StandardModel, model_path)
    mins, maxs = flux_variability_analysis_dict(
        model,
        Tulip.Optimizer;
        bounds = objective_bounds(0.99),
        modifications = [
            change_optimizer_attribute("IPM_IterationsLimit", 500),
            change_constraint("EX_glc__D_e", -10, -10),
            change_constraint("EX_o2_e", 0.0, 0.0),
        ],
    )

    @test isapprox(maxs["EX_ac_e"]["EX_ac_e"], 16.81521027177164, atol = TEST_TOLERANCE)
    @test isapprox(mins["EX_ac_e"]["EX_ac_e"], 1.553675043862341e-10, atol = TEST_TOLERANCE)
end


text:
using ExoplanetsSysSim
using ApproximateBayesianComputing
const ABC = ApproximateBayesianComputing
using SpecialFunctions
using Statistics
import ApproximateBayesianComputing.CompositeDistributions.CompositeDist
import ApproximateBayesianComputing.TransformedBetaDistributions.LinearTransformedBeta
#import EvalSysSimModel

# https://en.wikipedia.org/wiki/Trigamma_function
function trigamma_x_gr_4(x::T) where T<: Real
   1/x + 0.5/x^2 + 1/(6*x^3) - 1/(30*x^5) + 1/(42*x^7) - 1/(30*x^9) + 5/(66*x^11) - 691/(2730*x^13) + 7/(6*x^15)
end

function trigamma_x_lt_4(x::T) where T<: Real
  n = floor(Int64,5-x)
  z = x+n
  val = trigamma_x_gr_4(z)
  for i in 1:n
    z -= 1
    val += 1/z^2
  end
  val
end

function trigamma(x::T) where T<: Real
   x >= 4 ? trigamma_x_gr_4(x) : trigamma_x_lt_4(x)
end

function var_weighted(x::AbstractArray{Float64,1}, w::AbstractArray{Float64,1} )
  #println("# size(x) = ",size(x), " size(w) = ", size(w)); flush(stdout)
  @assert(length(x)==length(w) )
  sumw = sum(w)
  @assert( sumw > 0. )
  if(sumw!= 1.0)
     w /= sum(w)
     sumw = 1.0
  end
  sumw2 = sum(w.*w)
  xbar =  sum(x.*w)
  covar = sum((x.-xbar).*(x.-xbar) .* w) * sumw/(sumw*sumw-sumw2)
end

function mom_alpha(x_bar::T, v_bar::T) where T<: Real
    x_bar * (((x_bar * (1 - x_bar)) / v_bar) - 1)
end
function mom_beta(x_bar::T, v_bar::T) where T<: Real
    (1 - x_bar) * (((x_bar * (1 - x_bar)) / v_bar) - 1)
end

# For algorithm, see https://scholarsarchive.byu.edu/cgi/viewcontent.cgi?article=2613&context=etd
function fit_beta_mle(x::AbstractArray{T,1}; tol::T = 1e-6, max_it::Int64 = 10, init_guess::AbstractArray{T,1} = Array{T}(undef,0), w::AbstractArray{T,1} = Array{T}(undef,0), verbose::Bool = false ) where T<: Real
    lnxbar =   length(w)>1 ? Statistics.mean(log.(x),AnalyticWeights(w)) : Compat.Statistics.mean(log.(x))
    ln1mxbar = length(w)>1 ? Statistics.mean(log.(1.0.-x),AnalyticWeights(w)) : Compat.Statistics.mean(log.(1.0.-x))

    function itterate( mle_guess::Vector{T} ) where T<:Real
        (alpha, beta) = (mle_guess[1], mle_guess[2])
        dgab = digamma(alpha+beta)
        g1 = dgab - digamma(alpha) + lnxbar
        g2 = dgab - digamma(beta) + ln1mxbar
        tgab = trigamma(alpha+beta)
        G = [dgab-trigamma(alpha) tgab; tgab tgab-trigamma(beta)]
        mle_guess -= G \ [g1, g2]
    end

    local mle_new
    if length(init_guess) != 2
        xbar = length(w)>1 ? Compat.Statistics.mean(x,AnalyticWeights(w)) : Compat.Statistics.mean(x)
        vbar = length(w)>1 ? Compat.Statistics.varm(x,xbar,AnalyticWeights(w)) : Compat.Statistics.varm(x,xbar)
        mle_new = (vbar < xbar*(1.0-xbar)) ? [mom_alpha(xbar, vbar), mom_beta(xbar,vbar)] : ones(T,2)
    else
        mle_new = init_guess
    end
    if verbose
        println("it = 0: ", mle_new)
    end
    if any(mle_new.<=zero(T))
        println("# Warning: mean= ", xbar, " var= ",vbar," (alpha,beta)_init= ",mle_new," invalid, reinitializing to (1,1)")
        verbose = true
        mle_new = ones(T,2)
    end
    for i in 1:max_it
        mle_old = mle_new
        mle_new = itterate( mle_old )
        epsilon = max(abs.(mle_old.-mle_new))
        if verbose
            println("# it = ", i, ": ", mle_new, " max(Delta alpha, Delta beta)= ", epsilon)
        end
        if epsilon < tol
            break
        end
    end
    return mle_new
end

function make_beta(x::AbstractArray{T,1}, w::AbstractArray{T,1};
                   xmean::T = Compat.Statistics.mean(x,AnalyticWeights(w)),
                   xvar::T = Compat.Statistics.varm(x,xmean,AnalyticWeights(w)), tau_factor::T=one(T) ) where T<:Real
    alpha_beta = (xvar < xmean*(1.0-xmean)) ? [mom_alpha(xmean, xvar), mom_beta(xmean,xvar)] : ones(T,2)
    if any(alpha_beta.<=zero(T))
        alpha_beta = fit_beta_mle(x, w=w, init_guess=alpha_beta, verbose=true)
    end
    if any(alpha_beta.<=zero(T))
        alpha_beta = ones(T,2)
    else
        if minimum(alpha_beta)>1.5*tau_factor && sum(alpha_beta)>=20.0*tau_factor
            alpha_beta ./= tau_factor
        end
    end
    #println("Radius relative: a= ",alpha_beta[1], "  b= ",alpha_beta[2])
    Beta(alpha_beta[1], alpha_beta[2])
end

function make_beta_transformed(x::AbstractArray{T,1}, w::AbstractArray{T,1}; xmin::T=zero(T), xmax::T=one(T), xmean::T = Compat.Statistics.mean(x,AnalyticWeights(w)), xvar::T = Compat.Statistics.varm(x,xmean,AnalyticWeights(w)), tau_factor::T=one(T) ) where T<:Real
    alpha_beta = (xvar < xmean*(1.0-xmean)) ? [mom_alpha(xmean, xvar), mom_beta(xmean,xvar)] : ones(T,2)
    if any(alpha_beta.<=zero(T))
        alpha_beta = fit_beta_mle(x, w=w, init_guess=alpha_beta, verbose=true)
    end
    if any(alpha_beta.<=zero(T))
        alpha_beta = ones(T,2)
    else
        if minimum(alpha_beta)>1.5*tau_factor && sum(alpha_beta)>=20.0*tau_factor
            alpha_beta ./= tau_factor
        end
    end
    #println("Total: a= ",alpha_beta[1], "  b= ",alpha_beta[2])
    LinearTransformedBeta(alpha_beta[1], alpha_beta[2], xmin=xmin, xmax=xmax)
end

function make_proposal_dist_multidim_beta(theta::AbstractArray{Float64,2}, weights::AbstractArray{Float64,1},  tau_factor::Float64; verbose::Bool = false)
    global sim_param_closure
    local limitP::Array{Float64,1} = get_any(EvalSysSimModel.sim_param_closure, "p_lim_arr", Array{Float64,1})
    local p_dim = length(limitP)-1
    local r_dim = length(get_any(EvalSysSimModel.sim_param_closure, "r_lim_arr", Array{Float64,1}))-1

    theta_mean =  sum(theta.*weights',dims=2) # weighted mean for parameters
    theta_var = ABC.var_weighted(theta'.-theta_mean',weights)  # scaled, weighted covar for parameters
    tau_factor_indiv = fill(tau_factor,length(theta_var))

    dist_arr = ContinuousDistribution[]
    for j in 1:p_dim
        max_col_rate = 3*log(limitP[j+1]/limitP[j])/log(2)
        col_startidx = (j-1)*r_dim+1
        #tau_factor_indiv[col_startidx] = 2.0

        #=
        println("mean= ",theta_mean)
        println("var= ",theta_var)
        println("tau= ",tau_factor_indiv)
        for i in 1:length(theta_mean)
            println("a= ",alpha(theta_mean[i],tau_factor*theta_var[i]), "  b= ",beta(theta_mean[i],tau_factor*theta_var[i]))
        end
        =#

        dist_arr = vcat(dist_arr, ContinuousDistribution[make_beta_transformed(theta[i,:], weights, xmin=0.0, xmax=max_col_rate, xmean=theta_mean[i]/max_col_rate, xvar=theta_var[i]/max_col_rate^2, tau_factor=tau_factor_indiv[i]) for i in (col_startidx):(col_startidx+r_dim-1)])
    end

dist = CompositeDist(dist_arr)
end

function make_proposal_dist_multidim_beta(pop::abc_population_type, tau_factor::Float64; verbose::Bool = false)
    make_proposal_dist_multidim_beta(pop.theta, pop.weights, tau_factor, verbose=verbose)
end

function make_proposal_dist_multidim_beta_dirichlet(theta::AbstractArray{Float64,2}, weights::AbstractArray{Float64,1},  tau_factor::Float64; verbose::Bool = false)
    global sim_param_closure
    local limitP::Array{Float64,1} = get_any(EvalSysSimModel.sim_param_closure, "p_lim_arr", Array{Float64,1})
    local p_dim = length(limitP)-1
    local r_dim = length(get_any(EvalSysSimModel.sim_param_closure, "r_lim_arr", Array{Float64,1}))-1

    theta_mean =  sum(theta.*weights',dims=2) # weighted mean for parameters
    theta_var = ABC.var_weighted(theta'.-theta_mean',weights)  # scaled, weighted covar for parameters
    tau_factor_indiv = fill(tau_factor,length(theta_var))

    dist_arr = ContinuousDistribution[]
    for j in 1:p_dim
        max_col_rate = 3*log(limitP[j+1]/limitP[j])/log(2)
        col_startidx = (j-1)*(r_dim+1)+1
        #tau_factor_indiv[col_startidx] = 2.0

        # if verbose
        #     println("total: ",theta_mean[col_startidx]," ",theta_var[col_startidx])
        # end
        # for i in (col_startidx+1):(col_startidx+r_dim)
        #     mean_ratio = sum(theta[col_startidx,:].*theta[i,:].*weights) /(theta_mean[col_startidx]*theta_mean[i]) # weighted mean for parameters
        #     var_ratio = var_weighted(vec(theta[col_startidx,:].*theta[i,:]).-(theta_mean[col_startidx]*theta_mean[i]),weights)/(2 * theta_mean[col_startidx] * theta_var[i]) # scaled, weighted covar for parameters
        #     if verbose
        #         println("i=",i,": ",theta_mean[i]," ",theta_var[i]," ratios: ",mean_ratio, " ",var_ratio)
        #     end
        #     var_ratio  = var_ratio  >= one(var_ratio)  ? var_ratio  : one(var_ratio)
        #     tau_factor_indiv[i] = tau_factor*var_ratio
        # end
        # if verbose
        #     flush(stdout)
        # end

        #=
        println("mean= ",theta_mean)
        println("var= ",theta_var)
        println("tau= ",tau_factor_indiv)
        for i in 1:length(theta_mean)
            println("a= ",alpha(theta_mean[i],tau_factor*theta_var[i]), "  b= ",beta(theta_mean[i],tau_factor*theta_var[i]))
        end
        =#

        dist_arr = vcat(dist_arr, make_beta_transformed(theta[col_startidx,:], weights, xmin=0.0, xmax=max_col_rate, xmean=theta_mean[col_startidx]/max_col_rate, xvar=theta_var[col_startidx]/max_col_rate^2, tau_factor=tau_factor_indiv[col_startidx]), ContinuousDistribution[ make_beta(theta[i,:], weights, xmean=theta_mean[i], xvar=theta_var[i], tau_factor=tau_factor_indiv[i]) for i in (col_startidx+1):(col_startidx+r_dim)])
    end

dist = CompositeDist(dist_arr)
end

function make_proposal_dist_multidim_beta_dirichlet(pop::abc_population_type, tau_factor::Float64; verbose::Bool = false)
    make_proposal_dist_multidim_beta_dirichlet(pop.theta, pop.weights, tau_factor, verbose=verbose)
end


text:
@testset "sql" begin

sess = SparkSession()

# testing IO
mktempdir() do dir
    parquet_file = joinpath(dir, "people.parquet")
    ds = read_json(sess, joinpath(@__DIR__, "people.json"))
    write_parquet(ds, parquet_file)
    ds2 = read_parquet(sess, parquet_file)
    write_json(ds2, joinpath(dir, "people2.json"))
end

close(sess)

end


text:
using Distributed

@everywhere using Distributions, Random, LinearAlgebra, Statistics

#-------------------------------
# Notation for random variables
#-------------------------------
#
# f is a vector of length N
# K is the covariance matrix of size N×N
# α and β are vectors of length F, i.e. number of filters
# A is a N×N diagonal matrix whose elements take values in α
# b is N length vector whose elements take vales in β
# S is a N×N diagonal matrix that holds the squared standard errors of the measurements
# y are the observations collected from all filters, has length N

#-------------------------------
# Notation for posteriors
#-------------------------------
#
# q(α) = 𝓝(α|μα,Σα)
# q(β) = 𝓝(β|μβ,Σβ)
# q(f) = 𝓝(f|μf,Σf)

# In this file we code and numerically verify the expectations of the following terms:
#
# log𝓝(y|Af + b, S) + log𝓝(f|0,K) + log𝓝(a|0,v²ₐI) + log𝓝(b|0,v²ᵦ)
#        (1)               (2)             (3)             (4)
#

function check_functions_of_terms(seed=0)

    # set dummy dimensions and instantiate dummy values for the random variables
    N   = 6
    F   = 2
    rg  = MersenneTwister(seed)

    idx = ceil.(Int, rand(rg, N)*F)

    v²ₐ = 0.3
    μα = rand(rg, F)
    Σα = Diagonal(rand(rg, F))*0.2
    α  = rand(rg, MvNormal(μα, Σα))
    A  = Diagonal(α[idx])

    v²ᵦ = 0.4
    μβ = rand(rg, F)
    Σβ = Diagonal(rand(rg, F))*0.2
    β  = rand(rg, MvNormal(μβ, Σβ))
    b  = β[idx]

    K  = randn(rg,N,N)*0.2 + 0.1*I; K=K*K'
    μf = rand(rg,N)
    Σf = randn(rg,N,N)*0.2 + 0.1*I; Σf = Σf*Σf'
    f  = rand(rg, MvNormal(μf, Σf))

    S = Diagonal(rand(rg,N) .+ 0.1)
    y = rand(rg, MvNormal(f, S))

    #--------------------------
    # define terms
    #--------------------------

    t1 = logpdf(MvNormal(A*f + b, S), y)          # (1)

    t2 = logpdf(MvNormal(zeros(N), K), f)         # (2)

    t3 = logpdf(MvNormal(zeros(F), sqrt(v²ₐ)), α) # (3)

    t4 = logpdf(MvNormal(zeros(F), sqrt(v²ᵦ)), β) # (4)

    @assert(t1 ≈ term1(y, A, f, b, S) ≈ term1_broken_down(y, A, f, b, S))
    @assert(t2 ≈ term2(f, K))
    @assert(t3 ≈ term3(α, v²ₐ))
    @assert(t4 ≈ term4(β, v²ᵦ))

    return true

end

function numerical_expectation_inside_log_of_lower_bound(; y=y, idx=idx, S=S, μα=μα, Σα=Σα, v²ₐ=v²ₐ, μβ=μβ, Σβ=Σβ, v²ᵦ=v²ᵦ, μf=μf, Σf=Σf, K=K, numsamples=numsamples)

    qα = MvNormal(μα, Σα)
    qβ = MvNormal(μβ, Σβ)
    qf = MvNormal(μf, Σf)

    @distributed (+) for i=1:numsamples

        local α = rand(qα)
        local β = rand(qβ)
        local f = rand(qf)
        local A = Diagonal(α[idx])
        local b = β[idx]

        (term1(y, A, f, b, S) + term2(f, K) + term3(α, v²ₐ) + term4(β, v²ᵦ)) / numsamples

    end

end

function test_expectations(numsamples=1000; seed=0)

    # set dummy dimensions and instantiate dummy values for the random variables
    N   = 6
    F   = 2
    rg  = MersenneTwister(seed)

    idx = ceil.(Int, rand(rg, N)*F)

    v²ₐ = 0.3
    μα = rand(rg, F)
    Σα = Diagonal(rand(rg, F))*0.2
    α  = rand(rg, MvNormal(μα, Σα))
    A  = Diagonal(α[idx])

    v²ᵦ = 0.4
    μβ = rand(rg, F)
    Σβ = Diagonal(rand(rg, F))*0.2
    β  = rand(rg, MvNormal(μβ, Σβ))
    b  = β[idx]

    K  = randn(rg,N,N)*0.2 + 0.0*I; K=K*K'
    μf = rand(rg,N)
    Σf = randn(rg,N,N)*0.2 + 0.0*I; Σf = Σf*Σf'
    f  = rand(rg, MvNormal(μf, Σf))

    S = Diagonal(rand(rg,N) .+ 0.1)
    y = rand(rg, MvNormal(f, S))

    # #--------------------------
    # # Define posteriors
    # #--------------------------
    #
    # qα = MvNormal(μα, Σα)
    # qβ = MvNormal(μβ, Σβ)
    # qf = MvNormal(μf, Σf)

        exact_expectation_inside_log_of_lower_bound(; y=y, idx=idx, S=S, μα=μα, Σα=Σα, v²ₐ=v²ₐ, μβ=μβ, Σβ=Σβ, v²ᵦ=v²ᵦ, μf=μf, Σf=Σf, K=K),
    numerical_expectation_inside_log_of_lower_bound(; y=y, idx=idx, S=S, μα=μα, Σα=Σα, v²ₐ=v²ₐ, μβ=μβ, Σβ=Σβ, v²ᵦ=v²ᵦ, μf=μf, Σf=Σf, K=K, numsamples=numsamples)

end

function inside_log_of_lower_bound(numsamples=1000; seed=0)

    # set dummy dimensions and instantiate dummy values for the random variables
    N   = 6
    F   = 2
    rg  = MersenneTwister(seed)

    idx = ceil.(Int, rand(rg, N)*F)

    v²ₐ = 0.3
    μα = rand(rg, F)
    Σα = Diagonal(rand(rg, F))*0.2
    α  = rand(rg, MvNormal(μα, Σα))
    A  = Diagonal(α[idx])

    v²ᵦ = 0.4
    μβ = rand(rg, F)
    Σβ = Diagonal(rand(rg, F))*0.2
    β  = rand(rg, MvNormal(μβ, Σβ))
    b  = β[idx]

    K  = randn(rg,N,N)*0.2 + 0.0*I; K=K*K'
    μf = rand(rg,N)
    Σf = randn(rg,N,N)*0.2 + 0.0*I; Σf = Σf*Σf'
    f  = rand(rg, MvNormal(μf, Σf))

    S = Diagonal(rand(rg,N) .+ 0.1)
    y = rand(rg, MvNormal(f, S))

    #--------------------------
    # Define posteriors
    #--------------------------

    qα = MvNormal(μα, Σα)
    qβ = MvNormal(μβ, Σβ)
    qf = MvNormal(μf, Σf)

    #--------------------------------
    # verify expectation of term (1)
    #--------------------------------

    numexpectationterm1 = @distributed (+) for i=1:numsamples

        local α = rand(qα)
        local β = rand(qβ)
        local f = rand(qf)
        local A = Diagonal(α[idx])
        local b = β[idx]

        term1(y, A, f, b, S) / numsamples
    end

    @show expectationterm1(;y=y, μα=μα, Σα=Σα, μβ=μβ, Σβ=Σβ, μf=μf, Σf=Σf, S=S, idx=idx), numexpectationterm1 # ✅

    #--------------------------------
    # verify expectation of term (2)
    #--------------------------------

    numexpectationterm2 = @distributed (+) for i=1:numsamples

        term2(rand(qf), K) / numsamples

    end

    @show expectationterm2(μf=μf, Σf=Σf, K=K), numexpectationterm2 # ✅

    #--------------------------------
    # verify expectation of term (3)
    #--------------------------------

    numexpectationterm3 = @distributed (+) for i=1:numsamples

        term3(rand(qα), v²ₐ) / numsamples

    end

    @show expectationterm3(μα=μα, Σα=Σα, v²ₐ=v²ₐ), numexpectationterm3  # ✅

    #--------------------------------
    # verify expectation of term (4)
    #--------------------------------

    numexpectationterm4 = @distributed (+) for i=1:numsamples

        term4(rand(qβ), v²ᵦ) / numsamples

    end

    @show expectationterm4(μβ=μβ, Σβ=Σβ, v²ᵦ=v²ᵦ), numexpectationterm4

end

#------------------ Individual terms and their expectations ------------------#

##########
# TERM 1 #
##########

@everywhere function term1(y, A, f, b, S)

    N = length(y)

    -0.5*logdet(S) - 0.5*N*log(2π) - 0.5*(y - A*f - b)'*(S\(y - A*f - b))

end

@everywhere function term1_broken_down(y, A, f, b, S)

    N = length(y)

    subterm0 = -0.5*logdet(S) - 0.5*N*log(2π)

    subterm1 = -0.5*dot(y, S\y)

    subterm2 = -0.5*dot(y, -S\A*f) * 2

    subterm3 = -0.5*dot(y, -S\b) * 2

    subterm4 =  -0.5*dot(A*f,S\A*f)

    subterm5 = -0.5*dot(b, S\A*f) * 2

    subterm6 = -0.5*dot(b, S\b)

    (subterm0 + subterm1 + subterm2 + subterm3 + subterm4 + subterm5 + subterm6)

end

##########
# TERM 2 #
##########

@everywhere function term2(f, K)
    N = length(f)
    -0.5*N*log(2π) -0.5*logdet(K) - 0.5*f'*(K\f)
end

##########
# TERM 3 #
##########

@everywhere function term3(α, v²ₐ)
    F = length(α)
    - 0.5*F*log(2π) - 0.5*F*log(v²ₐ) - 0.5*dot(α,α)/v²ₐ
end

##########
# TERM 4 #
##########

@everywhere function term4(β, v²ᵦ)
    F = length(β)
    - 0.5*F*log(2π) - 0.5*F*log(v²ᵦ) - 0.5*dot(β,β)/v²ᵦ
end


text:
using CSV, DataFrames, DecisionTree, Distributions, Gurobi, JuMP, LightGraphs, Parameters, SparseArrays, Statistics, ArgParse
include("../oracles/portfolio_oracle.jl")
include("../solver/util.jl")
include("../solver/sgd.jl")
include("../solver/reformulation.jl")
include("../solver/random_forests_po.jl")
include("../solver/validation_set.jl")
include("../experiments/replication_functions.jl")

# Fixed parameter settings (these never change)
p_features = 5
num_assets = 50
num_trials = 50
n_test = 10000
num_factors = 4

num_lambda = 1
lambda_max = 10.0^(-6)
lambda_min_ratio = 1
holdout_percent = 0.25
different_validation_losses = false
data_type = :poly_kernel

# Fixed parameter sets (these are also the same for all experiments)
n_train_vec = [100; 1000]
n_sigmoid_polydegree_vec = [1; 4; 8; 16]
noise_multiplier_tau_vec = [1; 2]

# Set this based on expt_number (40 total)
rng_seed = 2223

# Run experiment and get results
# Note that Gurobi enviornments are set within this function call
expt_results = portfolio_multiple_replications(rng_seed, num_trials,
    num_assets, num_factors, n_train_vec, n_test, p_features,
    n_sigmoid_polydegree_vec, noise_multiplier_tau_vec;
    num_lambda = num_lambda, lambda_max = lambda_max, lambda_min_ratio = lambda_min_ratio, holdout_percent = holdout_percent,
    different_validation_losses = different_validation_losses, data_type = data_type)

csv_string = "portfolio_results.csv"
CSV.write(csv_string, expt_results)

max_stars_repo_name: paulgrigas/SmartPredictThenOptimize

max_issues_repo_name: paulgrigas/SmartPredictThenOptimize

max_forks_repo_name: paulgrigas/SmartPredictThenOptimize


text:
# Moran's I test of Global Spatial Autocorrelation
struct GlobalMoran <: AbstractGlobalSpatialAutocorrelation
    n::Int
    I::Float64
    EI::Float64
    p::Float64
    Iperms::Vector{Float64}
    Ipermsmean::Float64
    Ipermsstd::Float64
    z::Float64
end

"""
    moran(x, W)
Compute the  Moran's I test of global spatial autocorrelation.

# Optional Arguments
- `permutations=9999`: number of permutations for the randomization test.
- `rng=default_rng()`: random number generator for the randomization test.
"""
function moran(x::AbstractVector{T} where T<:Real, W::SpatialWeights; permutations::Int = 9999,
    rng::AbstractRNG = default_rng())::GlobalMoran

    n = length(x)
    z = x .- mean(x)

    S0 = sum(sum.(W.weights::Vector{Vector{Float64}}))

    # Auxiliar function to calculate Moran's I
    function moran_calc(z::AbstractVector, W::SpatialWeights)::Float64
        Wz = slag(W, z)
        return (sum(Wz .* z) / S0) / (sum(z .* z) / n)
    end

    # Morans's I
    I = moran_calc(z, W)
    EI = -1.0 / (n - 1)

    # Build permutations array
    Zi = Array{Float64}(undef, n, permutations)
    for i in 1:permutations
        Zi[:,i] = shuffle(rng, z)
    end

    # Calculate Moran's I for all the permutations
    Iperms = zeros(permutations)
    Threads.@threads for i in 1:permutations
        Iperms[i] = moran_calc(view(Zi, :, i), W)
    end

    larger = sum(Iperms .>= I)
    if (permutations - larger) < larger
        larger = permutations - larger
    end
    p = (larger + 1) / (permutations + 1)

    Ipermsstd = std(Iperms, corrected = false)
    Ipermsmean = mean(Iperms)
    z = (I - Ipermsmean)  / Ipermsstd

    return GlobalMoran(n, I, EI, p, Iperms, Ipermsmean, Ipermsstd, z)

end

score(x::GlobalMoran) = x.I;

scoreperms(x::GlobalMoran) = x.Iperms;

mean(x::GlobalMoran) = x.Ipermsmean;

expected(x::GlobalMoran) = x.EI;

std(x::GlobalMoran) = x.Ipermsstd;

zscore(x::GlobalMoran) = x.z;

pvalue(x::GlobalMoran) = x.p;

testname(x::GlobalMoran) = "Moran's I";


text:
using FunctionTabulations
using Unitful
using Test

@testset "warning" begin
    func_1d(x) = sin(x)

    # create interpolations and test them

    itp_1d_1 = create_tabulation_1D(
        func_1d,
        xmin = 0.0,
        xmax = 3.0,
        npoints = 100,
        x_scale = :linear,
        f_scale = :linear
    )

    func_1d(x) = sin(x)^2 + 3

    @test_logs (:warn,"The SHA for `func_1d` did not match the one of the stored tabulated function. Please check if the function definition has changed.") match_mode=:any create_tabulation_1D(
        func_1d,
        xmin = 0.0,
        xmax = 3.0,
        npoints = 100,
        x_scale = :linear,
        f_scale = :linear
    )

    rm("func_1d_data.jld2")

    func_2d(x, y) = sin(x) * sin(y)

    # create interpolations and test them

    itp_2d_1 = create_tabulation_2D(
        func_2d,
        xmin = 0.0,
        xmax = 1.0,
        ymin = 0.0,
        ymax = 2.0,
        npoints_x = 100,
        npoints_y = 200,
        x_scale = :linear,
        y_scale = :linear,
        f_scale = :linear
    )

    func_2d(x, y) = sin(x) * sin(y) * 3

    # create interpolations and test them

    @test_logs (:warn,"The SHA for `func_2d` did not match the one of the stored tabulated function. Please check if the function definition has changed.") match_mode=:any create_tabulation_2D(
        func_2d,
        xmin = 0.0,
        xmax = 1.0,
        ymin = 0.0,
        ymax = 2.0,
        npoints_x = 100,
        npoints_y = 200,
        x_scale = :linear,
        y_scale = :linear,
        f_scale = :linear
    )

    rm("func_2d_data.jld2")

    func_3d(x, y, z) = x * y + z

    # create interpolations and test them
    itp_3d_1 = create_tabulation_3D(
        func_3d,
        xmin = 0.0,
        xmax = 1.0,
        ymin = 0.0,
        ymax = 2.0,
        zmin = 0.0,
        zmax = 3.0,
        npoints_x = 100,
        npoints_y = 200,
        npoints_z = 200,
        x_scale = :linear,
        y_scale = :linear,
        z_scale = :linear,
        f_scale = :linear
    )

    func_3d(x, y, z) = x * y + z + 3

    # create interpolations and test them
    @test_logs (:warn,"The SHA for `func_3d` did not match the one of the stored tabulated function. Please check if the function definition has changed.") match_mode=:any create_tabulation_3D(
        func_3d,
        xmin = 0.0,
        xmax = 1.0,
        ymin = 0.0,
        ymax = 2.0,
        zmin = 0.0,
        zmax = 3.0,
        npoints_x = 100,
        npoints_y = 200,
        npoints_z = 200,
        x_scale = :linear,
        y_scale = :linear,
        z_scale = :linear,
        f_scale = :linear
    )

    rm("func_3d_data.jld2")
end

@testset "create_tabulation_1D_no_units" begin
    func_1d(x) = sin(x)

    # create interpolations and test them

    itp_1d_1 = create_tabulation_1D(
        func_1d,
        xmin = 0.0,
        xmax = 3.0,
        npoints = 100,
        x_scale = :linear,
        f_scale = :linear
    )

    itp_1d_2 = create_tabulation_1D(
        func_1d,
        jld_base_path = "interpolations",
        custom_name = "1d_2",
        xmin = 1e-1,
        xmax = 3.0,
        npoints = 100,
        interpolation_type = :cubic,
        x_scale = :log10,
        f_scale = :log10
    )

    itp_1d_3 = create_tabulation_1D(
        func_1d,
        jld_base_path = "interpolations",
        custom_name = "1d_3",
        collect(range(0,3,length=100)),
        x_scale = :linear,
        f_scale = :linear
    )

    itp_1d_4 = create_tabulation_1D(
        func_1d,
        10 .^ collect(range(log10(1e-1),log10(3),length=100)),
        jld_base_path = "interpolations",
        custom_name = "1d_4",
        x_scale = :log10,
        f_scale = :log10
    )

    @test isapprox(itp_1d_1(2.0), func_1d(2.0), rtol = 1e-3)
    @test isapprox(itp_1d_2(2.0), func_1d(2.0), rtol = 1e-3)
    @test isapprox(itp_1d_3(2.0), func_1d(2.0), rtol = 1e-3)
    @test isapprox(itp_1d_4(2.0), func_1d(2.0), rtol = 1e-3)

    @test isapprox(itp_1d_1(2.0), itp_1d_3(2.0), rtol = 1e-6)

    # load interpolations and test them

    itp_1d_1 = create_tabulation_1D(
        func_1d,
        xmin = 0.0,
        xmax = 3.0,
        npoints = 100,
        x_scale = :linear,
        f_scale = :linear
    )

    itp_1d_2 = create_tabulation_1D(
        func_1d,
        custom_name = "1d_2",
        jld_base_path = "interpolations",
        xmin = 1e-1,
        xmax = 3.0,
        npoints = 100,
        interpolation_type = :cubic,
        x_scale = :log10,
        f_scale = :log10
    )

    itp_1d_3 = create_tabulation_1D(
        func_1d,
        jld_base_path = "interpolations",
        custom_name = "1d_3",
        collect(range(0,3,length=100)),
        x_scale = :linear,
        f_scale = :linear
    )

    itp_1d_4 = create_tabulation_1D(
        func_1d,
        10 .^ collect(range(log10(1e-1),log10(3),length=100)),
        jld_base_path = "interpolations",
        custom_name = "1d_4",
        x_scale = :log10,
        f_scale = :log10
    )

    @test isapprox(itp_1d_1(2.0), func_1d(2.0), rtol = 1e-3)
    @test isapprox(itp_1d_2(2.0), func_1d(2.0), rtol = 1e-3)
    @test isapprox(itp_1d_3(2.0), func_1d(2.0), rtol = 1e-3)
    @test isapprox(itp_1d_4(2.0), func_1d(2.0), rtol = 1e-3)

    @test isapprox(itp_1d_1(2.0), itp_1d_3(2.0), rtol = 1e-6)

    rm("func_1d_data.jld2")
    rm("interpolations/1d_2_data.jld2")
    rm("interpolations/1d_3_data.jld2")
    rm("interpolations/1d_4_data.jld2")
    rm("interpolations")
end

@testset "create_tabulation_1D_with_units" begin
    func_1d(x) = x^2

    # create interpolations and test them

    itp_1d_1 = create_tabulation_1D(
        func_1d,
        custom_name = "1d_1",
        xmin = 0.0u"m",
        xmax = 3.0u"m",
        npoints = 100,
        x_scale = :linear,
        f_scale = :linear
    )

    itp_1d_2 = create_tabulation_1D(
        func_1d,
        custom_name = "1d_2",
        xmin = 1e-1u"m",
        xmax = 3.0u"m",
        npoints = 100,
        x_scale = :log10,
        f_scale = :log10
    )

    itp_1d_3 = create_tabulation_1D(
        func_1d,
        collect(range(0,3,length=100))*u"m",
        custom_name = "1d_3",
        x_scale = :linear,
        f_scale = :linear
    )

    itp_1d_4 = create_tabulation_1D(
        func_1d,
        10 .^ collect(range(log10(1e-1),log10(3),length=100))*u"m",
        custom_name = "1d_4",
        x_scale = :log10,
        f_scale = :log10
    )

    @test isapprox(itp_1d_1(2.0u"m"), func_1d(2.0u"m"), rtol = 1e-3)
    @test isapprox(itp_1d_2(2.0u"m"), func_1d(2.0u"m"), rtol = 1e-3)
    @test isapprox(itp_1d_3(2.0u"m"), func_1d(2.0u"m"), rtol = 1e-3)
    @test isapprox(itp_1d_4(2.0u"m"), func_1d(2.0u"m"), rtol = 1e-3)

    @test isapprox(itp_1d_1(2.0u"m"), itp_1d_3(2.0u"m"), rtol = 1e-6)

    # load interpolations and test them

    itp_1d_1 = create_tabulation_1D(
        func_1d,
        custom_name = "1d_1",
        xmin = 0.0u"m",
        xmax = 3.0u"m",
        npoints = 100,
        x_scale = :linear,
        f_scale = :linear
    )

    itp_1d_2 = create_tabulation_1D(
        func_1d,
        custom_name = "1d_2",
        xmin = 1e-1u"m",
        xmax = 3.0u"m",
        npoints = 100,
        x_scale = :log10,
        f_scale = :log10
    )

    itp_1d_3 = create_tabulation_1D(
        func_1d,
        collect(range(0,3,length=100))*u"m",
        custom_name = "1d_3",
        x_scale = :linear,
        f_scale = :linear
    )

    itp_1d_4 = create_tabulation_1D(
        func_1d,
        10 .^ collect(range(log10(1e-1),log10(3),length=100))*u"m",
        custom_name = "1d_4",
        x_scale = :log10,
        f_scale = :log10
    )

    @test isapprox(itp_1d_1(2.0u"m"), func_1d(2.0u"m"), rtol = 1e-3)
    @test isapprox(itp_1d_2(2.0u"m"), func_1d(2.0u"m"), rtol = 1e-3)
    @test isapprox(itp_1d_3(2.0u"m"), func_1d(2.0u"m"), rtol = 1e-3)
    @test isapprox(itp_1d_4(2.0u"m"), func_1d(2.0u"m"), rtol = 1e-3)

    @test isapprox(itp_1d_1(2.0u"m"), itp_1d_3(2.0u"m"), rtol = 1e-6)

    rm("1d_1_data.jld2")
    rm("1d_2_data.jld2")
    rm("1d_3_data.jld2")
    rm("1d_4_data.jld2")
end

@testset "create_tabulation_2D_no_units" begin
    func_2d(x, y) = sin(x) * sin(y)

    # create interpolations and test them

    itp_2d_1 = create_tabulation_2D(
        func_2d,
        xmin = 0.0,
        xmax = 1.0,
        ymin = 0.0,
        ymax = 2.0,
        npoints_x = 100,
        npoints_y = 200,
        x_scale = :linear,
        y_scale = :linear,
        f_scale = :linear
    )

    itp_2d_2 = create_tabulation_2D(
        func_2d,
        custom_name = "2d_2",
        jld_base_path = "interpolations",
        xmin = 1e-1,
        xmax = 1.0,
        ymin = 1e-1,
        ymax = 2.0,
        npoints_x = 100,
        npoints_y = 200,
        interpolation_type = :cubic,
        x_scale = :log10,
        y_scale = :log10,
        f_scale = :log10
    )

    itp_2d_3 = create_tabulation_2D(
        func_2d,
        collect(range(0,1,length=100)),
        collect(range(0,2,length=200)),
        custom_name = "2d_3",
        jld_base_path = "interpolations",
        x_scale = :linear,
        y_scale = :linear,
        f_scale = :linear
    )

    itp_2d_4 = create_tabulation_2D(
        func_2d,
        10 .^collect(range(log10(1e-1),log10(1),length=100)),
        10 .^collect(range(log10(1e-1),log10(2),length=200)),
        custom_name = "2d_4",
        jld_base_path = "interpolations",
        x_scale = :log10,
        y_scale = :log10,
        f_scale = :log10
    )

    @test isapprox(itp_2d_1(1.0, 1.3), func_2d(1.0, 1.3), rtol = 1e-3)
    @test isapprox(itp_2d_2(1.0, 1.3), func_2d(1.0, 1.3), rtol = 1e-3)
    @test isapprox(itp_2d_3(1.0, 1.3), func_2d(1.0, 1.3), rtol = 1e-3)
    @test isapprox(itp_2d_4(1.0, 1.3), func_2d(1.0, 1.3), rtol = 1e-3)

    @test isapprox(itp_2d_1(1.0, 1.3), itp_2d_3(1.0, 1.3), rtol = 1e-3)

    # load interpolations and test them

    itp_2d_1 = create_tabulation_2D(
        func_2d,
        xmin = 0.0,
        xmax = 1.0,
        ymin = 0.0,
        ymax = 2.0,
        npoints_x = 100,
        npoints_y = 200,
        x_scale = :linear,
        y_scale = :linear,
        f_scale = :linear
    )

    itp_2d_2 = create_tabulation_2D(
        func_2d,
        custom_name = "2d_2",
        jld_base_path = "interpolations",
        xmin = 1e-1,
        xmax = 1.0,
        ymin = 1e-1,
        ymax = 2.0,
        npoints_x = 100,
        npoints_y = 200,
        interpolation_type = :cubic,
        x_scale = :log10,
        y_scale = :log10,
        f_scale = :log10
    )

    itp_2d_3 = create_tabulation_2D(
        func_2d,
        collect(range(0,1,length=100)),
        collect(range(0,2,length=200)),
        custom_name = "2d_3",
        jld_base_path = "interpolations",
        x_scale = :linear,
        y_scale = :linear,
        f_scale = :linear
    )

    itp_2d_4 = create_tabulation_2D(
        func_2d,
        10 .^ collect(range(log10(1e-1),log10(1),length=100)),
        10 .^ collect(range(log10(1e-1),log10(2),length=200)),
        custom_name = "2d_4",
        jld_base_path = "interpolations",
        x_scale = :log10,
        y_scale = :log10,
        f_scale = :log10
    )

    @test isapprox(itp_2d_1(1.0, 1.3), func_2d(1.0, 1.3), rtol = 1e-3)
    @test isapprox(itp_2d_2(1.0, 1.3), func_2d(1.0, 1.3), rtol = 1e-3)
    @test isapprox(itp_2d_3(1.0, 1.3), func_2d(1.0, 1.3), rtol = 1e-3)
    @test isapprox(itp_2d_4(1.0, 1.3), func_2d(1.0, 1.3), rtol = 1e-3)

    @test isapprox(itp_2d_1(1.0, 1.3), itp_2d_3(1.0, 1.3), rtol = 1e-3)

    rm("func_2d_data.jld2")
    rm("interpolations/2d_2_data.jld2")
    rm("interpolations/2d_3_data.jld2")
    rm("interpolations/2d_4_data.jld2")
    rm("interpolations")
end

@testset "create_tabulation_2D_with_units" begin
    func_2d(x, y) = x^2 / y

    # create interpolations and test them

    itp_2d_1 = create_tabulation_2D(
        func_2d,
        custom_name = "2d_1",
        xmin = 0.0u"m",
        xmax = 1.0u"m",
        ymin = 0.0u"s",
        ymax = 2.0u"s",
        npoints_x = 100,
        npoints_y = 200,
        x_scale = :linear,
        y_scale = :linear,
        f_scale = :linear
    )

    itp_2d_2 = create_tabulation_2D(
        func_2d,
        custom_name = "2d_2",
        xmin = 1e-1u"m",
        xmax = 1.0u"m",
        ymin = 1e-1u"s",
        ymax = 2.0u"s",
        npoints_x = 100,
        npoints_y = 200,
        x_scale = :log10,
        y_scale = :log10,
        f_scale = :log10
    )

    itp_2d_3 = create_tabulation_2D(
        func_2d,
        collect(range(0,1,length=100))*u"m",
        collect(range(0,2,length=200))*u"s",
        custom_name = "2d_3",
        x_scale = :linear,
        y_scale = :linear,
        f_scale = :linear
    )

    itp_2d_4 = create_tabulation_2D(
        func_2d,
        custom_name = "2d_4",
        10 .^ collect(range(log10(1e-1),log10(1),length=100))*u"m",
        10 .^ collect(range(log10(1e-1),log10(2),length=200))*u"s",
        x_scale = :log10,
        y_scale = :log10,
        f_scale = :log10
    )
    @test isapprox(itp_2d_1(1.0u"m", 1.3u"s"), func_2d(1.0u"m", 1.3u"s"), rtol = 1e-3)
    @test isapprox(itp_2d_2(1.0u"m", 1.3u"s"), func_2d(1.0u"m", 1.3u"s"), rtol = 1e-3)
    @test isapprox(itp_2d_3(1.0u"m", 1.3u"s"), func_2d(1.0u"m", 1.3u"s"), rtol = 1e-3)
    @test isapprox(itp_2d_4(1.0u"m", 1.3u"s"), func_2d(1.0u"m", 1.3u"s"), rtol = 1e-3)

    @test isapprox(itp_2d_1(1.0u"m", 1.3u"s"), itp_2d_3(1.0u"m", 1.3u"s"), rtol = 1e-3)

    # load interpolations and test them

    itp_2d_1 = create_tabulation_2D(
        func_2d,
        custom_name = "2d_1",
        xmin = 0.0u"m",
        xmax = 1.0u"m",
        ymin = 0.0u"s",
        ymax = 2.0u"s",
        npoints_x = 100,
        npoints_y = 200,
        x_scale = :linear,
        y_scale = :linear,
        f_scale = :linear
    )

    itp_2d_2 = create_tabulation_2D(
        func_2d,
        custom_name = "2d_2",
        xmin = 1e-1u"m",
        xmax = 1.0u"m",
        ymin = 1e-1u"s",
        ymax = 2.0u"s",
        npoints_x = 100,
        npoints_y = 200,
        x_scale = :log10,
        y_scale = :log10,
        f_scale = :log10
    )
    itp_2d_3 = create_tabulation_2D(
        func_2d,
        collect(range(0,1,length=100))*u"m",
        collect(range(0,2,length=200))*u"s",
        custom_name = "2d_3",
        x_scale = :linear,
        y_scale = :linear,
        f_scale = :linear
    )

    itp_2d_4 = create_tabulation_2D(
        func_2d,
        custom_name = "2d_4",
        10 .^ collect(range(log10(1e-1),log10(1),length=100))*u"m",
        10 .^ collect(range(log10(1e-1),log10(2),length=200))*u"s",
        x_scale = :log10,
        y_scale = :log10,
        f_scale = :log10
    )
    @test isapprox(itp_2d_1(1.0u"m", 1.3u"s"), func_2d(1.0u"m", 1.3u"s"), rtol = 1e-3)
    @test isapprox(itp_2d_2(1.0u"m", 1.3u"s"), func_2d(1.0u"m", 1.3u"s"), rtol = 1e-3)
    @test isapprox(itp_2d_3(1.0u"m", 1.3u"s"), func_2d(1.0u"m", 1.3u"s"), rtol = 1e-3)
    @test isapprox(itp_2d_4(1.0u"m", 1.3u"s"), func_2d(1.0u"m", 1.3u"s"), rtol = 1e-3)

    @test isapprox(itp_2d_1(1.0u"m", 1.3u"s"), itp_2d_3(1.0u"m", 1.3u"s"), rtol = 1e-3)

    rm("2d_1_data.jld2")
    rm("2d_2_data.jld2")
    rm("2d_3_data.jld2")
    rm("2d_4_data.jld2")
end

@testset "create_tabulation_3D_no_units" begin
    func_3d(x, y, z) = x * y + z

    # create interpolations and test them
    itp_3d_1 = create_tabulation_3D(
        func_3d,
        xmin = 0.0,
        xmax = 1.0,
        ymin = 0.0,
        ymax = 2.0,
        zmin = 0.0,
        zmax = 3.0,
        npoints_x = 100,
        npoints_y = 200,
        npoints_z = 200,
        x_scale = :linear,
        y_scale = :linear,
        z_scale = :linear,
        f_scale = :linear
    )

    itp_3d_2 = create_tabulation_3D(
        func_3d,
        jld_base_path = "interpolations",
        custom_name = "3d_2",
        xmin = 1e-1,
        xmax = 1.0,
        ymin = 1e-1,
        ymax = 2.0,
        zmin = 1e-1,
        zmax = 3.0,
        interpolation_type = :cubic,
        npoints_x = 100,
        npoints_y = 200,
        npoints_z = 200,
        x_scale = :log10,
        y_scale = :log10,
        z_scale = :log10,
        f_scale = :log10
    )

    itp_3d_3 = create_tabulation_3D(
        func_3d,
        jld_base_path = "interpolations",
        custom_name = "3d_3",
        collect(range(0,1,length=100)),
        collect(range(0,2,length=200)),
        collect(range(0,3,length=200)),
        x_scale = :linear,
        y_scale = :linear,
        z_scale = :linear,
        f_scale = :linear
    )

    itp_3d_4 = create_tabulation_3D(
        func_3d,
        10 .^ collect(range(log10(1e-1),log10(1),length=100)),
        10 .^ collect(range(log10(1e-1),log10(2),length=200)),
        10 .^ collect(range(log10(1e-1),log10(3),length=200)),
        jld_base_path = "interpolations",
        custom_name = "3d_4",
        x_scale = :log10,
        y_scale = :log10,
        z_scale = :log10,
        f_scale = :log10
    )
    @test isapprox(itp_3d_1(1.0, 1.3, 2.5), func_3d(1.0, 1.3, 2.5), rtol = 1e-3)
    @test isapprox(itp_3d_2(1.0, 1.3, 2.5), func_3d(1.0, 1.3, 2.5), rtol = 1e-2)
    @test isapprox(itp_3d_3(1.0, 1.3, 2.5), func_3d(1.0, 1.3, 2.5), rtol = 1e-2)
    @test isapprox(itp_3d_4(1.0, 1.3, 2.5), func_3d(1.0, 1.3, 2.5), rtol = 1e-2)

    @test isapprox(itp_3d_1(1.0, 1.3, 2.5), itp_3d_3(1.0, 1.3, 2.5), rtol = 1e-2)

    # load interpolations and test them

    itp_3d_1 = create_tabulation_3D(
        func_3d,
        xmin = 0.0,
        xmax = 1.0,
        ymin = 0.0,
        ymax = 2.0,
        zmin = 0.0,
        zmax = 3.0,
        npoints_x = 100,
        npoints_y = 200,
        npoints_z = 200,
        x_scale = :linear,
        y_scale = :linear,
        z_scale = :linear,
        f_scale = :linear
    )

    itp_3d_2 = create_tabulation_3D(
        func_3d,
        jld_base_path = "interpolations",
        custom_name = "3d_2",
        xmin = 1e-1,
        xmax = 1.0,
        ymin = 1e-1,
        ymax = 2.0,
        zmin = 1e-1,
        zmax = 3.0,
        npoints_x = 100,
        npoints_y = 200,
        npoints_z = 200,
        interpolation_type = :cubic,
        x_scale = :log10,
        y_scale = :log10,
        z_scale = :log10,
        f_scale = :log10
    )
    itp_3d_3 = create_tabulation_3D(
        func_3d,
        jld_base_path = "interpolations",
        custom_name = "3d_3",
        collect(range(0,1,length=100)),
        collect(range(0,2,length=200)),
        collect(range(0,3,length=200)),
        x_scale = :linear,
        y_scale = :linear,
        z_scale = :linear,
        f_scale = :linear
    )

    itp_3d_4 = create_tabulation_3D(
        func_3d,
        10 .^ collect(range(log10(1e-1),log10(1),length=100)),
        10 .^ collect(range(log10(1e-1),log10(2),length=200)),
        10 .^ collect(range(log10(1e-1),log10(3),length=200)),
        jld_base_path = "interpolations",
        custom_name = "3d_4",
        x_scale = :log10,
        y_scale = :log10,
        z_scale = :log10,
        f_scale = :log10
    )
    @test isapprox(itp_3d_1(1.0, 1.3, 2.5), func_3d(1.0, 1.3, 2.5), rtol = 1e-3)
    @test isapprox(itp_3d_2(1.0, 1.3, 2.5), func_3d(1.0, 1.3, 2.5), rtol = 1e-2)
    @test isapprox(itp_3d_3(1.0, 1.3, 2.5), func_3d(1.0, 1.3, 2.5), rtol = 1e-2)
    @test isapprox(itp_3d_4(1.0, 1.3, 2.5), func_3d(1.0, 1.3, 2.5), rtol = 1e-2)

    @test isapprox(itp_3d_1(1.0, 1.3, 2.5), itp_3d_3(1.0, 1.3, 2.5), rtol = 1e-2)

    rm("func_3d_data.jld2")
    rm("interpolations/3d_2_data.jld2")
    rm("interpolations/3d_3_data.jld2")
    rm("interpolations/3d_4_data.jld2")
    rm("interpolations")
end

@testset "create_tabulation_3D_with_units" begin
    func_3d(x, y, z) = x * y + z

    # create interpolations and test them
    itp_3d_1 = create_tabulation_3D(
        func_3d,
        custom_name = "3d_1",
        xmin = 0.0u"m",
        xmax = 1.0u"m",
        ymin = 0.0u"s^-1",
        ymax = 2.0u"s^-1",
        zmin = 0.0u"m/s",
        zmax = 3.0u"m/s",
        npoints_x = 100,
        npoints_y = 200,
        npoints_z = 200,
        x_scale = :linear,
        y_scale = :linear,
        z_scale = :linear,
        f_scale = :linear
    )

    itp_3d_2 = create_tabulation_3D(
        func_3d,
        custom_name = "3d_2",
        xmin = 1e-1u"m",
        xmax = 1.0u"m",
        ymin = 1e-1u"s^-1",
        ymax = 2.0u"s^-1",
        zmin = 1e-1u"m/s",
        zmax = 3.0u"m/s",
        npoints_x = 100,
        npoints_y = 200,
        npoints_z = 200,
        x_scale = :log10,
        y_scale = :log10,
        z_scale = :log10,
        f_scale = :log10
    )

    itp_3d_3 = create_tabulation_3D(
        func_3d,
        collect(range(0,1,length=100))*u"m",
        collect(range(0,2,length=200))*u"s^-1",
        collect(range(0,3,length=200))*u"m/s",
        custom_name = "3d_3",
        x_scale = :linear,
        y_scale = :linear,
        z_scale = :linear,
        f_scale = :linear
    )
    itp_3d_4 = create_tabulation_3D(
        func_3d,
        10 .^ collect(range(log10(1e-1),log10(1),length=100))*u"m",
        10 .^ collect(range(log10(1e-1),log10(2),length=200))*u"s^-1",
        10 .^ collect(range(log10(1e-1),log10(3),length=200))*u"m/s",
        custom_name = "3d_4",
        x_scale = :log10,
        y_scale = :log10,
        z_scale = :log10,
        f_scale = :log10
    )
    @test isapprox(itp_3d_1(1.0u"m", 1.3u"s^-1", 2.5u"m/s"), func_3d(1.0u"m", 1.3u"s^-1", 2.5u"m/s"), rtol = 1e-3)
    @test isapprox(itp_3d_2(1.0u"m", 1.3u"s^-1", 2.5u"m/s"), func_3d(1.0u"m", 1.3u"s^-1", 2.5u"m/s"), rtol = 1e-2)
    @test isapprox(itp_3d_3(1.0u"m", 1.3u"s^-1", 2.5u"m/s"), func_3d(1.0u"m", 1.3u"s^-1", 2.5u"m/s"), rtol = 1e-2)
    @test isapprox(itp_3d_4(1.0u"m", 1.3u"s^-1", 2.5u"m/s"), func_3d(1.0u"m", 1.3u"s^-1", 2.5u"m/s"), rtol = 1e-2)

    @test isapprox(itp_3d_1(1.0u"m", 1.3u"s^-1", 2.5u"m/s"), itp_3d_3(1.0u"m", 1.3u"s^-1", 2.5u"m/s"), rtol = 1e-2)

    # load interpolations and test them

    itp_3d_1 = create_tabulation_3D(
        func_3d,
        custom_name = "3d_1",
        xmin = 0.0u"m",
        xmax = 1.0u"m",
        ymin = 0.0u"s^-1",
        ymax = 2.0u"s^-1",
        zmin = 0.0u"m/s",
        zmax = 3.0u"m/s",
        npoints_x = 100,
        npoints_y = 200,
        npoints_z = 200,
        x_scale = :linear,
        y_scale = :linear,
        z_scale = :linear,
        f_scale = :linear
    )

    itp_3d_2 = create_tabulation_3D(
        func_3d,
        custom_name = "3d_2",
        xmin = 1e-1u"m",
        xmax = 1.0u"m",
        ymin = 1e-1u"s^-1",
        ymax = 2.0u"s^-1",
        zmin = 1e-1u"m/s",
        zmax = 3.0u"m/s",
        npoints_x = 100,
        npoints_y = 200,
        npoints_z = 200,
        x_scale = :log10,
        y_scale = :log10,
        z_scale = :log10,
        f_scale = :log10
    )

    itp_3d_3 = create_tabulation_3D(
        func_3d,
        collect(range(0,1,length=100))*u"m",
        collect(range(0,2,length=200))*u"s^-1",
        collect(range(0,3,length=200))*u"m/s",
        custom_name = "3d_3",
        x_scale = :linear,
        y_scale = :linear,
        z_scale = :linear,
        f_scale = :linear
    )
    itp_3d_4 = create_tabulation_3D(
        func_3d,
        10 .^ collect(range(log10(1e-1),log10(1),length=100))*u"m",
        10 .^ collect(range(log10(1e-1),log10(2),length=200))*u"s^-1",
        10 .^ collect(range(log10(1e-1),log10(3),length=200))*u"m/s",
        custom_name = "3d_4",
        x_scale = :log10,
        y_scale = :log10,
        z_scale = :log10,
        f_scale = :log10
    )
    @test isapprox(itp_3d_1(1.0u"m", 1.3u"s^-1", 2.5u"m/s"), func_3d(1.0u"m", 1.3u"s^-1", 2.5u"m/s"), rtol = 1e-3)
    @test isapprox(itp_3d_2(1.0u"m", 1.3u"s^-1", 2.5u"m/s"), func_3d(1.0u"m", 1.3u"s^-1", 2.5u"m/s"), rtol = 1e-2)
    @test isapprox(itp_3d_3(1.0u"m", 1.3u"s^-1", 2.5u"m/s"), func_3d(1.0u"m", 1.3u"s^-1", 2.5u"m/s"), rtol = 1e-2)
    @test isapprox(itp_3d_4(1.0u"m", 1.3u"s^-1", 2.5u"m/s"), func_3d(1.0u"m", 1.3u"s^-1", 2.5u"m/s"), rtol = 1e-2)

    @test isapprox(itp_3d_1(1.0u"m", 1.3u"s^-1", 2.5u"m/s"), itp_3d_3(1.0u"m", 1.3u"s^-1", 2.5u"m/s"), rtol = 1e-2)

    rm("3d_1_data.jld2")
    rm("3d_2_data.jld2")
    rm("3d_3_data.jld2")
    rm("3d_4_data.jld2")
end

@testset "errors" begin
    func_1d(x) = sin(x)
    func_2d(x, y) = sin(x) * sin(y)
    func_3d(x, y, z) = x * y + z

    @test_throws ArgumentError FunctionTabulations.scaler(1.0, :nan)
    @test_throws ArgumentError FunctionTabulations.un_scaler(1.0, :nan)

    @test_throws ArgumentError create_tabulation_1D(
        func_1d,
        custom_name = "1d_1",
        xmin = 0.0,
        xmax = 3.0,
        npoints = 100,
        x_scale = :nan,
        f_scale = :linear
    )
    @test_throws ArgumentError create_tabulation_1D(
        func_1d,
        custom_name = "1d_1",
        xmin = 0.0,
        xmax = 3.0,
        npoints = 100,
        x_scale = :linear,
        f_scale = :nan
    )
    @test_throws ArgumentError create_tabulation_1D(
        func_1d,
        custom_name = "1d_1",
        xmin = 0.0,
        xmax = 3.0,
        npoints = 100,
        interpolation_type = :nan,
    )

    @test_throws ArgumentError create_tabulation_2D(
        func_2d,
        custom_name = "2d_1",
        xmin = 0.0,
        xmax = 1.0,
        ymin = 0.0,
        ymax = 2.0,
        npoints_x = 100,
        npoints_y = 200,
        x_scale = :nan,
        y_scale = :linear,
        f_scale = :linear
    )
    @test_throws ArgumentError create_tabulation_2D(
        func_2d,
        custom_name = "2d_1",
        xmin = 0.0,
        xmax = 1.0,
        ymin = 0.0,
        ymax = 2.0,
        npoints_x = 100,
        npoints_y = 200,
        x_scale = :linear,
        y_scale = :nan,
        f_scale = :linear
    )
    @test_throws ArgumentError create_tabulation_2D(
        func_2d,
        custom_name = "2d_1",
        xmin = 0.0,
        xmax = 1.0,
        ymin = 0.0,
        ymax = 2.0,
        npoints_x = 100,
        npoints_y = 200,
        x_scale = :linear,
        y_scale = :linear,
        f_scale = :nan
    )
    @test_throws ArgumentError create_tabulation_2D(
        func_2d,
        custom_name = "2d_1",
        xmin = 0.0,
        xmax = 1.0,
        ymin = 0.0,
        ymax = 2.0,
        npoints_x = 100,
        npoints_y = 200,
        interpolation_type = :nan,
    )

    @test_throws ArgumentError  create_tabulation_3D(
        func_3d,
        custom_name = "3d_1",
        xmin = 0.0,
        xmax = 1.0,
        ymin = 0.0,
        ymax = 2.0,
        zmin = 0.0,
        zmax = 3.0,
        npoints_x = 100,
        npoints_y = 200,
        npoints_z = 200,
        x_scale = :nan,
        y_scale = :linear,
        z_scale = :linear,
        f_scale = :linear
    )
    @test_throws ArgumentError  create_tabulation_3D(
        func_3d,
        custom_name = "3d_1",
        xmin = 0.0,
        xmax = 1.0,
        ymin = 0.0,
        ymax = 2.0,
        zmin = 0.0,
        zmax = 3.0,
        npoints_x = 100,
        npoints_y = 200,
        npoints_z = 200,
        x_scale = :linear,
        y_scale = :nan,
        z_scale = :linear,
        f_scale = :linear
    )
    @test_throws ArgumentError  create_tabulation_3D(
        func_3d,
        custom_name = "3d_1",
        xmin = 0.0,
        xmax = 1.0,
        ymin = 0.0,
        ymax = 2.0,
        zmin = 0.0,
        zmax = 3.0,
        npoints_x = 100,
        npoints_y = 200,
        npoints_z = 200,
        x_scale = :linear,
        y_scale = :linear,
        z_scale = :nan,
        f_scale = :linear
    )
    @test_throws ArgumentError  create_tabulation_3D(
        func_3d,
        custom_name = "3d_1",
        xmin = 0.0,
        xmax = 1.0,
        ymin = 0.0,
        ymax = 2.0,
        zmin = 0.0,
        zmax = 3.0,
        npoints_x = 100,
        npoints_y = 200,
        npoints_z = 200,
        x_scale = :linear,
        y_scale = :linear,
        z_scale = :linear,
        f_scale = :nan
    )
    @test_throws ArgumentError  create_tabulation_3D(
        func_3d,
        custom_name = "3d_1",
        xmin = 0.0,
        xmax = 1.0,
        ymin = 0.0,
        ymax = 2.0,
        zmin = 0.0,
        zmax = 3.0,
        npoints_x = 100,
        npoints_y = 200,
        npoints_z = 200,
        interpolation_type = :nan,
    )

    rm("1d_1_data.jld2")
    rm("2d_1_data.jld2")
    rm("3d_1_data.jld2")
end


text:
using ECOS
const ECOSOPT = with_optimizer(ECOS.Optimizer, verbose = 0)

@testset "ECOS conic Problems" begin
    @testset "ECOS Conic linear problems" begin
        list_of_soc_problems = [
            conic_linear1_test,
            conic_linear2_test,
            conic_linear3_test,
            conic_linear4_test
        ]
        test_strong_duality(list_of_soc_problems, ECOSOPT)
    end
    @testset "ECOS SOC problems" begin
        list_of_soc_problems = [
            soc1_test,
            soc2_test,
            soc3_test,
            soc4_test,
            soc5_test,
            soc6_test
        ]
        test_strong_duality(list_of_soc_problems, ECOSOPT)
    end
    @testset "ECOS RotatedSOC problems" begin
        list_of_rsoc_problems = [
            rsoc1_test,
            rsoc2_test,
            rsoc3_test,
            rsoc4_test
        ]
        test_strong_duality(list_of_rsoc_problems, ECOSOPT)
    end
end


text:
@doc raw"""
    Point2

An object `p` of the type [`Point2`](@ref) is a point in the two-dimensional
Euclidean plane ``Ε^2``.

# Example

The following declaration creates two points with Cartesian coordinates.

```julia-repl
julia> p, q = Point2(), Point2(1.0, 2.0)
(PointC2(0, 0), PointC2(1, 2))
```
"""
Point2

"""
    Point2(::Origin = ORIGIN)

Introduces a variable `p` with Cartesian coordinates ``(0, 0)``.
"""
Point2(::Origin = ORIGIN)

"""
    Point2(x::Real, y::Real)

Introduces a point `p` initialized to ``(x, y)``.
"""
Point2(x::Real, y::Real) = Point2(convert(FT, x), convert(FT, y))

@doc raw"""
    Point2(x::Real, y::Real, hw::Real = 1)

Introduces a point `p` initialized to `(hx/hw,hy/hw)`.

!!! info "Precondition"

    `hw` ``≠`` `0`
"""
Point2(x::Real, y::Real, hw::Real) =
    Point2(convert(RT, x), convert(RT, y), convert(RT, hw))

"""
    Point2(wp::WeightedPoint2)

Introduces a point from a weighted point.
"""
Point2(wp::WeightedPoint2)

"""
    ==(p::Point2, q::Point2)

Test for equality.

Two points are equal, iff their ``x`` and ``y`` coordinates are equal. The point
can be compared with [`ORIGIN`](@ref).
"""
==(p::Point2, q::Point2)

"""
    hx(p::Point2)

Returns the homogeneous ``x`` coordinate.
"""
hx(p::Point2)

"""
    hy(p::Point2)

Returns the homogeneous ``y`` coordinate.
"""
hy(p::Point2)

"""
    hw(p::Point2)

Returns the homogenizing coordinate.
"""
hw(p::Point2)

"""
    x(p::Point2)

Returns the Cartesian ``x`` coordinate, that is [`hx()`](@ref)/[`hw()`](@ref).
"""
x(p::Point2)

"""
    y(p::Point2)

Returns the Cartesian ``y`` coordinate, that is [`hy()`](@ref)/[`hw()`](@ref).
"""
y(p::Point2)

@doc raw"""
    homogeneous(p::Point2, i::Integer)

Returns the `i`ᵗʰ homogeneous coordinate of `p`.

!!! info "Precondition"

    ``0 ≤ i ≤ 2``.
"""
homogeneous(p::Point2, i::Integer)

@doc raw"""
    cartesian(p::Point2, i::Integer)

Returns the `i`ᵗʰ Cartesian coordinate of `p`.

!!! info "Precondition"

    ``0 ≤ i ≤ 1``.
"""
cartesian(p::Point2, i::Integer)

"""
    dimension(p::Point2)

Returns the dimension (the constant 2).
"""
dimension(p::Point2)

"""
    bbox(p::Point2)

Returns a bounding box containing `p`.
"""
bbox(p::Point2)

"""
    transform(p::Point2, t::AffTransformation2)

Returns the point obtained by applying `t` on `p`.
"""
transform(p::Point2, t::AffTransformation2)

"""
    <(p::Point2, q::Point2)

Return `true` iff `p` is lexicographically smaller than `q`, i.e. either if
`x(p) < x(q)` or if `x(p) == x(q)` and `y(p) < y(q)`.
"""
<(p::Point2, q::Point2)

"""
    >(p::Point2, q::Point2)

Return `true` iff `p` is lexicographically greater than `q`.
"""
>(p::Point2, q::Point2)

"""
    <=(p::Point2, q::Point2)

Return `true` iff `p` is lexicographically smaller or equal to `q`.
"""
<=(p::Point2, q::Point2)

"""
    >=(p::Point2, q::Point2)

Return `true` iff `p` is lexicographically greater or equal to `q`.
"""
>=(p::Point2, q::Point2)

"""
    -(p::Point2, q::Point2)

Returns the difference vector between `p` and `q`.

You can substitute [`ORIGIN`](@ref) for either `p` or `q`, but not for both.
"""
-(p::Point2, q::Point2)


text:
"""
    build_kernel_ex!(e, paramters)

Return an `Expr` to define a (tuple-argument) lambda whose body reflects the
structure of `e`.
"""
# function build_f_ex!(ex::Any, ds, index)::Tuple{Expr, Dict{Symbol, Vector{Symbol}}}
function build_f_ex!(srcs_used, ex::Any, index, primary)::Tuple{Expr, ArgsIndex}
    # ai is essentially a map from tokens (and hence sources) to sets of fields cited
    ai = ArgsIndex()
    body_ex = replace_symbols!(ai, ex, index, primary)
    # record which sources were used
    union!(srcs_used, [index[token] for token in keys(ai.index) ])
    args_ex = Expr(:tuple)
    _args = Vector{Symbol}(length(ai.index))
    for t in keys(ai.index)
        _args[1] = t
    end
    append!(args_ex.args, _args)
    return Expr(:->, args_ex, Expr(:block, body_ex)), ai
end