text:
using Flux
# set up data
# incomplete example pipeline with arbitrary cutoff points
...
X = imgs
Y = onehotbatch(labels, 0:4)
...
if checkmem() < 0.2
	@error("Memory free is less than 20%")
	callout("Stopping. Memory problem")
	#break
end
...
m = Chain(...)
function loss(...)
  ...
end
dataset = Flux.Data.DataLoader(X, Y, batchsize=32, shuffle=true)
...
# callback function
evalcb = () ->
begin
	@show(loss(X, Y)),
	if loss(X,Y) < 0.75
		callout("early stop")
		Flux.stop()
	end
	# cpu temperature
	cput = checkcpu()
  # in this example temperature of 70 C is the cutoff
	while cput > 70
		@warn("watching CPU... ($cput)")
    		callout("CPU is hot")
    		# rest the CPU for 60 seconds
		sleep(60)
		cput = checkcpu()
	end
end # callback

Flux.@epochs 7 Flux.train!(loss, params(m), dataset, opt, cb = throttle(evalcb,1))

callout("Training done")


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineScaleSetUpdateOSDisk <: SwaggerModel
    caching::Any # spec type: Union{ Nothing, Caching } # spec name: caching
    image::Any # spec type: Union{ Nothing, VirtualHardDisk } # spec name: image
    vhdContainers::Any # spec type: Union{ Nothing, Vector{String} } # spec name: vhdContainers
    managedDisk::Any # spec type: Union{ Nothing, VirtualMachineScaleSetManagedDiskParameters } # spec name: managedDisk

    function VirtualMachineScaleSetUpdateOSDisk(;caching=nothing, image=nothing, vhdContainers=nothing, managedDisk=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetUpdateOSDisk, Symbol("caching"), caching)
        setfield!(o, Symbol("caching"), caching)
        validate_property(VirtualMachineScaleSetUpdateOSDisk, Symbol("image"), image)
        setfield!(o, Symbol("image"), image)
        validate_property(VirtualMachineScaleSetUpdateOSDisk, Symbol("vhdContainers"), vhdContainers)
        setfield!(o, Symbol("vhdContainers"), vhdContainers)
        validate_property(VirtualMachineScaleSetUpdateOSDisk, Symbol("managedDisk"), managedDisk)
        setfield!(o, Symbol("managedDisk"), managedDisk)
        o
    end
end # type VirtualMachineScaleSetUpdateOSDisk

