text:
# OpenCL.Context

type Context <: CLObject
    id :: CL_context

    function Context(ctx_id::CL_context; retain=false)
        if retain
            @check api.clRetainContext(ctx_id)
        end
        ctx = new(ctx_id)
        finalizer(ctx, c -> begin
            retain || _deletecached!(c);
            if c.id != C_NULL
                @check api.clReleaseContext(c.id)
                c.id = C_NULL
            end
        end )
        return ctx
    end
end

Base.pointer(ctx::Context) = ctx.id

function Base.show(io::IO, ctx::Context)
    dev_strs = [replace(d[:name], r"\s+", " ") for d in devices(ctx)]
    devs_str = join(dev_strs, ",")
    ptr_val = convert(UInt, Base.pointer(ctx))
    ptr_address = "0x$(hex(ptr_val, Sys.WORD_SIZE>>2))"
    print(io, "OpenCL.Context(@$ptr_address on $devs_str)")
end

immutable _CtxErr
    handle :: Ptr{Void}
    err_info :: Ptr{Cchar}
    priv_info :: Ptr{Void}
    cb :: Csize_t
end

function ctx_notify_err(err_info::Ptr{Cchar}, priv_info::Ptr{Void},
                        cb::Csize_t, payload::Ptr{Void})
    ptr = convert(Ptr{_CtxErr}, payload)
    handle = unsafe_load(ptr, 1).handle

    val = _CtxErr(handle, err_info, priv_info, cb)
    unsafe_store!(ptr, val, 1)

    ccall(:uv_async_send, Void, (Ptr{Void},), handle)
    nothing
end

ctx_callback_ptr() = cfunction(ctx_notify_err, Void,
                              (Ptr{Cchar}, Ptr{Void}, Csize_t, Ptr{Void}))

function raise_context_error(error_info, private_info)
    throw(OpenCLException("OpenCL.Context error: $error_info"))
end

function Context(devs::Vector{Device};
                         properties=nothing,
                         callback::Union{Function, Void} = nothing)
    if isempty(devs)
        ArgumentError("No devices specified for context")
    end
    if properties !== nothing
        ctx_properties = _parse_properties(properties)
    else
        ctx_properties = C_NULL
    end

    n_devices = length(devs)
    device_ids = Array(CL_device_id, n_devices)
    for (i, d) in enumerate(devs)
        device_ids[i] = d.id
    end

    cb = Base.AsyncCondition()
    ctx_user_data = Ref(_CtxErr(Base.unsafe_convert(Ptr{Void}, cb), 0, 0, 0))

    err_code = Ref{CL_int}()
    ctx_id = api.clCreateContext(ctx_properties, n_devices, device_ids,
                                 ctx_callback_ptr(), ctx_user_data, err_code)
    if err_code[] != CL_SUCCESS
        throw(CLError(err_code[1]))
    end

    true_callback = callback == nothing ? raise_context_error : callback :: Function

    @async begin
        try
            Base.wait(cb)
            err = ctx_user_data[]
            error_info = unsafe_string(err.err_info)
            private_info = unsafe_string(err.priv_info)
            true_callback(error_info, private_info)
        catch
            rethrow()
        finally
            Base.close(cb)
        end
    end

    return Context(ctx_id)
end

Context(d::Device; properties=nothing, callback=nothing) =
        Context([d], properties=properties, callback=callback)

function Context(dev_type::CL_device_type;
                 properties=nothing, callback=nothing)
    if properties !== nothing
        ctx_properties = _parse_properties(properties)
    else
        ctx_properties = C_NULL
    end
    if callback !== nothing
        ctx_user_data = callback
    else
        ctx_user_data = raise_context_error
    end
    err_code = Array(CL_int, 1)
    ctx_id = api.clCreateContextFromType(ctx_properties, dev_type,
                                         ctx_callback_ptr(), ctx_user_data, err_code)
    if err_code[1] != CL_SUCCESS
        throw(CLError(err_code[1]))
    end
    return Context(ctx_id)
end

function Context(dev_type::Symbol;
                 properties=nothing, callback=nothing)
    Context(cl_device_type(dev_type),
            properties=properties, callback=callback)
end

function properties(ctx_id::CL_context)
    nbytes = Csize_t[0]
    @check api.clGetContextInfo(ctx_id, CL_CONTEXT_PROPERTIES, 0, C_NULL, nbytes)

    # Calculate length of storage array
    # At nbytes[1] the size of the properties array in bytes is stored
    # The length of the property array is then nbytes[1] / sizeof(CL_context_properties)
    # Note: nprops should be odd since it requires a C_NULL terminated array
    nprops = div(nbytes[1], sizeof(CL_context_properties))

    props = Array(CL_context_properties, nprops)
    @check api.clGetContextInfo(ctx_id, CL_CONTEXT_PROPERTIES,
                                nbytes[1], props, C_NULL)
    #properties array of [key,value..., C_NULL]
    result = Any[]
    for i in 1:2:nprops
        key = props[i]
        value = i < nprops ? props[i+1] : nothing

        if key == CL_CONTEXT_PLATFORM
            push!(result, (key, Platform(cl_platform_id(value))))
        elseif key == CL_GL_CONTEXT_KHR ||
           key == CL_EGL_DISPLAY_KHR ||
           key == CL_GLX_DISPLAY_KHR ||
           key == CL_WGL_HDC_KHR ||
           key == CL_CGL_SHAREGROUP_KHR
            push!(result, (key, value))
        elseif is_apple() ? (key == CL_CONTEXT_PROPERTY_USE_CGL_SHAREGROUP_APPLE) : false
            push!(result, (key, value))
        elseif key == 0
            if i != nprops
                warn("Encountered OpenCL.Context property key == 0 at position $i")
            end
            break
        else
            warn("Unknown OpenCL.Context property key encountered $key")
        end
    end
    return result
end

function properties(ctx::Context)
    properties(ctx.id)
end

#Note: properties list needs to be terminated with a NULL value!
function _parse_properties(props)
    if isempty(props)
        return C_NULL
    end
    cl_props = CL_context_properties[]
    for prop_tuple in props
        if length(prop_tuple) != 2
            throw(ArgumentError("Context property tuples must be of type (key, value)"))
        end
        prop, val = prop_tuple
        push!(cl_props, cl_context_properties(prop))
        if prop == CL_CONTEXT_PLATFORM
            isa(val, Platform) && (val = val.id)
            push!(cl_props, cl_context_properties(val))
        elseif prop == CL_WGL_HDC_KHR
            push!(cl_props, cl_context_properties(val))
        elseif is_apple() ? (prop == CL_CONTEXT_PROPERTY_USE_CGL_SHAREGROUP_APPLE) : false
            push!(cl_props, cl_context_properties(val))
        elseif prop == CL_GL_CONTEXT_KHR ||
            prop == CL_EGL_DISPLAY_KHR ||
            prop == CL_GLX_DISPLAY_KHR ||
            prop == CL_CGL_SHAREGROUP_KHR
            push!(cl_props, cl_context_properties(val))
        else
            throw(OpenCLException("Invalid OpenCL Context property"))
        end
    end
    push!(cl_props, cl_context_properties(C_NULL))
    return cl_props
end

function num_devices(ctx::Context)
    ndevices = Array(CL_uint, 1)
    @check api.clGetContextInfo(ctx.id, CL_CONTEXT_NUM_DEVICES,
                                sizeof(CL_uint), ndevices, C_NULL)
    return ndevices[1]
end

function devices(ctx::Context)
    n = num_devices(ctx)
    if n == 0
        return []
    end
    dev_ids = Array(CL_device_id, n)
    @check api.clGetContextInfo(ctx.id, CL_CONTEXT_DEVICES,
                                n * sizeof(CL_device_id), dev_ids, C_NULL)
    return [Device(id) for id in dev_ids]
end

function create_some_context()
    if isempty(platforms())
        throw(OpenCLException("No OpenCL.Platform available"))
    end
    gpu_devices = devices(:gpu)
    if !isempty(gpu_devices)
        for dev in gpu_devices
            local ctx::Context
            try
                ctx = Context(dev)
            catch
                continue
            end
            return ctx
        end
    end
    cpu_devices = devices(:cpu)
    if !isempty(cpu_devices)
        for dev in cpu_devices
            local ctx::Context
            try
                ctx = Context(dev)
            catch
                continue
            end
            return ctx
        end
    end
    if isempty(gpu_devices) && isempty(cpu_devices)
        throw(OpenCLException("Unable to create any OpenCL.Context, no available devices"))
    else
        throw(OpenCLException("Unable to create any OpenCL.Context, no devices worked"))
    end
end


text:
using ArgParse
using TMLEEpistasis

function parse_commandline()
    s = ArgParseSettings(
        description = "This program computes estimates of Genetic Epistatis from the UK-Biobank using TMLE."*
                      "Here is a list of the arguments that should be provided, you can also have a look at the "*
                      "test/data and test/config folders to see some examples.",
        commands_are_required = false,
        version = "0.2",
        add_version = true)

    @add_arg_table s begin
        "genotypes"
            help = "Path to genotypes .csv file"
            required = true
        "phenotypes"
            help = "A file (.csv format). The first row contains the column names with `eid` the sample ID"*
                   " and the rest of the columns are phenotypes of interest."
            required = true
        "confounders"
            help = "A file (.csv format) containing the confounding variables values and the sample ids associated"*
                   " with the participants. The first line of the file should contain the columns names and the sample ids "*
                   " column name should be: `SAMPLE_ID`."
            required = true
        "queries"
            help = "A file (.toml format) see: config/sample_query.toml for more information"
            required = true
        "estimator"
            help = "A file (.toml format) describing the tmle estimator to use, see config/sample_estimator.toml"*
                   " for a basic example."
            required = true
        "out"
            help = "Path where the ouput will be saved"
            required = true
        "--target-type", "-t"
            help = "The type of the target variable: Real or Bool"
            arg_type = String
            default = "Bool"
        "--phenotypes-list", "-p"
            help = "A file, one line for each phenotype, containing a restrictions of the phenotypes "*
                   "to consider for the analysis."
            required = false
            arg_type = String
        "--adaptive-cv", "-a"
            help = "Adaptively selects the number of folds used in cross validation and overrides the default used in the estimator file."
            action = :store_true
        "--save-full", "-f"
            help = "Also save the full TMLE estimators for each phenotype."
            action = :store_true
        "--verbosity", "-v"
            help = "Verbosity level"
            arg_type = Int
            default = 1
    end

    return parse_args(s)
end

parsed_args = parse_commandline()

UKBBVariantRun(parsed_args)

ext: jl


text:
# This file investigates MSE with respect to free variables

using FractalTools
using GeometryBasics
using Makie

# Construct interpolation data
f(x, y) = [
    x^2 + y^2 + 1,
    x^2 - y^2
    ]
# ngon = Triangle(
#     Point(BigFloat(0.), BigFloat(0.)),
#     Point(BigFloat(1.), BigFloat(0.)),
#     Point(BigFloat(0.5), BigFloat(1.)))

ngon = Triangle(
    Point(BigFloat(-5.), BigFloat(-5.)),
    Point(BigFloat(5), BigFloat(-5.)),
    Point(BigFloat(0), BigFloat(5.)))
npts = 100
pts = getdata(f, ngon, npts)

# Construct test data
tpts = getdata(ngon, npts)
ntpts = length(tpts)

# Compute errors
fvals = getindex.(map(pt -> f(pt...), tpts), 1)
freevars = 0.001 : 0.001 : 0.025
mse = map(freevars) do freevar
    interp = interpolate(pts, HInterp2D(fill(freevar, 2, 2)))
    ivals = getindex.(map(pt -> interp(pt...), tpts), 1)
    sum((fvals - ivals).^2) / ntpts
end

# Plot mse
fig = Figure()
ax = fig[1, 1] = Axis(fig, xlabel="Free Variable", ylabel="MSE", title="2D Hidden Interpolation MSE")
stem!(ax, freevars, mse, color=:black)
save(joinpath(@__DIR__, "hinterp2d_mse.png"), fig)
display(fig)


