text:
function regulargrid(mins, maxs, ns, dz)
	linearindex = (i1, i2)->i2 + ns[2] * (i1 - 1)
	coords = Array{Float64}(undef, 2, prod(ns))
	xs = range(mins[1]; stop=maxs[1], length=ns[1])
	ys = range(mins[2]; stop=maxs[2], length=ns[2])
	dx = xs[2] - xs[1]
	dy = ys[2] - ys[1]
	j = 1
	neighbors = Array{Pair{Int, Int}}(undef, 2 * prod(ns) - ns[1] - ns[2])
	areasoverlengths = Array{Float64}(undef, 2 * prod(ns) - ns[1] - ns[2])
	volumes = Array{Float64}(undef, 0)
	for i1 = 1:ns[1]
		for i2 = 1:ns[2]
			push!(volumes, dx * dy * dz)
			coords[1, linearindex(i1, i2)] = xs[i1]
			coords[2, linearindex(i1, i2)] = ys[i2]
			if i1 < ns[1]
				neighbors[j] = linearindex(i1, i2)=>linearindex(i1 + 1, i2)
				areasoverlengths[j] = dy * dz / dx
				j += 1
			end
			if i2 < ns[2]
				neighbors[j] = linearindex(i1, i2)=>linearindex(i1, i2 + 1)
				areasoverlengths[j] = dx * dz / dy
				j += 1
			end
		end
	end
	ρ = 997.0#kg/m^3
	g = 9.81#m/s^2
	α = 1e-8#Pa^-1 = s^2*m/kg -- pulled out of thin air, but wanting a contribution where the solid is more compressible than the water
	porosity = 0.25#unitless
	β = 45.8e-11#Pa^-1 = s^2*m/kg -- compressiblity of water is taken from http://hyperphysics.phy-astr.gsu.edu/hbase/Tables/compress.html
	specificstorage = fill(ρ * g * (α + porosity * β), size(coords, 2))#ρg(α+nβ) where ρ is fluid density, g, is gravitational constant, α is aquifer compressiblity, n is porosity, β is fluid compressibility
	return coords, neighbors, areasoverlengths, specificstorage, volumes
end

mins = [0, 0]
maxs = [100, 60]
#ns = [13, 7]
ns = [25, 15]
#ns = [51, 31]
#ns = [75, 45]
#ns = [101, 61]
#ns = [201, 121]
#ns = [401, 241]
#ns = [801, 481]
coords, neighbors, areasoverlengths, specificstorage, volumes = regulargrid(mins, maxs, ns, 1.0)
params = Dict(:clay=>(1.58e-4, 0.244, 1.09, 0.178947368), :claysilt=>(1e-2, 0.488, 1.37, 0.073913043))#(K, alpha, N, sr)
inclay(x, z) = x < 65 && z > 25 && z < 30
function setupparams(coords, neighbors)
	global params
	Ks = zeros(length(neighbors))
	#name	Kh	Kz	alpha	N	sr
	#claysilt	1.00E-02	1.00E-02	0.488	1.37	0.073913043
	#clay	1.58E-04	1.58E-04	0.244	1.09	0.178947368
	for (i, (node1, node2)) in enumerate(neighbors)
		x, z = 0.5 * (coords[:, node1] + coords[:, node2])
		if inclay(x, z)
			Ks[i], _, _ = params[:clay]
		else
			Ks[i], _, _ = params[:claysilt]
		end
	end
	dirichletnodes = Int[]
	Qs = zeros(size(coords, 2))
	for i = 1:size(coords, 2)
		if coords[2, i] == mins[2]
			push!(dirichletnodes, i)
		end
		if coords[2, i] == maxs[2]
			Qs[i] = 0.00055
		end
	end
	dirichleths = zeros(size(coords, 2))
	return Ks, dirichletnodes, dirichleths, Qs
end

Ks, dirichletnodes, dirichleths, Qs = setupparams(coords, neighbors)


text:
using DelimitedFiles

data= readdlm("input_day3", '\n')

horizontal = (size(data)[1])*3
lline = length(data[1])

#number of times to repet the pattern
rep = horizontal/lline

#println(horizontal)
#println(lline)
println(rep) #it is a 32, need to automatize it

#repeat the pattern the ammount needed
for i in eachindex(data)
   data[i] = data[i]^32
end

#horizontal index
h_i = 1

#number of trees
n_t = 0

for i in eachindex(data)
  if i < size(data)[1]
     line = data[i+1]
     global h_i += 3
     index = h_i

     if line[h_i] == '#'
        global n_t += 1
     end
  end
end

println(n_t)

ext: jl


text:
## using the PyQtGraph submodule
## Examples from pyqtgraphs plotting.py
## translated into julia style

using PySide
reload(Pkg.dir("PySide", "src", "pyqtgraph.jl"))
using PyQtGraph
using PyCall

w = Widget()
lyt = VBoxLayout(w)
setLayout(w, lyt)

## Graphics Layout Widget make grid of plot devices
win = GraphicsLayoutWidget(w)
addWidget(lyt, win)
set_size(w, 800, 600)

raise(w)

## addPlot adds a plot device in next column over and returns a plot object
p1 = addPlot(win, title="Basic array plotting")
p1.plot(randn(100))

p2 = addPlot(win, title="Multiple curves")
p2.plot(randn(100), pen=(255,0,0))
p2.plot(randn(100) .+ 5, pen=(0, 255, 0))
p2.plot(randn(100) .+ 10, pen=(0,0,255))

p3 = addPlot(win, title="Drawing with points")
p3.plot(randn(100), pen=(200, 200, 200), symbolBrush=(255,0,0), symbolPen="w")

nextRow(win)

p4 = addPlot(win, title="Parametric, grid enabled")
p4.plot(sin(linspace(0, 2pi, 1000)), cos(linspace(0, 2pi, 1000)))
p4.showGrid()                   # defaults x=true, y=true

p5 = addPlot(win, title="Scatter plot, axis labels, log scale")
x = 1e-5 * rand(1000)
y = x*1000 .+ 0.005 * randn(1000)
y -= minimum(y).-1.0
mask = x .> 1e-15
x = x[mask]
y = y[mask]
p5.plot(x, y, pen=nothing, symbol="t", symbolPen=nothing, symbolSize=10, symbolBrush=(100, 100, 255, 50))
p5.setLabel("left", text="Y Axis", units="A")
p5.setLabel("bottom", text="Y Axis", units="s")
p5.setLogMode(x=true, y=false)

p6 = addPlot(win, title="Updating plot")

curve = p6.plot(pen="y")
data = randn(10000, 10)
ptr = 0
function update()
    global curve, data, ptr, p6
    qinvoke(curve, :setData, 1:10_000, data[:, 1 + ptr]) # had to fix this
    if ptr == 0
        p6.enableAutoRange("xy", false)  ## stop auto-scaling after the first data set is plotted
    end
    ptr = (ptr + 1) % 10
end

timer = QtCore[:QTimer]()
qconnect(timer, :timeout, update)
qinvoke(timer, :start, 50)

nextRow(win)

p7 = addPlot(win, title="Filled plot, axis disabled")
y = sin(linspace(0.0, 10, 1000)) .+ randn(1000)*0.1
p7.plot(y, fillLevel=-0.3, brush=(50,50,200,100))
p7.showAxis("bottom", show=false)

p8 = addPlot(win, title="Region Selection")

x2 = linspace(-100, 100, 1000)
data2 = sin(x2) ./ x2
p8.plot(data2, pen=(255,255,255,200))

## no special functions for LinearRegionItem
lr = pyqtgraph.(:LinearRegionItem)([400,700])
qinvoke(lr, :setZValue, -10)
p8.addItem(lr)

p9 = addPlot(win, title="zoom")
p9.plot(data2)

function updatePlot()
    p9.setXRange(lr[:getRegion]()..., padding=0)
end
function updateRegion()
    view_range = qinvoke(p9, [:getViewBox, :viewRange]) # 2x2 array
    x_range = Int[view_range[1,j] for j in 1:2]         # a vector, not just view_range[1,:]'
    lr[:setRegion](x_range)
end
qconnect(lr, :sigRegionChanged, updatePlot)
#qconnect(p9, :sigXRangeChanged, updateRegion) ## this is giving issues. Find out why.
updatePlot()

 #raise(w)


text: ## WARNING: this library does a bad job at being functional-style. it was introduced to demonstrate the abilities of optimization and data structuring in Aetherlang.

max_stars_repo_name: Aetherlang/aetherlang

max_issues_repo_name: gojakuch/aetherlang

max_forks_repo_name: gojakuch/aetherlang


text:
if basename(pwd()) == "aoc"
    cd("2018/19")
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
module TestChaos

using Test

function run_chaos(filename; addenv = ())
    path = joinpath(@__DIR__, filename)
    script = """
    $(Base.load_path_setup_code())
    include($(repr(path)))
    """
    cmd = Base.julia_cmd()
    cmd = `$cmd -e $script`
    cmd = Base.addenv(cmd, addenv...)
    pipe = Pipe()
    cmd = pipeline(cmd; stdin = devnull, stdout = pipe, stderr = pipe)
    process = run(cmd; wait = false)
    close(pipe.in)
    output = read(pipe, String)
    wait(process)
    if !success(process)
        @info "Tests in subprocess failed" filename addenv Text(output)
    end
    @test success(process)
    return output
end

function test_recovery_takemanyto()
    output = run_chaos(
        "exec_chaos_recovery.jl";
        addenv = ("ASYNCFINALIZERSTESTS_CHAOS" => "takemanyto",),
    )
    @test occursin("Unexpected failure in finalizer queue", output)
end

function test_recovery_run_finalizers()
    output = run_chaos(
        "exec_chaos_recovery.jl";
        addenv = ("ASYNCFINALIZERSTESTS_CHAOS" => "run_finalizers",),
    )
    @test occursin("Error from async finalizer executor", output)
end

function test_fallback()
    output = run_chaos("exec_chaos_fallback.jl")
    @test occursin("Switching to fallback", output)
end

end  # module


text:
# Raw memory management

export Mem

module Mem

using ..CUDAdrv
using ..CUDAdrv: @enum_without_prefix, CUstream, CUdevice

using Base: @deprecate_binding

#
# untyped buffers
#

abstract type Buffer end

# expected interface:
# - similar()
# - ptr, bytesize and ctx fields
# - convert() to certain pointers

Base.sizeof(buf::Buffer) = buf.bytesize

# ccall integration
#
# taking the pointer of a buffer means returning the underlying pointer,
# and not the pointer of the buffer object itself.
Base.unsafe_convert(T::Type{<:Union{Ptr,CuPtr}}, buf::Buffer) = convert(T, buf)

## device buffer
##
## residing on the GPU

struct DeviceBuffer <: Buffer
    ptr::CuPtr{Cvoid}
    bytesize::Int
    ctx::CuContext
end

Base.similar(buf::DeviceBuffer, ptr::CuPtr{Cvoid}=buf.ptr,
             bytesize::Int=sizeof(buf), ctx::CuContext=buf.ctx) =
    DeviceBuffer(ptr, bytesize, ctx)

Base.convert(::Type{<:Ptr}, buf::DeviceBuffer) =
    throw(ArgumentError("cannot take the CPU address of a GPU buffer"))

Base.convert(::Type{CuPtr{T}}, buf::DeviceBuffer) where {T} =
    convert(CuPtr{T}, buf.ptr)

"""
    alloc(DeviceBuffer, bytesize::Integer)

Allocate `bytesize` bytes of memory on the device. This memory is only accessible on the
GPU, and requires explicit calls to `upload` and `download` for access on the CPU.
"""
function alloc(::Type{DeviceBuffer}, bytesize::Integer)
    bytesize == 0 && return DeviceBuffer(CU_NULL, 0, CuContext(C_NULL))

    ptr_ref = Ref{CUDAdrv.CUdeviceptr}()
    CUDAdrv.cuMemAlloc(ptr_ref, bytesize)

    return DeviceBuffer(reinterpret(CuPtr{Cvoid}, ptr_ref[]), bytesize, CuCurrentContext())
end

@deprecate_binding HOSTALLOC_DEFAULT 0 false
const HOSTALLOC_PORTABLE = CUDAdrv.CU_MEMHOSTALLOC_PORTABLE
const HOSTALLOC_DEVICEMAP = CUDAdrv.CU_MEMHOSTALLOC_DEVICEMAP
const HOSTALLOC_WRITECOMBINED = CUDAdrv.CU_MEMHOSTALLOC_WRITECOMBINED

function free(buf::DeviceBuffer)
    if buf.ptr != CU_NULL
        CUDAdrv.cuMemFree(buf)
    end
end

## host buffer
##
## pinned memory on the CPU, possibly accessible on the GPU

struct HostBuffer <: Buffer
    ptr::Ptr{Cvoid}
    bytesize::Int
    ctx::CuContext

    mapped::Bool
end

Base.similar(buf::HostBuffer, ptr::Ptr{Cvoid}=buf.ptr, bytesize::Int=sizeof(buf),
             ctx::CuContext=buf.ctx, mapped::Bool=buf.mapped) =
    HostBuffer(ptr, bytesize, ctx, mapped)

Base.convert(::Type{Ptr{T}}, buf::HostBuffer) where {T} =
    convert(Ptr{T}, buf.ptr)