const _property_map_VirtualMachineScaleSetUpdateOSDisk = Dict{Symbol,Symbol}(Symbol("caching")=>Symbol("caching"), Symbol("image")=>Symbol("image"), Symbol("vhdContainers")=>Symbol("vhdContainers"), Symbol("managedDisk")=>Symbol("managedDisk"))
const _property_types_VirtualMachineScaleSetUpdateOSDisk = Dict{Symbol,String}(Symbol("caching")=>"Caching", Symbol("image")=>"VirtualHardDisk", Symbol("vhdContainers")=>"Vector{String}", Symbol("managedDisk")=>"VirtualMachineScaleSetManagedDiskParameters")
Base.propertynames(::Type{ VirtualMachineScaleSetUpdateOSDisk }) = collect(keys(_property_map_VirtualMachineScaleSetUpdateOSDisk))
Swagger.property_type(::Type{ VirtualMachineScaleSetUpdateOSDisk }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualMachineScaleSetUpdateOSDisk[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetUpdateOSDisk }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetUpdateOSDisk[property_name]

function check_required(o::VirtualMachineScaleSetUpdateOSDisk)
    true
end

function validate_property(::Type{ VirtualMachineScaleSetUpdateOSDisk }, name::Symbol, val)
end

ext: jl


text:
using .C2BasisDefinition: C1, C1Norm

function opnormbound(N::Type{<:C1}, M, B)
    n = size(M, 2)

    est = 0

    for i in 1:n
        v = zeros(n)
        v[i] = 1
        z = C1Norm(B, v)
        est+= C1Norm(B, M*v)/z
    end
    return est.hi
end

using ProgressMeter
function norms_of_powers_basis(B, m::Integer, Q::DiscretizedOperator, f::AbstractArray;
    normv0::Real=-1., #used as "missing" value
    normQ::Real=-1.,
    normE::Real=-1.,
    normEF::Real=-1.,
    normIEF::Real=-1.,
    normN::Real=-1.)

    @assert eltype(f) <: Interval
    T = typeof(zero(eltype(Q.L)).hi) # gets "Float64" from Q.L
    n = size(Q.L, 1)
    M = mid.(Q.L)
    R = radius.(Q.L)
    #δ = opnormbound(N, R)
    #γz = gamma(T, max_nonzeros_per_row(Q.L))
    #γn = gamma(T, n+3) # not n+2 like in the paper, because we wish to allow for f to be the result of rounding
    #ϵ = zero(T)
    midf = mid.(f)

    # TODO: correct here
    norm(v) = C1Norm(B, v)

    norms = zeros(m)

    #S = zeros((n, m))
    #k = length(B.p)
    factor = rescaling_factor(B)

    @showprogress for (v, norm_0) in AverageZero(B)

        #@info v
        v/= norm_0
        #@info v
        #@info norm(v)
        for i in 1:m

            w = M*v
            v = w - Q.e * (midf*w)[1]
            #@info v

            #@info "norm_$i" norm(v)*factor
           # @info infnormoffunction(B, v)
           # @info infnormofderivative(B, v)
            norms[i] = max(norm(v).hi*factor, norms[i])
        end
        #@info norms
    end

    return norms

end


text:
function integrate(tree::Tree)
    state = zero(tree.state)
    for leaf ∈ cells(tree, filter=active)
        state += leaf.state * volume(leaf)
    end
    return state
end

# function polytope(tree::Tree{N}) where T
#     poly = Vector(undef, N)
#     poly[1] = cell.position[1] .+ [-1., 1., 1., -1., -1.] / (2 << (level(cell)))
#     poly[2] = cell.position[2] .+ [-1., -1., 1., 1., -1.] / (2 << (level(cell)))
# end

@inline volume(cell::Tree{N}) where N = 1. / (1 << (N*level(cell)))
@inline centroid(cell::Tree) = cell.position

@inline area(face::Face{N}) where N = 1. / (1 << ((N-1)*level(face)))
function face_area(cell::Tree{N}, face::Face{N}) where N
    if !at_refinement(face) || level(cell) != level(face)
        return area(face)
    else
        return area(face) / 1 << (N-1)
    end
end

function centroid(cells::Tuple{AbstractTree{N}, AbstractTree{N}}, direction) where N
    if initialized(cells[1]) && initialized(cells[2])
        idx = level(cells[1]) > level(cells[2]) ? 1 : 2
    else
        idx = initialized(cells[1]) ? 1 : 2
    end

    position = copy(cells[idx].position)
    if idx == 1
        position[direction] += 1. / (2<<level(cells[idx]))
    else
        position[direction] -= 1. / (2<<level(cells[idx]))
    end
    return position
end
@inline centroid(face::Face) = centroid(face.cells, face.face_direction)

function cell_distance(face::Face)
    dist = 0.0
    if initialized(face.cells[1])
        dist += 1. / (2<<level(face.cells[1]))
    end
    if initialized(face.cells[2])
        dist += 1. / (2<<level(face.cells[2]))
    end
    return dist
end

function cell_volume(face::Face{N}) where N
    if !at_refinement(face)
        return volume(face.cells[1])
    else
        return 1.5 / (1 << (N*level(face)))
    end
end


text:
# see math/arith.jl for these
#abs(x::DoubleFloat{T}) where {T<:IEEEFloat} = abs_db_db(x)
#(-)(x::DoubleFloat{T}) where {T<:IEEEFloat} = neg_db_db(x)
#negabs(x::DoubleFloat{T}) where {T<:IEEEFloat} = negabs_db_db(x)

inv(x::DoubleFloat{T}) where {T<:IEEEFloat} = inv_db_db(x)

square(x::DoubleFloat{T}) where {T<:IEEEFloat} = square_db_db(x)
cube(x::DoubleFloat{T}) where {T<:IEEEFloat} = cube_db_db(x)

sqrt(x::DoubleFloat{T}) where {T<:IEEEFloat} = sqrt_db_db(x)
cbrt(x::DoubleFloat{T}) where {T<:IEEEFloat} = cbrt_db_db(x)

(+)(x::T, y::DoubleFloat{T}) where {T<:IEEEFloat} = add_fpdb_db(x, y)
(+)(x::DoubleFloat{T}, y::T) where {T<:IEEEFloat} = add_dbfp_db(x, y)
(+)(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:IEEEFloat} = add_dbdb_db(x, y)

(-)(x::T, y::DoubleFloat{T}) where {T<:IEEEFloat} = sub_fpdb_db(x, y)
(-)(x::DoubleFloat{T}, y::T) where {T<:IEEEFloat} = sub_dbfp_db(x, y)
(-)(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:IEEEFloat} = sub_dbdb_db(x, y)

(*)(x::T, y::DoubleFloat{T}) where {T<:IEEEFloat} = mul_fpdb_db(x, y)
(*)(x::DoubleFloat{T}, y::T) where {T<:IEEEFloat} = mul_dbfp_db(x, y)
(*)(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:IEEEFloat} = mul_dbdb_db(x, y)

(/)(x::T, y::DoubleFloat{T}) where {T<:IEEEFloat} = dvi_fpdb_db(x, y)
(/)(x::DoubleFloat{T}, y::T) where {T<:IEEEFloat} = dvi_dbfp_db(x, y)
(/)(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:IEEEFloat} = dvi_dbdb_db(x, y)

# Also provide arithmetic to add two floating point numbers in a higher precision
# without first converting them
"""
    add2(x::T, y::T) where {T<:IEEEFloat}

Add `x` and `y` as if they are `DoubleFloat{T}` numbers. This is more
efficient as first converting to `DoubleFloat{T}` and then adding.
"""
add2(x::T, y::T) where {T<:IEEEFloat} = add_fpfp_db(x, y)

"""
    ⊕(x::T, y::T) where {T<:IEEEFloat}

Add `x` and `y` as if they are `DoubleFloat{T}` numbers. This is more
efficient as first converting to `DoubleFloat{T}` and then adding.
This is the inline version of [`add2`](@ref) and written as ``\\oplus``.
"""
⊕(x::T, y::T) where {T<:IEEEFloat} = add2(x, y)

"""
    sub2(x::T, y::T) where {T<:IEEEFloat}

Subtract `x` and `y` as if they are `DoubleFloat{T}` numbers. This is more
efficient as first converting to `DoubleFloat{T}` and then subtracting.
"""
sub2(x::T, y::T) where {T<:IEEEFloat} = sub_fpfp_db(x, y)

"""
    ⊖(x::T, y::T) where {T<:IEEEFloat}

Subtract `x` and `y` as if they are `DoubleFloat{T}` numbers. This is more
efficient as first converting to `DoubleFloat{T}` and then subtracting.
This is the inline version of [`sub2`](@ref) and written as `\\ominus`.
"""
⊖(x::T, y::T) where {T<:IEEEFloat} = sub2(x, y)

"""
    mul2(x::T, y::T) where {T<:IEEEFloat}

Multiply `x` and `y` as if they are `DoubleFloat{T}` numbers. This is more
efficient as first converting to `DoubleFloat{T}` and then multiplying.
"""
mul2(x::T, y::T) where {T<:IEEEFloat} = mul_fpfp_db(x, y)

"""
    ⊗(x::T, y::T) where {T<:IEEEFloat}

Multiply `x` and `y` as if they are `DoubleFloat{T}` numbers. This is more
efficient as first converting to `DoubleFloat{T}` and then multiplying.
This is the inline version of [`mul2`](@ref) and written as `\\otimes`.
"""
⊗(x::T, y::T) where {T<:IEEEFloat} = mul2(x, y)

"""
    div(x::T, y::T) where {T<:IEEEFloat}

Divide `x` and `y` as if they are `DoubleFloat{T}` numbers. This is more
efficient as first converting to `DoubleFloat{T}` and then dividing.
"""
div2(x::T, y::T) where {T<:IEEEFloat} = dvi_fpfp_db(x, y)

"""
    ⊘(x::T, y::T) where {T<:IEEEFloat}

Divide `x` and `y` as if they are `DoubleFloat{T}` numbers. This is more
efficient as first converting to `DoubleFloat{T}` and then dividing.
This is the inline version of [`div2`](@ref) and written as ``\\oslash``.
"""
⊘(x::T, y::T) where {T<:IEEEFloat} = div2(x, y)


text:
"""
Image with dimension `N` and memory type `M`.
"""
abstract type Image{N,M<:Memory} <: LavaAbstraction end

struct ImageMetaData
    dims::Vector{Int}
    format::Vk.Format
    samples::Vk.SampleCountFlag
    mip_levels::Int
    layers::Int
    usage::Vk.ImageUsageFlag
end

dim(::Type{<:Image{N}}) where {N} = N
dim(im) = dim(typeof(im))

memory_type(::Type{<:Image{N,M}}) where {N,M} = M

Vk.bind_image_memory(image::Image, memory::Memory) = Vk.bind_image_memory(device(image), image, memory, offset(memory))

struct ImageBlock{N,M} <: Image{N,M}
    handle::Vk.Image
    dims::NTuple{N,Int}
    format::Vk.Format
    samples::Vk.SampleCountFlag
    mip_levels::Int
    layers::Int
    usage::Vk.ImageUsageFlag
    queue_family_indices::Vector{Int8}
    sharing_mode::Vk.SharingMode
    is_linear::Bool
    layout::Ref{Vk.ImageLayout}
    memory::Ref{M}
end

dims(image::ImageBlock) = image.dims
format(image::ImageBlock) = image.format
image_layout(image::ImageBlock) = image.layout[]
image(image::ImageBlock) = image
memory(image::ImageBlock) = image.memory[]
isallocated(image::ImageBlock) = isdefined(image.memory, 1)
samples(image::ImageBlock) = image.samples
usage(image::ImageBlock) = image.usage

function Vk.Extent3D(image::Image)
    d = dims(image)
    Vk.Extent3D(d..., ntuple(Returns(1), 3 - length(d))...)
end

Vk.Offset3D(image::ImageBlock) = Vk.Offset3D(0, 0, 0)

vk_handle_type(::Type{ImageBlock}) = Vk.Image

function ImageBlock(device, dims, format, usage;
                    queue_family_indices = queue_family_indices(device),
                    sharing_mode = Vk.SHARING_MODE_EXCLUSIVE,
                    memory_type = MemoryBlock,
                    is_linear = false,
                    preinitialized = false,
                    mip_levels = 1,
                    array_layers = 1,
                    samples = Vk.SAMPLE_COUNT_1_BIT)
    N = length(dims)
    initial_layout = preinitialized ? Vk.IMAGE_LAYOUT_PREINITIALIZED : Vk.IMAGE_LAYOUT_UNDEFINED
    extent_dims = ones(Int, 3)
    extent_dims[1:N] .= dims
    info = Vk.ImageCreateInfo(
        flag(ImageBlock{N,memory_type}),
        format,
        Vk.Extent3D(extent_dims...),
        mip_levels,
        array_layers,
        samples,
        is_linear ? Vk.IMAGE_TILING_LINEAR : Vk.IMAGE_TILING_OPTIMAL,
        usage,
        sharing_mode,
        queue_family_indices,
        initial_layout,
    )
    code = Vk.get_physical_device_image_format_properties(physical_device(device), info.format, info.image_type, info.tiling, info.usage)
    if iserror(code) && unwrap_error(code).code == Vk.ERROR_FORMAT_NOT_SUPPORTED
        error("Format $format not supported for images with tiling $(info.tiling) and usage $(info.usage)")
    end
    handle = unwrap(create(ImageBlock, device, info))
    ImageBlock{N,memory_type}(handle, dims, format, samples, mip_levels, array_layers, usage, queue_family_indices, sharing_mode, is_linear, Ref(initial_layout), Ref{memory_type}())
end

function Base.similar(image::ImageBlock; memory_domain = nothing, usage = image.usage, is_linear = image.is_linear)
    similar = ImageBlock(device(image), dims(image), format(image), usage; image.queue_family_indices, image.sharing_mode, memory_type = memory_type(image), is_linear, image.mip_levels, array_layers = image.layers, image.samples)
    if isallocated(image)
        memory_domain = @something(memory_domain, memory(image).domain)
        allocate!(similar, memory_domain)
    end
    similar
end

function bind!(image::ImageBlock, memory::Memory)
    image.memory[] = memory
    unwrap(Vk.bind_image_memory(image, memory))
    image
end

"""
Allocate a `MemoryBlock` and bind it to provided image.
"""
function allocate!(image::ImageBlock, domain::MemoryDomain)
    _device = device(image)
    reqs = Vk.get_image_memory_requirements(_device, image)
    memory = MemoryBlock(_device, reqs.size, reqs.memory_type_bits, domain)
    bind!(image, memory)
end

"""
View of a resource, such as an image or buffer.
"""
abstract type View{O<:LavaAbstraction} <: LavaAbstraction end

struct ImageView{I<:Image} <: View{I}
    handle::Vk.ImageView
    image::I
    format::Vk.Format
    aspect::Vk.ImageAspectFlag
    mip_range::UnitRange
    layer_range::UnitRange
end

image(view::ImageView) = view.image

vk_handle_type(T::Type{<:ImageView}) = Vk.ImageView

image_type(::Type{<:ImageView{I}}) where {I} = I
image_type(view::ImageView) = image_type(typeof(view))

dim(T::Type{<:ImageView}) = dim(image_type(T))
memory_type(T::Type{<:ImageView}) = memory_type(image_type(T))
format(view::ImageView) = view.format

@forward ImageView.image samples, dims, image_layout, usage, Vk.Offset3D, Vk.Extent3D, isallocated, image

function flag(T::Type{<:Image})
    @match dim(T) begin
        1 => Vk.IMAGE_TYPE_1D
        2 => Vk.IMAGE_TYPE_2D
        3 => Vk.IMAGE_TYPE_3D
    end
end

function flag(T::Type{<:ImageView})
    @match dim(T) begin
        1 => Vk.IMAGE_VIEW_TYPE_1D
        2 => Vk.IMAGE_VIEW_TYPE_2D
        3 => Vk.IMAGE_VIEW_TYPE_3D
    end
end

View(image::Image, args...; kwargs...) = ImageView(image, args...; kwargs...)

const DEFAULT_ASPECT = Vk.IMAGE_ASPECT_COLOR_BIT

default_mip_range(image) = 0:image.mip_levels
default_layer_range(image) = 1:image.layers

function ImageView(
        image::I;
        view_type = flag(ImageView{I}),
        format = format(image),
        component_mapping = Vk.ComponentMapping(
            Vk.COMPONENT_SWIZZLE_IDENTITY,
            Vk.COMPONENT_SWIZZLE_IDENTITY,
            Vk.COMPONENT_SWIZZLE_IDENTITY,
            Vk.COMPONENT_SWIZZLE_IDENTITY
        ),
        aspect = DEFAULT_ASPECT,
        mip_range = default_mip_range(image),
        layer_range = default_layer_range(image),
    ) where {I<:Image}

    info = Vk.ImageViewCreateInfo(
        convert(Vk.Image, image),
        view_type,
        format,
        component_mapping,
        subresource_range(aspect, mip_range, layer_range),
    )
    handle = unwrap(create(ImageView, device(image), info))
    ImageView{I}(handle, image, format, aspect, mip_range, layer_range)
end

subresource_range(aspect::Vk.ImageAspectFlag, mip_range::UnitRange, layer_range::UnitRange) = Vk.ImageSubresourceRange(aspect, mip_range.start, mip_range.stop - mip_range.start, layer_range.start - 1, 1 + layer_range.stop - layer_range.start)
subresource_range(view::ImageView) = subresource_range(view.aspect, view.mip_range, view.layer_range)
subresource_range(image::Image) = subresource_range(DEFAULT_ASPECT, default_mip_range(image), default_layer_range(image))

subresource_layers(aspect::Vk.ImageAspectFlag = DEFAULT_ASPECT, mip_range::Integer = 0, layer_range::UnitRange = DEFAULT_LAYER_RANGE) = Vk.ImageSubresourceLayers(aspect, mip_range, layer_range.start - 1, 1 + layer_range.stop - layer_range.start)
subresource_layers(view::ImageView) = subresource_layers(view.aspect, first(view.mip_range), view.layer_range)
subresource_layers(image::Image) = subresource_layers(DEFAULT_ASPECT, first(default_mip_range(image)), default_layer_range(image))

"""
Opaque image that comes from the Window System Integration (WSI) as returned by `Vk.get_swapchain_images_khr`.
"""
struct ImageWSI <: Image{2,OpaqueMemory}
    handle::Vk.Image
end

function transfer(src::Union{<:Image, <:ImageView}, dst::Union{<:Image, <:ImageView}, device; command_buffer = request_command_buffer(device, Vk.QUEUE_TRANSFER_BIT), signal_fence = true, semaphore = nothing, free_src = false)

    @assert dims(src) == dims(dst)

    if image_layout(src) ≠ Vk.IMAGE_LAYOUT_TRANSFER_SRC_OPTIMAL
        transition_layout(command_buffer, device, src, Vk.IMAGE_LAYOUT_TRANSFER_SRC_OPTIMAL)
    end
    if image_layout(dst) ≠ Vk.IMAGE_LAYOUT_TRANSFER_DST_OPTIMAL
        transition_layout(command_buffer, device, dst, Vk.IMAGE_LAYOUT_TRANSFER_DST_OPTIMAL)
    end

    Vk.cmd_copy_image(command_buffer,
        image(src), image_layout(src),
        image(dst), image_layout(dst),
        [Vk.ImageCopy(subresource_layers(src), Vk.Offset3D(src), subresource_layers(dst), Vk.Offset3D(dst), Vk.Extent3D(src))]
    )
    signal_semaphores = []
    !isnothing(semaphore) && push!(signal_semaphores, semaphore)
    info = Vk.SubmitInfo2KHR([], [Vk.CommandBufferSubmitInfoKHR(command_buffer)], signal_semaphores)
    if free_src
        submit(device, command_buffer.queue_family_index, info; signal_fence, semaphore, free_after_completion = [Ref(src)])
    else
        submit(device, command_buffer.queue_family_index, info; signal_fence, semaphore, release_after_completion = [Ref(src)])
    end
end

function transition_layout_info(view_or_image::Union{<:Image, <:ImageView}, new_layout)
    Vk.ImageMemoryBarrier2KHR(image_layout(view_or_image), new_layout, 0, 0, handle(image(view_or_image)), subresource_range(view_or_image))
end

function transition_layout(command_buffer::CommandBuffer, device, view_or_image::Union{<:Image, <:ImageView}, new_layout)
    Vk.cmd_pipeline_barrier_2_khr(command_buffer,
        Vk.DependencyInfoKHR([], [], [transition_layout_info(view_or_image, new_layout)]),
    )
    image(view_or_image).layout[] = new_layout
end

function Base.collect(@nospecialize(T), image::ImageBlock, device)
    if image.is_linear && Vk.MEMORY_PROPERTY_HOST_VISIBLE_BIT in properties(memory(image))
        isbitstype(T) || error("Image type is not an `isbits` type.")
        bytes = collect(memory(image), prod(dims(image)) * sizeof(T), device)
        data = reinterpret(T, bytes)
        reshape(data, dims(image))
    else
        usage = Vk.IMAGE_USAGE_TRANSFER_DST_BIT
        dst = ImageBlock(device, dims(image), format(image), usage; is_linear = true)
        allocate!(dst, MEMORY_DOMAIN_HOST)
        wait(transfer(image, dst, device))
        collect(T, dst, device)
    end
end

function transfer(data::AbstractArray, image::Image, device; kwargs...)
    b = buffer(device, data; usage = Vk.BUFFER_USAGE_TRANSFER_SRC_BIT)
    transfer(b, image, device; kwargs...)
end

function transfer(buffer::Buffer, view_or_image::Union{<:Image, <:ImageView}, device; command_buffer = request_command_buffer(device, Vk.QUEUE_TRANSFER_BIT), free_src = false)
    transition_layout(command_buffer, device, view_or_image, Vk.IMAGE_LAYOUT_TRANSFER_DST_OPTIMAL)
    Vk.cmd_copy_buffer_to_image(command_buffer, buffer, image(view_or_image), Vk.IMAGE_LAYOUT_TRANSFER_DST_OPTIMAL,
        [Vk.BufferImageCopy(offset(buffer), dims(view_or_image)..., subresource_layers(view_or_image), Vk.Offset3D(view_or_image), Vk.Extent3D(view_or_image))])
    info = Vk.SubmitInfo2KHR([], [Vk.CommandBufferSubmitInfoKHR(command_buffer)], [])
    release_after_completion = Ref[Ref(view_or_image)]
    free_after_completion = Ref[]
    push!(free_src ? free_after_completion : release_after_completion, Ref(buffer))
    submit(device, command_buffer.queue_family_index, info; signal_fence = true, free_after_completion, release_after_completion)
end

function image(device, data::AbstractArray, format::Vk.Format; memory_domain = MEMORY_DOMAIN_DEVICE, optimal_tiling = true, usage = Vk.IMAGE_USAGE_SAMPLED_BIT, kwargs...)
    upload_usage = usage | Vk.IMAGE_USAGE_TRANSFER_DST_BIT
    optimal_tiling && (upload_usage |= Vk.IMAGE_USAGE_TRANSFER_SRC_BIT)
    img = ImageBlock(device, size(data), format, upload_usage; is_linear = !optimal_tiling)
    allocate!(img, memory_domain)
    wait(transfer(data, img, device; free_src = true))
    !optimal_tiling && return img

    dst = similar(img; is_linear = false, usage = usage | Vk.IMAGE_USAGE_TRANSFER_DST_BIT)
    wait(transfer(img, dst, device; free_src = true))
    dst
end


text:
facts("Reservoir data structures: rivers") do
  const name1 = "River_name 1"
  const name2 = "River_name*2"
  const a = 1.
  const a_discharge = a * 1.m^3/s
  const b = 4. # Value chosen for DivertedRiver, the function getMaximumAllowableFlow
  const b_discharge = b * 1.m^3/s
  const sc1_begin_year = 1999
  const sc1_n_weeks = 5
  const sc1_n_scenarios = 25
  const sc1 = gen_tses(sc1_n_weeks, sc1_begin_year, sc1_n_scenarios)
  const sc2_begin_year = 1804
  const sc2_n_weeks = 5
  const sc2_n_scenarios = 25
  const sc2 = gen_tses(sc2_n_weeks, sc2_begin_year, sc2_n_scenarios)
  const sc_wrong = [TimeArray(timestamp(sc1[1]), values(sc1[1]), colnames(sc1[1]), week)]
  const sc_nan = [TimeArray(timestamp(sc1[1]), values(sc1[1]) * NaN, colnames(sc1[1]), week)]
  const sc_empty = similar(sc1, 0)

  function naturalRiverTestShortVersions(r::NaturalRiver)
    @fact name(r) --> getName(r)
    @fact scenarios(r) --> getScenarios(r)
    @fact scenario(r, 1) --> getScenario(r, 1)
    @fact scenario(r, 1, 1) --> getScenario(r, 1, 1)
    @fact period(r) --> getPeriod(r)
    @fact countScenarios(r) --> ReservoirManagement.countScenariosWithConsistency(r)
    @fact countTimeSteps(r) --> ReservoirManagement.countTimeStepsWithConsistency(r)
    @fact getPeriod(r) --> ReservoirManagement.getPeriodWithConsistency(r)
    @fact typeof(getPeriod(r)) <: Base.Dates.Period --> true
  end

  context("NaturalRiver") do
    r = NaturalRiver(name=name1, scenarios=sc1)
    @fact isempty(r) --> false
    @fact getName(r) --> name1
    @fact getScenarios(r) --> sc1
    @fact countScenarios(r) --> sc1_n_scenarios
    @fact countTimeSteps(r) --> sc1_n_weeks
    @fact getPeriod(r) --> Week(1)
    @fact getScenario(r, 1) --> values(sc1[1])
    @fact getScenario(r, 1, 1) --> values(sc1[1])[1]
    naturalRiverTestShortVersions(r)

    r2 = copy(r, name=name2)
    @fact isempty(r2) --> false
    @fact getName(r2) --> name2
    @fact getScenarios(r2) --> sc1
    @fact countScenarios(r2) --> sc1_n_scenarios
    @fact countTimeSteps(r2) --> sc1_n_weeks
    @fact getPeriod(r2) --> Week(1)
    @fact getScenario(r2, 1) --> values(sc1[1])
    @fact getScenario(r2, 1, 1) --> values(sc1[1])[1]
    naturalRiverTestShortVersions(r2)

    r3 = copy(r2, scenarios=sc2)
    @fact isempty(r3) --> false
    @fact getName(r3) --> name2
    @fact getScenarios(r3) --> sc2
    @fact countScenarios(r3) --> sc2_n_scenarios
    @fact countTimeSteps(r3) --> sc2_n_weeks
    @fact getPeriod(r3) --> Week(1)
    @fact getScenario(r3, 1) --> values(sc2[1])
    @fact getScenario(r3, 1, 1) --> values(sc2[1])[1]
    naturalRiverTestShortVersions(r3)

    r4 = copy(r, name=name2, scenarios=sc2)
    @fact r3 --> exactly(r4)

    @fact_throws NaturalRiver(name=name1, scenarios=sc_wrong)
    @fact_throws NaturalRiver(name=name1, scenarios=sc_nan)
    @fact_throws NaturalRiver(name=name1, scenarios=sc_empty)
  end

  function divertedRiverTestShortVersions(r::DivertedRiver)
    @fact name(r) --> getName(r)
    @fact scenarios(r) --> getScenarios(r)
    @fact scenario(r, 1) --> getScenario(r, 1)
    @fact scenario(r, 1, 1) --> getScenario(r, 1, 1)
    @fact period(r) --> getPeriod(r)
    @fact countScenarios(r) --> ReservoirManagement.countScenariosWithConsistency(r)
    @fact countTimeSteps(r) --> ReservoirManagement.countTimeStepsWithConsistency(r)
    @fact getPeriod(r) --> ReservoirManagement.getPeriodWithConsistency(r)
    @fact typeof(getPeriod(r)) <: Base.Dates.Period --> true
    @fact environmentalFlow(r) --> getEnvironmentalFlow(r)
    @fact maxFlow(r) --> getMaximumFlow(r)
    @fact environmentalFlow(r, Week(1)) --> getEnvironmentalFlow(r, Week(1))
    @fact maxFlow(r, Week(1)) --> getMaximumFlow(r, Week(1))
    @fact maxAllowableFlow(r, 1, 1) --> getMaximumAllowableFlow(r, 1, 1)
  end

  context("DivertedRiver") do
    r = DivertedRiver(name=name1, scenarios=sc1, environmental_flow=a_discharge, maximum_flow=b_discharge)
    @fact isempty(r) --> false
    @fact getName(r) --> name1
    @fact getScenarios(r) --> sc1
    @fact countScenarios(r) --> sc1_n_scenarios
    @fact countTimeSteps(r) --> sc1_n_weeks
    @fact getPeriod(r) --> Week(1)
    @fact getScenario(r, 1) --> values(sc1[1])
    @fact getScenario(r, 1, 1) --> values(sc1[1])[1]
    @fact getEnvironmentalFlow(r) --> ReservoirManagement._from_unitful(a_discharge)
    @fact getMaximumFlow(r) --> ReservoirManagement._from_unitful(b_discharge)
    @fact getEnvironmentalFlow(r, Week(1)) --> ReservoirManagement._from_unitful(a_discharge) * 604800. # 10^6 m^3/s -> 10^6 m^3/week
    @fact getMaximumFlow(r, Week(1)) --> ReservoirManagement._from_unitful(b_discharge) * 604800. # 10^6 m^3/s -> 10^6 m^3/week
    @fact maxAllowableFlow(r, 1, 1) --> roughly(maximum([0., minimum([getMaximumFlow(r, Week(1)), getScenario(r, 1, 1) - getEnvironmentalFlow(r, Week(1))])]))
    divertedRiverTestShortVersions(r)

    r2 = copy(r, name=name2)
    @fact isempty(r2) --> false
    @fact getName(r2) --> name2
    @fact getScenarios(r2) --> sc1
    @fact countScenarios(r2) --> sc1_n_scenarios
    @fact countTimeSteps(r2) --> sc1_n_weeks
    @fact getPeriod(r2) --> Week(1)
    @fact getScenario(r2, 1) --> values(sc1[1])
    @fact getScenario(r2, 1, 1) --> values(sc1[1])[1]
    @fact getEnvironmentalFlow(r2) --> ReservoirManagement._from_unitful(a_discharge)
    @fact getMaximumFlow(r2) --> ReservoirManagement._from_unitful(b_discharge)
    @fact maxAllowableFlow(r2, 1, 1) --> roughly(maximum([0., minimum([getMaximumFlow(r2, Week(1)), getScenario(r2, 1, 1) - getEnvironmentalFlow(r2, Week(1))])]))
    divertedRiverTestShortVersions(r2)

    r3 = copy(r2, scenarios=sc2)
    @fact isempty(r3) --> false
    @fact getName(r3) --> name2
    @fact getScenarios(r3) --> sc2
    @fact countScenarios(r3) --> sc2_n_scenarios
    @fact countTimeSteps(r3) --> sc2_n_weeks
    @fact getPeriod(r3) --> Week(1)
    @fact getScenario(r3, 1) --> values(sc2[1])
    @fact getScenario(r3, 1, 1) --> values(sc2[1])[1]
    @fact getEnvironmentalFlow(r3) --> ReservoirManagement._from_unitful(a_discharge)
    @fact getMaximumFlow(r3) --> ReservoirManagement._from_unitful(b_discharge)
    @fact maxAllowableFlow(r3, 1, 1) --> roughly(maximum([0., minimum([getMaximumFlow(r3, Week(1)), getScenario(r3, 1, 1) - getEnvironmentalFlow(r3, Week(1))])]))
    divertedRiverTestShortVersions(r3)

    r4 = copy(r3, environmental_flow=b_discharge)
    @fact isempty(r4) --> false
    @fact getName(r4) --> name2
    @fact getScenarios(r4) --> sc2
    @fact countScenarios(r4) --> sc2_n_scenarios
    @fact countTimeSteps(r4) --> sc2_n_weeks
    @fact getPeriod(r4) --> Week(1)
    @fact getScenario(r4, 1) --> values(sc2[1])
    @fact getScenario(r4, 1, 1) --> values(sc2[1])[1]
    @fact getEnvironmentalFlow(r4) --> ReservoirManagement._from_unitful(b_discharge)
    @fact getMaximumFlow(r4) --> ReservoirManagement._from_unitful(b_discharge)
    @fact maxAllowableFlow(r4, 1, 1) --> roughly(maximum([0., minimum([getMaximumFlow(r4, Week(1)), getScenario(r4, 1, 1) - getEnvironmentalFlow(r4, Week(1))])]))
    divertedRiverTestShortVersions(r4)

    r5 = copy(r4, maximum_flow=a_discharge)
    @fact isempty(r5) --> false
    @fact getName(r5) --> name2
    @fact getScenarios(r5) --> sc2
    @fact countScenarios(r5) --> sc2_n_scenarios
    @fact countTimeSteps(r5) --> sc2_n_weeks
    @fact getPeriod(r5) --> Week(1)
    @fact getScenario(r5, 1) --> values(sc2[1])
    @fact getScenario(r5, 1, 1) --> values(sc2[1])[1]
    @fact getEnvironmentalFlow(r5) --> ReservoirManagement._from_unitful(b_discharge)
    @fact getMaximumFlow(r5) --> ReservoirManagement._from_unitful(a_discharge)
    @fact maxAllowableFlow(r5, 1, 1) --> roughly(maximum([0., minimum([getMaximumFlow(r5, Week(1)), getScenario(r5, 1, 1) - getEnvironmentalFlow(r5, Week(1))])]))
    divertedRiverTestShortVersions(r5)

    r6 = copy(r, name=name2, scenarios=sc2, environmental_flow=b_discharge, maximum_flow=a_discharge)
    @fact r6 --> exactly(r5)

    @fact_throws DivertedRiver(name=name1, scenarios=sc_wrong, environmental_flow=a_discharge, maximum_flow=b_discharge)
    @fact_throws DivertedRiver(name=name1, scenarios=sc_nan, environmental_flow=a_discharge, maximum_flow=b_discharge)
    @fact_throws DivertedRiver(name=name1, scenarios=sc_empty, environmental_flow=a_discharge, maximum_flow=b_discharge)
  end
end


text:
function post(𝒜::ASB07_decomposed,
              𝑃::InitialValueProblem{<:AbstractContinuousSystem},
              𝑂::Options)
    # =================================
    # Initialization and discretization
    # =================================

    𝑂 = merge(𝒜.options.defaults, 𝑂, 𝒜.options.specified)
    δ, T = 𝑂[:δ], 𝑂[:T]
    N = round(Int, T / δ)
    n = 𝑂[:n]
    partition = 𝑂[:partition]
    blocks = 𝑂[:blocks]
    max_order = 𝑂[:max_order]

    # compute and unrwap discretized system
    𝑃_discrete = discretize(𝑃, δ; algorithm="interval_matrix",
                            order=𝑂[:order_discretization],
                            set_operations=𝑂[:set_operations_discretization])
    Ω0, Φ = 𝑃_discrete.x0, 𝑃_discrete.s.A

    # decompose initial states
    # TODO this should use the concrete linear_map projection (see LazySets#1726)
    Ωhat0 = array(decompose(Ω0, partition, 𝑂[:block_options_init]))

    # ====================
    # Flowpipe computation
    # ====================

    # preallocate output
    Rsets = Vector{ReachSet{CartesianProductArray{Float64, LazySet{Float64}}}}(undef, N)

    info("Reachable States Computation...")
    @timing begin
    if inputdim(𝑃_discrete) == 0
        U = nothing
    else
        U = inputset(𝑃_discrete)
    end
    reach_ASB07_decomposed!(Rsets, Ωhat0, U, Φ, N, δ, max_order, n, partition,
                            blocks)
    end # timing

    Rsol = ReachSolution(Flowpipe(Rsets), 𝑂)

    # ==========
    # Projection
    # ==========

    if 𝑂[:project_reachset]
        info("Projection...")
        Rsol = @timing project(Rsol)
    end

    return Rsol
end


text:
"""
    abstract type Transformable
This is a bit of a weird name, but all scenes and plots are transformable,
so that's what they all have in common. This might be better expressed as traits.
"""
abstract type Transformable end

abstract type AbstractPlot{Typ} <: Transformable end
abstract type AbstractScene <: Transformable end
abstract type ScenePlot{Typ} <: AbstractPlot{Typ} end
abstract type AbstractScreen <: AbstractDisplay end

const SceneLike = Union{AbstractScene, ScenePlot}

"""
Main structure for holding attributes, for theming plots etc!
Will turn all values into nodes, so that they can be updated.
"""
struct Attributes
    attributes::Dict{Symbol, Observable}
end

struct Combined{Typ, T} <: ScenePlot{Typ}
    parent::SceneLike
    transformation::Transformable
    attributes::Attributes
    input_args::Tuple
    converted::Tuple
    plots::Vector{AbstractPlot}
end

function Base.show(io::IO, plot::Combined)
    print(io, typeof(plot))
end

Base.parent(x::AbstractPlot) = x.parent

struct Key{K} end
macro key_str(arg)
    :(Key{$(QuoteNode(Symbol(arg)))})
end
Base.broadcastable(x::Key) = (x,)

"""
Type to indicate that an attribute will get calculated automatically
"""
struct Automatic end

"""
Singleton instance to indicate that an attribute will get calculated automatically
"""
const automatic = Automatic()

abstract type Unit{T} <: Number end

"""
Unit in pixels on screen.
This one is a bit tricky, since it refers to a static attribute (pixels on screen don't change)
but since every visual is attached to a camera, the exact scale might change.
So in the end, this is just relative to some normed camera - the value on screen, depending on the camera,
will not actually sit on those pixels. Only camera that guarantees the correct mapping is the
`:pixel` camera type.
"""
struct Pixel{T} <: Unit{T}
    value::T
end

const px = Pixel(1)

"""
    Billboard([angle::Real])
    Billboard([angles::Vector{<: Real}])

Billboard attribute to always have a primitive face the camera.
Can be used for rotation.
"""
struct Billboard{T <: Union{Float32, Vector{Float32}}}
    rotation::T
end
Billboard() = Billboard(0f0)
Billboard(angle::Real) = Billboard(Float32(angle))
Billboard(angles::Vector) = Billboard(Float32.(angles))


text:
# This file was generated, do not modify it.

using MLJ
using PyPlot
ioff() # hide
MLJ.color_off() # hide
using StableRNGs

linear = @load LinearRegressor pkg=MLJLinearModels
knn = @load KNNRegressor; knn.K = 4
tree_booster = @load EvoTreeRegressor; tree_booster.nrounds = 100
forest = @load RandomForestRegressor pkg=DecisionTree; forest.n_trees = 500
svm = @load SVMRegressor;

X = source()
y = source()

model1 = linear
model2 = knn

m1 = machine(model1, X, y)
y1 = predict(m1, X)

m2 = machine(model2, X, y)
y2 = predict(m2, X)

yhat = 0.5*y1 + 0.5*y2

mach = machine(Deterministic(), X, y; predict=yhat)

@from_network mach begin
    mutable struct MyAverageTwo
        regressor1=model1
        regressor2=model2
    end
end

average_two = MyAverageTwo()

function print_performance(model, data...)
    e = evaluate(model, data...;
                 resampling=CV(rng=StableRNG(1234), nfolds=8),
                 measure=rms,
                 verbosity=0)
    μ = round(e.measurement[1], sigdigits=5)
    ste = round(std(e.per_fold[1])/sqrt(8), digits=5)
    println("$model = $μ ± $(2*ste)")
end;

X, y = @load_boston

print_performance(linear, X, y)
print_performance(knn, X, y)
print_performance(average_two, X, y)

folds(data, nfolds) =
    partition(1:nrows(data), (1/nfolds for i in 1:(nfolds-1))...);

f = folds(1:10, 3)

corestrict(string.(1:10), f, 2)

figure(figsize=(8,6))
steps(x) = x < -3/2 ? -1 : (x < 3/2 ? 0 : 1)
x = Float64[-4, -1, 2, -3, 0, 3, -2, 1, 4]
Xraw = (x = x, )
yraw = steps.(x);
idxsort = sortperm(x)
xsort = x[idxsort]
ysort = yraw[idxsort]
step(xsort, ysort, label="truth", where="mid")
plot(x, yraw, ls="none", marker="o", label="data")
xlim(-4.5, 4.5)
legend()

savefig(joinpath(@OUTPUT, "s1.svg")) # hide

model1 = linear
model2 = knn

judge = linear

X = source(Xraw)
y = source(yraw)

f = @node folds(X, 3)

f()

folds(X::AbstractNode, nfolds) = node(XX->folds(XX, nfolds), X)
f = folds(X, 3)
f()

MLJ.restrict(X::AbstractNode, f::AbstractNode, i) =
    node((XX, ff) -> restrict(XX, ff, i), X, f);
MLJ.corestrict(X::AbstractNode, f::AbstractNode, i) =
    node((XX, ff) -> corestrict(XX, ff, i), X, f);

m11 = machine(model1, corestrict(X, f, 1), corestrict(y, f, 1))
m12 = machine(model1, corestrict(X, f, 2), corestrict(y, f, 2))
m13 = machine(model1, corestrict(X, f, 3), corestrict(y, f, 3))

y11 = predict(m11, restrict(X, f, 1));
y12 = predict(m12, restrict(X, f, 2));
y13 = predict(m13, restrict(X, f, 3));

y1_oos = vcat(y11, y12, y13);

fit!(y1_oos, verbosity=0)

figure(figsize=(8,6))
step(xsort, ysort, label="truth", where="mid")
plot(x, y1_oos(), ls="none", marker="o", label="linear oos")
legend()

savefig(joinpath(@OUTPUT, "s2.svg")) # hide

m21 = machine(model2, corestrict(X, f, 1), corestrict(y, f, 1))
m22 = machine(model2, corestrict(X, f, 2), corestrict(y, f, 2))
m23 = machine(model2, corestrict(X, f, 3), corestrict(y, f, 3))
y21 = predict(m21, restrict(X, f, 1));
y22 = predict(m22, restrict(X, f, 2));
y23 = predict(m23, restrict(X, f, 3));

y2_oos = vcat(y21, y22, y23);
fit!(y2_oos, verbosity=0)

figure(figsize=(8,6))
step(xsort, ysort, label="truth", where="mid")
plot(x, y2_oos(), ls="none", marker="o", label="knn oos")
legend()

savefig(joinpath(@OUTPUT, "s3.svg")) # hide

X_oos = MLJ.table(hcat(y1_oos, y2_oos))
m_judge = machine(judge, X_oos, y)

m1 = machine(model1, X, y)
m2 = machine(model2, X, y)

y1 = predict(m1, X);
y2 = predict(m2, X);
X_judge = MLJ.table(hcat(y1, y2))
yhat = predict(m_judge, X_judge)

fit!(yhat, verbosity=0)

figure(figsize=(8,6))
step(xsort, ysort, label="truth", where="mid")
plot(x, yhat(), ls="none", marker="o", label="yhat")
legend()

savefig(joinpath(@OUTPUT, "s4.svg")) # hide

e1 = rms(y1(), y())
e2 = rms(y2(), y())
emean = rms(0.5*y1() + 0.5*y2(), y())
estack = rms(yhat(), y())
@show e1 e2 emean estack;

@from_network machine(Deterministic(), X, y; predict=yhat) begin
    mutable struct MyTwoModelStack
        regressor1=model1
        regressor2=model2
        judge=judge
    end
end

my_two_model_stack = MyTwoModelStack()

X, y = make_regression(1000, 20; sparse=0.75, noise=0.1, rng=123);

avg = MyAverageTwo(regressor1=tree_booster,
                   regressor2=svm)

stack = MyTwoModelStack(regressor1=tree_booster,
                        regressor2=svm,
                        judge=forest) # judge=linear

all_models = [tree_booster, svm, forest, avg, stack];

for model in all_models
    print_performance(model, X, y)
end

r = range(stack, :(regressor2.C), lower = 0.01, upper = 10, scale=:log)
tuned_stack = TunedModel(model=stack,
                         ranges=r,
                         tuning=Grid(shuffle=false),
                         measure=rms,
                         resampling=Holdout())

mach = fit!(machine(tuned_stack,  X, y), verbosity=0)
best_stack = fitted_params(mach).best_model
best_stack.regressor2.C

print_performance(best_stack, X, y)

PyPlot.close_figs() # hide


text:
# query_string = "DELETE FROM appuser WHERE id = \$1"
function delete_entity(deleted_object::IEntity,
                       data_type::DataType,
                       table_name::String,
                       columns_selection_and_mapping::Dict{Symbol,<:Any},
                       id_properties::Vector{Symbol},
                       dbconn::LibPQ.Connection;
                       editor::Union{IEntity,Missing} = missing)

    query_string = "DELETE FROM " * table_name

    # Check that we have a filter object
    if ismissing(deleted_object)
        throw(DomainError("You must specify an object to delete"))
    end

    # Loop over the properties of the object and build the appropriate where clause
    # props = PostgresORMUtil.getproperties_asdict(deleted_object,false)
    # props = util_remove_unmapped_properties_from_dict(props,
                                                      # columns_selection_and_mapping)

    props = util_get_entity_props_for_db_actions(deleted_object,
                                                 dbconn,
                                                 false # do not include missing
                                                       #  (although it doesn't make any difference anyway)
                                                 )

    # Compare versions if needed, by retrieving the current state of the object
    #   in the database
    orm_module = PostgresORM.get_orm(deleted_object)
    if util_get_track_changes(orm_module)

        # Create a string representation of the ID properties values.
        id_value_as_str = util_prettify_id_values_as_str(deleted_object)

        # Get the editor id if any
        editor_id = if ismissing(editor)
            missing
        else
            util_prettify_id_values_as_str(editor)
        end

        action_id = uuid4()

        modif = Modification(entity_type =  string(data_type),
                             entity_id = id_value_as_str,
                             appuser_id = editor_id,
                             creation_time = now(Dates.UTC),
                             action_type = CRUDType.delete,
                             action_id = action_id)
        create_entity!(modif,dbconn)

     end # if track_changes

     # Prepare the vectors used for creating the query
     properties_names = collect(keys(props))
     properties_values = collect(values(props))

     prepared_statement_args_values = Any[] # Holds the values of the prepared statement
     prepared_statement_args_counter = 0

     # Add the 'WHERE' clause:
     # id_columnname = columns_selection_and_mapping[id_property]
     # query_string *= " WHERE $id_columnname = \$1"
     # push!(prepared_statement_args_values,
     #      getproperty(deleted_object, id_property))

     # Add the 'WHERE' clause:
     # Loop over the ID properties
     query_string *= " WHERE "
     where_clause_elts = Vector{String}()
     # An object can have several ID properties
     for id_property in util_get_ids_props_names(deleted_object)
          # An ID property can have several columns (eg. when the property is a
          #   complex property that has several IDs)
          id_columnnames = tovector(columns_selection_and_mapping[id_property])
          for id_columnname in id_columnnames
              where_clause_elt =  "$id_columnname = \$$(prepared_statement_args_counter+=1)"
              push!(where_clause_elts, where_clause_elt)
              # If the ID property is a complex prop, then the value is stored
              #   under the name of the column name
              if isa(getproperty(deleted_object, id_property),IEntity)
                  push!(prepared_statement_args_values,
                        props[Symbol(id_columnname)])
              else
                  push!(prepared_statement_args_values,
                        props[id_property])
              end

          end
     end
     query_string *= join(where_clause_elts, " AND ")

     prepared_query = LibPQ.prepare(dbconn,
                                   query_string)

     # @info query_string
     # @info prepared_statement_args_values

     # Prepare the query aruments
     query_result = execute(prepared_query,
                           prepared_statement_args_values
                           ;throw_error=true)

      # Get the result as a dataframe
      result = LibPQ.num_affected_rows(query_result)

      # Close the result
      close(query_result)

      return result

end

function delete_entity(deleted_object::IEntity,
                       orm_module::Module,
                       dbconn::LibPQ.Connection;
                       editor::Union{IEntity,Missing} = missing)

    delete_entity(deleted_object,
                  orm_module.data_type,
                  util_get_table_with_schema_name(orm_module),
                  util_get_columns_selection_and_mapping(orm_module),
                  util_get_ids_props_names(deleted_object),
                  dbconn;
                  editor = editor)
end

function delete_entity(deleted_object::IEntity,
                       dbconn::LibPQ.Connection;
                       editor::Union{IEntity,Missing} = missing
                       )
    delete_entity(deleted_object,
                  PostgresORM.get_orm(deleted_object),
                  dbconn;
                  editor = editor)
end

# query_string = "DELETE FROM appuser WHERE id = \$1 AND attr1 = \2 AND ..."
function delete_entity_alike(filter_object::IEntity,
                             data_type::DataType,
                             table_name::String,
                             columns_selection_and_mapping::Dict{Symbol,<:Any},
                             dbconn::LibPQ.Connection;
                             editor::Union{IEntity,Missing} = missing)

    query_string = "DELETE FROM " * table_name

    # Check that we have a filter object
    if ismissing(filter_object)
        throw(DomainError("You must specify an object to delete"))
    end

    # Loop over the properties of the object and build the appropriate where clause
    props = PostgresORMUtil.getproperties_asdict(filter_object,false)
    props = util_remove_unmapped_properties_from_dict(props,
                                                      columns_selection_and_mapping)

      # Instantiate an empty instance of the data type just to get the module
      orm_module = PostgresORM.get_orm(data_type())

      # Loop over the properties of the object and build the appropriate where clause
      props = util_get_entity_props_for_db_actions(filter_object,
                                                   dbconn,
                                                   false # do not include missing
                                                   )
      if length(props) > 0
          n = 0
          query_string *= " WHERE "
          for (propertyname,value) in props
              n += 1
              # Add the 'AND' if it's not the first restriction
              if n > 1
                  query_string *= " AND "
              end

              # NOTE: The mapping may not have the properties because some of the
              #          properties are the columns names of the complex objects
              #          which are added at runtime
              columnname = if haskey(columns_selection_and_mapping, propertyname)
                              columns_selection_and_mapping[propertyname]
                           else
                               propertyname
                           end

              # Add the restriction
              query_string *= " $columnname = \$$n " # eg. 'a.login = \$1'

          end
      end

     # prepared_query = LibPQ.prepare(dbconn,
     #                                query_string)

     # Prepare the query aruments
     # NOTE: It would be more elegant to convert the keys in an array of String
     #         but I didn't find a way to do it
     query_args = collect(values(props))

     prepared_query = LibPQ.prepare(dbconn,
                                   query_string)

     # Prepare the query aruments
     query_result = execute(prepared_query,
                            query_args
                           ;throw_error=true)

      # Get the result as a dataframe
      result = LibPQ.num_affected_rows(query_result)

      # Close the result
      close(query_result)

      return result

end

function delete_entity_alike(filter_object::IEntity,
                       orm_module::Module,
                       dbconn::LibPQ.Connection;
                       editor::Union{IEntity,Missing} = missing)

    delete_entity_alike(filter_object,
                  orm_module.data_type,
                  util_get_table_with_schema_name(orm_module),
                  util_get_columns_selection_and_mapping(orm_module),
                  dbconn;
                  editor = editor)
end

function delete_entity_alike(filter_object::IEntity,
                       dbconn::LibPQ.Connection;
                       editor::Union{IEntity,Missing} = missing
                       )
    delete_entity_alike(filter_object,
                  PostgresORM.get_orm(filter_object),
                  dbconn;
                  editor = editor)
end


text:
# bandwidth(md::MinimumDivergenceProblem) = md.mf.sf.St
# κ₁(md::MinimumDivergenceProblem) = md.mf.sf.κ₁
# κ₂(md::MinimumDivergenceProblem) = md.mf.sf.κ₂

# bandwidth(mf::MomentFunction) = mf.sf.St
# κ₁(mf::MomentFunction) = mf.sf.κ₁
# κ₂(mf::MomentFunction) = mf.sf.κ₂

# nobs(md::MDP) = md.n
# ncond(md::MDP) = md.m

# npar(md::MinimumDivergenceProblem) = md.p
# npar(md::MinimumDivergenceProblemPlain) = 0.0

# lambda(md::MinimumDivergenceProblem) = md.lambda
# coef(md::MinimumDivergenceProblem) = md.fprob.x[nobs(md)+1:end]
# weights(md::MDP) = md.fprob.x[1:nobs(md)]
# obj_val(md::MDP) = md.fprob.obj_val     ## Should I multiply this by 2?

# function vcov!(md::MinimumDivergenceProblem)
#   global __p

#   if typeof(md.Ω)==Nothing
#     meat!(md)
#   end

#   if typeof(md.∇g_n)==Nothing
#     mfjacobian!(md)
#   end

#   w   = weights(md)
#   θ   = coef(md)
#   __p = ones(nobs(md))
#   uw_g_i = md.mf.g_i(θ)
#   __p = weights(md)
#   w_g_i = md.mf.g_i(θ)

#   nvar = length(md.Ω)

#   out = Array(PDMat, nvar)

#   for j=1:nvar
#     out[j] = inv(PDMat(Xt_invA_X(md.Ω[j], md.∇g_n[j].∇)))
#   end

#   md.Σ = out

# end

# vcov(md::MinimumDivergenceProblem) = vcov(md::MinimumDivergenceProblem, :weighted)  ## Default to weighted var
# vcov(md::MinimumDivergenceProblem, vtype::Symbol) = vtype==:unweighted ? md.Σ[2] : md.Σ[1]

# stderr(md::MinimumDivergenceProblem) = stderr(md::MinimumDivergenceProblem, :weighted)

# ## function stderr(md::MinimumDivergenceProblem, vtype=:Symbol)
# ##     if vtype==:unweighted
# ##         sqrt(diag(md.Σ[2]))
# ##     elseif vtype==:weighted
# ##         sqrt(diag(md.Σ[1]))
# ##     elseif vtype==:hessian
# ##         md.Σ[1]
# ## end

# function mfjacobian!(md::MinimumDivergenceProblem)
#   global __p
#   w   = weights(md)
#   θ   = coef(md)

#   out = Array(MomentFunctionJacobian, 2)

#   __p = ones(nobs(md))
#   out[1] = MomentFunctionJacobian(md.mf.Dg_n(θ))
#   __p = weights(md)
#   out[2] = MomentFunctionJacobian(md.mf.Dg_n(θ)./κ₁(md))
#   ## This is devided by κ₁
#   md.∇g_n = out
# end

# function meat!(md::MinimumDivergenceProblem)
#     global __p
#     if typeof(md.∇g_n)==Nothing
#       mfjacobian!(md)
#     end
#     ## Adjustment constant
#     γ = md.mf.sf.St/md.mf.sf.κ₂

#     θ   = coef(md)
#     __p = ones(nobs(md))
#     uw_g_i = md.mf.g_i(θ)
#     __p = weights(md)
#     w_g_i = md.mf.g_i(θ)
#     md.Ω =  [PDMat(γ.*uw_g_i'*uw_g_i), PDMat(γ.*w_g_i'*w_g_i)]
# end

# ## Jacobian mdproblem
# ## Calculate
# ##
# ## \partial obj(theta) / \partial \theta
# ##

# function gradient(md::MinimumDivergenceProblem)
#   global __p
#   __p = weights(md)
#   md.mf.Dg_n(coef(md))'*md.lambda[1:ncond(md)]
# end

# function gradient(mdobj::MinimumDivergenceProblem, theta_0::Array{Float64,1}, s::Symbol)
#   obj(theta) = obj_val(md(mdobj.mf.g_i(theta), mdobj.div))
#   Calculus.gradient(obj, theta_0)
# end

# function hessian!(mdobj::MinimumDivergenceProblem)
#   mdobj.H = hessian(mdobj, coef(mdobj))
# end

# function hessian(mdobj::MinimumDivergenceProblem)
#   if typeof(mdobj.H)==Nothing
#     hessian!(mdobj)
#   end
#   return mdobj.H
# end

# function hessian(mdobj::MinimumDivergenceProblem, theta_0::Array{Float64,1})
#   obj(theta) = obj_val(md(mdobj.mf.g_i(theta), mdobj.div))
#   PDMat(Calculus.hessian(obj, theta_0))
# end

# function show(io::IO, obj::MinimumDivergenceProblem)
#     println(io,"$(typeof(obj)):\n\nCoefficients:", "\n", coef(obj))
# end

# ## function coeftable(mm::MinimumDivergenceProblem)
# ##     cc = coef(mm)
# ## #    se = stderr(mm)
# ## #    tt = cc ./ se
# ##     ## CoefTable(hcat(cc,se,tt,ccdf(FDist(1, df_residual(mm)), abs2(tt))),
# ##     ##           ["Estimate","Std.Error","t value", "Pr(>|t|)"],
# ##     ##           ["θ\_$i" for i = 1:length(cc)], 4)
# ##     println(["θ\_$i: cc[$i]" for i = 1:length(cc)])
# ## end


text:
#using Debugger
using StructuralCausalModels, Test

ProjDir = @__DIR__
cd(ProjDir) #do

fname = scm_path("..", "examples", "Shipley", "fig2.6", "fig2.6.dot")
#Sys.isapple() && run(`open -a GraphViz.app $(fname)`)

# Read `=>` as `~` in regression models, or `<-` in causal models.
d = OrderedDict(
  :u => [:x, :v],
  :s1 => [:u],
  :w => [:v, :y],
  :s2 => [:w]
);

dag = DAG("fig2.6", d);

@testset "d_separation" begin

  @test to_ggm(dag) == "DAG(u ~ x + v, s1 ~ u, w ~ v + y, s2 ~ w)"
  @test d_separation(dag, :x, :v; debug=false) == true
  @test d_separation(dag, :x, [:v]; c=:u, debug=false) == false
  @test d_separation(dag, [:x], :v; c=:s1, debug=false) == false
  @test d_separation(dag, [:u], [:w]; debug=false) == false
  @test d_separation(dag, :u, [:w]; c=:v, debug=false) == true
  @test d_separation(dag, [:x], :y; debug=false) == true
  @test d_separation(dag, :x, :y; c=[:u, :w], debug=false) == false
  @test d_separation(dag, :x, :y; c=[:s1, :s2], debug=false) == false
  @test d_separation(dag, :x, :y; c=[:u, :v, :w], debug=false) == true
  @test d_separation(dag, :x, [:v]; c=:u, debug=false) == false
  @test d_separation(dag, [:x], :v; c=:s1, debug=false) == false

  include("test_descendents_02.jl")
  @test d_separation(dag, :k1, :v, c=[:x1]) == true
  @test d_separation(dag, :k1, :v, c=[:w]) == true
  @test d_separation(dag, :k1, :x1, c=[:w]) == true
  @test d_separation(dag, :k1, :x2, c=[:v]) == true
  @test d_separation(dag, :k1, :x2, c=[:x1]) == true
  @test d_separation(dag, :k1, :x2, c=[:w]) == true
  @test d_separation(dag, :k1, :x3, c=[:w]) == true
  @test d_separation(dag, :k1, :y, c=[:k2, :x2, :x3]) == true
  @test d_separation(dag, :k1, :y, c=[:k2, :v, :x3]) == true
  @test d_separation(dag, :k1, :y, c=[:k2, :x1, :x3]) == true
  @test d_separation(dag, :k1, :y, c=[:k2, :w]) == true
  @test d_separation(dag, :k2, :v, c=[:x1]) == true
  @test d_separation(dag, :k2, :v, c=[:w]) == true
  @test d_separation(dag, :k2, :v, c=[:k1]) == true
  @test d_separation(dag, :k2, :w, c=[:k1]) == true
  @test d_separation(dag, :k2, :x1, c=[:w]) == true
  @test d_separation(dag, :k2, :x1, c=[:k1]) == true
  @test d_separation(dag, :k2, :x2, c=[:v]) == true
  @test d_separation(dag, :k2, :x2, c=[:x1]) == true
  @test d_separation(dag, :k2, :x2, c=[:w]) == true
  @test d_separation(dag, :k2, :x2, c=[:k1]) == true
  @test d_separation(dag, :k2, :x3, c=[:w]) == true
  @test d_separation(dag, :k2, :x3, c=[:k1]) == true
  @test d_separation(dag, :v, :w; c=:x1) == true
  @test d_separation(dag, :v, :x3) == true
  @test d_separation(dag, :v, :y; c=[:k2, :x2, :x3]) == true
  @test d_separation(dag, :v, :y; c=[:k1, :x2, :x3]) == true
  @test d_separation(dag, :v, :y; c=[:w, :x2, :x3]) == true
  @test d_separation(dag, :v, :y; c=[:x2, :x1]) == true
  @test d_separation(dag, :w, :x2; c=[:v]) == true
  @test d_separation(dag, :w, :x2; c=[:x1]) == true
  @test d_separation(dag, :w, :y; c=[:k2, :x2, :x3]) == true
  @test d_separation(dag, :w, :y; c=[:k2, :v, :x3]) == true
  @test d_separation(dag, :w, :y; c=[:k1, :x2, :x3]) == true
  @test d_separation(dag, :w, :y; c=[:k1, :v, :x3]) == true
  @test d_separation(dag, :w, :y; c=[:k1, :x1, :x3]) == true
  @test d_separation(dag, :w, :y; c=[:k2, :x2, :x3]) == true
  @test d_separation(dag, :x1, :x2; c=[:v]) == true
  @test d_separation(dag, :x1, :x3) == true
  @test d_separation(dag, :x1, :y; c=[:k2, :x2, :x3]) == true
  @test d_separation(dag, :x1, :y; c=[:k1, :x2, :x3]) == true
  @test d_separation(dag, :x1, :y; c=[:w, :x2, :x3]) == true
  @test d_separation(dag, :x1, :y; c=[:k2, :v, :x3]) == true
  @test d_separation(dag, :x1, :y; c=[:k1, :v, :x3]) == true
  @test d_separation(dag, :x1, :y; c=[:v, :w, :x3]) == true
  @test d_separation(dag, :x2, :x3) == true

end


text:
""" Geometrical shapes. """

using LinearAlgebra, CoordinateTransformations

export Shape, Cylinder, Sphere, Cone, Plane, Empty, shapediff, inside,
    TransformedShape, mathematica

abstract type Shape; end

""" A Sphere centered at (`xc`, `yc`, `zc`) with radius `R`. """
struct Sphere <: Shape
    xc::Float64
    yc::Float64
    zc::Float64
    R::Float64
end

"""
A Cylinder between `bottom` < z < `top`, centered at (`xc`, `yc`) and with
radius `R`.
"""
struct Cylinder <: Shape
    bottom::Float64
    top::Float64
    xc::Float64
    yc::Float64
    R::Float64
end

"""
A Cone with a vertical axis and `bottom` < z < `top`, with vertex
at (`xv`, `yv`, `zv`) and slope m with
sqrt((x - xv)^2 + (y - yv)^2) = m (z - zv).
"""
struct Cone <: Shape
    bottom::Float64
    top::Float64

    xv::Float64
    yv::Float64
    zv::Float64
    m::Float64
end

"""
A plane at `z` that defines a shape containing points above
(if `up` is true) or below otherwise.
"""
struct Plane <: Shape
    z::Float64
    up::Bool
end

""" Union of several shapes. """
struct ShapeUnion{T<:Tuple} <: Shape
    shapes::T
end

Base.union(s::Shape...) = ShapeUnion(s)

Empty() = ShapeUnion(())

""" Intersection of several shapes. """
struct ShapeIntersect{T<:Tuple} <: Shape
    shapes::T
end

Base.intersect(s::Shape...) = ShapeIntersect(s)

""" Shape1 - Shape2. """
struct ShapeSubstract{S1,S2} <: Shape
    shape1::S1
    shape2::S2
end

shapediff(s1, s2) = ShapeSubstract(s1, s2)

# Shortcuts for unions and substractiong
Base.:-(s1::Shape, s2::Shape) = shapediff(s1, s2)
Base.:+(s1::Shape, rest::Shape...) = union(s1, rest...)

""" Transformed shape after a transformation or shift. """
struct TransformedShape{T<:AbstractAffineMap, S<:Shape} <: Shape
    trans::T
    inv::T
    shape::S
    function TransformedShape(trans::T, shape::S) where
        {T<:AbstractAffineMap, S<:Shape}

        inv_ = inv(trans)
        new{T, S}(trans, inv_, shape)
    end
end

# For some reason, dot for SVectors of size 3 is much slower than this
# explicit computation.
@inline @fastmath fdot(a::Point, b::Point) = a[1] * b[1] + a[2] * b[2] + a[3] * b[3]

# When we look for interface crossings sometimes we want a list of them,
# sometimes just the closest one.  To reuse the same function without
# creating too many closures we use accumulators for each possible use case
abstract type Accumulator; end

struct ListAccumulator{T <: AbstractVector} <: Accumulator
    list::T
end

@inline accum!(a::ListAccumulator, x) = push!(a.list, x)

mutable struct MinAccumulator{T <: Real} <: Accumulator
    value::T
end

@inline accum!(a::MinAccumulator, x) = (x < a.value) && (a.value = x)

# Just for debugging
struct ShowAccumulator <: Accumulator
end

accum!(a::ShowAccumulator, x) = @show x

"""
Compute the intersection points between a shape `c` and the line joining
points `a` and `b`.  Each intersection point s is used in the accumulator vwhere the intersection point is
p = (1-s)a + s b
"""
@fastmath function interfaces!(v::Accumulator, s::Sphere, a::Point, b::Point)
    (xc, yc, zc) = (s.xc, s.yc, s.zc)
    (xa, ya, za) = a .- @SVector [s.xc, s.yc, s.zc]
    (xb, yb, zb) = b .- @SVector [s.xc, s.yc, s.zc]

    a2 = xa^2 + ya^2 + za^2
    b2 = xb^2 + yb^2 + zb^2
    ab = xa * xb + ya * yb + za * zb

    # These are traditional a, b, c in the quadratic eq. but a, b are already
    # in use.
    α = a2 + b2 - 2 * ab
    β = -2*a2 + 2 * ab
    γ = a2 - s.R^2

    Δ = β^2 - 4 * α * γ

    (Δ >= 0) || return nothing

    s1 = (-β + sqrt(Δ)) / 2α
    s2 = (-β - sqrt(Δ)) / 2α

    (0 < s1 < 1) && accum!(v, s1)
    (0 < s2 < 1) && accum!(v, s2)

    nothing
end

@fastmath function interfaces!(v::Accumulator, c::Cylinder, a::Point, b::Point)
    # Aliases to simplify the formulas
    (xa, ya, za) = a
    (xb, yb, zb) = b
    (xc, yc, R) = (c.xc, c.yc, c.R)

    st = (c.top - za) / (zb - za)
    sb = (c.bottom - za) / (zb - za)

    h2 = (xa - xb)^2 + (ya - yb)^2

    # Vertical lines.
    if (h2 == 0)

        # Check if we are inside or out of the cylinder.
        ((xa - xc)^2 + (ya - yc)^2 > R^2) && return nothing

        # If we are inside, check if the caps are between a and b
        (0 < st < 1) && accum!(v, st)
        (0 < sb < 1) && accum!(v, sb)

        return nothing
    end

    Δ = R^2 * h2 - (xb * ya - xc * ya - xa * yb + xc * yb + xa * yc - xb * yc)^2
    m = xa^2 + xb * xc - xa * (xb + xc) + (ya - yb) * (ya - yc)

    Δ < 0 && return nothing

    s1 = (m - sqrt(Δ)) / h2
    s2 = (m + sqrt(Δ)) / h2

    z1 = (1 - s1) * a[3] + s1 * b[3]
    z2 = (1 - s2) * a[3] + s2 * b[3]

    (0 < st < 1) && (s1 < st < s2) && accum!(v, st)
    (0 < sb < 1) && (s1 < sb < s2) && accum!(v, sb)

    (0 < s1 < 1 && c.bottom < z1 < c.top) && accum!(v, s1)
    (0 < s2 < 1 && c.bottom < z2 < c.top) && accum!(v, s2)

    nothing
end

@fastmath function interfaces!(v::Accumulator, c::Cone, a::Point, b::Point)
    # Shift everything wrt the cone vertex
    (xa, ya, za) = a .- @SVector [c.xv, c.yv, c.zv]
    (xb, yb, zb) = b .- @SVector [c.xv, c.yv, c.zv]
    m = c.m
    m2 = m^2

    at2 = xa * xa + ya * ya
    bt2 = xb * xb + yb * yb
    adotb = xa * xb + ya * yb

    # with this we define a quadratic equation on s with coeffs
    α = at2 + bt2 - 2 * adotb - m2 * (za^2 + zb^2 - 2 * za * zb)
    β = -2 * at2 + 2 * adotb - m2 * (-2 * za^2 + 2 * za * zb)
    γ = at2 - m2 * za^2

    if α != 0
        Δ = sqrt(β^2 - 4 * α * γ)
        s1, s2 = (-β - Δ) / (2 * α), (-β + Δ) / (2 * α)
    else
        # Unfortunately we must allow the possibility that the line connecting
        # a and b passes though the vertex of the cone.  In that case
        # the line is either almost always inside or almost always outside
        # it.
        s1, s2 = -Inf, za / (za - zb)
    end

    # Make sure that s1 < s2
    (s1 < s2) || ((s1, s2) = (s2, s1))

    st = (c.top - za - c.zv) / (zb - za)
    sb = (c.bottom - za - c.zv) / (zb - za)

    (0 < st < 1) && (s1 < st < s2) && accum!(v, st)
    (0 < sb < 1) && (s1 < sb < s2) && accum!(v, sb)

    z1 = (1 - s1) * za + s1 * zb + c.zv
    z2 = (1 - s2) * za + s2 * zb + c.zv

    (0 < s1 < 1 && c.bottom < z1 < c.top) && accum!(v, s1)
    (0 < s2 < 1 && c.bottom < z2 < c.top) && accum!(v, s2)

    nothing
end

function interfaces!(v::Accumulator, p::Plane, a::Point, b::Point)
    # Aliases to simplify the formulas
    (xa, ya, za) = a
    (xb, yb, zb) = b
    zp = p.z

    s = (zp - za) / (zb - za)
    (0 < s < 1) && accum!(v, s)
    nothing
end

# Since a ShapeUnion or ShapeIntersect may be composed of shapes of different
# types, we use a generated function to unroll the loop and make the function
# type-stable.  This is cool.
@generated function interfaces!(v::Accumulator,
                                u::Union{ShapeUnion{T}, ShapeIntersect{T}},
                                a::Point, b::Point) where {T}
    L = fieldcount(T)
    out = quote end

    for i in 1:L
        push!(out.args,
              quote
              interfaces!(v, u.shapes[$i], a, b)
              end
              )
    end
    push!(out.args, :(return nothing))

    out
end

function interfaces!(v::Accumulator, s::ShapeSubstract, a::Point, b::Point)
    interfaces!(v, s.shape1, a, b)
    interfaces!(v, s.shape2, a, b)
    nothing
end

function interfaces!(v::Accumulator, s::TransformedShape,
                     a::Point, b::Point)
    at = s.inv(a)
    bt = s.inv(b)

    interfaces!(v, s.shape, at, bt)
end

"""
Checks whether point `a` is inside a given shape.
"""
function inside(s::Sphere, a::Point)
    (xa, ya, za) = a
    (xc, yc, zc, R) = (s.xc, s.yc, s.zc, s.R)
    ((xa - xc)^2 + (ya - yc)^2 + (za - zc)^2 < R^2)
end

function inside(c::Cylinder, a::Point)
    (xa, ya, za) = a
    (xc, yc, R, z1, z2) = (c.xc, c.yc, c.R, c.bottom, c.top)

    ((xa - xc)^2 + (ya - yc)^2 > R^2) && return false
    (z1 < za < z2)
end

function inside(c::Cone, a::Point)
    (c.bottom <= a[3] <= c.top) || return false

    (xa, ya, za) = a .- @SVector [c.xv, c.yv, c.zv]
    m = c.m

    (xa^2 + ya^2 > m^2 * za^2) && return false
    true
end

inside(p::Plane, a::Point) = (p.up == (a[3] > p.z))

# See above for why we use a generated function
@generated function inside(u::ShapeUnion{T}, a::Point) where {T}
    L = fieldcount(T)
    out = quote end

    for i in 1:L
        push!(out.args,
              quote
              inside(u.shapes[$i], a) && return true
              end
              )
    end
    push!(out.args, :(return false))

    out
end

@generated function inside(u::ShapeIntersect{T}, a::Point) where {T}
    L = fieldcount(T)
    out = quote end

    for i in 1:L
        push!(out.args,
              quote
              inside(u.shapes[$i], a) || return false
              end
              )
    end
    push!(out.args, :(return true))

    out
end

function inside(u::ShapeIntersect, a::Point)
    for item in u.shapes
        inside(item, a) || return false
    end
    true
end

inside(s::ShapeSubstract, a::Point) = inside(s.shape1, a) && !inside(s.shape2, a)

function inside(s::TransformedShape, a::Point)
    inside(s.shape, s.inv(a))
end

# To facilitate visualization we provide methods that allow the construction of
# a Mathematica expression that can be viewed with Mathematica versions above
# 11.2
mathlist(v) = "{" * join(v, ", ") * "}"
mathcall(func, v...) = func * "[" * join(v, ", ") * "]"

function mathematica(fig::Sphere)
    mathcall("Ball", mathlist([fig.xc, fig.yc, fig.zc]), fig.R)
end

function mathematica(fig::Cylinder)
    p1 = [fig.xc, fig.yc, fig.bottom]
    p2 = [fig.xc, fig.yc, fig.top]

    mathcall("Cylinder", mathlist([mathlist(p1), mathlist(p2)]),
             fig.R)
end

function mathematica(fig::Plane)
    mathcall("HalfSpace", mathlist([0, 0, fig.up ? -1 : 1]),
             mathlist([0, 0, fig.z]))
end

# Note that we are not currently supporting "double cones"
function mathematica(fig::Cone)
    rbottom = fig.m * abs(fig.bottom - fig.zv)
    rtop = fig.m * abs(fig.top - fig.zv)
    vertex = [fig.xv, fig.yv, fig.zv]

    if rbottom > rtop
        base = [fig.xv, fig.yv, fig.bottom]
        plane = Plane(fig.top, false)
    else
        base = [fig.xv, fig.yv, fig.top]
        plane = Plane(fig.top, true)
    end
    mathcall("RegionIntersection",
             mathcall("Cone", mathlist([mathlist(base), mathlist(vertex)]),
                      rbottom),
             mathematica(plane))
end

function mathematica(fig::ShapeUnion)
    mathcall("RegionUnion",
             [mathematica(x) for x in fig.shapes]...)
end

function mathematica(fig::ShapeIntersect)
    mathcall("RegionIntersection",
             [mathematica(x) for x in fig.shapes]...)
end

function mathematica(fig::ShapeSubstract)
    mathcall("RegionDifference",
             mathematica(fig.shape1), mathematica(fig.shape2))
end

function mathematica(fig::TransformedShape)
    mathcall(
        mathcall("AffineTransform",
                 mathlist(
                     [mathlist([mathlist(r) for r in eachrow(fig.trans.linear)]),
                      mathlist(fig.trans.translation)])),
        mathematica(fig.shape))
end


text:
## Step1Grid is a simple grid search to find optimal policy function

#module Step1Grid

using LinearAlgebra, Random, Distributions, Statistics, Plots
using BenchmarkTools, Interpolations

"""
    Performs a grid search to find the optimal policy function.
    Takes price q as given.
"""

# Define the function for Step1: iteration to find the optimal policy function

function step1(A, E, Π, q, β, σ; v_guess = 0.0, maxT = 600, tol = 0.01)

    """
    Function that finds the optimal policy function for a given price q.
    Inputs:
        - A: set of possible assets (endog. state var.)
        - E: set of exogenous shocks (exog. state var.)
        - Π: transition matrix
        - q: price
        - β: discount factor
        - σ: preferences parameter
    Output:
        - pol_func: optimal policy function
    """

    # Check that parameters and price make sense
    @assert β < 1
    @assert σ > 1
    @assert q > β

    a_size = length(A)
    e_size = size(E)[1]

    # Check that transition matrix is well defined
    for e = 1:e_size
        @assert sum(Π[e,:]) == 1
    end

    v0 = v_guess .* ones(a_size, e_size)
    v_update = zeros(a_size, e_size)
    pol_func = zeros(a_size, e_size)

    # Define the utility function (in the future, put it outside)
    u(c) = c.^(1-σ) ./ (1 - σ)

    # Iteration of value function
    normdiff = Inf
    t = 1

    while normdiff > tol && t <= maxT

        # Do linear interpolation of the value function on knots (A,E)
        v_interpol = LinearInterpolation((A, E), v0)

        # Replace new value function with solution
        for a = 1:a_size
            for e = 1:e_size

                # Generate value function V
                # If c <= 0, set value function to -∞
                V(x) = ((E[e] + A[a] - q*x) <= 0) ? -Inf :
                   u(E[e] + A[a] - q*x) + β*(v_interpol(x, E[1])* Π[e,1] +
                        v_interpol(x, E[2])* Π[e,2])

                # Find the maximum
                aux = findmax(V.(A))      # (maximum, index)
                v_update[a,e] = aux[1]    # actual maximum value
                pol_func[a,e] = A[aux[2]] # policy function a'

            end
        end

        t += 1
        normdiff = maximum(abs.(v_update - v0))

        v0 = copy(v_update)

    end

    return pol_func

end

#export step1

#end


text:
# Support Data for color maps
include("maps_data.jl")

# color scale generation
# ----------------------

@noinline function _generate_lab(l::Float32, c::Float32, h::Float32)
    lab = convert(Lab{Float32}, LCHab{Float32}(l, c, h))
    rgb = xyz_to_linear_rgb(convert(XYZ{Float32}, lab))
    convert(Lab{Float32}, linear_rgb_to_xyz(correct_gamut(rgb)))
end

"""
    colors = distinguishable_colors(n, seed=RGB{N0f8}[];
                                    dropseed=false,
                                    transform=identity,
                                    lchoices=range(0, stop=100, length=15),
                                    cchoices=range(0, stop=100, length=15),
                                    hchoices=range(0, stop=342, length=20))

Generate `n` maximally distinguishable colors.

This uses a greedy brute-force approach to choose `n` colors that are maximally
distinguishable. Given seed color(s), and a set of possible hue, chroma, and
lightness values (in LCHab space), it repeatedly chooses the next color as the
one that maximizes the minimum pairwise distance to any of the colors already
in the palette.

# Arguments

- `n`: Number of colors to generate.
- `seed`: Initial color(s) included in the palette.

# Keyword arguments

- `dropseed`: if true, the `seed` values will be dropped. This provides an easy
  mechanism to ensure that the chosen colors are distinguishable from the seed value(s).
  When true, `n` does not include the seed color(s).
- `transform`: Transform applied to colors before measuring distance. Default is `identity`;
  other choices include `deuteranopic` to simulate color-blindness.
- `lchoices`: Possible lightness values
- `cchoices`: Possible chroma values
- `hchoices`: Possible hue values

Returns a `Vector` of colors of length `n`, of the type specified in `seed`.
"""
function distinguishable_colors(n::Integer,
        @nospecialize(seed::AbstractVector{<:Color});
        dropseed::Bool = false,
        transform::Function = identity,
        lchoices::AbstractVector{<:Real} = range(0.0f0, stop=100.0f0, length=15),
        cchoices::AbstractVector{<:Real} = range(0.0f0, stop=100.0f0, length=15),
        hchoices::AbstractVector{<:Real} = range(0.0f0, stop=342.0f0, length=20))
    if n <= length(seed) && !dropseed
        return seed[1:n]
    end

    # Candidate colors
    N = length(lchoices) * length(cchoices) * length(hchoices)
    candidate = Vector{Lab{Float32}}(undef, N)
    j = 0
    for h in hchoices, c in cchoices, l in lchoices
        @inbounds candidate[j+=1] = _generate_lab(Float32(l), Float32(c), Float32(h))
    end

    _distinguishable_colors(Int(n), seed, dropseed, transform, candidate)
end

function _distinguishable_colors(n::Int,
        @nospecialize(seed::AbstractVector{<:Color}),
        dropseed::Bool,
        @nospecialize(transform::Function),
        candidate::Vector{Lab{Float32}})

    N = length(candidate)
    colors = Vector{eltype(seed)}(undef, n)

    # Transformed colors
    if transform === identity
        candidate_t = candidate
    else
        candidate_t = convert.(Lab{Float32}, transform.(candidate))::Vector{Lab{Float32}}
    end

    # Start with the seed colors, unless `dropseed`
    dropseed || copyto!(colors, seed)

    # Minimum distances of the current color to each previously selected color.
    ds = fill(Inf32, N)
    @inbounds for s in seed
        ts = convert(Lab{Float32}, transform(s))::Lab{Float32}
        for k = 1:N
            ds[k] = @fastmath min(ds[k], colordiff(ts, candidate_t[k]))
        end
    end

    n1 = dropseed ? 1 : length(seed) + 1
    @inbounds for i in n1:n
        j = argmax(ds)
        colors[i] = candidate[j]
        tc = candidate_t[j]
        ds[j] = 0.0f0
        for k = 1:N
            ds[k] == 0.0f0 && continue # already selected
            ds[k] = @fastmath min(ds[k], colordiff(tc, candidate_t[k]))
        end
    end

    return colors
end

distinguishable_colors(n::Integer, seed::Color; kwargs...) = distinguishable_colors(n, [seed]; kwargs...)
distinguishable_colors(n::Integer; kwargs...) = distinguishable_colors(n, Vector{RGB{N0f8}}(); kwargs...)

@deprecate distinguishable_colors(n::Integer,
                                transform::Function,
                                seed::Color,
                                ls::Vector{Float64},
                                cs::Vector{Float64},
                                hs::Vector{Float64})    distinguishable_colors(n, [seed], transform = transform, lchoices = ls, cchoices = cs, hchoices = hs)

"""
    sequential_palette(h, N::Int=100; <keyword arguments>)

Implements the color palette creation technique by
[Wijffelaars, M., et al. (2008)](https://dl.acm.org/doi/10.1111/j.1467-8659.2008.01203.x).

Colormaps are formed using Bézier curves in LCHuv colorspace
with some constant hue. In addition, start and end points can be given
that are then blended to the original hue smoothly.

# Arguments

- `N`        - number of colors
- `h`        - the main hue [0,360]
- `c`        - the overall lightness contrast [0,1]
- `s`        - saturation [0,1]
- `b`        - brightness [0,1]
- `w`        - cold/warm parameter, i.e. the strength of the starting color [0,1]
- `d`        - depth of the ending color [0,1]
- `wcolor`   - starting color (warmness)
- `dcolor`   - ending color (depth)
- `logscale` - `true`/`false` for toggling logspacing
"""
function sequential_palette(h,
                            N::Int=100;
                            w=0.15,
                            d=0.0,
                            c=0.88,
                            s=0.6,
                            b=0.75,
                            wcolor=RGB(1.0, 1.0, 0.0),
                            dcolor=RGB(0.0, 0.0, 1.0),
                            logscale::Bool=false)
    F = Float64
    return _sequential_palette(N, logscale, F(h), F(w), F(d), F(c), F(s), F(b),
                               RGB{F}(wcolor), RGB{F}(dcolor))
end

mix_hue(h0::Float64, h1::Float64, w) = h0 + w * (normalize_hue(180.0 + h1 - h0) - 180.0)
mix_linearly(a::Float64, b::Float64, w) = (1.0 - w) * a + w * b
function mix_linearly(a::LCHuv{Float64}, b::LCHuv{Float64}, w)
    l = mix_linearly(a.l, b.l, w)
    c = mix_linearly(a.c, b.c, w)
    h = mix_hue(a.h, b.h, w)
    LCHuv{Float64}(l, c, h) # without hue normalization
end

function _sequential_palette(N, logscale, h, w, d, c, s, b, wcolor, dcolor)
    function term_point(pb, l1, l2, weight)
        pl = mix_linearly(l1, l2, weight)
        ph = mix_hue(h, pb.h, weight)
        mc = find_maximum_chroma(LCHuv{Float64}(pl, 0.0, ph))
        pc = min(mc, weight * s * pb.c)
        LCHuv{Float64}(pl, pc, ph)
    end

    pstart0 = convert(LCHuv{Float64}, wcolor)
    pend0   = convert(LCHuv{Float64}, dcolor)
    # Modify the hue of gray
    pstart = LCHuv{Float64}(pstart0.l, pstart0.c, pstart0.c < 1e-4 ? h : pstart0.h)
    pend   = LCHuv{Float64}(pend0.l,   pend0.c,   pend0.c   < 1e-4 ? h : pend0.h)

    p2 = term_point(pstart, 100.0, pstart.l, w) # multi-hue start point
    p1 = MSC(h)
    p0 = term_point(pend, 0.0, 20.0, d) # multi-hue ending point

    q0 = mix_linearly(p0, p1, s)
    q2 = mix_linearly(p2, p1, s)
    q1 = mix_linearly(q0, q2, 0.5)

    pal = Vector{RGB{Float64}}(undef, N)

    step = 1.0 / (N - 1.0)
    for i = 0:N-1
        if logscale
            u = 1.0 - exp10(2.0 * i * step - 2.0)
        else
            u = 1.0 - i * step
        end

        # Change grid to favor light colors and to be uniform along the curve
        lt = 125.0 - 125.0 * 0.2^mix_linearly(b, u, c)
        tt = inv_bezier(lt, p0.l, p2.l, q0.l, q1.l, q2.l)

        # Get color components from Bezier curves
        ll = bezier(tt, p0.l, p2.l, q0.l, q1.l, q2.l)
        cc = bezier(tt, p0.c, p2.c, q0.c, q1.c, q2.c)
        hh = bezier(tt, p0.h, p2.h, q0.h, q1.h, q2.h)

        @inbounds pal[i + 1] = convert(RGB{Float64}, LCHuv{Float64}(ll, cc, hh))
    end
    pal
end

"""
    diverging_palette(h1, h2, N::Int=100; <keyword arguments>)

Create diverging palettes by combining 2 sequential palettes.

# Arguments

- `N`        - number of colors
- `mid`      - the position of the midpoint (0,1)
- `h1`       - the main hue of the first part [0,360]
- `h2`       - the main hue of the latter part [0,360]
- `c`        - the overall lightness contrast [0,1]
- `s`        - saturation [0,1]
- `b`        - brightness [0,1]
- `w`        - cold/warm parameter, i.e. the strength of the starting color [0,1]
- `d1`       - depth of the ending color in the first part [0,1]
- `d2`       - depth of the ending color in the latter part [0,1]
- `wcolor`   - starting color, i.e. the middle color (warmness)
- `dcolor1`  - ending color of the first part (depth)
- `dcolor2`  - ending color of the latter part (depth)
- `logscale` - `true`/`false` for toggling logspacing
"""
function diverging_palette(h1,
                           h2,
                           N::Int=100;
                           mid=0.5,
                           w=0.15,
                           d1=0.0,
                           d2=0.0,
                           c=0.88,
                           s=0.6,
                           b=0.75,
                           wcolor =RGB(1.0, 1.0, 0.0),
                           dcolor1=RGB(1.0, 0.0, 0.0),
                           dcolor2=RGB(0.0, 0.0, 1.0),
                           logscale::Bool=false)
    F = Float64
    return _diverging_palette(N, F(mid), logscale,
                              F(h1), F(h2), F(w), F(d1), F(d2), F(c), F(s), F(b),
                              RGB{F}(wcolor), RGB{F}(dcolor1), RGB{F}(dcolor2))
end
function _diverging_palette(N, mid, logscale, h1, h2, w, d1, d2, c, s, b, wcolor, dcolor1, dcolor2)
    n = isodd(N) ? (N - 1) : N
    N1 = max(ceil(Int, mid * n), 1)
    N2 = max(n - N1, 1)

    pal1 = _sequential_palette(N1 + 1, logscale, h1, w, d1, c, s, b, wcolor, dcolor1)
    pal2 = _sequential_palette(N2 + 1, logscale, h2, w, d2, c, s, b, wcolor, dcolor2)

    pal = Vector{RGB{Float64}}(undef, N1 + N2 + isodd(N))
    # concatenation without `vcat`, which is not precompilable
    @inbounds pal2[1] = weighted_color_mean(0.5, pal1[1], pal2[1])
    @inbounds for i in eachindex(pal)
        pal[i] = i <= N1 ? pal1[end - i + 1] : pal2[i - N1 + iseven(N)]
    end
    pal
end

# Colormap
# ----------------------
# Main function to handle different predefined colormaps
#
"""
    colormap(cname, N=100; logscale=false, <keyword arguments>)

Returns a predefined sequential or diverging colormap computed using
the algorithm by Wijffelaars, M., et al. (2008).

Sequential colormaps `cname` choices are:

- `"Blues"`
- `"Grays"`
- `"Greens"`
- `"Oranges"`
- `"Purples"`
- `"Reds"`

Diverging colormap choices are `"RdBu"`.

Optionally, you can specify the number of colors `N` (default 100).

Extra control is provided by keyword arguments.
- `mid` sets the position of the midpoint for diverging colormaps.
- `logscale=true` uses logarithmically-spaced steps in the colormap.
You can also use keyword argument names that match the argument names in
[`sequential_palette`](@ref) or [`diverging_palette`](@ref).
"""
function colormap(cname::AbstractString, N::Integer=100; kvs...)
    _colormap(String(cname), Int(N), kvs)
end

function _colormap(cname::String, N::Int, kvs)
    logscale = get(kvs, :logscale, false)::Bool
    lcname = lowercase(cname)
    F = Float64
    if haskey(colormaps_sequential, lcname)
        keys_s = (:h, :w, :d, :c, :s, :b, :wcolor, :dcolor)
        Ts = Tuple{F,  F,  F,  F,  F,  F,  RGB{F},  RGB{F}}
        pbs = colormaps_sequential[lcname]::Ts
        for k in keys(kvs)
            k === :logscale && continue
            k in keys_s || throw(ArgumentError("Unknown keyword argument: $k"))
        end
        ps(i) = oftype(pbs[i], get(kvs, keys_s[i], pbs[i]))
        return _sequential_palette(N, logscale, (ps(i) for i in eachindex(pbs))...)
    end
    if haskey(colormaps_diverging, lcname)
        keys_d = (:h1, :h2, :w, :d1, :d2, :c, :s, :b, :wcolor, :dcolor1, :dcolor2)
        Td = Tuple{F,   F,   F,  F,   F,   F,  F,  F,  RGB{F},  RGB{F},   RGB{F}}
        pbd = colormaps_diverging[lcname]::Td
        for k in keys(kvs)
            k === :logscale && continue
            k === :mid && continue
            k in keys_d || throw(ArgumentError("Unknown keyword argument: $k"))
        end
        mid = Float64(get(kvs, :mid, 0.5))
        pd(i) = oftype(pbd[i], get(kvs, keys_d[i], pbd[i]))
        return _diverging_palette(N, mid, logscale, (pd(i) for i in eachindex(pbd))...)
    end
    throw(ArgumentError(string("Unknown colormap: ", cname)))
end


text:
function Give_An_NPI_area(p :: AnArea_Regular ;
                          lon_west :: Float64 = 160.0, lon_east :: Float64 = 360-140.0,
                          lat_south :: Float64 = 30.0, lat_north :: Float64 = 65.0)
  NPI_area = A_selected_area(p, lon_west,lon_east, lat_south, lat_north )
end

function Give_An_PDO_area(p :: AnArea_Regular;
                          lon_west :: Float64 = 120.0, lon_east :: Float64 = 360-120.0,
                          lat_south :: Float64 = 30.0, lat_north :: Float64 = 65.0)
  PDO_area = A_selected_area(p, lon_west,lon_east, lat_south, lat_north )
end

export Give_An_NPI_area
export Give_An_PDO_area

function Give_An_Nino3_area( p :: AnArea_Regular;
                          lon_west :: Float64 = 360 - 150.0, lon_east :: Float64 = 360-90.0,
                          lat_south :: Float64 = -5.0, lat_north :: Float64 = 5.0)
  return  A_selected_area(p, lon_west,lon_east, lat_south, lat_north )
end

function Give_An_Nino4_area( p :: AnArea_Regular;
                          lon_west :: Float64 = 160.0, lon_east :: Float64 = 360-150.0,
                          lat_south :: Float64 = -5.0, lat_north :: Float64 = 5.0)
  return  A_selected_area(p, lon_west,lon_east, lat_south, lat_north )
end

function Give_An_Nino34_area( p :: AnArea_Regular;
                          lon_west :: Float64 = 170.0, lon_east :: Float64 = 360-120.0,
                          lat_south :: Float64 = -5.0, lat_north :: Float64 = 5.0)
  return  A_selected_area(p, lon_west,lon_east, lat_south, lat_north )
end

export Give_An_Nino3_area,Give_An_Nino4_area, Give_An_Nino34_area

function Give_An_ENSO_area(p :: AnArea_Regular;
                          lon_west :: Float64 = 120.0, lon_east :: Float64 = 300.0,
                          lat_south :: Float64 = -20.0, lat_north :: Float64 = 20.0)
  return  A_selected_area(p, lon_west,lon_east, lat_south, lat_north )
end
export Give_An_ENSO_area


text:
using DataArrays
using NullableArrays

srand(1)
M1 = rand(Bool, 5_000_000)
M2 = rand(Bool, 5_000_000, 2)
A1 = rand(5_000_000)
A2 = rand(Float64, 5_000_000, 2)
B1 = rand(Bool, 5_000_000)
B2 = rand(Bool, 5_000_000, 2)
C1 = rand(1:10, 5_000_000)
C2 = rand(Int, 5_000_000, 2)
L = Array(Float64, 5_000_000, 2)

U = NullableArray(Float64, 5_000_000, 2)
X1 = NullableArray(A1)
X2 = NullableArray(A2)
Y1 = NullableArray(A1, M1)
Y2 = NullableArray(A2, M2)
Z1 = NullableArray(B1)
Z2 = NullableArray(B2)
Q1 = NullableArray(C1)
Q2 = NullableArray(C2)
V1 = NullableArray(C1, M1)
V2 = NullableArray(C2, M2)

D = DataArray(Float64, 5_000_000, 2)
E1 = DataArray(A1)
E2 = DataArray(A2)
F1 = DataArray(A1, M1)
F2 = DataArray(A2, M2)
G1 = DataArray(B1)
G2 = DataArray(B2)
H1 = DataArray(C1)
H2 = DataArray(C2)
I1 = DataArray(C1, M1)
I2 = DataArray(C2, M2)

f(x, y) = x * y

function profile_broadcast(f, L, U, D, A1, A2, X1, X2, E1, E2, Y1, Y2, F1, F2)
    println("f(x, y) := x * y")
    println("Method: broadcast!(f, dest, A1, A2) (no empty entries):")
    broadcast!(f, L, A1, A2)
    print("  For Array{Float64}:          ")
    @time(broadcast!(f, L, A1, A2))

    broadcast!(f, U, X1, X2)
    print("  For NullableArray{Float64}:  ")
    @time(broadcast!(f, U, X1, X2))

    broadcast!(f, D, E1, E2)
    print("  For DataArray{Float64}:      ")
    @time(broadcast!(f, D, E1, E2))

    println()
    println("Method: broadcast!(f, dest, A1, A2) (~half empty entries):")
    broadcast!(f, U, Y1, Y2)
    print("  For NullableArray{Float64}:  ")
    @time(broadcast!(f, U, Y1, Y2))

    print("  For DataArray{Float64}:      ")
    broadcast!(f, D, F1, F2)
    @time(broadcast!(f, D, F1, F2))
    nothing
end

function profile_ops_nonulls(A1, A2, X1, X2, E1, E2)
    for op in (
        :(.+),
        :(.-),
        :(.*),
        :(./),
        :(.%),
        :(.^),
    )
        _op = symbol("$op")
        println("Method: $_op (no empty entries)")
        @eval begin
            $_op(A1, A2)
            $_op(X1, X2)
            $_op(E1, E2)
            print("  For Array{Float64}:          ")
            @time($_op(A1, A2))
            print("  For NullableArray{Float64}:  ")
            @time($_op(X1, X2))
            print("  For DataArray{Float64}:      ")
            @time($_op(E1, E2))
        end
    end

    for op in (
        :(.==),
        :(.!=),
        :(.<),
        :(.>),
        :(.<=),
        :(.>=),
    )
        _op = symbol("$op")
        println("Method: $_op (no empty entries)")
        @eval begin
            $_op(A1, A2)
            $_op(X1, X2)
            $_op(E1, E2)
            print("  For Array{Float64}:          ")
            @time($_op(A1, A2))
            print("  For NullableArray{Float64}:  ")
            @time($_op(X1, X2))
            print("  For DataArray{Float64}:      ")
            @time($_op(E1, E2))
        end
    end
    nothing
end

function profile_broadcasted_right(C1, C2, V1, V2, H1, H2, Q1, Q2)
    println("Method: .>> (no empty entries)")
    .>>(C2, C1)
    .>>(V2, V1)
    .>>(H2, H1)
    print("  For Array{Float64}:          ")
    @time(.>>(C2, C1))
    print("  For NullableArray{Float64}:  ")
    @time(.>>(V2, V1))
    print("  For DataArray{Float64}:      ")
    @time(.>>(H2, H1))

    println("Method: .>> (~half empty entries)")
    .>>(C2, C1)
    .>>(Q2, Q1)
    # .>>(I2, I1)
    print("  For Array{Float64}:          ")
    @time(.>>(C2, C1))
    print("  For NullableArray{Float64}:  ")
    @time(.>>(Q2, Q1))
    # print("For DataArray{Float64}:      ")
    # @time(.>>(I2, I1))
    nothing
end

function profile_ops_halfnulls(A1, A2, Y1, Y2, F1, F2)
    for op in (
        :(.+),
        :(.-),
        :(.*),
        :(./),
        :(.%),
        :(.^),
    )
        _op = symbol("$op")
        println("Method: $_op (~half empty entries)")
        @eval begin
            $_op(A1, A2)
            $_op(Y1, Y2)
            $_op(F1, F2)
            print("  For Array{Float64}:          ")
            @time($_op(A1, A2))
            print("  For NullableArray{Float64}:  ")
            @time($_op(Y1, Y2))
            print("  For DataArray{Float64}:      ")
            @time($_op(F1, F2))
        end
    end

    for op in (
        :(.==),
        :(.!=),
        :(.<),
        :(.>),
        :(.<=),
        :(.>=),
    )
        _op = symbol("$op")
        println("Method: $_op (~half empty entries)")
        @eval begin
            $_op(A1, A2)
            $_op(Y1, Y2)
            $_op(F1, F2)
            print("  For Array{Float64}:          ")
            @time($_op(A1, A2))
            print("  For NullableArray{Float64}:  ")
            @time($_op(Y1, Y2))
            print("  For DataArray{Float64}:      ")
            @time($_op(F1, F2))
        end
    end
    nothing
end

function profile_broadcast_all()
    profile_broadcast(f, L, U, D, A1, A2, X1, X2, E1, E2, Y1, Y2, F1, F2)
    profile_ops_nonulls(A1, A2, X1, X2, E1, E2)
    profile_ops_halfnulls(A1, A2, Y1, Y2, F1, F2)
    profile_broadcasted_right(C1, C2, V1, V2, H1, H2, Q1, Q2)
end


text:
module HartreeFockModule
export computeEnergyHartreeFock, evaluateSCF
using ..BaseModule
using ..BasisModule
using ..GeometryModule
using ..SpecialMatricesModule
using ..IntegralsModule
using ..ShellModule
using ..LibInt2Module

function computeEnergyHartreeFock(
  density::Matrix,
  kinetic::Matrix,
  nuclearAttraction::Matrix,
  coulomb::Matrix,
  exchange::Matrix,
  interatomicRepulsion::Float64;
  info::Bool=true)

  P = density
  info && println("=================================")
  eT = 2*trace(kinetic*P)
  info && println("Kinetic Energy: $eT")
  eV = 2*trace(nuclearAttraction*P)
  info && println("Nuclear Attraction Energy: $eV")
  info && println("Sum 1-Electron Energy: $(eT+eV)")
  info && println("---------------------------------")
  # G = 2J-K
  eJ = 2*trace(coulomb*P)
  info && println("Coulomb Energy: $eJ")
  eK = -trace(exchange*P)
  info && println("Exchange Energy: $eK")
  info && println("Total Electronic Energy: $(eT+eV+eJ+eK)")
  info && println("Total Energy: $(eT+eV+eJ+eK+interatomicRepulsion)")
  info && println("=================================")

  return eT+eV+eJ+eK
end

function computeEnergyHartreeFock(
  basis::GaussianBasis,
  geometry::Geometry,
  density::Matrix;
  info::Bool=true)

  P = density
  T = computeMatrixKinetic(basis)
  V = computeMatrixNuclearAttraction(basis,geometry)
  J = computeMatrixCoulomb(basis,density)
  K = computeMatrixExchange(basis,density)
  Vnn = computeEnergyInteratomicRepulsion(geometry)
  computeEnergyHartreeFock(P,T,V,J,K,Vnn;info=info)
end

function computeEnergyHartreeFock(
  density::Matrix,
  matrixKinetic::Matrix,
  matrixNuclearAttraction::Matrix,
  electronRepulsionIntegrals::Array{Float64,4},
  energyInteratomicRepulsion::Float64;
  info::Bool=true,
  # inceptions
  electronRepulsionIntegralsExchange::Array{Float64,4}=electronRepulsionIntegrals)

  J = computeMatrixCoulomb(electronRepulsionIntegrals,density)
  K = computeMatrixExchange(electronRepulsionIntegralsExchange,density)
  computeEnergyHartreeFock(density, matrixKinetic, matrixNuclearAttraction, J, K, energyInteratomicRepulsion;info=info)
end

function evaluateSCFStep(
  initialGuessDensity::Matrix,
  fock::Matrix,
  overlap::Matrix,
  electronNumber::Integer)

  eig = eigfact(Symmetric(fock),Symmetric(overlap)) # because Symmetric is faster and allows for sorted eigenvalues
  C = eig.vectors
  moEnergies = eig.values
  Cocc = eig.vectors[:,1:electronNumber]
  P = Cocc * Cocc.'

  return (moEnergies,P)
end

function evaluateSCF(
  initialGuessDensity::Matrix,
  overlap::Matrix,
  kinetic::Matrix,
  nuclearAttraction::Matrix,
  coulomb::Function,
  exchange::Function,
  interatomicRepulsion::Float64,
  electronNumber::Integer;
  detailedinfo::Bool=true,
  info::Bool=true,
  energyConvergenceCriterion::Float64 = 1e-8)

  P = initialGuessDensity
  energyConverged = false
  local totalEnergy = computeEnergyHartreeFock(P,kinetic,nuclearAttraction,coulomb(P),exchange(P),interatomicRepulsion;info=false)
  local energies::Array{Float64,1}
  local P::Matrix
  i = 0
  while (!energyConverged)
    i+=1
    J = coulomb(P)
    K = exchange(P)
    F = kinetic+nuclearAttraction+2J-K
    energies,P = evaluateSCFStep(P,F,overlap,electronNumber)
    oldEnergy = totalEnergy
    totalEnergy = computeEnergyHartreeFock(P,kinetic,nuclearAttraction,J,K,interatomicRepulsion,info=detailedinfo)
    info && println("E[$i]: $totalEnergy")
    if (abs(totalEnergy-oldEnergy) < energyConvergenceCriterion)
      info && println("Converged!")
      energyConverged = true
    end
  end
  return energies, totalEnergy, P
end

function evaluateSCF(
  basis::GaussianBasis,
  geometry::Geometry,
  initialGuessDensity::Matrix,
  electronNumber::Integer;
  detailedinfo::Bool=true,
  info::Bool=true,
  energyConvergenceCriterion::Float64 = 1e-8,
  computeTensorElectronRepulsionIntegralsCoulomb = computeTensorElectronRepulsionIntegrals,
  computeTensorElectronRepulsionIntegralsExchange = computeTensorElectronRepulsionIntegrals)

  S = computeMatrixOverlap(basis)
  T = computeMatrixKinetic(basis)
  V = computeMatrixNuclearAttraction(basis,geometry)
  Vnn = computeEnergyInteratomicRepulsion(geometry)
  ERIsCoulomb = computeTensorElectronRepulsionIntegralsCoulomb(basis)
  ERIsExchange = computeTensorElectronRepulsionIntegralsCoulomb==computeTensorElectronRepulsionIntegralsExchange ? ERIsCoulomb : computeTensorElectronRepulsionIntegralsExchange(basis)

  evaluateSCF(initialGuessDensity, S, T, V, P->computeMatrixCoulomb(ERIsCoulomb,P), P->computeMatrixExchange(ERIsExchange,P), Vnn, electronNumber; energyConvergenceCriterion=energyConvergenceCriterion, info=info, detailedinfo=detailedinfo)
end

function evaluateSCF(
  shells::Vector{Shell},
  geometry::Geometry,
  initialGuessDensity::Matrix,
  electronNumber::Integer;
  detailedinfo::Bool=true,
  info::Bool=true,
  energyConvergenceCriterion::Float64 = 1e-8)

  S = computeMatrixOverlap(shells)
  T = computeMatrixKinetic(shells)
  V = computeMatrixNuclearAttraction(shells,geometry)
  Vnn = computeEnergyInteratomicRepulsion(geometry)

  #too slow because trivial ERI recomputation in every step:
  evaluateSCF(initialGuessDensity, S, T, V, P->computeMatrixCoulomb(shells,P), P->computeMatrixExchange(shells,P), Vnn, electronNumber; energyConvergenceCriterion=energyConvergenceCriterion, info=info, detailedinfo=detailedinfo)
  #this would be faster but doesn't allow any shell-structure reuse as we might later want (with sparsity):
  #bas = Basis([])
  #for sh in shells
  #  append!(bas,expandShell(sh))
  #end
  #ERIs = computeTensorElectronRepulsionIntegrals(bas)
  #evaluateSCF(initialGuessDensity, S, T, V, P->computeMatrixCoulomb(ERIs,P), P->computeMatrixExchange(ERIs,P), Vnn, electronNumber; energyConvergenceCriterion=energyConvergenceCriterion)
end

function evaluateSCF(
  shells::Vector{LibInt2Shell},
  geometry::Geometry,
  initialGuessDensity::Matrix,
  electronNumber::Integer;
  detailedinfo::Bool=true,
  info::Bool=true,
  energyConvergenceCriterion::Float64 = 1e-8)

  S = computeMatrixOverlap(shells)
  T = computeMatrixKinetic(shells)
  V = computeMatrixNuclearAttraction(shells,geometry)
  Vnn = computeEnergyInteratomicRepulsion(geometry)

  evaluateSCF(initialGuessDensity, S, T, V, P->computeMatrixCoulomb(shells,P), P->computeMatrixExchange(shells,P), Vnn, electronNumber; energyConvergenceCriterion=energyConvergenceCriterion, info=info, detailedinfo=detailedinfo)
end

end # module


text:
function ftest()
   unpack_expr = Meta.parse(string(join(string.(p), ","), ", = a"))
   G = :(a1*t + a2)
   J = :(a2*t + a1)
   @eval return(TT(
      ((t, a) -> ($unpack_expr; $G)) |> FunctionWrapper{Float64, Tuple{Float64, Vector{Float64}}},
      ((t, a) -> ($unpack_expr; $J)) |> FunctionWrapper{Float64, Tuple{Float64, Vector{Float64}}},

      ((ta, a) -> ($unpack_expr; [$G for t in ta])) |> FunctionWrapper{Vector{Float64}, Tuple{Vector{Float64}, Vector{Float64}}},
      ((ta, a) -> ($unpack_expr; [$J for t in ta])) |> FunctionWrapper{Vector{Float64}, Tuple{Vector{Float64}, Vector{Float64}}},
      ))
end

macro mtest()
   # test base case
   p=(:a1, :a2)
   unpack_expr = Meta.parse(string(join(string.(p), ","), ", = a"))
   G = :(a1*t + a2)
   J = :(a2*t + a1)
   return quote TT(
      ((t, a) -> ($unpack_expr; $G)) |> FunctionWrapper{Float64, Tuple{Float64, Vector{Float64}}},
      ((t, a) -> ($unpack_expr; $J)) |> FunctionWrapper{Float64, Tuple{Float64, Vector{Float64}}},

      ((ta, a) -> ($unpack_expr; [$G for t in ta])) |> FunctionWrapper{Vector{Float64}, Tuple{Vector{Float64}, Vector{Float64}}},
      ((ta, a) -> ($unpack_expr; [$J for t in ta])) |> FunctionWrapper{Vector{Float64}, Tuple{Vector{Float64}, Vector{Float64}}},
     ) end
end

ftest()
@mtest

max_stars_repo_path: 20190606 macrotest.jl

max_stars_repo_name: JuliaRheology/RheoBenchAndTest

max_issues_repo_path: 20190606 macrotest.jl

max_issues_repo_name: JuliaRheology/RheoBenchAndTest

max_forks_repo_path: 20190606 macrotest.jl

max_forks_repo_name: JuliaRheology/RheoBenchAndTest


text:
using HorizonSideRobots
include("Functions.jl")

function paint_inside(robot::Robot)
    w = move_with_path!(robot, West)
    n = move_with_path!(robot, Nord)

    putmarker!(robot)

    while !isborder(robot, Sud)
        side = HorizonSide((isborder(robot, West) * 3) + (isborder(robot, Ost) * 1))
        put_markers_to_border!(robot, side)
        move!(robot, Sud)
        putmarker!(robot)
    end
    side = HorizonSide((isborder(robot, West) * 3) + (isborder(robot, Ost) * 1))
    put_markers_to_border!(robot, side)

    _ = move_with_path!(robot, West)
    __ = move_with_path!(robot, Nord)
    movements!(robot, Ost, w)
    movements!(robot, Sud, n)
end

max_stars_repo_name: Mananhina/mirea-progs

max_issues_repo_name: Mananhina/mirea-progs

max_forks_repo_name: Mananhina/mirea-progs


text:
### FixedDAryTree

"""
A fixed full-d-ary tree, i.e., in a tree each node has either d or 0 children and each level is
full.
"""
struct FullDAryTree <: FiniteFixedTree
    spec::TreeSpec
    height::Int
end

"""
    FullDAryTree(height::Int, d::Int)

Construct a `FullDAryTree` of height `height` and maximum degree `d`.
"""
function FullDAryTree(height::Int, d::Int)
    if d > 1
        return FullDAryTree(TreeSpec((d^(height+1)-1)//(d-1), "Full $d-arry tree (height $height)", d), height)
    else
        return FullDAryTree(TreeSpec(height+1, "Full $d-arry tree (height $height)", d), height)
    end
end

# FixedDAryTree

degrees(tree::FullDAryTree) = tree
eltype(tree::FullDAryTree) = Int
length(tree::FullDAryTree) = size(tree)

function iterate(tree::FullDAryTree)
    stack = Stack{Tuple{Int, Int}}()
    if tree.height == 0
        return (0, stack)
    else
        push!(stack, (0, maxdegree(tree)))
        return (maxdegree(tree), stack)
    end
end

# We have to do a DFS to work out the degrees
function iterate(tree::FullDAryTree, stack)
    if isempty(stack)
        return nothing
    end

    d = maxdegree(tree)

    parent_height, child_left = pop!(stack)

    if child_left > 1
        push!(stack, (parent_height, child_left - 1))
    end

    current_height = parent_height += 1
    current_degree = 0

    if current_height < tree.height
        current_degree = d
        push!(stack, (current_height, d))
    end

    return (current_degree, stack)
end


text:
module aI_Calculations

include("scratchPadFile.jl")

export RMS, Decibels

end