text:
using Documenter, ModelingToolkitStandardLibrary
using ModelingToolkitStandardLibrary.Blocks
using ModelingToolkitStandardLibrary.Mechanical
using ModelingToolkitStandardLibrary.Mechanical.Rotational
using ModelingToolkitStandardLibrary.Magnetic
using ModelingToolkitStandardLibrary.Magnetic.FluxTubes
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Thermal

include("pages.jl")

makedocs(
    sitename="ModelingToolkitStandardLibrary.jl",
    authors="Julia Computing",
    clean=true,
    doctest=false,
    modules=[ModelingToolkitStandardLibrary,
             ModelingToolkitStandardLibrary.Blocks,
             ModelingToolkitStandardLibrary.Mechanical,
             ModelingToolkitStandardLibrary.Mechanical.Rotational,
             ModelingToolkitStandardLibrary.Magnetic,
             ModelingToolkitStandardLibrary.Magnetic.FluxTubes,
             ModelingToolkitStandardLibrary.Electrical,
             ModelingToolkitStandardLibrary.Thermal],

    format=Documenter.HTML(assets=["assets/favicon.ico"],
                           canonical="https://mtkstdlib.sciml.ai/stable/"),

    pages=pages
)

deploydocs(
    repo="github.com/SciML/ModelingToolkitStandardLibrary.jl";
    push_preview=true
)


text:
# BioSymbols.jl
# =============
#
# Types representing Biological Symbols.
#
# This file is a part of the BioJulia ecosystem.
# License is MIT: https://github.com/BioJulia/NucleicAcids.jl/blob/master/LICENSE.md

__precompile__()

module BioSymbols

export
    # Types
    BioSymbol,
    NucleicAcid,
    DNA,
    RNA,
    AminoAcid,

    # Constants
    DNA_A,
    DNA_C,
    DNA_G,
    DNA_T,
    DNA_M,
    DNA_R,
    DNA_W,
    DNA_S,
    DNA_Y,
    DNA_K,
    DNA_V,
    DNA_H,
    DNA_D,
    DNA_B,
    DNA_N,
    DNA_Gap,
    ACGT,
    ACGTN,
    RNA_A,
    RNA_C,
    RNA_G,
    RNA_U,
    RNA_M,
    RNA_R,
    RNA_W,
    RNA_S,
    RNA_Y,
    RNA_K,
    RNA_V,
    RNA_H,
    RNA_D,
    RNA_B,
    RNA_N,
    RNA_Gap,
    ACGU,
    ACGUN,
    AA_A,
    AA_R,
    AA_N,
    AA_D,
    AA_C,
    AA_Q,
    AA_E,
    AA_G,
    AA_H,
    AA_I,
    AA_L,
    AA_K,
    AA_M,
    AA_F,
    AA_P,
    AA_S,
    AA_T,
    AA_W,
    AA_Y,
    AA_V,
    AA_O,
    AA_U,
    AA_B,
    AA_J,
    AA_Z,
    AA_X,
    AA_Term,
    AA_Gap,

    # Bit Operations
    gap,
    isGC,
    ispurine,
    ispyrimidine,
    isambiguous,
    iscertain,
    isgap,
    complement,
    iscompatible,
    compatbits,
    alphabet

import Automa
import Automa.RegExp: @re_str

"""
The BioSymbol type is an abstract type that represents
any kind of biological symbol that may appear in data
such as biological sequences, SNP datasets and more.

Every abstract or concrete subtype of BioSymbol is
expected to have the following methods defined:

`isterm`
`bytemask`
`prefix`
`type_text`

"""
abstract type BioSymbol end

Base.length(::BioSymbol) = 1
Base.iterate(sym::BioSymbol) = (sym, nothing)
Base.iterate(sym::BioSymbol, state) = nothing

include("nucleicacid.jl")
include("aminoacid.jl")

"""
    isgap(symbol::BioSymbol)

Test if `symbol` is a gap.
"""
isgap(symbol::BioSymbol) = symbol == gap(typeof(symbol))

# Arithmetic and Order
# --------------------

# These methods are necessary when deriving some algorithims
# like iteration, sort, comparison, and so on.

Base.:-(x::S, y::S) where S <: BioSymbol = convert(Int, x) - convert(Int, y)

Base.:+(x::Integer, y::BioSymbol) = y + x

Base.isless(x::S, y::S) where S <: BioSymbol = isless(convert(Integer, x), convert(Integer, y))

function Base.:~(x::BioSymbol)
    return reinterpret(typeof(x), ~reinterpret(UInt8, x) & bytemask(x))
end

function Base.:&(x::S, y::S) where S <: BioSymbol
    return reinterpret(S, convert(Integer, x) & convert(Integer, y))
end

@inline function Base.count_ones(symbol::BioSymbol)
    return count_ones(convert(Integer, symbol))
end

@inline function Base.trailing_zeros(symbol::BioSymbol)
    return trailing_zeros(convert(Integer, symbol))
end

# Printing BioSymbols
# -------------------

function suffix(symbol::BioSymbol)
    if isterm(symbol)
        return "_Term"
    end
    if isgap(symbol)
        return "_Gap"
    end
    return string('_', convert(Char, symbol))
end

function Base.show(io::IO, symbol::BioSymbol)
    if isvalid(symbol)
        write(io, prefix(symbol), suffix(symbol))
    else
        write(io, "Invalid ", type_text(symbol))
    end
    return nothing
end

function Base.print(io::IO, symbol::BioSymbol)
    if !isvalid(symbol)
        throw(ArgumentError("Invalid $(type_text(symbol))"))
    end
    write(io, convert(Char, symbol))
    return nothing
end

Base.write(io::IO, symbol::BioSymbol) = write(io, convert(Integer, symbol))
Base.read(io::IO, ::Type{T}) where T<:BioSymbol = reinterpret(T, read(io, UInt8))

"""
    iscompatible(x::S, y::S) where S <: BioSymbol

Test if `x` and `y` are compatible with each other.

Examples
--------

```jldoctest
julia> iscompatible(AA_A, AA_R)
false

julia> iscompatible(AA_A, AA_X)
true

julia> iscompatible(DNA_A, DNA_A)
true

julia> iscompatible(DNA_C, DNA_N)  # DNA_N can be DNA_C
true

julia> iscompatible(DNA_C, DNA_R)  # DNA_R (A or G) cannot be DNA_C
false

```
"""
@inline function iscompatible(x::S, y::S) where S <: BioSymbol
    return compatbits(x) & compatbits(y) != 0
end

end


text:
using Test: include
@testset "31.next-permutation.jl" begin
    @test next_permutation([1,1,5]) == [1, 5, 1]
    @test next_permutation([1,1,5]) == [1, 5, 1]
    @test next_permutation([1, 1, 1]) == [1, 1, 1]
    @test next_permutation([1]) == [1]
    @test next_permutation([1,2,3]) == [1, 3, 2]
    @test next_permutation(Int[]) == Int[]
end


text:
export ishom, ishet, _ishom, _ishet, isbiallelic, isbinary

"""
    isbiallelic(data::GenoArray)
Returns `true` if the `GenoArray` is biallelic, `false` if not.
"""
@inline function isbiallelic(data::T) where T<:GenoArray
    length(unique(Base.Iterators.flatten(skipmissing(data)))) == 2
end

# method for using just a genodata input
function isbiallelic(data::DataFrame)
    tmp = sort(data, [:locus, :name])
    prt = Base.Iterators.partition(tmp.genotype, length(data.name.pool))
    bi = findfirst(!isbiallelic, collect(prt))
    isnothing(bi) ? true : false
end

"""
    isbiallelic(data::PopData)
Returns `true` all the loci in the `PopData` are biallelic, `false` if not.
"""
isbiallelic(data::PopData) = data.metadata.biallelic

#TODO how to treat haploids?
"""
```
ishom(locus::T) where T <: GenoArray
ishom(locus::T) where T<:Base.SkipMissing
ishom(locus::Genotype)
ishom(locus::Missing)
```
A series of methods to test if a locus or loci are homozygous and return `true` if
it is, `false` if it isn't (or missing). For calculations, we recommend using `_ishom()`,
which returns `missing` if the genotype is `missing`. The vector version
simply maps the function over the elements.
"""
ishom(locus::Genotype) = all(@inbounds first(locus) .== locus)

# public facing method
ishom(locus::Missing) = false
ishom(locus::T) where T<:GenoArray = @inbounds map(ishom, locus)
ishom(locus::T) where T<:Base.SkipMissing = @inbounds map(ishom, locus)

# API computational methods
_ishom(locus::Genotype) = all(@inbounds first(locus) .== locus)
_ishom(locus::T) where T<:GenoArray = @inbounds map(_ishom, locus)
_ishom(locus::T) where T<:Base.SkipMissing = @inbounds map(_ishom, locus)
_ishom(locus::Missing) = missing

#= scales for size, which isn't super necessary yet
function ishom(geno::Genotype)
    x = true
    y = 2
    while x
        x = geno[1] == geno[y]
        y += 1
    end
    x
end
=#2
"""
    ishom(locus::Genotype, allele::Signed)
    ishom(loci::GenoArray, allele::Signed)
Returns `true` if the `locus`/`loci` is/are homozygous for the specified `allele`.
"""
function ishom(geno::T, allele::U) where T<:Genotype where U<:Integer
    ∈(allele, geno) & ishom(geno) ? true : false
end

ishom(geno::T, allele::U) where T<:GenoArray where U<:Integer = map(i -> ishom(i, allele), geno)

# public facing method
ishom(geno::Missing, allele::U) where U<:Integer = false

# API computational method
_ishom(geno::T, allele::U) where T<:Genotype where U<:Integer = ∈(allele, geno) & _ishom(geno) ? true : false
_ishom(geno::T, allele::U) where T<:GenoArray where U<:Integer = map(i -> _ishom(i, allele), geno)
_ishom(geno::Missing, allele::U) where U<:Integer = missing

# public facing method
"""
```
ishet(locus::T) where T <: GenoArray
ishet(locus::Genotype)
ishet(locus::Missing)
```
A series of methods to test if a locus or loci are heterozygous and return `true` if
it is, `false` if it isn't (or missing). For calculations, we recommend using `_ishet()`,
which returns `missing` if the genotype is `missing`. The vector version
simply maps the function over the elements.
"""
ishet(locus::Genotype) = !ishom(locus)
ishet(locus::Missing) = false
ishet(locus::T) where T<:GenoArray = @inbounds map(ishet, locus)
ishet(locus::T) where T<:Base.SkipMissing = @inbounds map(ishet, locus)

_ishet(locus::Genotype) = !_ishom(locus)
_ishet(locus::Missing) = missing
_ishet(locus::T) where T<:GenoArray = @inbounds map(_ishet, locus)
_ishet(locus::T) where T<:Base.SkipMissing = @inbounds map(_ishet, locus)

# public facing methods
"""
    ishet(locus::Genotype, allele::Signed)
    ishet(loci::GenoArray, allele::Signed)
Returns `true` if the `locus`/`loci` is/are heterozygous for the specified `allele`.
"""
function ishet(geno::T, allele::U) where T<:Genotype where U<:Integer
    ∈(allele, geno) & !ishom(geno) ? true : false
end

ishet(geno::T, allele::U) where T<:GenoArray where U<:Integer = map(i -> ishet(i, allele), geno)
ishet(geno::Missing, allele::U) where U<:Integer = false

_ishet(geno::T, allele::U) where T<:Genotype where U<:Integer = ∈(allele, geno) & !ishom(geno) ? true : false
_ishet(geno::T, allele::U) where T<:GenoArray where U<:Integer = map(i -> ishet(i, allele), geno)
_ishet(geno::Missing, allele::U) where U<:Integer = missing