function Base.convert(::Type{CuPtr{T}}, buf::HostBuffer) where {T}
    if buf.mapped
        buf.ptr == C_NULL && return convert(CuPtr{T}, CU_NULL)
        ptr_ref = Ref{CuPtr{Cvoid}}()
        CUDAdrv.cuMemHostGetDevicePointer(ptr_ref, buf.ptr, #=flags=# 0)
        convert(CuPtr{T}, ptr_ref[])
    else
        throw(ArgumentError("cannot take the GPU address of a pinned but not mapped CPU buffer"))
    end
end

"""
    alloc(HostBuffer, bytesize::Integer, [flags])

Allocate `bytesize` bytes of page-locked memory on the host. This memory is accessible from
the CPU, and makes it possible to perform faster memory copies to the GPU. Furthermore, if
`flags` is set to `HOSTALLOC_DEVICEMAP` the memory is also accessible from the GPU.
These accesses are direct, and go through the PCI bus.
"""
function alloc(::Type{HostBuffer}, bytesize::Integer, flags=0)
    bytesize == 0 && return HostBuffer(C_NULL, 0, CuContext(C_NULL), false)

    ptr_ref = Ref{Ptr{Cvoid}}()
    CUDAdrv.cuMemHostAlloc(ptr_ref, bytesize, flags)

    mapped = (flags & HOSTALLOC_DEVICEMAP) != 0
    return HostBuffer(ptr_ref[], bytesize, CuCurrentContext(), mapped)
end

@enum_without_prefix CUDAdrv.CUmemAttach_flags CU_MEM_

"""
    register(HostBuffer, ptr::Ptr, bytesize::Integer, [flags])

Page-lock the host memory pointed to by `ptr`. Subsequent transfers to and from devices will
be faster, and can be executed asynchronously. If the `MEMHOSTREGISTER_DEVICEMAP` flag is
specified, the buffer will also be accessible directly from the GPU. These accesses are
direct, and go through the PCI bus.
"""
function register(::Type{HostBuffer}, ptr::Ptr, bytesize::Integer, flags=0)
    bytesize == 0 && throw(ArgumentError())

    CUDAdrv.cuMemHostRegister(ptr, bytesize, flags)

    mapped = (flags & HOSTREGISTER_DEVICEMAP) != 0
    return HostBuffer(ptr, bytesize, CuCurrentContext(), mapped)
end

function unregister(buf::HostBuffer)
    CUDAdrv.cuMemHostUnregister(buf)
end

function free(buf::HostBuffer)
    if buf.ptr != CU_NULL
        CUDAdrv.cuMemFreeHost(buf)
    end
end

## unified buffer
##
## managed buffer that is accessible on both the CPU and GPU

struct UnifiedBuffer <: Buffer
    ptr::CuPtr{Cvoid}
    bytesize::Int
    ctx::CuContext
end

Base.similar(buf::UnifiedBuffer, ptr::CuPtr{Cvoid}=buf.ptr,
             bytesize::Int=sizeof(buf), ctx::CuContext=buf.ctx) =
    UnifiedBuffer(ptr, bytesize, ctx)

Base.convert(::Type{Ptr{T}}, buf::UnifiedBuffer) where {T} =
    convert(Ptr{T}, reinterpret(Ptr{Cvoid}, buf.ptr))

Base.convert(::Type{CuPtr{T}}, buf::UnifiedBuffer) where {T} =
    convert(CuPtr{T}, buf.ptr)

"""
    alloc(UnifiedBuffer, bytesize::Integer, [flags::CUmemAttach_flags])

Allocate `bytesize` bytes of unified memory. This memory is accessible from both the CPU and
GPU, with the CUDA driver automatically copying upon first access.
"""
function alloc(::Type{UnifiedBuffer}, bytesize::Integer,
              flags::CUDAdrv.CUmemAttach_flags=ATTACH_GLOBAL)
    bytesize == 0 && return UnifiedBuffer(CU_NULL, 0, CuContext(C_NULL))

    ptr_ref = Ref{CuPtr{Cvoid}}()
    CUDAdrv.cuMemAllocManaged(ptr_ref, bytesize, flags)

    return UnifiedBuffer(ptr_ref[], bytesize, CuCurrentContext())
end

function free(buf::UnifiedBuffer)
    if buf.ptr != CU_NULL
        CUDAdrv.cuMemFree(buf)
    end
end

const HOSTREGISTER_PORTABLE = CUDAdrv.CU_MEMHOSTREGISTER_PORTABLE
const HOSTREGISTER_DEVICEMAP = CUDAdrv.CU_MEMHOSTREGISTER_DEVICEMAP
const HOSTREGISTER_IOMEMORY = CUDAdrv.CU_MEMHOSTREGISTER_IOMEMORY

function prefetch(buf::UnifiedBuffer, bytes=sizeof(buf);
                  device::CuDevice=device(), stream::CuStream=CuDefaultStream())
    bytes > sizeof(buf) && throw(BoundsError(buf, bytes))
    CUDAdrv.cuMemPrefetchAsync(buf, bytes, device, stream)
end

@enum_without_prefix CUDAdrv.CUmem_advise CU_MEM_

function advise(buf::UnifiedBuffer, advice::CUDAdrv.CUmem_advise, bytes=sizeof(buf),
                device=device(buf.ctx))
    bytes > sizeof(buf) && throw(BoundsError(buf, bytes))
    CUDAdrv.cuMemAdvise(buf, bytes, advice, device)
end

## convenience aliases

const Device  = DeviceBuffer
const Host    = HostBuffer
const Unified = UnifiedBuffer

#
# typed pointers
#

## initialization

"""
    set!(buf::CuPtr, value::Union{UInt8,UInt16,UInt32}, len::Integer;
         async::Bool=false, stream::CuStream)

Initialize device memory by copying `val` for `len` times. Executed asynchronously if
`async` is true, in which case a valid `stream` is required.
"""
set!

for T in [UInt8, UInt16, UInt32]
    bits = 8*sizeof(T)
    fn_sync = Symbol("cuMemsetD$(bits)")
    fn_async = Symbol("cuMemsetD$(bits)Async")
    @eval function set!(ptr::CuPtr{$T}, value::$T, len::Integer;
                        async::Bool=false, stream::Union{Nothing,CuStream}=nothing)
        if async
          stream===nothing &&
              throw(ArgumentError("Asynchronous memory operations require a stream."))
            $(getproperty(CUDAdrv, fn_async))(ptr, value, len, stream)
        else
          stream===nothing ||
              throw(ArgumentError("Synchronous memory operations cannot be issues on a stream."))
            $(getproperty(CUDAdrv, fn_sync))(ptr, value, len)
        end
    end
end

## copy operations

for (f, srcPtrTy, dstPtrTy) in (("cuMemcpyDtoH", CuPtr, Ptr),
                                ("cuMemcpyHtoD", Ptr,   CuPtr),
                                ("cuMemcpyDtoD", CuPtr, CuPtr),
                               )
    @eval function Base.unsafe_copyto!(dst::$dstPtrTy{T}, src::$srcPtrTy{T}, N::Integer;
                                       stream::Union{Nothing,CuStream}=nothing,
                                       async::Bool=false) where T
        if async
            stream===nothing &&
                throw(ArgumentError("Asynchronous memory operations require a stream."))
            $(getproperty(CUDAdrv, Symbol(f * "Async")))(dst, src, N*sizeof(T), stream)
        else
            stream===nothing ||
                throw(ArgumentError("Synchronous memory operations cannot be issued on a stream."))
            $(getproperty(CUDAdrv, Symbol(f)))(dst, src, N*sizeof(T))
        end
        return dst
    end
end

#
# utilities
#

## memory info

function info()
    free_ref = Ref{Csize_t}()
    total_ref = Ref{Csize_t}()
    CUDAdrv.cuMemGetInfo(free_ref, total_ref)
    return convert(Int, free_ref[]), convert(Int, total_ref[])
end

end

"""
    available_memory()

Returns the available_memory amount of memory (in bytes), available for allocation by the CUDA context.
"""
available_memory() = Mem.info()[1]

"""
    total_memory()

Returns the total amount of memory (in bytes), available for allocation by the CUDA context.
"""
total_memory() = Mem.info()[2]


text:
using TensorOperations

include("RCCSDHelper.jl")

function RCCSD(alg::RCCSDa)
    aoints = IntegralHelper{Float64}()
    rhf = RHF(aoints)

    if typeof(aoints.eri_type) === JKFIT || Options.get("precision") == "single"
        aoints = IntegralHelper(eri_type=RIFIT())
    end
    moints = IntegralHelper(orbitals=rhf.orbitals)
    RCCSD(moints, aoints, alg)
end

function RCCSD(moints::IntegralHelper{T,<:AbstractERI,<:AbstractRestrictedOrbitals}, aoints::IntegralHelper{T,<:AbstractERI,AtomicOrbitals}, alg::RCCSDa) where T<:AbstractFloat

    Fermi.Integrals.compute!(moints, aoints, "F")
    Fermi.Integrals.compute!(moints, aoints, "OOOO")
    Fermi.Integrals.compute!(moints, aoints, "OOOV")
    Fermi.Integrals.compute!(moints, aoints, "OOVV")
    Fermi.Integrals.compute!(moints, aoints, "OVOV")
    Fermi.Integrals.compute!(moints, aoints, "OVVV")
    Fermi.Integrals.compute!(moints, aoints, "VVVV")
    RCCSD(moints, alg)
end

function RCCSD(moints::IntegralHelper{T,E1,<:AbstractRestrictedOrbitals}, aoints::IntegralHelper{T,E2,AtomicOrbitals}, alg::RCCSDa) where {T<:AbstractFloat,T2<:AbstractFloat,
                                                                                E1<:AbstractDFERI,E2<:AbstractDFERI}
    Fermi.Integrals.compute!(moints, aoints, "F")
    Fermi.Integrals.compute!(moints, aoints, "BOO")
    Fermi.Integrals.compute!(moints, aoints, "BOV")
    Fermi.Integrals.compute!(moints, aoints, "BVV")
    RCCSD(moints, alg)
end

function RCCSD(moints::IntegralHelper{T,E,O}, alg::RCCSDa) where {T<:AbstractFloat,E<:AbstractERI,O<:AbstractRestrictedOrbitals}

    # Create zeroed guesses for amplitudes

    o = moints.molecule.Nα - Options.get("drop_occ")
    v = size(moints.orbitals.C,1) - Options.get("drop_vir") - moints.molecule.Nα

    output("Using MP2 guess")
    T1guess = deepcopy(moints["Fia"])
    T2guess = permutedims(moints["OVOV"], (1,3,2,4))

    # Orbital energies line
    if haskey(moints.cache, "D1")
        d = moints["D1"]
    else
        Fd = moints["Fd"]
        ndocc = moints.molecule.Nα
        frozen = Options.get("drop_occ")
        inac = Options.get("drop_vir")
        ϵo = Fd[(1+frozen:ndocc)]
        ϵv = Fd[(1+ndocc):end-inac]

        d = [ϵo[i]-ϵv[a] for i=eachindex(ϵo), a=eachindex(ϵv)]
        moints["D1"] = d
    end

    if haskey(moints.cache, "D2")
        D = moints["D2"]
    else
        Fd = moints["Fd"]
        ndocc = moints.molecule.Nα
        frozen = Options.get("drop_occ")
        inac = Options.get("drop_vir")
        ϵo = Fd[(1+frozen:ndocc)]
        ϵv = Fd[(1+ndocc):end-inac]

        D = [ϵo[i]+ϵo[j]-ϵv[a]-ϵv[b] for i=eachindex(ϵo), j=eachindex(ϵo), a=eachindex(ϵv), b=eachindex(ϵv)]
        moints["D2"] = D
    end

    T1guess ./= d
    T2guess ./= D

    RCCSD(moints, T1guess, T2guess, alg)
end

function RCCSD(moints::IntegralHelper{T,E,O}, newT1::AbstractArray{T,2}, newT2::AbstractArray{T,4},
                    alg::RCCSDa) where {T<:AbstractFloat, E<:AbstractERI, O<:AbstractRestrictedOrbitals}

    # Print intro
    cc_header()
    Eref = moints.orbitals.sd_energy
    ndocc = moints.molecule.Nα
    nvir = size(moints.orbitals.C,1) - ndocc
    core = Options.get("drop_occ")
    inac = Options.get("drop_vir")

    if core ≥ ndocc
        throw(FermiException("invalid number of frozen orbitals ($core) for $ndocc doubly occupied orbitals"))
    end
    if inac ≥ nvir
        throw(FermiException("invalid number of inactive virtual orbitals ($inac) for $nvir total virtual orbitals"))
    end

    # Get relevant Options
    cc_max_iter = Options.get("cc_max_iter")
    cc_e_conv = Options.get("cc_e_conv")
    cc_max_rms = Options.get("cc_max_rms")
    precision_override = Options.get("precision_override")

    dp = Options.get("cc_damp_ratio")
    do_diis = Options.get("cc_diis")
    ndiis = Options.get("cc_ndiis")
    diis_start = Options.get("diis_start")
    diis_relax = Options.get("cc_diis_relax")

    diis_prec = Options.get("diis_prec")
    diis_prec = diis_prec == "double" ? Float64 : diis_prec == "single" ? Float32 : Float16
    do_diis ? DM_T1 = Fermi.DIIS.DIISManager{T, diis_prec}(size=ndiis) : nothing
    do_diis ? DM_T2 = Fermi.DIIS.DIISManager{T, diis_prec}(size=ndiis) : nothing

    output("\tDropped Occupied Orbitals →  {:3.0f}", core)
    output("\tDropped Virtual Orbitals  →  {:3.0f}", inac)
    output("\n"*repeat("-",80))
    output("\tOptions:")
    output("\tPrecision           {}", T)
    output("\tDIIS                {}", do_diis)
    output("\tDIIS Vectors        {}", (do_diis ? ndiis : 0))
    output("\tDIIS Precision      {}", diis_prec)
    output("\tDamping percentage  {}", dp)
    output("\t\tcc_max_iter →  {:3.0d}", cc_max_iter)
    output("\t\tcc_e_conv   →  {:2.0e}", cc_e_conv)
    output("\t\tcc_max_rms  →  {:2.0e}", cc_max_rms)
    if (cc_e_conv < eps(T) || cc_max_rms < eps(T)) && !(precision_override)
        output("\t⚠ WARNING ⚠  cc_e_conv set to less than ϵ ({}) for precision {}", eps(T), T)
        output("\t             CCSD is unlikely to converge to your standards.")
        output("\t             OVERRIDING set convergence criterion.")
        output("\t             Use @set `precision_override` true to perform the computation as entered.")
        cc_e_conv = eps(T)
        cc_max_rms = eps(T)
    end

    # Initialize Loop parameters
    r1 = 1
    r2 = 1
    dE = 1
    rms = 1
    ite = 1
    T1 = deepcopy(newT1)
    T2 = deepcopy(newT2)

    output(repeat("-", 80))

    # Compute Guess Energy
    Ecc = cc_update_energy(newT1, newT2, moints, alg)
    Eguess = Ecc + Eref

    output("\tGuess Correlation Energy:   {:15.10f}", Ecc)
    output("\tGuess Total Energy:         {:15.10f}\n", Eguess)
    output("    Starting CC Iterations\n")

    dE = 1
    rms = 1
    ite = 1

    if do_diis
        DM_T1 = Fermi.DIIS.DIISManager{T, diis_prec}(size=ndiis)
        DM_T2 = Fermi.DIIS.DIISManager{T, diis_prec}(size=ndiis)
    end

    main_time = 0
    output("{:10s}    {: 15s}    {: 12s}    {:12s}    {:10s}    {:s5}","Iteration","CC Energy","ΔE","Max RMS","time (s)", "diis")
    while (abs(dE) > cc_e_conv || rms > cc_max_rms)
        if ite > cc_max_iter
            output("\n⚠  CC Equations did not converge in {:1.0d} iterations.", cc_max_iter)
            break
        end
        t = @elapsed begin

            T1 .= newT1
            T2 .= newT2
            oldE = Ecc

            #debug
            update_amp!(newT1, newT2, T1, T2, moints, alg)

            # Compute residues
            r1 = sqrt(sum((newT1 .- T1).^2)/length(T1))
            r2 = sqrt(sum((newT2 .- T2).^2)/length(T2))

            diis = false
            if do_diis
                # Update DIIS vector
                e1 = (newT1 - T1)
                e2 = (newT2 - T2)
                push!(DM_T1,newT1,e1)
                push!(DM_T2,newT2,e2)
                # Using DIIS at every iteration is BAD! So do it every three steps
                if ite ≥ diis_start && ite % diis_relax == 0
                    # Get new extrapolated amplitudes from DIIS
                    newT2 = Fermi.DIIS.extrapolate(DM_T2;add_res=true)
                    newT1 = Fermi.DIIS.extrapolate(DM_T1;add_res=true)
                    diis = true
                end
            end

            # Apply damping
            if dp > 0
                newT1 .= (1-dp)*newT1 .+ dp*T1
                newT2 .= (1-dp)*newT2 .+ dp*T2
            end

            Ecc = cc_update_energy(newT1, newT2, moints, alg)
        end
        rms = max(r1,r2)
        dE = Ecc - oldE
        main_time += t
        output("    {:<5.0d}    {:<15.10f}    {:>+12.10f}    {:<12.10f}    {:<10.5f}    {:s5}", ite, Ecc, dE, rms, t, diis)
        ite += 1
    end
    output("\nMain CCSD iterations done in {:5.5f} s", main_time)
    output("Average time per iteration {:5.5f}", main_time/(ite-1))

    # Converged?
    conv = false
    if abs(dE) < cc_e_conv && rms < cc_max_rms
        output("\n 🍾 Equations Converged!")
        conv = true
    end

    output(" @Final CCSD Correlation Energy:     {:15.10f}", Ecc)
    output(" @Final CCSD Energy:                 {:15.10f}", Ecc+Eref)
    output(repeat("-",80))

    if Options.get("return_ints")
        return RCCSD(Eguess, Ecc, Ecc+Eref, newT1, newT2, dE, rms), moints
    else
        return RCCSD(Eguess, Ecc, Ecc+Eref, newT1, newT2, dE, rms)
    end
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
function indices_z(s::Simulation)
    model = s.model
    nq = model.nq
    nc = model.nc
    q = collect(1:nq)
    γ = collect(index_γ1(s.model, s.env))
    sγ = collect(1:0)
    ψ = collect(1:0)
    b = collect(index_b1(s.model, s.env))
    sψ = collect(1:0)
    sb = collect(1:0)
    IndicesZ(q, γ, sγ, ψ, b, sψ, sb)
end

function simulator(s, T;
    h=0.01,
    policy=empty_policy(s.model),
    dist=empty_disturbances(s.model),
    f=friction_coefficients(s.model),
    residual=s.res.r!,
    jacobian_z=s.res.rz!,
    jacobian_θ=s.res.rθ!,
    diff_sol=false,
    solver_opts=InteriorPointOptions(
        undercut=Inf,
        γ_reg=0.1,
        r_tol=1e-8,
        κ_tol=1e-8,
        max_ls=25,
        ϵ_min=0.25,
        diff_sol=diff_sol,
        verbose=false),
    stats=SimulatorStatistics(T),
    sim_opts=SimulatorOptions()
    )

    idx_z = indices_z(s)
    idx_θ = indices_θ(s.model, nf=length(f))
    idx_opt = IndicesOptimization(s.model, s.env)
    nz = num_var(s.model, s.env)
    nθ = num_data(s.model)
    z0 = zeros(nz)
    θ0 = zeros(nθ)
    q0 = zeros(s.model.nq)
    initialize_z!(z0, s.model, idx_z, q0)
    initialize_θ!(θ0, s.model, idx_θ, q0, q0, zeros(s.model.nu), zeros(s.model.nw), f, h)

    ip = interior_point(
         z0,
         θ0,
         idx=idx_opt,
         r! = residual,
         rz! = jacobian_z,
         rθ! = jacobian_θ,
         rz=zeros(nz, nz),
         rθ=zeros(nz, nθ),
         opts=solver_opts)

    traj = Trajectory(s.model, T, nb=s.model.nc * friction_dim(s.env))
    grad = GradientTrajectory(s.model, T, nb=s.model.nc * friction_dim(s.env))

    return Simulator(s.model, policy, dist, traj, grad, ip, idx_z, idx_θ, f, h, stats, sim_opts)
end

"""
    control saturation
"""
control_saturation(u, uL, uU) = min.(max.(uL, u), uU)


text:
using Mocking

create_http_namespace_1_patch = @patch function create_http_namespace(a...; b...)
    return Dict{String, Any}("OperationId" => "httpvoqozuhfet5kzxoxg-a-response-example")
end

create_private_dns_namespace_1_patch = @patch function create_private_dns_namespace(a...; b...)
    return Dict{String, Any}("OperationId" => "gv4g5meo7ndmeh4fqskygvk23d2fijwa-k9302yzd")
end

create_public_dns_namespace_1_patch = @patch function create_public_dns_namespace(a...; b...)
    return Dict{String, Any}("OperationId" => "dns2voqozuhfet5kzxoxg-a-response-example")
end

create_service_1_patch = @patch function create_service(a...; b...)
    return Dict{String, Any}("Service" => Dict{String, Any}("Arn" => "arn:aws:servicediscovery:us-west-2:123456789012:service/srv-p5zdwlg5uvvzjita", "Id" => "srv-p5zdwlg5uvvzjita", "DnsConfig" => Dict{String, Any}("NamespaceId" => "ns-ylexjili4cdxy3xm", "RoutingPolicy" => "MULTIVALUE", "DnsRecords" => Any[Dict{String, Any}("TTL" => 60, "Type" => "A")]), "NamespaceId" => "ns-ylexjili4cdxy3xm", "CreateDate" => 1.587081768334e9, "Name" => "myservice", "CreatorRequestId" => "567c1193-6b00-4308-bd57-ad38a8822d25"))
end

delete_namespace_1_patch = @patch function delete_namespace(a...; b...)
    return Dict{String, Any}("OperationId" => "gv4g5meo7ndmeh4fqskygvk23d2fijwa-k98y6drk")
end

delete_service_1_patch = @patch function delete_service(a...; b...)
    return Dict{String, Any}()
end

deregister_instance_1_patch = @patch function deregister_instance(a...; b...)
    return Dict{String, Any}("OperationId" => "4yejorelbukcjzpnr6tlmrghsjwpngf4-k98rnaiq")
end

discover_instances_1_patch = @patch function discover_instances(a...; b...)
    return Dict{String, Any}("Instances" => Any[Dict{String, Any}("NamespaceName" => "example.com", "HealthStatus" => "UNKNOWN", "Attributes" => Dict{String, Any}("AWS_INSTANCE_IPV4" => "172.2.1.3", "AWS_INSTANCE_PORT" => "808"), "InstanceId" => "myservice-53", "ServiceName" => "myservice")])
end

get_instance_1_patch = @patch function get_instance(a...; b...)
    return Dict{String, Any}("Instance" => Dict{String, Any}("Id" => "i-abcd1234", "Attributes" => Dict{String, Any}("AWS_INSTANCE_IPV4" => "192.0.2.44", "color" => "green", "stage" => "beta", "region" => "us-west-2", "AWS_INSTANCE_PORT" => "80")))
end

get_instances_health_status_1_patch = @patch function get_instances_health_status(a...; b...)
    return Dict{String, Any}("Status" => Dict{String, Any}("i-abcd1234" => "HEALTHY", "i-abcd1235" => "UNHEALTHY"))
end

get_namespace_1_patch = @patch function get_namespace(a...; b...)
    return Dict{String, Any}("Namespace" => Dict{String, Any}("Arn" => "arn:aws:servicediscovery:us-west-2: 123456789120:namespace/ns-e1tpmexample0001", "Id" => "ns-e1tpmexample0001", "Description" => "Example.com AWS Cloud Map HTTP Namespace", "Properties" => Dict{String, Any}("HttpProperties" => Dict{String, Any}("HttpName" => "example-http.com"), "DnsProperties" => Dict{String, Any}()), "Type" => "HTTP", "CreateDate" => "20181118T211712Z", "Name" => "example-http.com", "CreatorRequestId" => "example-creator-request-id-0001"))
end

get_operation_1_patch = @patch function get_operation(a...; b...)
    return Dict{String, Any}("Operation" => Dict{String, Any}("Id" => "gv4g5meo7ndmeh4fqskygvk23d2fijwa-k9302yzd", "Status" => "SUCCESS", "Type" => "CREATE_NAMESPACE", "UpdateDate" => 1.587055900469e9, "CreateDate" => 1.587055860121e9, "Targets" => Dict{String, Any}("NAMESPACE" => "ns-ylexjili4cdxy3xm")))
end

get_service_1_patch = @patch function get_service(a...; b...)
    return Dict{String, Any}("Service" => Dict{String, Any}("Arn" => "arn:aws:servicediscovery:us-west-2:123456789120:service/srv-e4anhexample0004", "HealthCheckConfig" => Dict{String, Any}("ResourcePath" => "/", "Type" => "HTTPS", "FailureThreshold" => 3), "Description" => "Example.com AWS Cloud Map HTTP Service", "Id" => "srv-e4anhexample0004", "NamespaceId" => "ns-e4anhexample0004", "CreateDate" => "20181118T211707Z", "Name" => "example-http-service", "CreatorRequestId" => "example-creator-request-id-0004"))
end

list_instances_1_patch = @patch function list_instances(a...; b...)
    return Dict{String, Any}("Instances" => Any[Dict{String, Any}("Id" => "i-06bdabbae60f65a4e", "Attributes" => Dict{String, Any}("AWS_INSTANCE_IPV4" => "172.2.1.3", "AWS_INSTANCE_PORT" => "808"))])
end

list_namespaces_1_patch = @patch function list_namespaces(a...; b...)
    return Dict{String, Any}("Namespaces" => Any[Dict{String, Any}("Arn" => "arn:aws:servicediscovery:us-west-2:123456789012:namespace/ns-a3ccy2e7e3a7rile", "Id" => "ns-a3ccy2e7e3a7rile", "Type" => "DNS_PRIVATE", "CreateDate" => 1.585354387357e9, "Name" => "local", "Properties" => Dict{String, Any}("HttpProperties" => Dict{String, Any}("HttpName" => "local"), "DnsProperties" => Dict{String, Any}("HostedZoneId" => "Z06752353VBUDTC32S84S"))), Dict{String, Any}("Arn" => "arn:aws:servicediscovery:us-west-2:123456789012:namespace/ns-pocfyjtrsmwtvcxx", "Id" => "ns-pocfyjtrsmwtvcxx", "Description" => "My second namespace", "Type" => "HTTP", "CreateDate" => 1.586468974698e9, "Name" => "My-second-namespace", "Properties" => Dict{String, Any}("HttpProperties" => Dict{String, Any}("HttpName" => "My-second-namespace"), "DnsProperties" => Dict{String, Any}())), Dict{String, Any}("Arn" => "arn:aws:servicediscovery:us-west-2:123456789012:namespace/ns-ylexjili4cdxy3xm", "Id" => "ns-ylexjili4cdxy3xm", "Type" => "DNS_PRIVATE", "CreateDate" => 1.587055896798e9, "Name" => "example.com", "Properties" => Dict{String, Any}("HttpProperties" => Dict{String, Any}("HttpName" => "example.com"), "DnsProperties" => Dict{String, Any}("HostedZoneId" => "Z09983722P0QME1B3KC8I")))])
end

list_operations_1_patch = @patch function list_operations(a...; b...)
    return Dict{String, Any}("Operations" => Any[Dict{String, Any}("Id" => "76yy8ovhpdz0plmjzbsnqgnrqvpv2qdt-kexample", "Status" => "SUCCESS"), Dict{String, Any}("Id" => "prysnyzpji3u2ciy45nke83x2zanl7yk-dexample", "Status" => "SUCCESS"), Dict{String, Any}("Id" => "ko4ekftir7kzlbechsh7xvcdgcpk66gh-7example", "Status" => "PENDING")])
end

list_services_1_patch = @patch function list_services(a...; b...)
    return Dict{String, Any}("Services" => Any[Dict{String, Any}("Arn" => "arn:aws:servicediscovery:us-west-2:123456789012:service/srv-p5zdwlg5uvvzjita", "Id" => "srv-p5zdwlg5uvvzjita", "DnsConfig" => Dict{String, Any}("RoutingPolicy" => "MULTIVALUE", "DnsRecords" => Any[Dict{String, Any}("TTL" => 60, "Type" => "A")]), "CreateDate" => 1.587081768334e9, "Name" => "myservice")])
end

list_tags_for_resource_1_patch = @patch function list_tags_for_resource(a...; b...)
    return Dict{String, Any}("Tags" => Any[Dict{String, Any}("Value" => "Zeta", "Key" => "Project"), Dict{String, Any}("Value" => "Engineering", "Key" => "Department")])
end

register_instance_1_patch = @patch function register_instance(a...; b...)
    return Dict{String, Any}("OperationId" => "4yejorelbukcjzpnr6tlmrghsjwpngf4-k95yg2u7")
end

tag_resource_1_patch = @patch function tag_resource(a...; b...)
    return Dict{String, Any}()
end

untag_resource_1_patch = @patch function untag_resource(a...; b...)
    return Dict{String, Any}()
end

update_service_1_patch = @patch function update_service(a...; b...)
    return Dict{String, Any}("OperationId" => "m35hsdrkxwjffm3xef4bxyy6vc3ewakx-jdn3y5g5")
end


text:
using Test
using Random
using ProximalOperators

Random.seed!(0)

@testset "CubeNormL2" begin

for R in [Float16, Float32, Float64]
    for T in [R, Complex{R}]
        for shape in [(5,), (3, 5), (3, 4, 5)]
            lambda = R(0.1) + 5*rand(R)
            f = CubeNormL2(lambda)
            predicates_test(f)
            x = randn(T, shape)
            call_test(f, x)
            gamma = R(0.5)+rand(R)
            y, f_y = prox_test(f, x, gamma)
            grad_f_y, f_y = gradient(f, y)
            @test grad_f_y ≈ (x - y)/gamma
        end
    end
end

end


text:
@testset "Lists" begin

@test filter(iseven, List(1, 2, 3)) == List(2)
@test filter(!iseven, List(1, 2, 3)) == List(1, 3)
@test List(1, 2, 3) ++ List(4, 5, 6) == List(1, 2, 3, 4, 5, 6)
@test append(List(1), List(2), List(3)) == List(1, 2, 3)

@testset "repr" begin
    @test repr(list(1, 2, 3)) == "(1 2 3)"
    @test repr(list(1, true, nothing)) == "(1 #t #<void>)"
end

end


text:
mutable struct Texture
    _image::Any
end

image(texture::Texture) = texture._image

function Texture(path::String)
    data = open(read, path)
    return Texture(data)
end
function TextureData(data::Vector{UInt8})
    return Texture(data)
end

#-----------------------------------------------------------------
mutable struct Material <: AbstractMaterial
    _color::RGBA{Float64}
    _line_width::Float64
    _wireframe::Bool
    _wireframe_line_width::Float64
    _size::Float64
    _texture::Union{Nothing, Texture}
end

function Material(;
    color = Colors.color_names["magenta"],
    line_width::Float64  = 1.0,
    wireframe::Bool = false,
    wireframe_line_width::Float64  = 1.0,
    size::Float64 = 0.004,
    texture::Union{Nothing, Texture} = nothing,
    kwargs...)
    return Material(to_color(color), line_width, wireframe, wireframe_line_width, size, texture)
end

color(m::Material) = m._color
function color!(m::Material, color::RGBA{Float64})
    m._color = color
end

texture(m::Material) = m._texture
function texture!(m::Material, texture::Union{Nothing, Texture})
    m._texture = texture
end

lineWidth(m::Material) = m._line_width
wireframe(m::Material) = m._wireframe
wireframeLineWidth(m::Material) = m._wireframe_line_width
Base.size(m::Material) = m._size

function material(m::Material)
    tex = (texture(m) !== nothing) ? MeshCat.Texture(wrap= (300, 300), repeat= (1, 1), image=MeshCat.PngImage(Glimmer.image(Glimmer.texture(m)))) : nothing

    mat = MeshCat.MeshPhongMaterial(
        color=color(m),
        linewidth = lineWidth(m),
        wireframe = wireframe(m),
        wireframeLinewidth = wireframeLineWidth(m),
        map=tex
    )

    return mat
end

function lines_material(m::Material)
    mat = MeshCat.LineBasicMaterial(
        color=color(m),
        linewidth = lineWidth(m),
        wireframe = wireframe(m),
        wireframeLinewidth = wireframeLineWidth(m),
    )
    return mat
end

function points_material(m::Material)
    mat = MeshCat.PointsMaterial(
        color=color(m),
        size = size(m),
    )
    return mat
end

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


text:
module HouseHolders

  function QRDecompose(R)
    println("computing I")
    I=eye(size(R)[1])
    println("phew")
    Q,R= decompose(R)
    return Q',R
  end
  function decompose(R)
    println(size(R))
    if 0 in size(R)

      return eye(size(R)[1]),1
    end
    a1 = R[:,1]
    lena1 = norm(a1)
    d11 = a1[1]>0?-1*lena1:lena1
    w11 = a1[1]-d11
    f1 = sqrt(-2*w11*d11)
    v = [w11; a1[2:end]]/f1
    #println(a1,lena1,d11,w11,f1,v)
    #A = 2*v*v'
    println("computing H...")
    H=eye(size(R)[1])
    BLAS.gemm!('N','T',-2.,v,v,1.0,H)
  #  println(H)

    println("did i compute H?")
    for i=1:size(R)[1]
      R[i,1]=0.
    end

    println("did i clear column1?")
    for i=2:size(R)[2]
      ai =R[:,i]
      fi = (2*v'*ai)[1]
      R[:,i] = ai-fi*v
    end
    # R=H*R
    R[1,1]=d11
  #  println(R)
  println("is it the matrix multiply??")

    newH,newR = decompose(R[2:end,2:end])
    newHH = eye(size(H)[1])
    # println(newHH)
    # println(newH)
    newHH[2:end,2:end]=newH
    println(newHH)

    R[2:end,2:end]=newR

    println("heh?")
    return BLAS.gemm('N','N',newHH,H),R
  end

  function householdDecompose(R,y)
    return decompose(R,y,1)
  end

  function decompose(R,y,n)
    if n-1 in size(R)
      return
    end
    #computes our needed values
    a1 = R[n:end,n]
    lena1 = norm(a1)
    d11 = a1[1]>0?-1*lena1:lena1
    w11 = a1[1]-d11
    f1 = sqrt(-2*w11*d11)
    v = [w11; a1[2:end]]/f1
    #clear out that first column
    for i=n:size(R)[1]
      R[i,n]=0.
    end
    #iterative approach to do Hai.
    for i=n+1:size(R)[2]
      ai =R[n:end,i]
      fi = (2*v'*ai)[1]
      R[n:end,i] = ai-fi*v
    end
    R[n,n]=d11

    #Optimized way to update y without matrices!
    y[n:end]-=2*(y[n:end]'*v)[1]*v

    return decompose(R,y,n+1)
  end
end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
abstract type Prior end

# make priors iterable
Base.iterate(design::Prior, state = 0) = state > 0 ? nothing : (design, state + 1)
Base.length(design::Prior) = 1

Base.show(io::IO, prior::Prior) = print(io, string(prior))
Base.show(io::IO, ::MIME"application/prs.juno.inline", prior::Prior) = print(io, string(prior))

bounds(prior::Prior) = [prior.low, prior.high]

pdf(p::Real, prior::Prior) = error("not implemented")

cdf(p::Real, prior::Prior) = error("not implemented")

expectation(f::Function, prior::Prior) = error("not implemented")

mean(prior::Prior) = error("not implemented")

condition(prior::Prior; low::Real = prior.low, high::Real = prior.high) = error("not implemented")
function |(prior::Prior, interval::Vector{T}) where {T<:Real}

    @assert (length(interval) == 2) @sprintf "must condition on interval (vector of length 2, is %i)" length(interval)
    @assert (0 <= interval[1] < interval[2] <= 1)  @sprintf "invalid interval specification must be 0<=a<b<=1, is [%.2f, %.2f]" interval[1] interval[2]
    condition(prior; low = interval[1], high = interval[2])
end
condition(prior::Prior, point::Real) = (pdf(point, prior) > 0) ? PointMass(point) : error("conditioning on point with pdf = 0 is not well-defined")
|(prior::Prior, point::Real) = condition(prior, point)

<=(low::T, prior::Prior) where {T<:Real} = condition(prior; low = low)
<=(prior::Prior, high::T) where {T<:Real} = condition(prior; high = high)
>=(prior::Prior, low::T) where {T<:Real} = condition(prior; low = low)
>=(high::T, prior::Prior) where {T<:Real} = condition(prior; high = high)

update(prior::Prior, x::Integer, n::Integer) = error("not implemented")

string(prior::Prior) = error("not implemented")


text:
"""
    graph_preferential(n::Integer, m::Integer)

Scale-free network model of Barabási and Albert (1999).

# Arguments
* `n`: number of nodes,
* `m`: number of edges added with each subsequent node.

# References
Barabási, A.-L., & Albert, R. (1999).
Emergence of Scaling in Random Networks. Science, 286(5439), 509.
"""
function graph_preferential(n::Integer, m::Integer; kwargs...)
    g = Graph(n; kwargs...)
    p = DiscreteSampler(Int[])
    for i = 1:m         # initial complete subgraph on `m` nodes
        for j = 1:i - 1
            addedge!(g, i, j)
        end
        push!(p, max(1, m - 1))
    end
    for i = m + 1:n     # preferential attachment for remaining nodes
        while outdegree(g, i) < m
            j = randd(p)
            hasedge(g, i, j) && continue
            addedge!(g, i, j)
            inc_index!(p, j)
        end
        push!(p, m)
    end
    return g
end


text:
import Base: precision

"""
    precision(a::Padic)

Return the number of digits of precision in `a`
"""
precision(a::Padic) = a.x.parent.prec_max

import Nemo: valuation, base

export base
"""
    base(a::Padic)

Return the prime number `p` for this `p`-adic number.
"""
base(a::Padic{P}) where {P} = P

export valuation
valuation(a::Padic) = valuation(a.x)

import Base: sqrt, exp, log, abs

sqrt(a::Padic) = Padic(sqrt(a.x))
exp(a::Padic) = Padic(exp(a.x))
log(a::Padic) = Padic(log(a.x))

function abs(a::Padic{P}) where {P}
    if a == 0
        return 0.0
    end
    return Float64(P)^(-valuation(a))
end


text:
# This file is a part of TextClassification.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

using IterTools
using Random

function filtered_power_set(set, lowersize=0, uppersize=5)
    lst = collect(subsets(set))
    filter(x -> lowersize <= length(x) <= uppersize, lst)
end

@with_kw struct TextConfigSpace <: AbstractSolutionSpace
    del_diac = [true]
    del_dup = [false]
    del_punc = [false]

    group_num = [true]
    group_url = [true]
    group_usr = [true]
    group_emo = [false]

    lc = [true]
    qlist = filtered_power_set([2, 3, 4, 5], 1, 3)
    nlist = filtered_power_set([1, 2, 3], 0, 2)
    slist = [[Skipgram(2, 1), Skipgram(2, 2)], [Skipgram(2, 1)], [Skipgram(2, 2)], Skipgram[]]

    qlist_space = 1:9
    nlist_space = 1:5
    slist_space = [Skipgram(2, 1), Skipgram(2, 2), Skipgram(3, 1), Skipgram(3, 2)]
end

Base.eltype(::TextConfigSpace) = TextConfig

function Base.rand(space::TextConfigSpace)
    qlist = isempty(space.qlist) ? Int[] : rand(space.qlist)
    nlist = isempty(space.nlist) ? Int[] : rand(space.nlist)
    slist = isempty(space.slist) ? Skipgram[] : rand(space.slist)

    if length(qlist) + length(nlist) + length(slist) == 0
        nlist = [1]
    end

    TextConfig(
        rand(space.del_diac),
        rand(space.del_dup),
        rand(space.del_punc),

        rand(space.group_num),
        rand(space.group_url),
        rand(space.group_usr),
        rand(space.group_emo),

        rand(space.lc),
        qlist,
        nlist,
        slist
    )
end

function combine(a::TextConfig, b::TextConfig)
    L = [a, b]

    qlist = rand() < 0.5 ? union(a.qlist, b.qlist) : intersect(a.qlist, b.qlist)
    nlist = rand() < 0.5 ? union(a.nlist, b.nlist) : intersect(a.nlist, b.nlist)
    slist = rand() < 0.5 ? union(a.slist, b.slist) : intersect(a.slist, b.slist)

    if length(qlist) + length(nlist) + length(slist) == 0
        nlist = [1]
    end

    TextConfig(
        rand(L).del_diac,
        rand(L).del_dup,
        rand(L).del_punc,
        rand(L).group_num,
        rand(L).group_url,
        rand(L).group_usr,
        rand(L).group_emo,
        rand(L).lc,
        sort!(qlist),
        sort!(nlist),
        sort!(slist)
    )
end

function mutate_token_list(lst, L; p1=0.5, p2=0.5)
    if length(L) == 0 || rand() < p1
        return lst
    end

    if rand() < p2
        if length(lst) > 0
            lst = copy(lst)
            shuffle!(lst)
            pop!(lst)
            sort!(lst)
        end
    else
        lst = union(lst, [rand(L)])
    end

    lst
end

function mutate(space::TextConfigSpace, c::TextConfig, iter)
    qlist = mutate_token_list(c.qlist, space.qlist_space)
    nlist = mutate_token_list(c.nlist, space.nlist_space)
    slist = mutate_token_list(c.slist, space.slist_space)
    if length(qlist) + length(nlist) + length(slist) == 0
        nlist = [1]
    end

    TextConfig(
        SearchModels.change(c.del_diac, space.del_diac),
        SearchModels.change(c.del_dup, space.del_dup),
        SearchModels.change(c.del_punc, space.del_punc),
        SearchModels.change(c.group_num, space.group_num),
        SearchModels.change(c.group_url, space.group_url),
        SearchModels.change(c.group_usr, space.group_usr),
        SearchModels.change(c.group_emo, space.group_emo),
        SearchModels.change(c.lc, space.lc),
        qlist,
        nlist,
        slist
    )
end


text:
# TODO
# - Document extent option.

# extent should be a vector of integers of length 6.
function vtk_grid{T<:AbstractFloat}(
        filename_noext::AbstractString,
        x::AbstractVector{T}, y::AbstractVector{T}, z::AbstractVector{T};
        compress::Bool=true, append::Bool=true, extent=nothing)
    xvtk = XMLDocument()

    Ni, Nj, Nk = length(x), length(y), length(z)
    Npts = Ni*Nj*Nk
    Ncls = num_cells_structured(Ni, Nj, Nk)

    ext = extent_attribute(Ni, Nj, Nk, extent)

    vtk = DatasetFile(xvtk, filename_noext*".vtr", "RectilinearGrid",
                      Npts, Ncls, compress, append)

    # VTKFile node
    xroot = vtk_xml_write_header(vtk)

    # RectilinearGrid node
    xGrid = new_child(xroot, vtk.gridType_str)
    set_attribute(xGrid, "WholeExtent", ext)

    # Piece node
    xPiece = new_child(xGrid, "Piece")
    set_attribute(xPiece, "Extent", ext)

    # Coordinates node
    xPoints = new_child(xPiece, "Coordinates")

    # DataArray node
    data_to_xml(vtk, xPoints, x, "x")
    data_to_xml(vtk, xPoints, y, "y")
    data_to_xml(vtk, xPoints, z, "z")

    return vtk::DatasetFile
end

# 2D version
vtk_grid{T<:AbstractFloat}(
        filename_noext::AbstractString,
        x::AbstractVector{T}, y::AbstractVector{T};
        compress::Bool=true, append::Bool=true, extent=nothing) =
    vtk_grid(filename_noext, x, y, zeros(T,1), compress=compress, append=append, extent=extent)


text:
# # Saved Experiments
#
#md # [`notebook`](@__NBVIEWER_ROOT_URL__notebooks/two_problems.ipynb)
#
using WaterWaves1D,Plots,FFTW,Statistics;gr()
#include("../src/dependencies.jl")

#---- Cold-Hot states transition for Serre-Green-Naghdi, discussed in the paper of Sergey Gavrilyuk, Boniface Nkonga, Keh-Ming Shyue and Lev Truskinovsky
function ColdHot(;scenario)
	if scenario !=1 && scenario !=2
		error("the scenario must be 1 or 2")
	end
	param = ( μ  = 1, ϵ  = 1, N  = 2^13,
			h₀=1.0962,h₁=1.1,h₂=1.2)
	@time (η,u,v,mesh,par)=CnoidalWaveWhithamGreenNaghdi(param;P=500,SGN=true)
	x=mesh.x
	h=1 .+ param.ϵ*η
	m=-sqrt(param.h₀*param.h₁*param.h₂)
	if scenario == 1
		meanη=mean(h)-1
		meanu=m*mean(1 ./h)
	else
		meanη=1.09808-1
		meanu=m/1.09808
	end

	M=-m*mean(1 ./ h)/sqrt(mean(h)) #Froude number, close to 1
	L=250*par.λ

	u.-=mean(u)
	u.+=m*mean(1 ./h)

	interp(x,δ)=sech.(x/δ).^2
	η₀(x,δ)=(meanη.+(0.2-meanη)*interp(x,δ))
	u₀(x,δ)=(meanu .+(u[Int(end/2)]-meanu)*interp(x,δ))

	Dx(v)=real.(ifft(1im*mesh.k.*fft(v)))

	η1=η.*(x.>=-L).*(x.<=L)+η₀(x.+L,0.1).*(x.<-L)+η₀(x.-L,0.1).*(x.>L)
	u1=u.*(x.>=-L).*(x.<=L)+u₀(x.+L,0.1).*(x.<-L)+u₀(x.-L,0.1).*(x.>L)
	v1=u1-param.μ/3*Dx(h.^3 .*Dx(u1))./h

	η2=η.*(x.>=-L).*(x.<=L)+η₀(x.+L,2).*(x.<-L)+η₀(x.-L,2).*(x.>L)
	u2=u.*(x.>=-L).*(x.<=L)+u₀(x.+L,2).*(x.<-L)+u₀(x.-L,2).*(x.>L)
	v2=u2-param.μ/3*Dx(h.^3 .*Dx(u2))./h

	η3=η.*(x.>=-L).*(x.<=L)+η₀(x.+L,3).*(x.<-L)+η₀(x.-L,3).*(x.>L)
	u3=u.*(x.>=-L).*(x.<=L)+u₀(x.+L,3).*(x.<-L)+u₀(x.-L,3).*(x.>L)
	v3=u3-param.μ/3*Dx(h.^3 .*Dx(u3))./h

	init1=Init(mesh,η1,v1)
	init2=Init(mesh,η2,v2)
	init3=Init(mesh,η3,v3)

	param2=merge(param,(L=-mesh.xmin,T = 8000, dt = 0.05, ns=1000))
	@time model=WhithamGreenNaghdi(param2;SGN=true,dealias=0,precond=false, gtol=1e-12)

	problem1 = Problem(model, init1, param2)
	problem2 = Problem(model, init2, param2)
	problem3 = Problem(model, init3, param2)
	solve!( problem1 )
	solve!( problem2 )
	solve!( problem3 )

	p = plot(layout=(2,1))
	plot_solution!(p,problem1,t=2000)
	plot_solution!(p,problem2,t=2000)
	plot_solution!(p,problem3,t=2000)
	savefig("ColdHot2000.pdf")
	xlims!(p[1,1],-L-200,-L+50)
	savefig("ColdHot2000Xzoom.pdf")
	ylims!(p[1,1],0.095,0.1)
	savefig("ColdHot2000Yzoom.pdf")

	p = plot(layout=(2,1))
	plot_solution!(p,problem1,t=4000)
	plot_solution!(p,problem2,t=4000)
	plot_solution!(p,problem3,t=4000)
	savefig("ColdHot4000.pdf")
	xlims!(p[1,1],-L-200,-L+50)
	savefig("ColdHot4000Xzoom.pdf")
	ylims!(p[1,1],0.095,0.1)
	savefig("ColdHot4000Yzoom.pdf")

	p = plot(layout=(2,1))
	plot_solution!(p,problem1,t=8000)
	plot_solution!(p,problem2,t=8000)
	plot_solution!(p,problem3,t=8000)
	savefig("ColdHot8000.pdf")
	xlims!(p[1,1],-L-200,-L+50)
	savefig("ColdHot8000Xzoom.pdf")
	ylims!(p[1,1],0.095,0.1)
	savefig("ColdHot8000Yzoom.pdf")

	display(p)

	save(problem1,"ColdHot1")
	save(problem2,"ColdHot2")
	save(problem3,"ColdHot3")
end

#---- Wavebreaking for the quasilinear Boussinesq-Whitham equation
function Wavebreaking1()
	param = ( μ  = 1/2,
			ϵ  = 1,
        	N  = 2^18,
            L  = 3,
            T  = .04,
            dt = 0.00001)
	η = x -> 2 .^(-abs.((x).^4))
	v = x -> -100*x.*η(x)
	init = Init(η,v)
	model = WhithamBoussinesq(merge(param,(α=1/2,)))
	problem = Problem(model, init, param)
	@time solve!( problem )
	p = plot(layout=(2,1))
	plot_solution!( p, problem, .03898)
	plot!(xlims=[(0.0018,0.0025) (-100000,100000 )])
	display(p)
end

#---- Wavebreaking for the quasilinear Boussinesq-Whitham equation ?
function Wavebreaking2()
	param = ( μ  = 1/2,
				ϵ  = 1,
	        	N  = 2^17, # you need to add more modes
	            L  = 3,
	            T  = 1,
	            dt = 0.0001) # you need to lower this
	η = x -> 2 .^(-abs.((x).^4))
	v = x -> -10*x.*η(x)
	init = Init(η,v)
	model = WhithamBoussinesq(merge(param,(α=1/2,)))
	problem = Problem(model, init, param)
	@time solve!( problem )
	p = plot(layout=(2,1))
	plot_solution!( p, problem,t= .92)
	plot!(xlims=[(0.275,0.285) (-100000,100000 )])
	display(p)
end


text:
include("utils.jl")

build_so("006")

immutable Foo
   i::Cint
end

x = Foo(13)

y = ccall(
    (:bar, "binary/lib006.so"),
    Cint,
    (Foo, ),
    x
)

@assert y == x.i


text:
using Test
@time @testset "Acoustics" begin include("test_acoustics.jl") end
true


text:
using Distributions
using RmathDist

distlist = [
            Geometric(0.5),
            Geometric(0.999999999),
            Geometric(1e-10)
            ]

for d in distlist
    for x in [0.0:20.0,linspace(21.0,maxintfloat()/4,1_000_000)]
        x = iround(x)

        test_rmath(pdf,d,x)
        test_rmath(logpdf,d,x)

        test_rmath(cdf,d,x)
        test_rmath(ccdf,d,x)
        test_rmath(logcdf,d,x)
        test_rmath(logccdf,d,x)
    end

    for p in logspace(-300,0,1_000_000)
        test_rmath(quantile,d,p)
        test_rmath(cquantile,d,p)
    end

    for lp in -logspace(-300,300,1_000_000)
        test_rmath(invlogcdf,d,lp)
        test_rmath(invlogccdf,d,lp)
    end
end


text:
@defcomp PermafrostTotal begin
    perm_sib_ce_co2 = Parameter(index=[time], unit="Mtonne CO2")
    perm_jul_ce_co2 = Parameter(index=[time], unit="Mtonne CO2")
    perm_sib_e_co2 = Parameter(index=[time], unit="MtonCO2/yr")
    perm_jul_e_co2 = Parameter(index=[time], unit="MtonCO2/yr")
    perm_sib_ce_ch4 = Parameter(index=[time], unit="Mtonne CH4")
    perm_jul_ce_ch4 = Parameter(index=[time], unit="Mtonne CH4")

    perm_initial_c_stock_uncertainty = Parameter(default=0.0)

    perm_tot_ce_co2 = Variable(index=[time], unit="Mtonne")
    perm_tot_e_co2 = Variable(index=[time], unit="Mtonne")
    perm_tot_ce_ch4 = Variable(index=[time], unit="Mtonne")

    function run_timestep(p, v, d, tt)
        # Cumulative emissions CO2 perm
        v.perm_tot_ce_co2[tt] = 0.5 *(p.perm_sib_ce_co2[tt] + p.perm_jul_ce_co2[tt]) * (1 + p.perm_initial_c_stock_uncertainty/100)
        # Annual emissions CO2 perm
        v.perm_tot_e_co2[tt] = 0.5 *(p.perm_sib_e_co2[tt] + p.perm_jul_e_co2[tt]) * (1 + p.perm_initial_c_stock_uncertainty/100)
        # Cumulative emissions CH4 perm
        v.perm_tot_ce_ch4[tt] = 0.5 *(p.perm_sib_ce_ch4[tt] + p.perm_jul_ce_ch4[tt]) * (1 + p.perm_initial_c_stock_uncertainty/100)
    end
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
# using Dists: ValueSupport, VariateForm

# struct MixedSupport <: ValueSupport end
# struct MixedVariate <: VariateForm end

"""
    AbstractModel{A,B}

Gives an abstract type for all Tilde models

Type variables ending in T are type-level representations used to reconstruct

N gives the Names of arguments (each a Symbol)
B gives the Body, as an Expr
M gives the Module where the model is defined
"""
abstract type AbstractModel{A,B,M} <: AbstractKleisli end

abstract type AbstractConditionalModel{M, Args, Obs} <: AbstractMeasure end

argstype(::AbstractModel{A,B,M}) where {A,B,M} = A

bodytype(::AbstractModel{A,B,M}) where {A,B,M} = B

getmodule(::Type{AMF}) where {A,B,M, AMF<:AbstractModel{A,B,M}} = from_type(M)
getmodule(::AbstractModel{A,B,M}) where {A,B,M} = from_type(M)

# getmoduletypencoding(::Type{AbstractModel{A,B}}) where  {M,A,O,AM<:AbstractModel{A,B}} = M
# getmoduletypencoding(::AbstractModel{A,B}) where  {M,A,O,AM<:AbstractModel{A,B}} = M

argvalstype(::AbstractModel{A}) where {A} = A
argvalstype(::Type{AM}) where {A,AM<:AbstractModel{A}} = A

obstype(::AbstractModel) = NamedTuple{(), Tuple{}}
obstype(::Type{<:AbstractModel}) = NamedTuple{(), Tuple{}}

(m::AbstractModel)(;argvals...)= m((;argvals...))

(m::AbstractModel{A})(args...) where {A} = m(A(args))

body(::AbstractModel{A,B}) where {A,B} = from_type(B)
body(::Type{AM}) where {A,B,AM<:AbstractModel{A,B}} = from_type(B)


text:
function plot_resource_utilization(fig, gp, ax, x, y, smoothing,ticks)
    xlow = fill(0, length(unique(x)))
    ylow = zeros(length(xlow))
    xhigh = zero(xlow)
    yhigh = zero(ylow)
    yall = zero(ylow)
    xall = zero(xlow)
    last = first(x)
    idx = 1
    for (new, ut) in zip(x, y)
        idx += new !== last
        last = new
        if ut > 1
            yhigh[idx] += ut
            xhigh[idx] += 1
        else
            ylow[idx] += ut
            xlow[idx] += 1
        end
        yall[idx] += ut
        xall[idx] += 1
    end
    ylow ./=  [n == 0 ? 1 : n for n in xlow]
    yhigh ./=  [n == 0 ? 1 : n for n in xhigh]
    yall ./= xall
    perhigh = xhigh ./(xlow .+ xhigh)
    if smoothing
        ylow = accumulate!((o,n) -> o*0.8+n*0.2, ylow,ylow)
        yhigh = accumulate!((o,n) -> o*0.8+n*0.2, yhigh,yhigh)
        yall = accumulate!((o,n) -> o*0.8+n*0.2, yall,yall)
    end
    b = barplot!(ax,ticks, perhigh, color = :gray80)
    l1 = lines!(ax,ticks, yhigh, color = :red)
    l2 = lines!(ax,ticks, ylow, color = :green)
    l3 = lines!(ax,ticks, yall, color = :blue)
    gp[2,1] = Legend(fig, [b,l1,l2,l3], ["% over capacity", "mean over capacity","mean under capacity", "mean utilization"], tellwidth = false, tellhight = false, valign = :top , halign = :right)
    nothing
end
function plot_ressource_setups(fig,gp,ax,x,y, smoothing,ticks)
    yint = Int.(y)
    ys = sort(unique(yint)) #grps
    counts = zeros(length(ys), maximum(x))
    xcounts = fill(0, 1, maximum(x))
    last = first(x)
    idx = 1
    for (new, b) in zip(x, yint)
        idx += new !== last
        last = new
        xcounts[idx] += 1
        counts[setups, idx] += 1
    end
    xs = repeat(unique(x), inner = length(ys))
    counts ./= xcounts
    grp = repeat(ys, outer = maximum(x))
    pal = ax.attributes[:palette][:color][]
    barplot!(ax,ticks,xs, vec(counts), stack=grp, color = [pal[g] for g in grp])
    labels = string.(ys)
    gp[2,1] = Legend(fig, [PolyElement(color = i) for i in [pal[j] for j in ys]], labels, "# setups", tellwidth = false, tellhight = false, valign = :top , halign = :right)
    nothing
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: jl:percent
#     text_representation:
#       extension: .jl
#       format_name: percent
#       format_version: '1.3'
#       jupytext_version: 1.4.2
#   kernelspec:
#     display_name: Julia 1.4.1
#     language: julia
#     name: julia-1.4
# ---

# %%
# Ensure your working directory is: https://github.com/henriquebecker91/phd/tree/master/latex/revised_PPG2KP
include("notebook_setup.jl")

# %%
# Read the data, show nothing.
raw_vudf_path = "./data/vu_experiment.csv"
raw_vudf = DataFrame(CSV.File(raw_vudf_path))
@show collect(zip(names(raw_vudf), typeof.(eachcol(raw_vudf))))
#showtable(raw_vudf)
nothing

# %%
# Just assert some basic things, and show the solved runs.
let raw_vudf = deepcopy(raw_vudf)
    @show size(raw_vudf)
    #describe(raw_vudf)
    @assert @with(raw_vudf, all(:had_timeout .== :run_ended_by_exception))
    #@assert @with(raw_vudf, .!((isnan.(:solution_profit) .& isnan.(:obj_value)) .|
    #    (:solution_profit .≈ :obj_value)))
    @show 160 - sum(raw_vudf[!, :had_timeout])
    # Should be an empty table. If it is not, then the solution assembled from model variables has
    # a different profit value from the objective function.
    @assert isempty(@where(raw_vudf, .!((isnan.(:solution_profit) .& isnan.(:obj_value)) .|
        (:solution_profit .≈ :obj_value))))
    solved = @linq raw_vudf |> select(:instance_name, :pricing_method, :run_total_time, :solution_profit, :had_timeout) |>
        where(.!:had_timeout)
end

# %%
# Clean the data: replace sentinel values by missing; replace full paths by only the filename.
# No need to create config id (because there is only one parameter that changes between them);
# no need to find the right number of vars/plates (because for both is just the after_purge).
vudf = let vudf = deepcopy(raw_vudf)
    # Keep only the instance name (not the path).
    @with(vudf, :instance_name .= basename.(:instance_name))
    @with(vudf, :datafile .= basename.(:datafile))
    colnames = names(vudf)
    for name in colnames
        column = vudf[!, name]
        if eltype(column) <: AbstractFloat
            vudf[!, name] = ifelse.(isnan.(column), missing, column)
        elseif eltype(column) <: Integer
            vudf[!, name] = ifelse.(column .< 0, missing, column)
        end
    end

    vudf
end
showtable(vudf)
#nothing

# %%
@show names(vudf)
@show unique(vudf[!, :build_stop_reason])
@show @where(vudf, :had_timeout .!= :run_ended_by_exception)

# %%
vu2019_path = "./data/vu_2019.csv"
vu2019 = DataFrame(CSV.File(vu2019_path))
#nothing
showtable(vu2019)

# %%
# Get some statistics regarding the heuristic.
import Statistics: mean
stats_heuristic = let vudf = deepcopy(vudf)
    # Normalizes how some information appears in the restricted and unrestricted steps.
    vudf[!, :restricted_had_timeout] = occursin.("TIME_LIMIT", vudf[!, :restricted_stop_reason]) .|
        occursin.("TIME_LIMIT", vudf[!, :restricted_lp_stop_reason])
    vudf[!, :built_restricted_model] = .!occursin.("TIME_LIMIT", vudf[!, :restricted_lp_stop_reason])
    vudf[!, :built_model] = vudf[!, :build_stop_reason] .== "BUILT_MODEL"
    # Sets of columns used extensively below.
    shared_columns = [:instance_name, :pricing_method]
    restricted_columns = [
        :heuristic_lb, :restricted_obj_value, :restricted_had_timeout, :built_restricted_model
    ]
    unrestricted_columns = [
        :heuristic_lb, :obj_value, :had_timeout, :built_model
    ]
    # Clean the dataframe.
    used_columns = unique(vcat(shared_columns, restricted_columns, unrestricted_columns))
    select!(vudf, used_columns) # Get only the used_columns
    # Create a restricted dataframe, with the same column names as the unrestricted counterpart.
    restricted = select(vudf, vcat(shared_columns, restricted_columns))
    filter!("pricing_method" => isequal("furini"), restricted)
    select!(restricted, vcat(shared_columns, restricted_columns .=> unrestricted_columns)...)
    restricted[!, :pricing_method] .= "restricted"
    # Create a single dataframe with restricted data as extra rows, not extra columns in
    # priced rows.
    select!(vudf, vcat(shared_columns, unrestricted_columns))
    vudf = vcat(vudf, restricted)
    filter!(:built_model => identity, vudf)
    #@show vudf
    select!(vudf, Not(:pricing_method), :pricing_method => :variant)
    vudf[!, :heuristic_gap] = @with vudf (100.0 .- (:heuristic_lb ./ :obj_value) .* 100.0)
    @show(@where(vudf, iszero.(:heuristic_gap)))
    #filter!([:instance_name, :variant] => !((i, v) -> i == "P3_250_250_25_3" && v == "none"), vudf)
    @show maximum(vudf[!, :heuristic_gap])
    @show minimum(vudf[!, :heuristic_gap])
    @show mean(vudf[!, :heuristic_gap])
    @linq select(vudf, :heuristic_gap, :variant) |> groupby(:variant) |>
        based_on(; mean = mean(:heuristic_gap), max = maximum(:heuristic_gap), min = minimum(:heuristic_gap))
    #sort!(vudf, [:instance_name, :variant])
end
#nothing
showtable(stats_heuristic)

# %%
# Summary dataframe. Dataframe with summarized time, #variables, etc...
function bitmask_of_rows_with_missing_values(df)
    return reduce(broadwrap(|), (map(broadwrap(ismissing), eachcol(df))))
end
sdf = let vudf = deepcopy(vudf)
    # Normalizes how some information appears in the restricted and unrestricted steps.
    vudf[!, :restricted_had_timeout] = occursin.("TIME_LIMIT", vudf[!, :restricted_stop_reason]) .|
        occursin.("TIME_LIMIT", vudf[!, :restricted_lp_stop_reason])
    vudf[!, :restricted_early_opt] = .!vudf[!, :solved_priced_restricted_model] .&
        .!occursin.("TIME_LIMIT", vudf[!, :restricted_stop_reason]) .&
        .!occursin.("TIME_LIMIT", vudf[!, :restricted_lp_stop_reason])
    vudf[!, :early_opt] = vudf[!, :build_stop_reason] .== "FOUND_OPTIMUM"

    # Sets of columns used extensively below.
    shared_columns = [:instance_name, :pricing_method]
    restricted_columns = [
        :qt_pevars_priced_restricted,
        :qt_cmvars_priced_restricted,
        :qt_plates_priced_restricted,
        :restricted_pricing_time,
        :restricted_had_timeout,
        :restricted_early_opt,
    ]
    unrestricted_columns = [
        :qt_pevars_after_purge,
        :qt_cmvars_after_purge,
        :qt_plates_after_purge,
        :build_and_solve_time,
        :had_timeout,
        :early_opt,
    ]
    # Clean the dataframe.
    used_columns = vcat(shared_columns, restricted_columns, unrestricted_columns)
    select!(vudf, used_columns) # Get only the used_columns
    # Create a restricted dataframe, with the same column names as the unrestricted counterpart.
    restricted = select(vudf, vcat(shared_columns, restricted_columns))
    filter!("pricing_method" => isequal("furini"), restricted)
    select!(restricted, vcat(shared_columns, restricted_columns .=> unrestricted_columns)...)
    restricted[!, :pricing_method] .= "restricted"
    # Create a single dataframe with restricted data as extra rows, not extra columns in
    # priced rows.
    select!(vudf, vcat(shared_columns, unrestricted_columns))
    vudf = vcat(vudf, restricted)
    # Create the instance groups and remove/rename some columns.
    vudf[!, :group] = (m -> parse(Int, m.captures[1])).(
        match.(r"^P([1-4]).*", vudf[!, :instance_name])
    )
    select!(vudf, Not(:instance_name))
    select!(vudf, Not(:pricing_method), :pricing_method => :variant)
    # If no instance makes the pricing return early, then we can exclude this from the comparison.
    @assert iszero(sum(vudf[!, :early_opt]))
    select!(vudf, Not(:early_opt))
    @assert all(ifelse.(
        bitmask_of_rows_with_missing_values(vudf), # If a row has some missing value
        vudf[!, :had_timeout], # then the same row has an indication of timeout
        true # if there is not missing value, then the run can have hit the time limit or not
    ))
    #@show maximum(skipmissing(vudf[!, :build_and_solve_time]))
    #showtable(@where(vudf, (:group .== 4) .& (:variant .== "furini")))
    vudf = @linq vudf |> groupby([:group, :variant]) |>
        based_on(;
            solved = sum(.!:had_timeout),
            variables = sum(skipmissing(:qt_pevars_after_purge)) + sum(skipmissing(:qt_cmvars_after_purge)),
            plates = sum(skipmissing(:qt_plates_after_purge)),
            built = sum(.!ismissing.(:qt_cmvars_after_purge)),
            total_time = sum(replace(
                :build_and_solve_time, missing => (3 * 60 * 60)
            ))
        ) |> sort([:group, :variant])
    vudf
end
#nothing
#showtable(@where(sdf, .!reduce(broadwrap(|), (map(broadwrap(ismissing), eachcol(sdf))))))

# %%
st = let sdf = deepcopy(sdf)
    variant_order = Dict{String, Int}(
        "none"       => 1,
        "restricted" => 2,
        "furini"     => 3,
    )
    sdf[!, :variant_order] = (e -> variant_order[e]).(sdf[!, :variant])
    sort!(sdf, [:group, :variant_order])
    select!(sdf, Not(:variant_order))
    pretty_names = Dict{String, String}(
        "none"       => "Not Priced",
        "restricted" => "Restricted Priced",
        "furini"     => "Priced",
    )
    replace!((e -> pretty_names[e]), sdf[!, :variant])
    # Transform some values in averages
    sdf[!, :variables] ./= sdf[!, :built]
    sdf[!, :plates] ./= sdf[!, :built]
    sdf[!, :solved_time] = sdf[!, :total_time] .- ((20 .- sdf[!, :solved]) .* (3 * 60 * 60))
    sdf[!, :solved_time] ./= sdf[!, :solved]
    sdf[!, :total_time] = trunc.(Int, sdf[!, :total_time])
    # define header names and column order
    select!(sdf,
        :group => "G.",
        :variant => "Variant",
        :built => "#m",
        :variables => #=(c -> c ./ sdf[!, :built]) =>=# "Avg. #v",
        :plates => #=(c -> c ./ sdf[!, :built]) =>=# "Avg. #p",
        :total_time => "T. T.",
        :solved => "#s",
        :solved_time => #=(c -> c ./ sdf[!, :solved]) =>=# "Avg. S. T.",
    )
    select!(sdf, names(sdf) .=> esc_latex.(names(sdf)))
    for column in names(sdf) # latexify all numeric columns
        if nonmissingtype(eltype(sdf[!, column])) <: Number
            sdf[!, column] = number2latex.(sdf[!, column]; enclose = false)
        end
    end
    sdf
end

pretty_table(
    st; backend = :latex, nosubheader = true, alignment = vcat([:l], repeat([:r], ncol(st) - 1))
)

# %%
# New Results DataFrame. Dataframe detailing new results obtained.
nrdf = let vudf = deepcopy(vudf), vu2019 = select(vu2019, :instance_name, :lb => :vu_lb, :ub => :vu_ub)
    used_columns = [
        :instance_name, :pricing_method, :heuristic_lb, :restricted_obj_value, :restricted_obj_bound,
        :obj_value, :obj_bound
    ]
    select!(vudf, used_columns) # Get only the used_columns
    for column in names(vudf) # Integralize all numerical columns
        column in ("instance_name", "pricing_method") && continue
        #@assert all(v -> abs(round(v, RoundNearest) - v) < 1e-6, filter(!ismissing, vudf[!, column]))
        vudf[!, column] = trunc.(Union{Int,Missing}, vudf[!, column] .+ 1e-6)
    end
    # Separate the priced and non-priced rows to after join them.
    priced = filter!(:pricing_method => isequal("furini"), deepcopy(vudf))
    select!(priced,
        :instance_name, :heuristic_lb => :ph_lb, :restricted_obj_value => :r_lb,
        :restricted_obj_bound => :r_ub, :obj_value => :p_lb, :obj_bound => :p_ub
    )
    non_priced = filter!(:pricing_method => isequal("none"), deepcopy(vudf))
    select!(non_priced,
        :instance_name, :heuristic_lb => :nph_lb, :obj_value => :np_lb, :obj_bound => :np_ub
    )
    vudf = innerjoin(priced, non_priced, vu2019; on = :instance_name)
    @assert @with(vudf, all(:ph_lb .== :nph_lb))
    select!(vudf, Not(:nph_lb), :ph_lb => :h_lb)
    @show names(vudf)
    vudf = @linq vudf |> where(
            (.!ismissing.(:p_lb) .& ((:p_lb .== :p_ub) .| (:p_lb .> :vu_lb) .| (:p_ub .< :vu_ub))) .|
            (.!ismissing.(:r_lb) .& ((:r_lb .== :r_ub) .| (:r_lb .> :vu_lb))) .|
            (.!ismissing.(:np_lb) .& ((:np_lb .== :np_ub) .| (:np_lb .> :vu_lb) .| (:np_ub .< :vu_ub)))
        )

    p_h_gap = @with vudf skipmissing(100.0 .- (:h_lb ./ :p_lb) .* 100.0)
    @show minimum(p_h_gap)
    @show maximum(p_h_gap)
    @show mean(p_h_gap)

    r_h_gap = @with vudf skipmissing(100.0 .- (:h_lb ./ :r_lb) .* 100.0)
    @show minimum(r_h_gap)
    @show maximum(r_h_gap)
    @show mean(r_h_gap)

    np_h_gap = @with vudf skipmissing(100.0 .- (:h_lb ./ :np_lb) .* 100.0)
    @show minimum(np_h_gap)
    @show maximum(np_h_gap)
    @show mean(np_h_gap)

    vudf
end

showtable(nrdf)

# %%
nrt = let nrdf = deepcopy(nrdf)
    select!(nrdf, Not(:h_lb)) # For now, let us disconsider the heuristic column.
    nrdf[!, :instance_name] = esc_latex.(nrdf[!, :instance_name])
    # Sometimes the restricted priced was stopped too early and gave terrible upper bounds,
    # this breaks the column layout and bring no useful information, better to display a dash instead.
    nrdf[!, :r_ub] .= replace(x -> (ismissing(x) ? missing : x > 100000 ? missing : x), nrdf[!, :r_ub])
    for column in names(nrdf) # latexify all numeric columns
        column == "instance_name" && continue
        nrdf[!, column] = number2latex.(nrdf[!, column]; enclose = false)
    end
    # First let us hightlight (with textbf), the best values among the LB, and the best values
    # among the UB.
    highlight_best_values!(nrdf;
        columns = [#=:h_lb,=# :r_lb, :p_lb, :np_lb, :vu_lb], chooser = find_index_allmax,
        cleaner = dirt2number
    )
    highlight_best_values!(nrdf;
        columns = [:p_ub, :np_ub, :vu_ub], chooser = find_index_allmin, cleaner = dirt2number
    )
    # Then, let us italicize the UBs and LBs from the same method, if they have the same value.
    select_all_if_equal(a) = all(e -> e == a[1], a) ? collect(Int, keys(a)) : Int[]
    latex_wrap(command, s) = "\\$(command){$(s)}"
    struct LatexWrapper; command :: String; end
    function (w::LatexWrapper)(s); return "\\$(w.command){$(s)}"; end
    latex_wrap(command) = LatexWrapper(command)
    highlight_best_values!(nrdf;
        columns = [:r_lb, :r_ub], chooser = select_all_if_equal,
        cleaner = dirt2number, changer = latex_wrap("underline")
    )
    highlight_best_values!(nrdf;
        columns = [:p_lb, :p_ub], chooser = select_all_if_equal,
        cleaner = dirt2number, changer = latex_wrap("underline")
    )
    highlight_best_values!(nrdf;
        columns = [:np_lb, :np_ub], chooser = select_all_if_equal,
        cleaner = dirt2number, changer = latex_wrap("underline")
    )
    highlight_best_values!(nrdf;
        columns = [:vu_lb, :vu_ub], chooser = select_all_if_equal,
        cleaner = dirt2number, changer = latex_wrap("underline")
    )
    sort!(nrdf, :instance_name)
    select!(nrdf,
        :instance_name => "Instance",
        #:h_lb => "H LB",
        :r_lb => "R LB",
        :p_lb => "P LB",
        :np_lb => "NP LB",
        :vu_lb => "VU LB",
        :r_ub => "R UB",
        :p_ub => "P UB",
        :np_ub => "NP UB",
        :vu_ub => "V UB",
    )
    select!(nrdf, names(nrdf) .=> esc_latex.(names(nrdf)))
    nrdf
end

pretty_table(
    nrt; backend = :latex, nosubheader = true, alignment = vcat([:l], repeat([:r], ncol(nrt) - 1))
)

# %%


text:
export computeMisfit

"""
function jInv.InverseSolve.computeMisfit(...)

Computes misfit for PDE parameter estimation problem. There are several
ways to use computeMisfit, some of which are parallelized.

Inputs:

    sig  - current model
    pMis - description of misfit term (MisfitParam, Array{MisfitParam}, or Array{Future})

Optional arguments:

    doDerivative - flag for computing derivatives (default: true)
    doClear      - flag to clear memory after computation (default: false)

Output:

   Dc    - current data
   F     - current misfit
   dF    - gradient
   d2F   - Hessian of misfit
   pMis  - modified misfit param (e.g., with PDE factorizations)
   times - some runtime statistics
"""
function computeMisfit(sig,
                       pMis::MisfitParam,doDerivative::Bool=true, doClear::Bool=false;
                       printProgress=false)
    if printProgress
        error("Print progress only works with multiple pFors")
    end
#=
 computeMisfit for a single forward problem. Everything is stored in memory on the node executing this function.
=#

    times = zeros(4)
    sigma,dsigma = pMis.modelfun(sig)
    times[1] = @elapsed   sigmaloc = interpGlobalToLocal(sigma,pMis.gloc.PForInv,pMis.gloc.sigmaBackground);
    times[2] = @elapsed   Dc,pMis.pFor  = getData(sigmaloc,pMis.pFor)      # fwd model to get predicted data
    times[3] = @elapsed   F,dF,d2F = pMis.misfit(Dc,pMis.dobs,pMis.Wd)
    if doDerivative
        times[4] = @elapsed dF = dsigma'*interpLocalToGlobal(getSensTMatVec(dF,sigmaloc,pMis.pFor),pMis.gloc.PForInv)
    end

    if doClear; clear!(pMis.pFor.Ainv); end
    return Dc,F,dF,d2F,pMis,times
end

function computeMisfit(sigmaRef::RemoteChannel,
                        pMisRef::RemoteChannel,
				      dFRef::RemoteChannel,
                  doDerivative,doClear::Bool=false)
#=
 computeMisfit for single forward problem

 Note: model (including interpolation matrix) and forward problems are RemoteRefs
=#

    rrlocs = [ pMisRef.where  dFRef.where]
    if !all(rrlocs .== myid())
        warn("computeMisfit: Problem on worker $(myid()) not all remote refs are stored here, but rrlocs=$rrlocs")
    end

    sigma = fetch(sigmaRef)
    pMis  = take!(pMisRef)

    Dc,F,dFi,d2F,pMis,times = computeMisfit(sigma,pMis,doDerivative,doClear)

    put!(pMisRef,pMis)
    # add to gradient
    if doDerivative
        dF = take!(dFRef)
        put!(dFRef,dF += dFi)
    end
    # put predicted data and d2F into remote refs (no need to communicate them)
    Dc  = remotecall(identity,myid(),Dc)
    d2F = remotecall(identity,myid(),d2F)

    return Dc,F,d2F,times
end

function computeMisfit(sigma,
	pMisRefs::Array{RemoteChannel,1},
	doDerivative::Bool=true;
	indCredit::AbstractVector=1:length(pMisRefs),
    printProgress::Bool=false)
#=
computeMisfit for multiple forward problems

This method runs in parallel (iff nworkers()> 1 )

Note: ForwardProblems and Mesh-2-Mesh Interpolation are RemoteRefs
    (i.e. they are stored in memory of a particular worker).
=#

    n = 1

	F   = 0.0
	dF  = (doDerivative) ? zeros(length(sigma)) : []
	d2F = Array{Any}(undef, length(pMisRefs));
	Dc  = Array{Future}(undef,size(pMisRefs))

	indDebit = []
	updateRes(Fi,idx) = (F+=Fi;push!(indDebit,idx))
	updateDF(x) = (dF+=x)

    workerList = []
    for k=indCredit
        push!(workerList,pMisRefs[k].where)
    end
    workerList = unique(workerList)
    sigRef = Array{RemoteChannel}(undef,maximum(workers()))
	dFiRef = Array{RemoteChannel}(undef,maximum(workers()))

	times = zeros(4);
	updateTimes(tt) = (times+=tt)

	@sync begin
		for p=workerList
			@async begin
				# communicate model and allocate RemoteRef for gradient
				sigRef[p] = initRemoteChannel(identity,p,sigma)   # send conductivity to workers
				dFiRef[p] = initRemoteChannel(zeros,p,length(sigma)) # get remote Ref to part of gradient
				# solve forward problems
				for idx in indCredit
					if pMisRefs[idx].where==p
						Dc[idx],Fi,d2F[idx],tt = remotecall_fetch(computeMisfit,p,sigRef[p],pMisRefs[idx],dFiRef[p],doDerivative)
						updateRes(Fi,idx)
						updateTimes(tt)
                        if printProgress && ((length(indDebit)/length(indCredit)) > n*0.1)
                            if doDerivative
                                println("Misfit and gradients computed for $(10*n)% of forward problems")
                            else
                                println("Misfit and gradients computed for $(10*n)% of forward problems")
                            end
                            n += 1
                        end
					end
				end

				# sum up gradients
				if doDerivative
					updateDF(fetch(dFiRef[p]))
				end
			end
		end
	end
	return Dc,F,dF,d2F,pMisRefs,times,indDebit
end

function computeMisfit(sigma,pMis::Array,doDerivative::Bool=true,indCredit=collect(1:length(pMis));
                       printProgress=false)
	#
	#	computeMisfit for multiple forward problems
	#
	#	This method runs in parallel (iff nworkers()> 1 )
	#
	#	Note: ForwardProblems and Mesh-2-Mesh Interpolation are stored on the main processor
	#		  and then sent to a particular worker, which returns an updated pFor.
	#
	numFor   = length(pMis)
 	F        = 0.0
    dF       = (doDerivative) ? zeros(length(sigma)) : []
 	d2F      = Array{Any}(undef,numFor)
 	Dc       = Array{Any}(undef,numFor)
	indDebit = []

	# draw next problem to be solved
	nextidx() = (idx = (isempty(indCredit)) ? -1 : pop!(indCredit))

 	updateRes(Fi,dFi,idx) = (F+=Fi; dF= (doDerivative) ? dF+dFi : []; push!(indDebit,idx))

	times = zeros(4);
	updateTimes(tt) = (times+=tt)

 	@sync begin
 		for p = workers()
 				@async begin
 					while true
 						idx = nextidx()
 						if idx == -1
 							break
 						end
 							Dc[idx],Fi,dFi,d2F[idx],pMis[idx],tt = remotecall_fetch(computeMisfit,p,sigma,pMis[idx],doDerivative)
 							updateRes(Fi,dFi,idx)
							updateTimes(tt)
 					end
 				end
 		end
 	end

 	return Dc,F,dF,d2F,pMis,times,indDebit
 end


text:
using Test
using StaticArrays
import ClimaCore.DataLayouts: IJFH, VF
import ClimaCore:
    Geometry, Fields, Domains, Topologies, Meshes, Spaces, Operators
using LinearAlgebra, IntervalSets

FT = Float64
domain = Domains.RectangleDomain(
    Geometry.XPoint{FT}(-pi)..Geometry.XPoint{FT}(pi),
    Geometry.YPoint{FT}(-pi)..Geometry.YPoint{FT}(pi);
    x1periodic = true,
    x2periodic = true,
)

Nq = 5
quad = Spaces.Quadratures.GLL{Nq}()

grid_mesh = Meshes.EquispacedRectangleMesh(domain, 17, 16)
grid_topology = Topologies.GridTopology(grid_mesh)
grid_space = Spaces.SpectralElementSpace2D(grid_topology, quad)
grid_coords = Geometry.Cartesian12Point.(Fields.coordinate_field(grid_space))

ts_mesh = Meshes.TensorProductMesh(domain, 17, 16)
ts_topology = Topologies.GridTopology(ts_mesh)
ts_space = Spaces.SpectralElementSpace2D(ts_topology, quad)
ts_coords = Geometry.Cartesian12Point.(Fields.coordinate_field(ts_space))

grid_test_setup = (grid_topology, grid_space, grid_coords)
ts_test_setup = (ts_topology, ts_space, ts_coords)

@testset "interpolate / restrict" begin

    for (topology, space, coords) in (grid_test_setup, ts_test_setup)
        INq = 9
        Iquad = Spaces.Quadratures.GLL{INq}()
        Ispace = Spaces.SpectralElementSpace2D(topology, Iquad)

        I = Operators.Interpolate(Ispace)
        R = Operators.Restrict(space)

        f = sin.(coords.x1 .+ 2 .* coords.x2)

        interpolated_field = I.(f)
        Spaces.weighted_dss!(interpolated_field)

        @test axes(interpolated_field).quadrature_style == Iquad
        @test axes(interpolated_field).topology == topology

        restrict_field = R.(f)
        Spaces.weighted_dss!(restrict_field)

        @test axes(restrict_field).quadrature_style == quad
        @test axes(restrict_field).topology == topology

        interp_restrict_field = R.(I.(f))
        Spaces.weighted_dss!(interp_restrict_field)

        @test axes(interp_restrict_field).quadrature_style == quad
        @test axes(interp_restrict_field).topology == topology

        @test norm(interp_restrict_field .- f) ≤ 3.0e-4
    end
end

@testset "gradient" begin

    for (topology, space, coords) in (grid_test_setup, ts_test_setup)
        f = sin.(coords.x1 .+ 2 .* coords.x2)

        grad = Operators.Gradient()
        gradf = grad.(f)
        Spaces.weighted_dss!(gradf)

        @test gradf ≈
              Geometry.Covariant12Vector.(
            Geometry.Cartesian12Vector.(
                cos.(coords.x1 .+ 2 .* coords.x2),
                2 .* cos.(coords.x1 .+ 2 .* coords.x2),
            ),
        ) rtol = 1e-2

        fv =
            Geometry.Cartesian12Vector.(
                sin.(coords.x1 .+ 2 .* coords.x2),
                cos.(coords.x1 .+ 2 .* coords.x2),
            )
        gradfv = grad.(fv)
        Spaces.weighted_dss!(gradfv)
        @test eltype(gradfv) <: Geometry.Axis2Tensor
    end
end

@testset "weak gradient" begin
    for (topology, space, coords) in (grid_test_setup, ts_test_setup)
        f = sin.(coords.x1 .+ 2 .* coords.x2)

        wgrad = Operators.WeakGradient()
        gradf = wgrad.(f)
        Spaces.weighted_dss!(gradf)

        @test Geometry.Cartesian12Vector.(gradf) ≈
              Geometry.Cartesian12Vector.(
            cos.(coords.x1 .+ 2 .* coords.x2),
            2 .* cos.(coords.x1 .+ 2 .* coords.x2),
        ) rtol = 1e-2
    end
end

@testset "curl" begin
    for (topology, space, coords) in (grid_test_setup, ts_test_setup)
        v =
            Geometry.Cartesian12Vector.(
                sin.(coords.x1 .+ 2 .* coords.x2),
                cos.(3 .* coords.x1 .+ 4 .* coords.x2),
            )

        curl = Operators.Curl()
        curlv = curl.(v)
        Spaces.weighted_dss!(curlv)
        curlv_ref =
            Geometry.Contravariant3Vector.(
                .-3 .* sin.(3 .* coords.x1 .+ 4 .* coords.x2) .-
                2 .* cos.(coords.x1 .+ 2 .* coords.x2),
            )

        @test curlv ≈ curlv_ref rtol = 1e-2
    end
end

@testset "curl-curl" begin
    for (topology, space, coords) in (grid_test_setup, ts_test_setup)
        v =
            Geometry.Cartesian12Vector.(
                sin.(coords.x1 .+ 2 .* coords.x2),
                cos.(3 .* coords.x1 .+ 2 .* coords.x2),
            )
        curlv_ref =
            .-3 .* sin.(3 .* coords.x1 .+ 2 .* coords.x2) .-
            2 .* cos.(coords.x1 .+ 2 .* coords.x2)
        curlcurlv_ref1 =
            .-6 .* cos.(3 .* coords.x1 .+ 2 .* coords.x2) .+
            4 .* sin.(coords.x1 .+ 2 .* coords.x2)
        curlcurlv_ref2 =
            9 .* cos.(3 .* coords.x1 .+ 2 .* coords.x2) .-
            2 .* sin.(coords.x1 .+ 2 .* coords.x2)

        curl = Operators.Curl()
        curlcurlv =
            Geometry.Cartesian12Vector.(
                curl.(Geometry.Covariant3Vector.(curl.(v))),
            )
        Spaces.weighted_dss!(curlcurlv)

        @test curlcurlv ≈
              Geometry.Cartesian12Vector.(curlcurlv_ref1, curlcurlv_ref2) rtol =
            4e-2
    end
end

@testset "weak curl-strong curl" begin
    for (topology, space, coords) in (grid_test_setup, ts_test_setup)
        v =
            Geometry.Cartesian12Vector.(
                sin.(coords.x1 .+ 2 .* coords.x2),
                cos.(3 .* coords.x1 .+ 2 .* coords.x2),
            )
        curlv_ref =
            .-3 .* sin.(3 .* coords.x1 .+ 2 .* coords.x2) .-
            2 .* cos.(coords.x1 .+ 2 .* coords.x2)
        curlcurlv_ref1 =
            .-6 .* cos.(3 .* coords.x1 .+ 2 .* coords.x2) .+
            4 .* sin.(coords.x1 .+ 2 .* coords.x2)
        curlcurlv_ref2 =
            9 .* cos.(3 .* coords.x1 .+ 2 .* coords.x2) .-
            2 .* sin.(coords.x1 .+ 2 .* coords.x2)

        curl = Operators.Curl()
        wcurl = Operators.WeakCurl()
        curlcurlv =
            Geometry.Cartesian12Vector.(
                wcurl.(Geometry.Covariant3Vector.(curl.(v))),
            )
        Spaces.weighted_dss!(curlcurlv)

        @test curlcurlv ≈
              Geometry.Cartesian12Vector.(curlcurlv_ref1, curlcurlv_ref2) rtol =
            4e-2
    end
end

@testset "weak curl" begin
    for (topology, space, coords) in (grid_test_setup, ts_test_setup)
        v =
            Geometry.Cartesian12Vector.(
                sin.(coords.x1 .+ 2 .* coords.x2),
                cos.(3 .* coords.x1 .+ 4 .* coords.x2),
            )

        wcurl = Operators.WeakCurl()
        curlv = wcurl.(v)
        Spaces.weighted_dss!(curlv)
        curlv_ref =
            .-3 .* sin.(3 .* coords.x1 .+ 4 .* coords.x2) .-
            2 .* cos.(coords.x1 .+ 2 .* coords.x2)

        @test curlv ≈ Geometry.Contravariant3Vector.(curlv_ref) rtol = 1e-2
    end
end

@testset "div" begin
    for (topology, space, coords) in (grid_test_setup, ts_test_setup)
        v =
            Geometry.Cartesian12Vector.(
                sin.(coords.x1 .+ 2 .* coords.x2),
                cos.(3 .* coords.x1 .+ 2 .* coords.x2),
            )

        div = Operators.Divergence()
        divv = div.(v)
        Spaces.weighted_dss!(divv)
        divv_ref =
            cos.(coords.x1 .+ 2 .* coords.x2) .-
            2 .* sin.(3 .* coords.x1 .+ 2 .* coords.x2)

        @test divv ≈ divv_ref rtol = 1e-2
    end
end

@testset "weak div" begin
    for (topology, space, coords) in (grid_test_setup, ts_test_setup)
        v =
            Geometry.Cartesian12Vector.(
                sin.(coords.x1 .+ 2 .* coords.x2),
                cos.(3 .* coords.x1 .+ 2 .* coords.x2),
            )

        wdiv = Operators.WeakDivergence()
        divv = wdiv.(v)
        Spaces.weighted_dss!(divv)
        divv_ref =
            cos.(coords.x1 .+ 2 .* coords.x2) .-
            2 .* sin.(3 .* coords.x1 .+ 2 .* coords.x2)

        @test divv ≈ divv_ref rtol = 1e-2
    end
end

@testset "annhilator property: curl-grad" begin
    for (topology, space, coords) in (grid_test_setup, ts_test_setup)
        f = sin.(coords.x1 .+ 2 .* coords.x2)

        grad = Operators.Gradient()
        gradf = grad.(f)
        Spaces.weighted_dss!(gradf)

        curl = Operators.Curl()
        curlgradf = curl.(gradf)
        Spaces.weighted_dss!(curlgradf)

        @test norm(curlgradf) < 1e-12
    end
end

@testset "annhilator property: div-curl" begin
    for (topology, space, coords) in (grid_test_setup, ts_test_setup)
        v = Geometry.Covariant3Vector.(sin.(coords.x1 .+ 2 .* coords.x2))
        curl = Operators.Curl()
        curlv = curl.(v)
        Spaces.weighted_dss!(curlv)

        div = Operators.Divergence()
        divcurlv = div.(curlv)
        Spaces.weighted_dss!(divcurlv)

        @test norm(divcurlv) < 1e-12
    end
end

@testset "scalar hyperdiffusion" begin
    for (topology, space, coords) in (grid_test_setup, ts_test_setup)
        k = 2
        l = 3
        y = @. sin(k * coords.x1 + l * coords.x2)
        ∇⁴y_ref = @. (k^2 + l^2)^2 * sin(k * coords.x1 + l * coords.x2)

        wdiv = Operators.WeakDivergence()
        grad = Operators.Gradient()
        χ = Spaces.weighted_dss!(@. wdiv(grad(y)))
        ∇⁴y = Spaces.weighted_dss!(@. wdiv(grad(χ)))

        @test ∇⁴y_ref ≈ ∇⁴y rtol = 2e-2
    end
end

@testset "vector hyperdiffusion" begin
    for (topology, space, coords) in (grid_test_setup, ts_test_setup)
        k = 2
        l = 3
        y = @. Geometry.Cartesian12Vector(
            sin(k * coords.x1 + l * coords.x2),
            0.0,
        )
        ∇⁴y_ref = @. Geometry.Cartesian12Vector(
            (k^2 + l^2)^2 * sin(k * coords.x1 + l * coords.x2),
            0.0,
        )
        curl = Operators.Curl()
        wcurl = Operators.WeakCurl()

        sdiv = Operators.Divergence()
        wgrad = Operators.WeakGradient()

        χ = Spaces.weighted_dss!(
            @. Geometry.Cartesian12Vector(wgrad(sdiv(y))) -
               Geometry.Cartesian12Vector(
                wcurl(Geometry.Covariant3Vector(curl(y))),
            )
        )
        ∇⁴y = Spaces.weighted_dss!(
            @. Geometry.Cartesian12Vector(wgrad(sdiv(χ))) -
               Geometry.Cartesian12Vector(
                wcurl(Geometry.Covariant3Vector(curl(χ))),
            )
        )

        @test ∇⁴y_ref ≈ ∇⁴y rtol = 2e-2
    end
end


text:
using ChaosTools
using Test, StaticArrays

println("\nTesting orbit diagrams...")

@testset "Orbit Diagram" begin
    @testset "Discrete 1D" begin
        ds = Systems.logistic()
        i = 1
        parameter = 1
        pvalues = 2:0.01:4
        ics = [rand() for m in 1:10]
        n = 50
        Ttr = 5000
        output = orbitdiagram(ds, i, parameter, pvalues; n = n, Ttr = Ttr)
        @test output[1][1] ≈ 0.5
        @test output[2][1] == 0.5024875621890547
        @test output[end][1] != output[end][2] != output[end][3]
    end

    @testset "Discrete 2D" begin
        ds = Systems.standardmap()
        i = 2
        parameter = 1
        pvalues = 0:0.005:2
        ics = [0.001rand(2) for m in 1:10]
        n = 50
        Ttr = 5000
        output = orbitdiagram(ds, i, parameter, pvalues; n = n, Ttr = Ttr)
        @test length(output[1]) == n
    end

    @testset "Discrete IIP" begin
        ds = Systems.henon_iip()
        i = 2
        parameter = 1
        pvalues = [1.4, 1.5]
        ics = [rand(2), rand(2)]
        n = 50
        Ttr = 5000
        output = orbitdiagram(ds, i, parameter, pvalues; n = n, Ttr = Ttr)
        @test length(output[1]) == n
    end
end


text:
#
println("Da) Test of the PhotoEmission module with ASF from an internally generated initial- and final-state multiplet.")
#
# setDefaults("unit: rate", "a.u.")
setDefaults("unit: rate", "1/s")
setDefaults("print summary: open", "zzz-radiative.sum")

if  true

    grid = Radial.Grid(true)
    setDefaults("standard grid", grid)
    defaultsSettings = PhotoEmission.Settings()
    photoSettings = PhotoEmission.Settings(defaultsSettings, multipoles=[E1, M1, E2, M2], gauges=[UseCoulomb, UseBabushkin], printBefore=true)

    comp = Atomic.Computation(Atomic.Computation(), name="Energies and Einstein coefficients for the spectrum Fe X",
              grid=grid, nuclearModel=Nuclear.Model(26.);
              initialConfigs = [Configuration("[Ne] 3s 3p^6"), Configuration("[Ne] 3s^2 3p^4 3d")],
              finalConfigs   = [Configuration("[Ne] 3s^2 3p^5")],
              process = Radiative(), processSettings = photoSettings );
    @show comp
    perform(comp)

elseif false
    #== wa = Atomic.Computation("xx",  Nuclear.Model(36.);
                            initialConfigs=[Configuration("1s 2s^2"), Configuration("1s 2s 2p"), Configuration("1s 2p^2")],
                            finalConfigs  =[Configuration("1s 2s^2"), Configuration("1s 2s 2p"), Configuration("1s 2p^2")],
                            process = Radiative(),
                            processSettings=PhotoEmission.Settings([E1, M1, E2, M2], [JAC.UseCoulomb, JAC.UseBabushkin], true, true,
                                true, Tuple{Int64,Int64}[(5,0), (7,0), (10,0), (11,0), (12,0), (13,0), (14,0), (15,0), (16,0)], 0., 0., 10000. ) )  ==#
elseif false
    wa = Atomic.Computation(Atomic.Computation(), name="Energies and Einstein coefficients for the spectrum Fe X",
                            grid=JAC.Radial.Grid(true), nuclearModel=Nuclear.Model(26.),
                            initialConfigs = [Configuration("[Ne] 3s 3p^6"), Configuration("[Ne] 3s^2 3p^4 3d")],
                            finalConfigs   = [Configuration("[Ne] 3s^2 3p^5")],
                            process = Radiative(),
                            processSettings = PhotoEmission.Settings([E1, M1, E2, M2], [JAC.UseCoulomb, JAC.UseBabushkin],
                                              false, true, false, Tuple{Int64,Int64}[], 0., 0., 10000. ) );
    wb = @time( perform(wa) )
elseif true
    # Test for dielectronic computations
    wa = Atomic.Computation(Atomic.Computation(), name="Energies and Einstein rates for DR configurations",
                            grid=JAC.Radial.Grid(true), nuclearModel=Nuclear.Model(26.),
                            initialConfigs = [Configuration("1s 2s^2 2p")],
                            finalConfigs   = [Configuration("1s^2 2s^2")],
                            process = Radiative(),
                            processSettings = PhotoEmission.Settings([E1, M1], [JAC.UseCoulomb, JAC.UseBabushkin],
                                              false, true, true, Tuple{Int64,Int64}[(2,0)], 0., 0., 10000. ) );
    wb = @time( perform(wa) )
end
setDefaults("print summary: close", "")

## , Configuration("[Ne] 3s^2 3p^4 3d")  Configuration("[Ne] 3s 3p^6")