text:
# ---
# title: 697. Degree of an Array
# id: problem697
# author: Indigo
# date: 2021-02-20
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/degree-of-an-array/description/>
# hidden: true
# ---
#
# Given a non-empty array of non-negative integers `nums`, the **degree** of
# this array is defined as the maximum frequency of any one of its elements.
#
# Your task is to find the smallest possible length of a (contiguous) subarray
# of `nums`, that has the same degree as `nums`.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,2,2,3,1]
#     Output: 2
#     Explanation:
#     The input array has a degree of 2 because both elements 1 and 2 appear twice.
#     Of the subarrays that have the same degree:
#     [1, 2, 2, 3, 1], [1, 2, 2, 3], [2, 2, 3, 1], [1, 2, 2], [2, 2, 3], [2, 2]
#     The shortest length is 2. So return 2.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1,2,2,3,1,4,2]
#     Output: 6
#     Explanation:
#     The degree is 3 because the element 2 is repeated 3 times.
#     So [2,2,3,1,4,2] is the shortest subarray, therefore returning 6.
#
#
#
#
# **Constraints:**
#
#   * `nums.length` will be between 1 and 50,000.
#   * `nums[i]` will be an integer between 0 and 49,999.
#
#
## @lc code=start
using LeetCode

function find_shortest_sub_array(nums::Vector{Int})
    mp = Dict{Int,Vector{Int}}()

    for (idx, num) in enumerate(nums)
        if num in keys(mp)
            mp[num][1] += 1
            mp[num][3] = idx
        else
            mp[num] = [1, idx, idx]
        end
    end
    maxn = minlen = 0
    for (cnt, l, r) in values(mp)
        if maxn < cnt
            maxn = cnt
            minlen = r - l + 1
        elseif maxn == cnt
            minlen = min(minlen, r - l + 1)
        end
    end
    minlen
end
## @lc code=end


text:
### A Pluto.jl notebook ###
# v0.14.0

using Markdown
using InteractiveUtils

# ╔═╡ 086c95a0-279f-11ec-3378-6924b1d7777e
begin

# Import packages used

using Distributed
using DelimitedFiles
using DataFrames
using Optim
using CSV
using LsqFit
using ForwardDiff
using LinearAlgebra
using StatFiles

end

# ╔═╡ aa3922f0-279e-11ec-33b6-6bb73cf7c16a
# Replication of DellaVigna and Pope, 2018, "What Motivates Effort? Evidence and Expert Forecasts"

# The code in this notebook estimates the parameters of column (2) (4) panel A and columns (3) (6) panel B of "Table 5: Estimates of behavioural parameters I: Mturkers actual effort and expert beliefs" and table 6 panel A: "Estimates of reference-dependent parameters: Mturker actual effort and expert beliefs" using Non-Linear-Leat-Squares

# ╔═╡ 08a0ec12-279f-11ec-12a3-f15a088873f4
# Authors:
# Massimiliano Pozzi (Bocconi University, pozzi.massimiliano@studbocconi.it)
# Salvatore Nunnari (Bocconi University, https://snunnari.github.io, salvatore.nunnari@unibocconi.it)

# ╔═╡ 6783e0c0-5904-11ec-2a41-5784eec17590
# This notebook works with both Neptune and Pluto; the "begin" and "end" tags in each cell are necessary for Pluto notebooks but do not play any role in Neptune

# ╔═╡ 903cd5d0-5904-11ec-0a95-75a9f8821a44
# This notebook was tested with the following packages versions:

# Pozzi: julia 1.6.1, Pluto 0.14.7, Neptune 0.14.0, DataFrames 1.2.2, CSV 0.9.5

# ╔═╡ 086c6e90-279f-11ec-2ce7-ef05700c68b3
begin

# 1. Data Cleaning and Data Preparation

# We import the relevant dataset containing data on the number of buttonpresses in the different treatments and for different piece rates wage that the participants received to complete the task. We then create a series of variables that are needed for estimation.

dt = DataFrame(load(raw"../input/mturk_clean_data_short.dta"))

# Create new variables needed for estimation:

# Create payoff per 100 button presses (p)

dt.payoff_per_100 = 0.0
dt[(dt[:treatment] .== "1.1"), :payoff_per_100] = 0.01
dt[(dt[:treatment] .== "1.2"), :payoff_per_100] = 0.1
dt[(dt[:treatment] .== "1.3"), :payoff_per_100] = 0.0
dt[(dt[:treatment] .== "2"),   :payoff_per_100] = 0.001
dt[(dt[:treatment] .== "1.4"), :payoff_per_100] = 0.04
dt[(dt[:treatment] .== "4.1"), :payoff_per_100] = 0.01
dt[(dt[:treatment] .== "4.2"), :payoff_per_100] = 0.01
dt[(dt[:treatment] .== "6.2"), :payoff_per_100] = 0.02
dt[(dt[:treatment] .== "6.1"), :payoff_per_100] = 1.0

# (ALPHA) create payoff per 100 to charity and dummy charity

dt.payoff_charity_per_100 = 0.0
dt[(dt[:treatment] .== "3.1"), :payoff_charity_per_100] = 0.01
dt[(dt[:treatment] .== "3.2"), :payoff_charity_per_100] = 0.1
dt.dummy_charity = 0
dt[(dt[:treatment] .== "3.1"), :dummy_charity] = 1
dt[(dt[:treatment] .== "3.2"), :dummy_charity] = 1

# (BETA/DELTA) create payoff per 100 delayed by 2 weeks and dummy delay

dt.delay_wks = 0
dt[(dt[:treatment] .== "4.1"), :delay_wks] = 2
dt[(dt[:treatment] .== "4.2"), :delay_wks] = 4
dt.delay_dummy = 0
dt[(dt[:treatment] .== "4.1"), :delay_dummy] = 1
dt[(dt[:treatment] .== "4.2"), :delay_dummy] = 1

# probability weights to back out curvature and dummy

dt.prob = 1.0
dt[(dt[:treatment] .== "6.2"), :prob] = 0.5
dt[(dt[:treatment] .== "6.1"), :prob] = 0.01
dt.weight_dummy = 0
dt[(dt[:treatment] .== "6.1"), :weight_dummy] = 1

# dummy for gift exchange

dt.gift_dummy = 0
dt[(dt[:treatment] .== "10"), :gift_dummy] = 1

# generating effort and log effort. authors round buttonpressed to nearest 100 value. If 0 set it to 25. julia rounds 50 to 0, while stata to 100. by adding a small value to the column 'buttonpresses' we avoid this mismatch

dt.buttonpresses = dt.buttonpresses .+ 0.1
dt.buttonpresses_nearest_100 = round.(dt.buttonpresses, digits = -2)
dt[(dt[:buttonpresses_nearest_100] .== 0.0), :buttonpresses_nearest_100] = 25
dt.logbuttonpresses_nearest_100 = log.(dt.buttonpresses_nearest_100)

end;

# ╔═╡ 082f8ca0-279f-11ec-07c3-838f8c01603e
# 2. Define the Model and the estimation technique

# For a more detailed explanation of the model we refer the reader to our python notebook or to section 2 in the paper.

# The model is one of costly effort, where an agent needs to choose the optimal effort (in this case the number of buttons pressed in a 10 minute session) to solve a simple tradeoff problem between disutility of effort and consumption utility derived from the consequent payment. On top of this simple problem the authors use 18 different treatments to examine the effects of standard monetary incentives, behavioral factors like social preferences and reference dependence and non-monetary incentives.

# ╔═╡ 531e546e-27a0-11ec-091f-ab9fb76190b2
begin

# 3. Point estimates and standard errors

# Define the benchmark sample by creating dummies equal to one if in treatment 1.1, 1.2, 1.3

dt.t11 = dt.treatment .== "1.1"
dt.t12 = dt.treatment .== "1.2"
dt.t13 = dt.treatment .== "1.3"
dt.dummy1 = dt.t11 .+ dt.t12 .+ dt.t13

# Set the initial values for the optimization procedure and scalers for k and s in the exp cost function case

gamma_init_exp, k_init_exp, s_init_exp =  0.015645717, 1.69443, 3.69198
st_values_exp = [gamma_init_exp, k_init_exp, s_init_exp]
k_scaler_exp, s_scaler_exp = 1e+16,1e+6

end;

# ╔═╡ 7b307e00-27a7-11ec-2070-9f6d80983e32
begin

# Since there are many different specifications (5 columns for the power cost function and 5 for the exponential cost function) we preferred to write each function separately instead of writing a single function with many if statements. Hopefully this will make each specification more clear.

# Estimate procedure for s, k, gamma in benchmark case with Exp cost function

# The non-linear least square procedure will choose the parameters that minimize the sum of squared residuals (y-f(x,θ))^2. Here we are computing f(x,θ) while y is the number of rounded buttonpresses.

# pay100 is the column we created containing the piece rate for different treatments
# g, k, s are the parameters to estimate (our θ vector). g stands for gamma.

function benchmarkExp(pay100, params)

	g, k, s = params

    check1= k / k_scaler_exp  # 'first' component to compute f(x,θ). We call                                           it check1 since it will enter a log, so we                                             need to be careful with its value being > 0
   	check2= s / s_scaler_exp .+ pay100

    f_x = (-1 /g * log.(check1) .+1 /g * log.(check2))

    return f_x
end

# Find the solution to the problem by non-linear least squares

sol = curve_fit(benchmarkExp,
				dt[(dt[:dummy1] .== 1.0), :payoff_per_100],
				dt[(dt[:dummy1] .== 1.0), :buttonpresses_nearest_100],
				st_values_exp,
		        autodiff=:forwarddiff)
be54 = sol.param     # Point estimates
se54 = stderror(sol) # Standard errors. They are the square root of the diagonal                              elements of the estimate for the var-cov matrix
end;

# ╔═╡ d153ea60-590c-11ec-09a4-d7d6bffb1d5a
begin

# Estimate procedure for s, k, gamma in benchmark case with power cost function

gamma_init_power, k_init_power, s_init_power =  19.8117987, 1.66306e-10, 7.74996
st_values_power = [gamma_init_power, k_init_power, s_init_power]
k_scaler_power, s_scaler_power = 1e+57, 1e+6

# Define f(x,θ) in the power case

function benchmarkPower(pay100, params)

	g, k, s = params

    check1= max(k, 1e-50) / k_scaler_power
   	check2= max(s, 1e-10) / s_scaler_power .+ pay100

    f_x = (-1 /g .* log.(check1) .+1 /g .* log.(check2))

    return f_x
end

# Find the solution to the problem by non-linear least squares

sol52 = curve_fit(benchmarkPower,
				dt[(dt[:dummy1] .== 1.0), :payoff_per_100],
				dt[(dt[:dummy1] .== 1.0), :logbuttonpresses_nearest_100],
				st_values_power,
		        autodiff=:forwarddiff)
bp52 = sol52.param
sp52 = stderror(sol52)

end;

# ╔═╡ 05e3eb00-7527-11ec-3380-cdcafa241aaa
# We find slightly different estimates compared to the ones obtained by the authors, but we see also see that different minimization algorithms implemented with the same programming language (julia) result in different estimates and/or values of the objective function, so it is not surprising that there are small discrepancies between our estimates and the authors' estimates (the authors use the Gauss-Newton minimization algorithm implemented in Stata). At the same time, the differences are small in absolute value (and limited to the NLSS estimation method, there are no discrepancies when using GMM) and the estimated values of k and s are always statistically indistinguishable from 0. More importantly, the economic implications of the estimated parameters and the qualitative conclusions on what motivates effort in the experiment are unaffected by the choice of programming language and minimization algorithm. We report the results we obtained with the curve_fit function since this also returns an estimate for the variance-covariance matrix for the parameters

# ╔═╡ 66dd6fb0-2b32-11ec-2a0c-5593db218da8
begin

# Allnoweight Exp. Create dummies for this specification

dt.t31 = dt.treatment .== "3.1"
dt.t32 = dt.treatment .== "3.2"
dt.t41 = dt.treatment .== "4.1"
dt.t42 = dt.treatment .== "4.2"
dt.t10 = dt.treatment .== "10"
dt.samplenw = dt.dummy1 .+ dt.t31 .+ dt.t32 .+ dt.t41 .+ dt.t42 .+ dt.t10

# Define the initial guesses for all parameters

alpha_init, a_init, beta_init, delta_init, gift_init = 0.003, 0.13, 1.16, 0.75, 5e-6
stvale_spec = [alpha_init, a_init, gift_init, beta_init, delta_init]
st_valuesnoweight_exp = vcat(st_values_exp,stvale_spec)

end;

# ╔═╡ 27a92130-2c64-11ec-2c22-c55b0094de30
begin

# Define the f(x,θ) to estimate all parameters but the probability weight in the exp case

# xdata is the vector containing the explanatory variables x:

# gd is gift dummy
# dd is delay dummy
# dw is delay weeks
# paychar is pay in charity treatment
# dc is dummy charity

# parameters:

# g, k, s are the same parameters from before
# alpha is the pure altruism coefficient
# a is the warm glow coefficient
# gift is the gift exchange coefficient Δs
# beta is the present bias paramater
# delta is the (weekly) discount factor

function noweightExp(pay100, params)

	# We know that passing variables from outside the function instead of using them         as function arguments is bad practice and it slows down the code, but curve_fit       will return an error if the first argument for the independent variable is a           vector of vectors. The problem lies in the fact that the package checks for           quality of the data using isinf() on the independent variables, but fails for a       vector of vectors since the code does not check element-wise. Writing isinf.()         in the package code and then check for any true should solve the problem.

    gd = dt[(dt[:samplenw] .== 1.0), :gift_dummy]
    dd = dt[(dt[:samplenw] .== 1.0), :delay_dummy]
    dw = dt[(dt[:samplenw] .== 1.0), :delay_wks]
    paychar = dt[(dt[:samplenw] .== 1.0), :payoff_charity_per_100]
    dc = dt[(dt[:samplenw] .== 1.0), :dummy_charity]
	g, k, s, alpha, a, gift, beta, delta = params

    check1 = k/k_scaler_exp
   	check2 = s/s_scaler_exp .+ gift.*0.4.*gd .+ (beta.^dd).*(delta.^dw).*pay100 .+ alpha.*paychar .+ a.*0.01.*dc

    f_x = (-1/g .* log.(check1) .+ 1/g .* log.(check2))

    return f_x
end

# Find the solution to the problem by non-linear least squares

sol56 = curve_fit(noweightExp,
				dt[(dt[:samplenw] .== 1.0), :payoff_per_100],
				dt[(dt[:samplenw] .== 1.0), :buttonpresses_nearest_100],
				st_valuesnoweight_exp,
		        autodiff=:forwarddiff)
be56 = sol56.param
se56 = stderror(sol56)

end;

# ╔═╡ 592baf2e-599f-11ec-3f47-1782a76c828d
begin

# Define the f(x,θ) to estimate all parameters but the probability weight in the power case

function noweightPower(pay100, params)

    gd = dt[(dt[:samplenw] .== 1.0), :gift_dummy]
    dd = dt[(dt[:samplenw] .== 1.0), :delay_dummy]
    dw = dt[(dt[:samplenw] .== 1.0), :delay_wks]
    paychar = dt[(dt[:samplenw] .== 1.0), :payoff_charity_per_100]
    dc = dt[(dt[:samplenw] .== 1.0), :dummy_charity]
	g, k, s, alpha, a, gift, beta, delta = params

    check1 = max(k,1e-50)/k_scaler_power
   	check2 = max(s,1e-10)/s_scaler_power .+ gift.*0.4.*gd .+ (beta.^dd).*(delta.^dw).*pay100 .+ alpha.*paychar .+ a.*0.01.*dc

    f_x = (-1/g .* log.(check1) .+ 1/g .* log.(check2))

    return f_x
end

# Find the solution to the problem by non-linear least squares

st_valuesnoweight_power = vcat(st_values_power, stvale_spec)
sol53 = curve_fit(noweightPower,
				  dt[(dt[:samplenw] .== 1.0), :payoff_per_100],
				  dt[(dt[:samplenw] .== 1.0), :logbuttonpresses_nearest_100],
				  st_valuesnoweight_power,
		          autodiff=:forwarddiff)
bp53 = sol53.param
sp53 = stderror(sol53)

end;

# ╔═╡ eff34620-59a0-11ec-37b4-f372da908fac
begin

# Create the arrays that will be the columns for table 5. Create the dataframe for table 5 and save it as a csv file

# Point estimates for power case do not coincide precisely as explained above. Standard errors do not coincide precisely because of the differences in the point estimates and because we leave here non-robust standard errors provided by curve_fit. To see an implementation of the formula for robust standard errors please refer to the python or julia notebooks for table_1 of augenblick-rabin or table_1 of bruhin-fehr-schunk. The formula is the same as in the cited notebooks without considering the clustering at the individual level.

params_name = ["Curvature γ of cost function", "Level k of cost of effort", "Intrinsic motivation s","Social preferences α","Warm glow coefficient a","Gift exchange Δs", "Present bias β","(Weekly) discount factor δ"]

be5 = [be54[1], be54[2]/1e+16, be54[3]/1e+6, be56[4], be56[5], be56[6], be56[7], be56[8]]
se5 = [se54[1], se54[2]/1e+16, se54[3]/1e+6, se56[4], se56[5], se56[6], se56[7], se56[8]]
bp5 = [bp52[1], bp52[2]/1e+57, bp52[3]/1e+6, bp53[4], bp53[5], bp53[6], bp53[7], bp53[8]]
sp5 = [sp52[1], sp52[2]/1e+57, sp52[3]/1e+6, sp53[4], sp53[5], sp53[6], sp53[7], sp53[8]]
t5 = DataFrame()
t5.parameters_name = params_name
t5.power_est=bp5
t5.power_se=sp5
t5.exp_est=be5
t5.exp_se=se5
CSV.write(raw"../output/table5NLS_julia.csv", t5)
t5

end

# ╔═╡ 8e8624e0-59a3-11ec-1934-6bf58626ec8a
# We now look at panel A of Table 6

# ╔═╡ 8ca6c9f0-2c66-11ec-1cb8-bf7d2915a69d
begin

# Create the sample used for Table 6 panel A

dt.t61 = dt.treatment .== "6.1"
dt.t62 = dt.treatment .== "6.2"
dt.samplepr = dt.dummy1 .+ dt.t61 .+ dt.t62

end;

# ╔═╡ 2f9046b0-2c6b-11ec-3d18-51475bf3a66c
begin

# Define f(x,θ) for the exponential cost function. Here we assume curvature of utility over piece rate = 1, (Column 4)
# wd is the weight_dummy
# prob is the prob_dummy
# g, k and s are the same parameters as before
# p_weight is the probability weighting coefficient under the assumption of linear value function in this case
# curv is the curvature of the value function. Here curv = 1

function probweight4Exp(pay100, params)

	wd = dt[(dt[:samplepr] .== 1.0), :weight_dummy]
    prob = dt[(dt[:samplepr] .== 1.0), :prob]
	g, k, s, p_weight = params

    check1= k / k_scaler_exp
   	check2= s / s_scaler_exp .+ p_weight.^wd.*prob.*pay100

    f_x = (-1 /g * log.(check1) .+1 /g * log.(check2))

    return f_x
end

# Find the solution to the problem by non-linear least squares

prob_weight_init = [0.2]
st_valuesprobweight_exp = vcat(st_values_exp, prob_weight_init)

sol64 = curve_fit(probweight4Exp,
				dt[(dt[:samplepr] .== 1.0), :payoff_per_100],
				dt[(dt[:samplepr] .== 1.0), :buttonpresses_nearest_100],
				st_valuesprobweight_exp,
		        autodiff=:forwarddiff)
be64 = sol64.param
se64 = stderror(sol64)

# we rescale the variables

be64 = [be64[1];be64[2]/1e+16;be64[3]/1e+6;be64[4]]
se64 = [se64[1];se64[2]/1e+16;se64[3]/1e+6;se64[4]]

end;

# ╔═╡ 3dd7ee00-59a5-11ec-021d-7fb70238f5a9
begin

# Define f(x,θ). Here we assume curvature of utility over piece rate = 0.88, Column (5)

function probweight5Exp(pay100, params)

	wd = dt[(dt[:samplepr] .== 1.0), :weight_dummy]
    prob = dt[(dt[:samplepr] .== 1.0), :prob]
	g, k, s, p_weight = params

    check1= k / k_scaler_exp
   	check2= s / s_scaler_exp .+ p_weight.^wd.*prob.*pay100.^0.88

    f_x = (-1 /g * log.(check1) .+1 /g * log.(check2))

    return f_x
end

# Find the solution to the problem by non-linear least squares

sol65 = curve_fit(probweight5Exp,
				dt[(dt[:samplepr] .== 1.0), :payoff_per_100],
				dt[(dt[:samplepr] .== 1.0), :buttonpresses_nearest_100],
				st_valuesprobweight_exp,
		        autodiff=:forwarddiff)
be65 = sol65.param
se65 = stderror(sol65)

be65 = [be65[1];be65[2]/1e+16;be65[3]/1e+6;be65[4]]
se65 = [se65[1];se65[2]/1e+16;se65[3]/1e+6;se65[4]]

# Define f(x,θ). Here we assume curvature of utility over piece rate = 0.88, Column (5)

function probweight6Exp(pay100, params)

	wd = dt[(dt[:samplepr] .== 1.0), :weight_dummy]
    prob = dt[(dt[:samplepr] .== 1.0), :prob]
	g, k, s, p_weight, curv = params

    check1= k / k_scaler_exp
   	check2= s / s_scaler_exp .+ p_weight.^wd.*prob.*pay100.^curv

    f_x = (-1 /g * log.(check1) .+1 /g * log.(check2))

    return f_x
end

# Find the solution to the problem by non-linear least squares

curv_init = [0.5]
st_valuesprobweight6_exp = vcat(st_valuesprobweight_exp, curv_init)

sol66 = curve_fit(probweight6Exp,
				dt[(dt[:samplepr] .== 1.0), :payoff_per_100],
				dt[(dt[:samplepr] .== 1.0), :buttonpresses_nearest_100],
				st_valuesprobweight6_exp,
		        autodiff=:forwarddiff)
be66 = sol66.param
se66 = stderror(sol66)

be66 = [be66[1];be66[2]/1e+16;be66[3]/1e+6;be66[4:5]]
se66 = [se66[1];se66[2]/1e+16;se66[3]/1e+6;se66[4:5]]

end;

# ╔═╡ b5b94c80-2c6d-11ec-0ab0-71387355858d
begin

# We do the same for the power cost function specification

# For column 4

function probweight4Power(pay100, params)

	wd = dt[(dt[:samplepr] .== 1.0), :weight_dummy]
    prob = dt[(dt[:samplepr] .== 1.0), :prob]
	g, k, s, p_weight = params

    check1= max(k,1e-50) / k_scaler_power
   	check2= max(s,1e-10) / s_scaler_power .+ p_weight.^wd.*prob.*pay100

    f_x = (-1 /g * log.(check1) .+1 /g * log.(check2))

    return f_x
end

# Find the solution to the problem by non-linear least squares

st_valuesprobweight_power = vcat(st_values_power, prob_weight_init)

sol61 = curve_fit(probweight4Power,
				dt[(dt[:samplepr] .== 1.0), :payoff_per_100],
				dt[(dt[:samplepr] .== 1.0), :logbuttonpresses_nearest_100],
				st_valuesprobweight_power,
		        autodiff=:forwarddiff)
bp61 = sol61.param
sp61 = stderror(sol61)

bp61 = [bp61[1];bp61[2]/1e+57;bp61[3]/1e+6;bp61[4]]
sp61 = [sp61[1];sp61[2]/1e+57;sp61[3]/1e+6;sp61[4]]

# For column 5

function probweight5Power(pay100, params)

	wd = dt[(dt[:samplepr] .== 1.0), :weight_dummy]
    prob = dt[(dt[:samplepr] .== 1.0), :prob]
	g, k, s, p_weight = params

    check1= max(k,1e-50) / k_scaler_power
   	check2= max(s,1e-10) / s_scaler_power .+ p_weight.^wd.*prob.*pay100.^0.88

    f_x = (-1 /g * log.(check1) .+1 /g * log.(check2))

    return f_x
end

# Find the solution to the problem by non-linear least squares

sol62 = curve_fit(probweight5Power,
				dt[(dt[:samplepr] .== 1.0), :payoff_per_100],
				dt[(dt[:samplepr] .== 1.0), :logbuttonpresses_nearest_100],
				st_valuesprobweight_power,
		        autodiff=:forwarddiff)
bp62 = sol62.param
sp62 = stderror(sol62)

bp62 = [bp62[1];bp62[2]/1e+57;bp62[3]/1e+6;bp62[4]]
sp62 = [sp62[1];sp62[2]/1e+57;sp62[3]/1e+6;sp62[4]]

# For column 6

function probweight6Power(pay100, params)

	wd = dt[(dt[:samplepr] .== 1.0), :weight_dummy]
    prob = dt[(dt[:samplepr] .== 1.0), :prob]
	g, k, s, p_weight, curv = params

    check1= max(k,1e-50) / k_scaler_power
   	check2= max(s,1e-10) / s_scaler_power .+ p_weight.^wd.*prob.*pay100.^curv

    f_x = (-1 /g * log.(check1) .+1 /g * log.(check2))

    return f_x
end

# Find the solution to the problem by non-linear least squares

st_valuesprobweight6_power = vcat(st_valuesprobweight_power, curv_init)

sol63 = curve_fit(probweight6Power,
				dt[(dt[:samplepr] .== 1.0), :payoff_per_100],
				dt[(dt[:samplepr] .== 1.0), :logbuttonpresses_nearest_100],
				st_valuesprobweight6_power,
		        autodiff=:forwarddiff)
bp63 = sol63.param
sp63 = stderror(sol63)

bp63 = [bp63[1];bp63[2]/1e+57;bp63[3]/1e+6;bp63[4:5]]
sp63 = [sp63[1];sp63[2]/1e+57;sp63[3]/1e+6;sp63[4:5]]

# columns must be of the same length
abp61 = append!(bp61,1)
asp61 = append!(sp61,0)
abp62 = append!(bp62,0.88)
asp62 = append!(sp62,0)
abe64 = append!(be64,1)
ase64 = append!(se64,0)
abe65 = append!(be65,0.88)
ase65 = append!(se65,0)

end;

# ╔═╡ d7f36f90-2c6f-11ec-155a-db4f5d543604
begin

# Create the dataframe relative to table 6 and save it as a csv file

# There are some differences in the standard errors since we leave here non robust standard errors provided by curve_fit. Point estimates for the power cost function are again a little different from the authors', while they are the same for the exponential cost function

t6 = DataFrame()

pnames = ["Curvature γ of cost function", "Level k of cost of effort", "Intrinsic motivation s", "Probability weighting π (1%) (in %)", "Curvature of utility over piece rate"]

t6.parameters=pnames
t6.p_est1 = abp61
t6.p_se1  = asp61
t6.p_est2 = abp62
t6.p_se2  = asp62
t6.p_est3 = bp63
t6.p_se3  = sp63
t6.e_est4 = abe64
t6.e_se4  = ase64
t6.e_est5 = abe65
t6.e_se5  = ase65
t6.e_est6 = be66
t6.e_se6  = se66

CSV.write(raw"../output/table6_julia.csv", t6)
t6

end

# ╔═╡ Cell order:
# ╠═aa3922f0-279e-11ec-33b6-6bb73cf7c16a
# ╠═08a0ec12-279f-11ec-12a3-f15a088873f4
# ╠═6783e0c0-5904-11ec-2a41-5784eec17590
# ╠═903cd5d0-5904-11ec-0a95-75a9f8821a44
# ╠═086c95a0-279f-11ec-3378-6924b1d7777e
# ╠═086c6e90-279f-11ec-2ce7-ef05700c68b3
# ╠═082f8ca0-279f-11ec-07c3-838f8c01603e
# ╠═531e546e-27a0-11ec-091f-ab9fb76190b2
# ╠═7b307e00-27a7-11ec-2070-9f6d80983e32
# ╠═d153ea60-590c-11ec-09a4-d7d6bffb1d5a
# ╠═05e3eb00-7527-11ec-3380-cdcafa241aaa
# ╠═66dd6fb0-2b32-11ec-2a0c-5593db218da8
# ╠═27a92130-2c64-11ec-2c22-c55b0094de30
# ╠═592baf2e-599f-11ec-3f47-1782a76c828d
# ╠═eff34620-59a0-11ec-37b4-f372da908fac
# ╠═8e8624e0-59a3-11ec-1934-6bf58626ec8a
# ╠═8ca6c9f0-2c66-11ec-1cb8-bf7d2915a69d
# ╠═2f9046b0-2c6b-11ec-3d18-51475bf3a66c
# ╠═3dd7ee00-59a5-11ec-021d-7fb70238f5a9
# ╠═b5b94c80-2c6d-11ec-0ab0-71387355858d
# ╠═d7f36f90-2c6f-11ec-155a-db4f5d543604


text:
@testset "Get All Files Function Tests" begin

  @test isdefined(Julz, :get_all_files) == true

  cd("dummy") do

    @test length(Julz.get_all_files("src")) == 1
    @test length(Julz.get_all_files("test")) == 0

    Julz.generate("type", "foo")

    @test length(Julz.get_all_files("src")) == 2
    @test length(Julz.get_all_files("test")) == 1

    Julz.destroy("type", "foo")

    @test length(Julz.get_all_files("src")) == 1
    @test length(Julz.get_all_files("test")) == 0

  end

  main_file = Julz.get_all_files("dummy/src")[1]

  @test main_file == "dummy/src/JulzDummy.jl"

end


text:
py"""
import pickle
builtin_type_samples = {
  'str': 'Julia!',
  'int': 42,
  'bool': {True: False, False: True},
  'float': 3.1415926,
  'bytes': b'1234',
  'tuple': (1, 2.0, '3', b'4'),
  'set': {1, 2, 3, 12, 21},
  'bigint': 1234567890987654321012345678909876543210,
  'list': ['February', 14, 2012]
}

def pyload(f):
  with open(f, "rb") as io:
    return pickle.load(io)

def pyloads(s):
  return pickle.loads(s)

def pystore(f, x, p):
  with open(f, "wb+") as io:
    return pickle.dump(x, io, protocol=p)

def pystores(x, p):
  return pickle.dumps(x, protocol=p)

def check_bts(f):
  global builtin_type_samples
  return builtin_type_samples == pyload(f)

def scheck_bts(f):
  global builtin_type_samples
  with open(f, "rb") as io:
    s = io.read()
  return builtin_type_samples == pyloads(s)

def to_nparray(x):
  return x.toarray()

# pre-built pickle files for testing `load`
# for i in range(5):
#   pickle.dump(x, open(f"./test_pkl/builtin_type_p{i}.pkl", "wb+"), protocol=i)
"""
pyload = py"pyload"
pyloads = py"pyloads"
pystore = py"pystore"
pystores = py"pystores"
check_bts = py"check_bts"
scheck_bts = py"scheck_bts"
to_nparray = py"to_nparray"


text:
using SecondQuantization
using Documenter

DocMeta.setdocmeta!(SecondQuantization, :DocTestSetup, :(using SecondQuantization); recursive=true)

makedocs(;
    modules=[SecondQuantization],
    authors="Alberto Mercurio",
    repo="https://github.com/albertomercurio/SecondQuantization.jl/blob/{commit}{path}#{line}",
    sitename="SecondQuantization.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://albertomercurio.github.io/SecondQuantization.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/albertomercurio/SecondQuantization.jl",
    devbranch="main",
)


text:
const Ctypes = Dict{String, DataType}(
    "void"      => Cvoid,
    "int"       => Cint,
    "slong"     => Clong,
    "ulong"     => Culong,
    "double"    => Cdouble,
    "arf_t"     => Arf,
    "arb_t"     => Arb,
    "acb_t"     => Acb,
    "mag_t"     => Mag,
    "arf_rnd_t" => arb_rnd
)

struct Carg{ArgT}
    name::String
    isconst::Bool
end

function Carg(str)
    m = match(r"(?<const>const)?\s*(?<type>\w+)\s+(?<name>\w+)", str)
    isnothing(m) && throw(ArgumentError("string doesn't match c-argument pattern"))

    return Carg{Ctypes[m[:type]]}(m[:name], !isnothing(m[:const]))
end

name(ca::Carg) = ca.name
isconst(ca::Carg) = ca.isconst
jltype(::Carg{ArgT}) where ArgT = ArgT
ctype(ca::Carg) = jltype(ca)
ctype(::Carg{ArgT}) where ArgT <: Union{Arf, Arb, Acb, Mag}  = Ref{ArgT}

struct Arbfunction{ReturnT}
    fname::String
    args::Vector{Carg}
end

function Arbfunction(str)
    m = match(r"(?<returntype>\w+)\s+(?<arbfunction>[\w_]+)\((?<args>.*)\)",
        str)
    isnothing(m) && throw(ArgumentError("string doesn't match arblib function signature pattern"))

    args = Carg.(strip.(split(m[:args], ",")))

    return Arbfunction{Ctypes[m[:returntype]]}(m[:arbfunction], args)
end

function jlfname(arbfname,
        prefixes=("arf", "arb", "acb", "mag"),
        suffixes=("si", "ui", "d", "arf", "arb");
        inplace=false)
    strs = split(arbfname, "_")
    k = findfirst(s->s ∉ prefixes, strs)
    l = findfirst(s->s ∉ suffixes, reverse(strs))
    fname = join(strs[k:end-l+1], "_")
    return inplace ? Symbol(fname, "!") : Symbol(fname)
end

arbfname(af::Arbfunction) = af.fname
returntype(af::Arbfunction{ReturnT}) where ReturnT = ReturnT
arguments(af::Arbfunction) = af.args

function inplace(af::Arbfunction)
    firstarg = first(arguments(af))
    return !isconst(firstarg) && ctype(firstarg) <: Ref
end

function jlfname(af::Arbfunction,
        prefixes=("arf", "arb", "acb", "mag"),
        suffixes=("si", "ui", "d", "arf", "arb");
        inplace=inplace(af))
    return jlfname(arbfname(af), prefixes, suffixes, inplace=inplace)
end

function jlcode(af::Arbfunction, jl_fname=jlfname(af))
    returnT = returntype(af)
    args = arguments(af)
    c_types = ctype.(args)

    arg_names = Symbol.(name.(args))
    jl_types = jltype.(args)

    kwargs = Expr[]

    k = findfirst(==(:prec), arg_names)
    if !isnothing(k)
        @assert jl_types[k] == Int64
        p = :prec
        a = first(args)
        default = if jltype(a) ∈ (Arf, Arb, Acb)
            :(precision($(Symbol(name(a)))))
        else
            :(DEFAULT_PRECISION[])
        end
        push!(kwargs, Expr(:kw, :($p::Integer), default))
        deleteat!(arg_names, k)
        deleteat!(jl_types, k)
    end

    k = findfirst(==(:rnd), arg_names)
    if !isnothing(k)
        @assert jl_types[k] == arb_rnd
        r = :rnd
        push!(kwargs, Expr(:kw, :($r::Union{arb_rnd, RoundingMode}), :(RoundNearest)))
        deleteat!(arg_names, k)
        deleteat!(jl_types, k)
    end

    jl_args = [:($a::$T) for (a, T) in zip(arg_names, jl_types)]

    res = first(arg_names)

    return :(
        function $jl_fname($(jl_args...); $(kwargs...))
            ccall(Arblib.@libarb($(arbfname(af))),
            $returnT,
            $(Expr(:tuple, c_types...)),
            $(Symbol.(name.(args))...))
            return $res
        end
        )
end

macro arbcall_str(str)
    af = Arbfunction(str)
    return jlcode(af)
end


text:
################################################################################
#
#  Compute unit group from class group context
#
################################################################################

function _unit_group(O::NfOrd, c::ClassGrpCtx)
  u = UnitGrpCtx{FacElem{nf_elem, AnticNumberField}}(O)
  _unit_group_find_units(u, c)
  return u
end

function find_candidates(x::ClassGrpCtx, u::UnitGrpCtx, add::Int = 0)
  time_kernel = 0.0
  add_units = Int[]
  rel = SMat{fmpz}()
  K = nf(x)
  r1, r2 = signature(K)
  nrel = max(10, r1+r2-1)
  if nrel > nrows(x.M.rel_gens)
    nrel = nrows(x.M.rel_gens)
  end
  while length(add_units) < nrel
    xj = rand(1:nrows(x.M.rel_gens))
    if length(u.relations_used) != nrows(x.M.rel_gens)
      if xj in add_units || xj in u.relations_used
        continue
      end
      push!(u.relations_used, xj)
    else
      if xj in add_units
        continue
      end
    end
    push!(add_units, xj)
    push!(rel, x.M.rel_gens[xj])
  end
  for i = 1:add
    xj = rand(1:nrows(x.M.rel_gens))
    if length(u.relations_used) != nrows(x.M.rel_gens)
      if xj in add_units || xj in u.relations_used
        continue
      end
      push!(u.relations_used, xj)
    else
      if xj in add_units
        continue
      end
    end
    push!(add_units, xj)
    push!(rel, x.M.rel_gens[xj])
  end
  time_kernel += @elapsed k, d = solve_dixon_sf(x.M.bas_gens, rel)
  @vprint :UnitGroup 1 "Saturating the kernel\n"
  @vtime_add_elapsed :UnitGroup 1 x :saturate_time s = saturate(hcat(k, (-d)*identity_matrix(SMat, FlintZZ, k.r)))
  @vprint :UnitGroup 1 "Done\n"
  s1 = matrix(s)
  lll!(s1)
  @vprint :UnitGroup 1 "Kernel time: $time_kernel\n"
  @vtime_add :UnitGroup 1 x :unit_hnf_time time_kernel
  return k, add_units, s1
end

function _unit_group_find_units(u::UnitGrpCtx, x::ClassGrpCtx; add_orbit::Bool = true, expected_reg::arb = ArbField(32, cached = false)(-1), add::Int = 0)
  add_orbit = false
  @vprint :UnitGroup 1 "Processing ClassGrpCtx to find units ... (using orbits: $add_orbit)\n"
  @vprint :UnitGroup 1 "Relation module $(x.M)\n"

  O = order(u)

  K = nf(order(x.FB.ideals[1]))
  r = unit_group_rank(O)

  if r == 0
    Ar = ArbField(u.indep_prec, cached = false)
    u.tentative_regulator = Ar(1)
    u.regulator = Ar(1)
    u.regulator_precision = u.indep_prec
    u.full_rank = true
    return 1, 0
  end

  # I am not allowed to do this before the other block
  if nrows(x.M.rel_gens) == 0
    @vprint :UnitGroup 1 "No additional relations. Going back ...\n"
    return 0, 0
  end

  time_indep = 0.0
  time_add_dep_unit = 0.0
  time_torsion = 0.0
  not_larger = 0

  @vprint :UnitGroup 1 "Enlarging unit group by adding kernel elements ...\n"

  starting_full_rank = has_full_rank(u)
  if starting_full_rank
    starting_idx = _validate_class_unit_group(x, u)[1]
  else
    starting_idx = 0
  end

  not_larger_bound = min(20, nrows(x.M.rel_gens), r)

  first = true
  if has_full_rank(u)
    first = false
    add += 2 #+ div(nrows(x.M.rel_gens), r)
  end
  if !isdefined(u, :relations_used)
    u.relations_used = Vector{Int}()
  end
  finished = false

  if add_orbit
    aut = automorphisms(u)
    gens_aut = small_generating_set(aut)
    indices_aut = Int[]
    for s = 1:length(gens_aut)
      for j = 1:length(aut)
        if aut[j] == gens_aut[s]
          push!(indices_aut, j)
          break
        end
      end
    end
  end
  new_add = 0
  while not_larger < not_larger_bound
    add += new_add
    new_add = 2
    k, add_units, s1 = find_candidates(x, u, add)

    ge = vcat(x.R_gen[1:k.c], x.R_rel[add_units])
    elements = Vector{FacElem{nf_elem, AnticNumberField}}(undef, nrows(s1))
    for i = 1:nrows(s1)
      elements[i] = FacElem(ge, fmpz[s1[i, j] for j = 1:ncols(s1)])
    end
    p = 32
    if has_full_rank(u)
      elements = reduce_mod_units(elements, u)
    end
    #I order the elements by the maximum of the conjugate log.
    m_conjs = Vector{fmpz}(undef, length(elements))
    for i = 1:length(m_conjs)
      m_conjs[i] = maximum(fmpz[abs_upper_bound(fmpz, x) for x in conjugates_arb_log(elements[i], p)])
    end
    p_elements = sortperm(m_conjs)
    elements = elements[p_elements]
    done = falses(length(elements))
    for i = 1:length(elements)
      @vprint :UnitGroup 1 "Element $(i) / $(length(elements))\n"
      y = elements[i]
      @hassert :UnitGroup 2 _isunit(y)
      @vprint :UnitGroup 1 "(It really is a unit.)\n"
      @hassert :UnitGroup 9000 denominator(minpoly(evaluate(y))) == 1

      if has_full_rank(u)
        @vprint :UnitGroup 2 "have full rank, can reduce unit first...\n"
        y = reduce_mod_units(FacElem{nf_elem, AnticNumberField}[y], u)[1]
      else
        @vprint :UnitGroup 2 "no full rank, cannot reduce unit first...\n"
      end

      @vprint :UnitGroup 1 "Exponents are of bit size $(isempty(y.fac) ? 0 : maximum([ nbits(o) for o in values(y.fac)]))\n"
      @vprint :UnitGroup 1 "Test if kernel element yields torsion unit ... \n"

      @v_do :UnitGroup 2 pushindent()
      time_torsion += @elapsed is_tors, p1 = is_torsion_unit(y, false, u.tors_prec)
      @v_do :UnitGroup 2 popindent()

      u.tors_prec = max(p1, u.tors_prec)
      p = max(p, u.tors_prec)
      if is_tors
        @vprint :UnitGroup 1 "Element is torsion unit\n"
        not_larger += 1
        # We do break out of the for loop if not_larger > not_larger_bound,
        # because otherwise we do not check all kernel elements
        done[i] = true
        continue
      end

      @vprint :UnitGroup 1 "Element is non-torsion unit\n"

      @v_do :UnitGroup 2 pushindent()

      m = add_unit!(u, y)
      if m
        new_add = 0
        done[i] = true
        not_larger = 0
        if has_full_rank(u)
          @vprint :UnitGroup 1 "improved reg, reg is $(tentative_regulator(u))\n"
          if first
            idx, expected_reg = _validate_class_unit_group(x, u)
            first = false
          end
          if expected_reg > divexact(abs(tentative_regulator(u)), 2)
            done = trues(length(elements))
            not_larger = not_larger_bound + 1
            finished = true
            break
          end
          if add_orbit
            for j  in indices_aut
              phiu = apply_automorphism(u, j, y)
              if add_unit!(u, phiu)
                @vprint :UnitGroup 1 "New regulator: $(tentative_regulator(u))\n"
                not_larger = 0
              end
            end
          end
        else
          @vprint :UnitGroup 1 "Increased rank by 1 (now $(rank(u)))\n"
        end
      else
        if has_full_rank(u)
          done[i] = true
        end
        not_larger = not_larger + 1
        # We do break out of the for loop if not_larger > not_larger_bound,
        # because otherwise we do not check all kernel elements
      end

      @v_do :UnitGroup 2 popindent()
    end

    if has_full_rank(u)
      add_done = false
      for i = 1:length(elements)
        if !done[i]
          time_torsion += @elapsed is_tors, p1 = is_torsion_unit(elements[i], false, u.tors_prec)
          @v_do :UnitGroup 2 popindent()

          u.tors_prec = max(p1, u.tors_prec)
          p = max(p, u.tors_prec)
          if is_tors
            @vprint :UnitGroup 1 "Element is torsion unit\n"
            done[i] = true
            continue
          end
          add_done = add_unit!(u, elements[i]) || add_done
          if expected_reg > divexact(abs(tentative_regulator(u)), 2)
            done = trues(length(elements))
            not_larger = not_larger_bound + 1
            finished = true
            break
          end
        end
      end
      u.units = reduce(u.units, u.tors_prec)
    end
    if finished
      break
    end
    if not_larger > not_larger_bound && length(u.relations_used) != nrows(x.M.rel_gens)
      not_larger = 0
    end
  end
  @vprint :UnitGroup 1 "Finished processing\n"
  if has_full_rank(u)
    u.tentative_regulator = regulator(u.units, 64)
    @vprint :UnitGroup 1 "Regulator of current unit group is $(u.tentative_regulator)\n"
  else
    @vprint :UnitGroup 1 "current rank is $(length(u.units)), need $r\n"
  end
  @vprint :UnitGroup 1 "-"^80 * "\n"
  @vprint :UnitGroup 1 "Independent unit time: $time_indep\n"
  @vprint :UnitGroup 1 "Adding dependent unit time: $time_add_dep_unit\n"
  @vprint :UnitGroup 1 "Torsion test time: $time_torsion\n"

  if starting_full_rank
    return 1, div(starting_idx, _validate_class_unit_group(x, u)[1])
  elseif has_full_rank(u)
    return 1, 0
  else
    return 0, 0
  end
end

function compute_galois_closure!(U::UnitGrpCtx, c::ClassGrpCtx)
  @vprint :UnitGroup 1 "Computing Galois closure \n"
  aut = automorphisms(U)
  gens_aut = small_generating_set(aut)
  indices_aut = Int[]
  for s = 1:length(gens_aut)
    for j = 1:length(aut)
      if aut[j] == gens_aut[s]
        push!(indices_aut, j)
        break
      end
    end
  end
  found_new = true
  while found_new
    found_new = false
    for i = 1:length(indices_aut)
      for j = 1:length(U.units)
        @vprint :UnitGroup 1 "Applying auto \n"
        uphi = apply_automorphism(U, indices_aut[i], U.units[j])
        @vprint :UnitGroup 1 "Adding unit \n"
        fl = add_unit!(U, uphi)
        if fl
          @vprint :UnitGroup 1 "Found new unit! \n"
          found_new = true
          idx = _validate_class_unit_group(c, U)[1]
          if isone(idx)
            return nothing
          end
        end
      end
    end
  end
  return nothing

end


text:
using DbgMacro

@test (@qn foo(1, x, :y, $z)) == :(foo(1, x, :y, $(Expr(:$, :z))))


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct OpDeprecation <: ProtoType
    version::Int32
    explanation::AbstractString
    OpDeprecation(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct OpDeprecation

mutable struct OpDef_ArgDef <: ProtoType
    name::AbstractString
    description::AbstractString
    _type::Int32
    type_attr::AbstractString
    number_attr::AbstractString
    type_list_attr::AbstractString
    is_ref::Bool
    OpDef_ArgDef(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct OpDef_ArgDef
const __fnum_OpDef_ArgDef = Int[1,2,3,4,5,6,16]
meta(t::Type{OpDef_ArgDef}) = meta(t, ProtoBuf.DEF_REQ, __fnum_OpDef_ArgDef, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct OpDef_AttrDef <: ProtoType
    name::AbstractString
    _type::AbstractString
    default_value::AttrValue
    description::AbstractString
    has_minimum::Bool
    minimum::Int64
    allowed_values::AttrValue
    OpDef_AttrDef(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct OpDef_AttrDef

mutable struct OpDef <: ProtoType
    name::AbstractString
    input_arg::Base.Vector{OpDef_ArgDef}
    output_arg::Base.Vector{OpDef_ArgDef}
    control_output::Base.Vector{AbstractString}
    attr::Base.Vector{OpDef_AttrDef}
    deprecation::OpDeprecation
    summary::AbstractString
    description::AbstractString
    is_commutative::Bool
    is_aggregate::Bool
    is_stateful::Bool
    allows_uninitialized_input::Bool
    OpDef(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct OpDef
const __fnum_OpDef = Int[1,2,3,20,4,8,5,6,18,16,17,19]
meta(t::Type{OpDef}) = meta(t, ProtoBuf.DEF_REQ, __fnum_OpDef, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct OpList <: ProtoType
    op::Base.Vector{OpDef}
    OpList(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct OpList

export OpDef_ArgDef, OpDef_AttrDef, OpDef, OpDeprecation, OpList


text:
module FakePkg

using Test
using Comonicon

"""
fake add

# Args

- `package`: package to add.
"""
@cast add(package) = @test package == "ABC"

"""
fake rm

# Args

- `package`: package to add.
"""
@cast rm(package) = @test package == "ABC"

"""
fake activate

# Args

- `env`: environment to activate.

# Flags

- `-s, --shared`: fake flag share.
"""
@cast function activate(env; shared::Bool = false)
    @test env == "fake"
    @test shared == true
end

include("registry.jl")

@main

end


text:
using .CuArrays

# Discard `effort` argument for CuArrays
plan_flows_fft(a::CuArray, effort) = plan_fft(a)
plan_flows_rfft(a::CuArray, effort) = plan_rfft(a)

OneDGrid(dev::GPU, args...; kwargs...) = OneDGrid(args...; ArrayType=CuArray, kwargs...)
TwoDGrid(dev::GPU, args...; kwargs...) = TwoDGrid(args...; ArrayType=CuArray, kwargs...)

function Base.zeros(::GPU, T, dims)
    a = CuArray{T}(undef, dims...)
    a .= 0
    return a
end

ArrayType(::GPU, T, dim) = CuArray{T, dim}
ArrayType(::GPU) = CuArray

supersize(a::CuArray) = size(a)

getetdcoeffs(dt, L::CuArray; kwargs...) =
    (CuArray(ζ) for ζ in getetdcoeffs(dt, Array(L); kwargs...))

makefilter(K::CuArray; kwargs...) = CuArray(makefilter(Array(K); kwargs...))

function makefilter(g::AbstractGrid{Tg, <:CuArray}, T, sz; kwargs...) where Tg
    CuArray(ones(T, sz)) .* makefilter(g; realvars=sz[1]==g.nkr, kwargs...)
end


text:
module TestManual

using Test
using Trixi

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

# Run various manual (= non-parameter-file-triggered tests)
@testset "Manual tests" begin
  @testset "Tree" begin
    @testset "constructors" begin
      @test_nowarn Trixi.Tree(Val(1), 10, 0.0, 1.0)
    end

    @testset "helper functions" begin
      t = Trixi.Tree(Val(1), 10, 0.0, 1.0)
      @test_nowarn show(t)
      @test Trixi.ndims(t) == 1
      @test Trixi.ndims(Trixi.Tree{1}) == 1
      @test Trixi.has_any_neighbor(t, 1, 1) == true
      @test Trixi.isperiodic(t, 1) == true
      @test Trixi.n_children_per_cell(t) == 2
      @test Trixi.n_children_per_cell(2) == 4
      @test Trixi.n_directions(t) == 2
    end

    @testset "refine!/coarsen!" begin
      t = Trixi.Tree(Val(1), 10, 0.0, 1.0)
      @test Trixi.refine!(t) == [1]
      @test Trixi.coarsen!(t) == [1]
      @test Trixi.refine!(t) == [1]
      @test Trixi.coarsen!(t, 1) == [1]
      @test Trixi.coarsen!(t) == Int[] # Coarsen twice to check degenerate case of single-cell tree
      @test Trixi.refine!(t) == [1]
      @test Trixi.refine!(t) == [2,3]
      @test Trixi.coarsen_box!(t, [-0.5], [0.0]) == [2]
      @test Trixi.coarsen_box!(t, 0.0, 0.5) == [3]
      @test isnothing(Trixi.reset_data_structures!(t))
    end
  end

  @testset "interpolation" begin
    @testset "nodes and weights" begin
      @test Trixi.gauss_nodes_weights(1) == ([0.0], [2.0])
    end

    @testset "multiply_dimensionwise" begin
      nodes_in  = [0.0, 0.5, 1.0]
      nodes_out = [0.0, 1/3, 2/3, 1.0]
      matrix = Trixi.polynomial_interpolation_matrix(nodes_in, nodes_out)
      data_in = [3.0 4.5 6.0]
      @test isapprox(Trixi.multiply_dimensionwise(matrix, data_in), [3.0 4.0 5.0 6.0])

      n_vars   = 3
      size_in  = 4
      size_out = 5
      matrix   = randn(size_out, size_in)
      # 1D
      data_in  = randn(n_vars, size_in)
      data_out = Trixi.multiply_dimensionwise_naive(matrix, data_in)
      @test isapprox(data_out, Trixi.multiply_dimensionwise(matrix, data_in))
      # 2D
      data_in  = randn(n_vars, size_in, size_in)
      data_out = Trixi.multiply_dimensionwise_naive(matrix, data_in)
      @test isapprox(data_out, Trixi.multiply_dimensionwise(matrix, data_in))
      # 3D
      data_in  = randn(n_vars, size_in, size_in, size_in)
      data_out = Trixi.multiply_dimensionwise_naive(matrix, data_in)
      @test isapprox(data_out, Trixi.multiply_dimensionwise(matrix, data_in))
    end
  end

  @testset "containers" begin
    # Set up mock container
    mutable struct MyContainer <: Trixi.AbstractContainer
      data::Vector{Int}
      capacity::Int
      length::Int
      dummy::Int
    end
    function MyContainer(data, capacity)
      c = MyContainer(Vector{Int}(undef, capacity+1), capacity, length(data), capacity+1)
      c.data[1:length(data)] .= data
      return c
    end
    MyContainer(data::AbstractArray) = MyContainer(data, length(data))
    Trixi.invalidate!(c::MyContainer, first, last) = (c.data[first:last] .= 0; c)
    function Trixi.raw_copy!(target::MyContainer, source::MyContainer, first, last, destination)
      Trixi.copy_data!(target.data, source.data, first, last, destination)
      return target
    end
    Trixi.move_connectivity!(c::MyContainer, first, last, destination) = c
    Trixi.delete_connectivity!(c::MyContainer, first, last) = c
    Trixi.reset_data_structures!(c::MyContainer) = (c.data = Vector{Int}(undef, c.capacity+1); c)
    function Base.:(==)(c1::MyContainer, c2::MyContainer)
      return (c1.capacity == c2.capacity &&
              c1.length == c2.length &&
              c1.dummy == c2.dummy &&
              c1.data[1:c1.length] == c2.data[1:c2.length])
    end

    @testset "size" begin
      c = MyContainer([1, 2, 3])
      @test size(c) == (3,)
    end

    @testset "resize!" begin
      c = MyContainer([1, 2, 3])
      @test length(resize!(c, 2)) == 2
    end

    @testset "copy!" begin
      c1 = MyContainer([1, 2, 3])
      c2 = MyContainer([4, 5])
      @test Trixi.copy!(c1, c2, 2, 1, 2) == MyContainer([1, 2, 3]) # no-op

      c1 = MyContainer([1, 2, 3])
      c2 = MyContainer([4, 5])
      @test Trixi.copy!(c1, c2, 1, 2, 2) == MyContainer([1, 4, 5])

      c1 = MyContainer([1, 2, 3])
      @test Trixi.copy!(c1, c2, 1, 2) == MyContainer([1, 4, 3])

      c1 = MyContainer([1, 2, 3])
      @test Trixi.copy!(c1, 2, 3, 1) == MyContainer([2, 3, 3])

      c1 = MyContainer([1, 2, 3])
      @test Trixi.copy!(c1, 1, 3) == MyContainer([1, 2, 1])
    end

    @testset "move!" begin
      c = MyContainer([1, 2, 3])
      @test Trixi.move!(c, 1, 1) == MyContainer([1, 2, 3]) # no-op

      c = MyContainer([1, 2, 3])
      @test Trixi.move!(c, 1, 2) == MyContainer([0, 1, 3])
    end

    @testset "swap!" begin
      c = MyContainer([1,2])
      @test Trixi.swap!(c, 1, 1) == MyContainer([1, 2]) # no-op

      c = MyContainer([1,2])
      @test Trixi.swap!(c, 1, 2) == MyContainer([2,1])
    end

    @testset "erase!" begin
      c = MyContainer([1, 2])
      @test Trixi.erase!(c, 2, 1) == MyContainer([1, 2]) # no-op

      c = MyContainer([1, 2])
      @test Trixi.erase!(c, 1) == MyContainer([0, 2])
    end

    @testset "remove_shift!" begin
      c = MyContainer([1, 2, 3, 4])
      @test Trixi.remove_shift!(c, 2, 1) == MyContainer([1, 2, 3, 4]) # no-op

      c = MyContainer([1, 2, 3, 4])
      @test Trixi.remove_shift!(c, 2, 2) == MyContainer([1, 3, 4], 4)

      c = MyContainer([1, 2, 3, 4])
      @test Trixi.remove_shift!(c, 2) == MyContainer([1, 3, 4], 4)
    end

    @testset "remove_fill!" begin
      c = MyContainer([1, 2, 3, 4])
      @test Trixi.remove_fill!(c, 2, 1) == MyContainer([1, 2, 3, 4]) # no-op

      c = MyContainer([1, 2, 3, 4])
      @show "jo"
      @test Trixi.remove_fill!(c, 2, 2) == MyContainer([1, 4, 3], 4)
    end

    @testset "reset!" begin
      c = MyContainer([1, 2, 3])
      @test Trixi.reset!(c, 2) == MyContainer(Int[], 2)
    end
  end

  @testset "example parameters" begin
    @test basename(examples_dir()) == "examples"
    @test !isempty(get_examples())
    @test endswith(default_example(), "parameters.toml")
  end

  @testset "DG L2 mortar container debug output" begin
    c2d = Trixi.L2MortarContainer2D{1, 1}(1)
    @test isnothing(show(c2d))
    c3d = Trixi.L2MortarContainer3D{1, 1}(1)
    @test isnothing(show(c3d))
  end
end

end #module


text:
### A Pluto.jl notebook ###
# v0.12.18

using Markdown
using InteractiveUtils

# ╔═╡ 82deace6-507d-11eb-19f1-8711faa73635
using PlutoUI, Plots, LaTeXStrings

# ╔═╡ 9a10a468-507c-11eb-0fbb-0f453554e499
md"""
# Week 1: Secant and Tangent lines

**Author:** Marlin Figgins. **Date:** Jan 06, 2021.

The goal of this notebook is to elucidate the relationship between the secant and the tangent line to a curve. We'll also reframe this relationship in terms of average and instanteous rates of change.
"""

# ╔═╡ b0b6027a-507f-11eb-2e7a-27b3ca9569e5
md"""
## Secant Line

Simply put, a secant line is any straight line which joins two points on a curve.
"""

# ╔═╡ 1a0bb026-5080-11eb-25e4-57bac77b033d
begin
	f_1(x) = x^3 + x^2 - x - 1

	a1 = -1.5
	b1 = 2.25

	plot(f_1,
		-3, 3,
		label = L"f(x)",
		linewidth = 1.5,
		color = "black",
		xticks = false,
		yticks = false)

	plot!( x -> (f_1(b1) - f_1(a1))/(b1 - a1) * (x - a1) + f_1(a1),
		-3, 3,
		color = "orange",
		label = "Secant line")

	scatter!( [a1], [f_1(a1)],
		color = "white",
		markersize = 6,
		label = false)

	scatter!([b1], [f_1(b1)],
		color = "white",
		markersize = 6,
		label = false)
end

# ╔═╡ 371ed380-5082-11eb-3081-3f922e11d1d0
md"""
Since we've primarily interested in the functional forms of such curves, we'll re-cast the definition of the secant curve in terms of a function $f$. Let's say that we want to find secant line between two points $a$ and $b$ on the function $f$. We can write the points on the function itself as $(a, f(a))$ and $(b, f(b))$ which we visualize below.
"""

# ╔═╡ 8aa2e226-5082-11eb-1863-ad4351d89db2
begin
	plot(f_1,
		-3, 3,
		label = L"f(x)",
		linewidth = 1.5,
		color = "black",
		xticks = false,
		yticks = false)

	plot!( x -> (f_1(b1) - f_1(a1))/(b1 - a1) * (x - a1) + f_1(a1),
		-3, 3,
		color = "orange",
		label = "Secant line")

	scatter!( [a1], [f_1(a1)],
		color = "white",
		markersize = 6,
		label = false)

	annotate!( [a1], [f_1(a1) - 4.5], L"(a, f(a))", textsize = 5)

	scatter!([b1], [f_1(b1)],
		color = "white",
		markersize = 6,
		label = false)

	annotate!( [b1 + 0.3], [f_1(b1) - 4], L"(b, f(b))", textsize = 5)

end

# ╔═╡ 4c9acdfa-5083-11eb-3ef1-1d3bf8880bef
md"""
Since we know the position of two points on the secant line, we can write the slope of our secant line as

```math
m = \frac{\text{Rise}}{\text{Run}} = \frac{f(b) - f(a)}{b - a}
```

We can think of  this slope $m$ as describing the *average* change in the function $f$ between $a$ and $b$. Let's do a quick example to see what's going on here.
"""

# ╔═╡ b48e1642-5084-11eb-385e-8b50c677cb92
md"""

#### Example 1: Finding the secant line.

**Problem:** Compute the formula for the secant line between $a = 0$ and $b = 2$ for $f(x) = x^2$. In order to find the formula of the secant line, we need a slope and a point. We'll begin by finding the slope using the formula we derived above, so that

```math
m = \frac{f(2) - f(0)}{2 - 0} = \frac{2^2 - 0^2}{2 - 0} = 2.
```
Additionally, we know that the point $(a, f(a))$ lies on the curve of the function $f$, so we can solve for the equation of the tangent line as

```math
y - f(0) = m ( x - 0) \implies y = 2x.
```
"""

# ╔═╡ 9876359e-508f-11eb-20af-7be789d66b2c
begin
	f_2(x) = x^2
	a2 = 0
	b2 = 2

	plot(f_2,
		-3, 3,
		label = L"x^2",
		linewidth = 1.5,
		color = "black")

	plot!( x -> (f_2(b2) - f_2(a2))/(b2 - a2) * (x - a2) + f_2(a2),
		-3, 3,
		color = "orange",
		label = L"y = 2x")

	scatter!( [a2], [f_2(a2)],
		color = "white",
		markersize = 6,
		label = false)

	scatter!([b2], [f_2(b2)],
		color = "white",
		markersize = 6,
		label = false)

end

# ╔═╡ ce774fac-5085-11eb-139c-878de40d2658
md"""
## Tangent lines

In previous sections, we've shown how to find the tangent line to a circle using geometry and other properties of lines. We'll now expand this notion to more general curves. As opposed to the secant line, the tangent line is a line which touches a curve at a point and has the same direction of that curve at the point of contact.
"""

# ╔═╡ 67d2ebc0-5086-11eb-24c6-99f48646bcaa
begin

	f_1_prime(x) = 3*x^2 + 2*x - 1

	plot(f_1,
			-3, 3,
			label = L"f(x)",
			linewidth = 1.5,
			color = "black",
			xticks = false,
			yticks = false)

		plot!( x ->  f_1_prime(a1) * (x - a1) + f_1(a1),
			-3, 3,
			color = "green",
			label = "Tangent line")

		scatter!( [a1], [f_1(a1)],
			color = "white",
			markersize = 6,
			label = false)
end

# ╔═╡ b959d364-5086-11eb-0296-5d9363d44711
md"""
In many simple cases, we can find the tangent line by approximating it using secant lines. We'll use the above as an example. Suppose that we want to find the tangent line to $f$ at $x = a$. We can try and approximate this tangent line using the secant line between $a$ and $b$ for $b\approx a$. As $b$ gets closer to $a$, the slope of the secant line or the average rate of change of $f$ between $a$ and $b$ gets closer and closer to the slope of the tangent line or the *instanteous* rate of change of $f$ at $a$. We can see this process of convergence by picking $b$ sufficiently close to $a$, which we have done below.
"""

# ╔═╡ 27f6d934-5087-11eb-338a-dd3d28fe06c1
begin
	h1 = 0.05

		plot(f_1,
			-3, 3,
			label = L"f(x)",
			linewidth = 1.5,
			color = "black",
			xticks = false,
			yticks = false)

		plot!( x ->  f_1_prime(a1) * (x - a1) + f_1(a1),
			-3, 3,
			color = "green",
			label = "Tangent line")

		plot!( x -> (f_1(a1 + h1) - f_1(a1))/(h1) * (x - a1) + f_1(a1),
		-3, 3,
		color = "orange",
		label = "Secant line")

		scatter!( [a1], [f_1(a1)],
			color = "white",
			markersize = 6,
			label = false)

		scatter!( [a1 + h1], [f_1(a1 + h1)],
			color = "white",
			markersize = 6,
			label = false)
end

# ╔═╡ 75a962a0-5087-11eb-2d7b-6d00fa8c005c
md"""
As you can see, these two lines are exceedingly close. For several key functions, we can actually directly solve for the slopes of their tangent lines using approximating secant lines.
"""

# ╔═╡ ca762d9a-5087-11eb-1d3d-85bf5ecfe5ef
md"""

#### Example 2: Estimating the slope of the tangent line.
**Problem**: Estimate the slope of tangent line of $x^2$ at the point $x=0$.

We'll begin by finding the secant line between points $x=0$ and $x=a$. Following the equation derived above, we can compute the slope of this secant line as

```math
m = \frac{a^2 - 0^2}{a - 0} = \frac{a^2}{a}  = a.
```

Notice that this holds for any $a$, but since we're interested in values $a \approx 0$, we estimate that the slope of the tangent line is $m = 0$.
"""

# ╔═╡ 3d0c72fc-5091-11eb-3865-71342b0c55cc
begin
	plot(f_2,
		-3, 3,
		label = L"x^2",
		linewidth = 1.5,
		color = "black")

	plot!( x -> 0,
		-3, 3,
		color = "green",
		label = L"y = 0")

	scatter!( [a2], [f_2(a2)],
		color = "white",
		markersize = 6,
		label = false)
end

# ╔═╡ Cell order:
# ╟─9a10a468-507c-11eb-0fbb-0f453554e499
# ╟─b0b6027a-507f-11eb-2e7a-27b3ca9569e5
# ╟─1a0bb026-5080-11eb-25e4-57bac77b033d
# ╟─371ed380-5082-11eb-3081-3f922e11d1d0
# ╟─8aa2e226-5082-11eb-1863-ad4351d89db2
# ╟─4c9acdfa-5083-11eb-3ef1-1d3bf8880bef
# ╟─b48e1642-5084-11eb-385e-8b50c677cb92
# ╟─9876359e-508f-11eb-20af-7be789d66b2c
# ╟─ce774fac-5085-11eb-139c-878de40d2658
# ╟─67d2ebc0-5086-11eb-24c6-99f48646bcaa
# ╟─b959d364-5086-11eb-0296-5d9363d44711
# ╟─27f6d934-5087-11eb-338a-dd3d28fe06c1
# ╟─75a962a0-5087-11eb-2d7b-6d00fa8c005c
# ╟─ca762d9a-5087-11eb-1d3d-85bf5ecfe5ef
# ╟─3d0c72fc-5091-11eb-3865-71342b0c55cc
# ╟─82deace6-507d-11eb-19f1-8711faa73635


text:
module Courgette

using SparseArrays
using JuMP
using LinearAlgebra: dot, I, Diagonal
import Random

include("cut_generating_cone.jl")
include("normalizing.jl")

end # module


text:
const Simplex{M, N, T} = Union{MVector{M, SVector{N, T}}, SVector{M, SVector{N, T}}}
const SimplexType{M, N, T} = Union{Type{MVector{M, SVector{N, T}}}, Type{SVector{M, SVector{N, T}}}}

dimension(::SimplexType{M, N, T}) where {M,N,T} = Val(N)

any_inside(simplex::Simplex{M, N, T}) where {M,N,T} = Tagged(simplex[1])

function support_vector_max(simplex::Simplex{M, N, T}, direction, initial_guess::Tagged) where {M,N,T}
    best_pt, score = gt.argmax(p -> dot(p, direction), simplex)
    Tagged(best_pt)
end

@generated function simplex_face(simplex::Simplex{M, N, T}, i::Integer) where {N,M,T}
    simplex_face_impl(simplex, i)
end

function simplex_face_impl(simplex::SimplexType{M, N, T}, i) where {N,M,T}
    Expr(:call, :(SVector),
        Expr(:tuple, [:(i > $j ? simplex[$j] : simplex[$(j+1)]) for j in 1:(M - 1)]...))
end


text:
"""
"""
function QPointCellField(value::Number,cell_map::AbstractArray{<:Field},quad::CellQuadrature)
  q = get_coordinates(quad)
  v_q = [ fill(value,size(qi)) for qi in q ]
  array = ArrayOfEvaluatedFields(v_q,q)
  GenericCellField(array, cell_map)
end

"""
"""
function CellField(value::Number,cell_map::AbstractArray{<:Field},quad::CellQuadrature)
  QPointCellField(value,cell_map,quad)
end

struct EvaluatedField{A<:AbstractArray,P<:AbstractArray} <:Field
  array::A
  points::P
end

function field_cache(f::EvaluatedField,x)
  nothing
end

function evaluate_field!(cache,f::EvaluatedField,x)
  @assert length(x) == length(f.array)
  @assert x === f.points || x == f.points
  f.array
end

struct ArrayOfEvaluatedFields{T,N,A,B} <: AbstractArray{EvaluatedField{T},N}
  array::A
  points::B
  function ArrayOfEvaluatedFields(array::AbstractArray{T,N},points::AbstractArray) where {T<:AbstractArray,N}
    A = typeof(array)
    B = typeof(points)
    new{T,N,A,B}(array,points)
  end
end

Base.size(a::ArrayOfEvaluatedFields) = size(a.array)

Base.IndexStyle(::Type{<:ArrayOfEvaluatedFields{T,N,A}}) where {T,N,A} = IndexStyle(A)

@inline function Base.getindex(a::ArrayOfEvaluatedFields,i::Integer)
  EvaluatedField(a.array[i],a.points[i])
end

@inline function Base.getindex(a::ArrayOfEvaluatedFields{T,N},i::Vararg{Int,N}) where {T,N}
  EvaluatedField(a.array[i...],a.points[i...])
end

function array_cache(a::ArrayOfEvaluatedFields)
  ca = array_cache(a.array)
  cp = array_cache(a.points)
  (ca,cp)
end

@inline function getindex!(cache,a::ArrayOfEvaluatedFields,i...)
  ca, cp = cache
  array = getindex!(ca,a.array,i...)
  points = getindex!(ca,a.points,i...)
  EvaluatedField(array,points)
end

function evaluate_field_array(a::ArrayOfEvaluatedFields,x::AbstractArray)
  @assert a.points === x || a.points == x
  @assert length(a) == length(x)
  a.array
end

function update_state_variables!(updater::Function,quad::CellQuadrature,f::CellField...)
  x = get_coordinates(quad)
  update_state_variables!(updater,x,f...)
end

function update_state_variables!(updater::Function,x::AbstractArray,f::CellField...)
  fx = map(i->evaluate(i,x),f)
  caches = array_caches(fx...)
  cache_x = array_cache(x)
  _update_state_variables!(updater,caches,fx,cache_x,x)
end

function update_state_variables!(quad::CellQuadrature,updater::Function,f::CellField...)
  msg =
  """
  The method
      update_state_variables!(quad::CellQuadrature,updater::Function,f::CellField...)
  has been removed. Use
      update_state_variables!(updater::Function,quad::CellQuadrature,f::CellField...)
  instead
  """
  error(msg)
end

@noinline function  _update_state_variables!(updater,caches,fx,cache_x,x)
  ncells = length(x)
  for cell in 1:ncells
    fxi = getitems!(caches,fx,cell)
    xi = getindex!(cache_x,x,cell)
    for q in 1:length(xi)
      fxiq = getitems(fxi,q)
      r = updater(fxiq...)
      need_to_update, states = Arrays._split(r...)
      if need_to_update
        _update_states!(fxi,q,states,Val{length(states)}())
      end
    end
  end
end

@inline function _update_states!(b,q,states,::Val{i}) where i
  _update_state!(b,q,states,Val{i}())
  _update_states!(b,q,states,Val{i-1}())
  nothing
end

@inline function _update_states!(b,q,states,::Val{0})
  nothing
end

@inline function _update_state!(b,q,states,::Val{i}) where i
  m = length(b)
  n = length(states)
  o = m-n
  b[i+o][q] = states[i]
  nothing
end