text:
abstract type GlobalFunctionApproximator end

"""
    fit!(gfa::GlobalFunctionApproximator, dataset_input::AbstractMatrix, dataset_output::AbstractVector)
Fit the global function approximator to the dataset using some optimization method and a chosen
loss function.
"""
function fit! end

"""
    compute_value(gfa::GlobalFunctionApproximator, v::AbstractVector)
Return the value of the function at some query point v, based on the global function approximator

    compute_value(gfa::GlobalFunctionApproximator, v_list::AbstractVector{V}) where V <: AbstractVector{Float64}
Return the value of the function for a list of query points, based on the global function approximator
"""
function compute_value end


text:
import Dates

s = Dates.Date(1901, 1, 1)
c = s
e = Dates.Date(2000, 12, 31)
d = Dates.Day(1)

counter = 0
while c != e
  global c, counter
  if Dates.day(c) == 1 && Dates.dayofweek(c) == 7
    counter += 1
  end
  c += d
end
println(counter)

ext: jl

max_stars_repo_name: sujimodern/project-euler-jl

max_issues_repo_name: sujimodern/project-euler-jl

max_forks_repo_name: sujimodern/project-euler-jl


text:
# Example to demonstrate the use of a codemirror.

using ThreeJS
using LightGraphs
using NetworkViz

g = Graph(10)

main(window) =  begin
    push!(window.assets,("ThreeJS","threejs"))
    push!(window.assets,"widgets")
    push!(window.assets, "codemirror")
    push!(window.assets, "layout2")
    default = "g = CompleteGraph(10)"
    inp = Signal(Dict{Any, Any}(:name=>default))
    s = sampler() # A thing that lets you watch widgets/behaviors upon updates to other behaviors
    editor = watch!(s, :code, codemirror(default))
    code_cell = trigger!(s, :submit, keypress("ctrl+enter shift+enter", editor))
    t, plots = wire(
                        tabs(["3D View";"2D View";]),
                        pages(
                        [
                            map(inp) do f
                                fn = get(f,:code,default)
                                eval(parse(fn))
                                drawGraph(g,z=1)
                            end;
                            map(inp) do f
                                fn = get(f,:code,default)
                                eval(parse(fn))
                                drawGraph(g,z=0)
                            end;
                        ]
                        ),
                        :tab_channel,
                        :selected
                    )
        plugsampler(s,
        vbox(

                md"""Enter a valid graph function.
                    `ctrl+enter` or `shift+enter` to redraw the plot.
                    Use the mouse the drag, zoom and pan.
                    Try resizing the browser if you cant see a codebox.
                    """,
                code_cell,
                vskip(2em),
                t, plots
            ) |> pad(6em)
   ) >>> inp
end


text:
using BaseOnnx
using Test

@testset "BaseOnnx.jl" begin
    include("readwrite.jl")
end


text:
"""
    IOExtras

This module defines extensions to the `Base.IO` interface to support:
 - `startwrite`, `closewrite`, `startread` and `closeread` for streams
    with transactional semantics.
"""
module IOExtras

using ..Sockets
using MbedTLS: MbedException

export bytes, isbytes, nbytes, ByteView, nobytes, IOError, isioerror,
       startwrite, closewrite, startread, closeread,
       tcpsocket, localport, safe_getpeername

"""
    bytes(x)

If `x` is "castable" to an `AbstractVector{UInt8}`, then an
`AbstractVector{UInt8}` is returned; otherwise `x` is returned.
"""
function bytes end
bytes(s::AbstractVector{UInt8}) = s
bytes(s::AbstractString) = codeunits(s)
bytes(x) = x

"""whether `x` is "castable" to an `AbstractVector{UInt8}`; i.e. you can call `bytes(x)` if `isbytes(x)` === true"""
isbytes(x) = x isa AbstractVector{UInt8} || x isa AbstractString

"""
    nbytes(x) -> Int

Length in bytes of `x` if `x` is `isbytes(x)`.
"""
function nbytes end
nbytes(x) = nothing
nbytes(x::AbstractVector{UInt8}) = length(x)
nbytes(x::AbstractString) = sizeof(x)
nbytes(x::Vector{T}) where T <: AbstractString = sum(sizeof, x)
nbytes(x::Vector{T}) where T <: AbstractVector{UInt8} = sum(length, x)
nbytes(x::IOBuffer) = bytesavailable(x)
nbytes(x::Vector{IOBuffer}) = sum(bytesavailable, x)

"""
    isioerror(exception)

Is `exception` caused by a possibly recoverable IO error.
"""
isioerror(e) = false
isioerror(::Base.EOFError) = true
isioerror(::Base.IOError) = true
isioerror(e::ArgumentError) = e.msg == "stream is closed or unusable"
isioerror(::MbedException) = true

"""
    IOError <: Exception

The request terminated with due to an IO-related error.

Fields:
 - `e`, the error.
"""
struct IOError <: Exception
    e
    message
end

Base.show(io::IO, e::IOError) = print(io, "IOError(", e.e, " ", e.message, ")\n")

_doc = """
    startwrite(::IO)
    closewrite(::IO)
    startread(::IO)
    closeread(::IO)

Signal start/end of write or read operations.
"""
if isdefined(Base, :startwrite)
    "$_doc"
    Base.startwrite(io) = nothing
else
    "$_doc"
    startwrite(io) = nothing
end

if isdefined(Base, :closewrite)
    "$_doc"
    Base.closewrite(io) = nothing
else
    "$_doc"
    closewrite(io) = nothing
end

if isdefined(Base, :startread)
    "$_doc"
    Base.startread(io) = nothing
else
    "$_doc"
    startread(io) = nothing
end

if isdefined(Base, :closeread)
    "$_doc"
    Base.closeread(io) = nothing
else
    "$_doc"
    closeread(io) = nothing
end

using MbedTLS: SSLContext
tcpsocket(io::SSLContext)::TCPSocket = io.bio
tcpsocket(io::TCPSocket)::TCPSocket = io

localport(io) = try !isopen(tcpsocket(io)) ? 0 :
                    Sockets.getsockname(tcpsocket(io))[2]
                catch
                    0
                end

function safe_getpeername(io)
    try
        if isopen(tcpsocket(io))
            return Sockets.getpeername(tcpsocket(io))
        end
    catch
    end
    return IPv4(0), UInt16(0)
end

const ByteView = typeof(view(UInt8[], 1:0))
const nobytes = view(UInt8[], 1:0)

"""
Read from an `IO` stream until `find_delimiter(bytes)` returns non-zero.
Return view of bytes up to the delimiter.
"""
function Base.readuntil(buf::IOBuffer,
                    find_delimiter::Function #= Vector{UInt8} -> Int =#
                   )::ByteView

    l = find_delimiter(view(buf.data, buf.ptr:buf.size))
    if l == 0
        return nobytes
    end
    bytes = view(buf.data, buf.ptr:buf.ptr + l - 1)
    buf.ptr += l
    return bytes
end

end


text:
function ned2enu(p)
    if length(p) == 3
        return [p[2], p[1], -p[3]]
    else
        return error("Invalid position")
    end
end

function enu2ned(p)
    ned2enu(p)  # actually the same
end


text:
# # Polarizability by linear response using ForwardDiff

using DFTK
using LinearAlgebra
using ForwardDiff

a = 10.
lattice = a * I(3)  # cube of ``a`` bohrs
He = ElementPsp(:He, psp=load_psp("hgh/lda/He-q2"))
atoms = [He => [[1/2; 1/2; 1/2]]]  # Helium at the center of the box

kgrid = [1, 1, 1]  # no kpoint sampling for an isolated system
Ecut = 30
tol = 1e-8

## dipole moment of a given density (assuming the current geometry)
function dipole(basis, ρ)
    rr = [a * (r[1] - 1/2) for r in r_vectors(basis)]
    d = sum(rr .* ρ) * basis.dvol
end

function dipole(ε)
    basis_ε = make_basis(ε)
    res_ε = self_consistent_field(basis_ε, tol=tol)
    dipole(basis_ε, res_ε.ρ)
end

function make_basis(ε, lattice)
    model_ε = model_DFT(lattice, atoms, [:lda_x, :lda_c_vwn];
                        extra_terms=[ExternalFromReal(r -> -ε * (r[1] - a/2))],
                        symmetries=false)
    PlaneWaveBasis(model_ε; Ecut, kgrid)
end
make_basis(ε) = make_basis(ε, lattice)

# workaround: enforce the Model to promote its type from Float64 to Dual, by promoting the lattice too
make_basis(ε::T) where T <: ForwardDiff.Dual = make_basis(ε, T.(lattice))

# # Forward-mode implicit differentiation of SCF

# We keep both a non-dual basis, and a basis including Duals for easy bookkeeping (but redundant computation)
function self_consistent_field_dual(basis::PlaneWaveBasis, basis_dual::PlaneWaveBasis{T}; kwargs...) where T <: ForwardDiff.Dual
    scfres = self_consistent_field(basis; kwargs...)
    ψ = DFTK.select_occupied_orbitals(basis, scfres.ψ)
    filled_occ = DFTK.filled_occupation(basis.model)
    n_spin = basis.model.n_spin_components
    n_bands = div(basis.model.n_electrons, n_spin * filled_occ)
    # number of kpoints and occupation
    Nk = length(basis.kpoints)
    occupation = [filled_occ * ones(n_bands) for ik = 1:Nk]

    # promote everything eagerly to Dual numbers
    occupation_dual = [T.(occupation[1])]
    ψ_dual = [Complex.(T.(real(ψ[1])), T.(imag(ψ[1])))]
    ρ_dual = compute_density(basis_dual, ψ_dual, occupation_dual)

    _, δH = energy_hamiltonian(basis_dual, ψ_dual, occupation_dual; ρ=ρ_dual)
    δHψ = δH * ψ_dual
    δHψ = [ForwardDiff.partials.(δHψ[1], 1)]
    δψ = DFTK.solve_ΩplusK(basis, ψ, -δHψ, occupation)
    δρ = DFTK.compute_δρ(basis, ψ, δψ, occupation)
    ρ = ForwardDiff.value.(ρ_dual)
    ψ, ρ, δψ, δρ
end

function dipole(ε::ForwardDiff.Dual)
    T = ForwardDiff.tagtype(ε)
    basis = make_basis(ForwardDiff.value(ε))
    basis_dual = make_basis(ε)
    ψ, ρ, δψ, δρ = self_consistent_field_dual(basis, basis_dual; tol=tol)
    ρ_dual = ForwardDiff.Dual{T}.(ρ, δρ)
    dipole(basis_dual, ρ_dual)
end

polarizability = ForwardDiff.derivative(dipole, 0.0)
polarizability_ref = let
    ε = 1e-4
    (dipole(ε) - dipole(0.0)) / ε
end

println("Polarizability via ForwardDiff: $polarizability")
println("Polarizability via finite difference: $polarizability_ref")


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

const mintrials = 5
const mintime = 2000.0
print_output = isempty(ARGS)
codespeed = length(ARGS) > 0 && ARGS[1] == "codespeed"

if codespeed
    using JSON
    using HTTPClient.HTTPC

    # Ensure that we've got the environment variables we want:
    if !haskey(ENV, "JULIA_FLAVOR")
        error( "You must provide the JULIA_FLAVOR environment variable identifying this julia build!" )
    end

    # Setup codespeed data dict for submissions to codespeed's JSON endpoint.  These parameters
    # are constant across all benchmarks, so we'll just let them sit here for now
    csdata = Dict()
    csdata["commitid"] = Base.GIT_VERSION_INFO.commit
    csdata["project"] = "Julia"
    csdata["branch"] = Base.GIT_VERSION_INFO.branch
    csdata["executable"] = ENV["JULIA_FLAVOR"]
    csdata["environment"] = chomp(readstring(`hostname`))
    csdata["result_date"] = join( split(Base.GIT_VERSION_INFO.date_string)[1:2], " " )    #Cut the timezone out
end

# Takes in the raw array of values in vals, along with the benchmark name, description, unit and whether less is better
function submit_to_codespeed(vals,name,desc,unit,test_group,lessisbetter=true)
    # Points to the server
    codespeed_host = "julia-codespeed.csail.mit.edu"

    csdata["benchmark"] = name
    csdata["description"] = desc
    csdata["result_value"] = mean(vals)
    csdata["std_dev"] = std(vals)
    csdata["min"] = minimum(vals)
    csdata["max"] = maximum(vals)
    csdata["units"] = unit
    csdata["units_title"] = test_group
    csdata["lessisbetter"] = lessisbetter

    println( "$name: $(mean(vals))" )
    ret = post( "http://$codespeed_host/result/add/json/", Dict("json" => json([csdata])) )
    println( json([csdata]) )
    if ret.http_code != 200 && ret.http_code != 202
        error("Error submitting $name [HTTP code $(ret.http_code)], dumping headers and text: $(ret.headers)\n$(String(ret.body))\n\n")
        return false
    end
    return true
end

macro output_timings(t,name,desc,group)
    t = esc(t)
    name = esc(name)
    desc = esc(desc)
    group = esc(group)
    quote
        # If we weren't given anything for the test group, infer off of file path!
        test_group = length($group) == 0 ? basename(dirname(Base.source_path())) : $group[1]
        if codespeed
            submit_to_codespeed( $t, $name, $desc, "seconds", test_group )
        elseif print_output
            @printf "julia,%s,%f,%f,%f,%f\n" $name minimum($t) maximum($t) mean($t) std($t)
        end
        gc()
    end
end

macro timeit(ex,name,desc,group...)
    quote
        t = Float64[]
        tot = 0.0
        i = 0
        while i < mintrials || tot < mintime
            e = 1000*(@elapsed $(esc(ex)))
            tot += e
            if i > 0
                # warm up on first iteration
                push!(t, e)
            end
            i += 1
        end
        @output_timings t $(esc(name)) $(esc(desc)) $(esc(group))
    end
end

macro timeit_init(ex,init,name,desc,group...)
    quote
        t = zeros(mintrials)
        for i=0:mintrials
            $(esc(init))
            e = 1000*(@elapsed $(esc(ex)))
            if i > 0
                # warm up on first iteration
                t[i] = e
            end
        end
        @output_timings t $name $desc $group
    end
end

function maxrss(name)
    # FIXME: call uv_getrusage instead here
    @static if is_linux()
        rus = Array{Int64}(div(144,8))
        fill!(rus, 0x0)
        res = ccall(:getrusage, Int32, (Int32, Ptr{Void}), 0, rus)
        if res == 0
            mx = rus[5]/1024
            @printf "julia,%s.mem,%f,%f,%f,%f\n" name mx mx mx 0
        end
    end
end

# seed rng for more consistent timings
srand(1776)


text:
module Pricing

using Instruments

function blackformula()

end

end #module


text:
"""
Validation PSSE/HYGOV:
This case study defines a three bus system with an infinite bus, GENROU+SEXS+HYGOV and a load.
The fault drop the line connecting the infinite bus and GENROU
"""

##################################################
############### SOLVE PROBLEM ####################
##################################################

raw_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/HYGOV/ThreeBusMulti.raw")
dyr_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/HYGOV/ThreeBus_HYGOV.dyr")
csv_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/HYGOV/HYGOV_RESULTS.csv")

@testset "Test 31 HYGOV ResidualModel" begin
    path = (joinpath(pwd(), "test-psse-gast"))
    !isdir(path) && mkdir(path)
    try
        sys = System(raw_file, dyr_file)

        # Define Simulation Problem
        sim = Simulation!(
            ResidualModel,
            sys, #system
            path,
            (0.0, 20.0), #time span
            BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
        )

        # Test Initial Condition
        diff = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in test31_hygov_x0_init
            diff[1] += LinearAlgebra.norm(res[k] - v)
        end

        @test (diff[1] < 1e-3)

        # Obtain small signal results for initial conditions
        small_sig = small_signal_analysis(sim)
        eigs = small_sig.eigenvalues
        @test small_sig.stable

        # Test Eigenvalues
        @test LinearAlgebra.norm(eigs - test31_eigvals) < 1e-3

        # Solve problem
        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
              PSID.SIMULATION_FINALIZED
        results = read_results(sim)

        # Obtain results
        t, δ = get_state_series(results, ("generator-102-1", :δ))
        _, Vt = get_voltage_magnitude_series(results, 102)
        _, ω = get_state_series(results, ("generator-102-1", :ω))

        # Obtain PSSE results
        M = get_csv_data(csv_file)
        t_psse = M[:, 1]
        Vt_psse = M[:, 2]
        δ_psse = M[:, 3]
        ω_psse = M[:, 4] .+ 1.0

        # Test Transient Simulation Results
        # PSSE results are in Degrees
        @test LinearAlgebra.norm(δ - (δ_psse .* pi / 180), Inf) <= 1e-2
        @test LinearAlgebra.norm(t - round.(t_psse, digits = 3)) == 0.0
        @test LinearAlgebra.norm(Vt - Vt_psse, Inf) <= 1e-3
        @test LinearAlgebra.norm(ω - ω_psse, Inf) <= 1e-3
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end

@testset "Test 31 HYGOV MassMatrixModel" begin
    path = (joinpath(pwd(), "test-psse-gast"))
    !isdir(path) && mkdir(path)
    try
        sys = System(raw_file, dyr_file)

        # Define Simulation Problem
        sim = Simulation!(
            MassMatrixModel,
            sys, #system
            path,
            (0.0, 20.0), #time span
            BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
        )

        # Test Initial Condition
        diff = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in test31_hygov_x0_init
            diff[1] += LinearAlgebra.norm(res[k] - v)
        end

        @test (diff[1] < 1e-3)

        # Obtain small signal results for initial conditions
        small_sig = small_signal_analysis(sim)
        eigs = small_sig.eigenvalues
        @test small_sig.stable

        # Test Eigenvalues
        @test LinearAlgebra.norm(eigs - test31_eigvals) < 1e-3

        # Solve problem
        @test execute!(sim, Rodas4(), dtmax = 0.005, saveat = 0.005) ==
              PSID.SIMULATION_FINALIZED
        results = read_results(sim)

        # Obtain results
        t, δ = get_state_series(results, ("generator-102-1", :δ))
        _, Vt = get_voltage_magnitude_series(results, 102)
        _, ω = get_state_series(results, ("generator-102-1", :ω))

        # Obtain PSSE results
        M = get_csv_data(csv_file)
        t_psse = M[:, 1]
        Vt_psse = M[:, 2]
        δ_psse = M[:, 3]
        ω_psse = M[:, 4] .+ 1.0

        # Test Transient Simulation Results
        # PSSE results are in Degrees
        @test LinearAlgebra.norm(δ - (δ_psse .* pi / 180), Inf) <= 1e-2
        @test LinearAlgebra.norm(t - round.(t_psse, digits = 3)) == 0.0
        @test LinearAlgebra.norm(Vt - Vt_psse, Inf) <= 1e-3
        @test LinearAlgebra.norm(ω - ω_psse, Inf) <= 1e-3
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end


text:
@testset "DRO models" begin
    D = 10

    financial_data = get_financial_data(D=D)

    μ = mean(Matrix(financial_data[!, 1:D]), dims=1)[:]
    Σ = cov(Matrix(financial_data[!, 1:D]), dims=1)
    sqrtΣ = Matrix(cholesky(Σ).U)

    lower_bound = 0.0
    upper_bound = 25.0
    wealth = 100.0
    max_risk = 500.0

    solver = optimizer_with_attributes(
        SCS.Optimizer,
        "eps" => 1.0e-5,
        "max_iters" => 20000,
        "verbose" => false
    )

    function BasicWeightConvexSet(D, wealth, lower_bound, upper_bound, max_risk, μ, sqrtΣ)
        model = Model()

        @variable(model, w[i=1:D])

        @constraint(model, sum(w) <= wealth)
        @constraint(model, lower_bound .<= w .<= upper_bound)
        @constraint(model, [max_risk; sqrtΣ*w] in JuMP.SecondOrderCone())
        return model
    end

    basic_model = BasicWeightConvexSet(
        D,
        wealth,
        lower_bound,
        upper_bound,
        max_risk,
        μ,
        sqrtΣ
    )

    weights, constraints = DiffBasicWeightConvexSet(
        D,
        wealth,
        lower_bound,
        upper_bound,
        max_risk,
        μ,
        sqrtΣ
    )

    @testset "Bertsimas model" begin
        Δ = collect(0.1:0.1:1.0)
        Γ = 0.8

        formulation = RobustBertsimas(
            predicted_mean=μ,
            predicted_covariance=Σ,
            uncertainty_delta=Δ,
            bertsimas_budget=Γ,
        )
        model = copy(basic_model)
        @objective(model, Max, portfolio_return!(model, model[:w], formulation))
        w_po_opt = compute_solution(model, solver) * wealth

        w_diff_dro = DiffBertsimas(D, weights, constraints, μ, sqrtΣ).get_w(vcat(Δ, Γ)).data.numpy()

        @test w_po_opt ≈ w_diff_dro atol=0.001
    end

    @testset "BenTal model" begin
        δ = 10.0

        formulation = RobustBenTal(
            predicted_mean=μ,
            predicted_covariance=Σ,
            uncertainty_delta=δ,
        )
        model = copy(basic_model)
        @objective(model, Max, portfolio_return!(model, model[:w], formulation))
        w_po_opt = compute_solution(model, solver) * wealth

        w_diff_dro = DiffBenTal(D, weights, constraints, μ, sqrtΣ).get_w([δ]).data.numpy()

        @test w_po_opt ≈ w_diff_dro atol=0.001
    end

    @testset "Delage model" begin
        γ1 = 1.5
        γ2 = 1.1

        formulation = RobustDelague(
            predicted_mean=μ,
            predicted_covariance=Σ,
            γ1=γ1,
            γ2=γ2,
            utility_coeficients=[1.0],
            utility_intercepts=[0.0]
        )
        model = copy(basic_model)
        model = po_max_utility_return(formulation, current_wealth=wealth, model=model)
        w_po_opt = compute_solution(model, solver) * wealth

        w_diff_dro = DiffDelage(D, weights, constraints, μ, sqrtΣ).get_w([γ1, γ2]).data.numpy()

        @test w_po_opt ≈ w_diff_dro atol=0.01
    end
end

max_stars_repo_name: molet/DiffDRO

max_issues_repo_name: molet/DiffDRO

max_forks_repo_name: molet/DiffDRO


text:
module QuartzImageIO

using ImageCore, ColorTypes, ColorVectorSpace, FixedPointNumbers, Libdl
import FileIO: DataFormat, @format_str, File, Stream, filename, stream

const CFURLRef = Ptr{Cvoid}
const CFStringRef = Ptr{UInt8}
const CFDictionaryRef = Ptr{Cvoid}
const CGImageDestinationRef = Ptr{Cvoid}
const CGImageRef = Ptr{Cvoid}
const CGColorSpaceRef = Ptr{Cvoid}
const CGContextRef = Ptr{Cvoid}

load(imagefile::File{T}, args...; key_args...) where T <: DataFormat =
    load_(filename(imagefile), args...; key_args...)
load(filename::AbstractString, args...; key_args...) =
    load_(filename, args...; key_args...)
load(imgstream::Stream{T}, args...; key_args...) where T <: DataFormat =
    load_(read(stream(imgstream)), args...; key_args...)
load(imgstream::IO, args...; key_args...) =
    load_(read(imgstream), args...; key_args...)

save(imagefile::File{T}, args...; key_args...) where T <: DataFormat =
    save_(imagefile, args...; key_args...)
save(filename::AbstractString, args...; key_args...) =
    save_(filename, args...; key_args...)
save(imgstream::Stream{T}, args...; key_args...) where T <: DataFormat =
    save_(imgstream, args...; key_args...)

function load_(b::Array{UInt8, 1})
    data = CFDataCreate(b)
    imgsrc = CGImageSourceCreateWithData(data)
    CFRelease(data)
    read_and_release_imgsrc(imgsrc)
end

function load_(filename::String)
    myURL = CFURLCreateWithFileSystemPath(abspath(filename))
    imgsrc = CGImageSourceCreateWithURL(myURL)
    CFRelease(myURL)
    read_and_release_imgsrc(imgsrc)
end

## core, internal function
function read_and_release_imgsrc(imgsrc)
    if imgsrc == C_NULL
        @warn "QuartzImageIO created no image source"
        return nothing
    end
    # Get image information
    imframes = convert(Int, CGImageSourceGetCount(imgsrc))
    if imframes == 0
        # Bail out to ImageMagick
        @warn "QuartzImageIO found no frames"
        CFRelease(imgsrc)
        return nothing
    end
    dict = CGImageSourceCopyPropertiesAtIndex(imgsrc, 0)
    imheight = CFNumberGetValue(CFDictionaryGetValue(dict, "PixelHeight"), Int16)
    imwidth = CFNumberGetValue(CFDictionaryGetValue(dict, "PixelWidth"), Int16)
    isindexed = CFBooleanGetValue(CFDictionaryGetValue(dict, "IsIndexed"))
    if isindexed
        # Bail out to ImageMagick
        @warn "QuartzImageIO: indexed color images not implemented"
        CFRelease(imgsrc)
        return nothing
    end
    hasalpha = CFBooleanGetValue(CFDictionaryGetValue(dict, "HasAlpha"))

    pixeldepth = CFNumberGetValue(CFDictionaryGetValue(dict, "Depth"), Int16)
    # Colormodel is one of: "RGB", "Gray", "CMYK", "Lab"
    colormodel = CFStringGetCString(CFDictionaryGetValue(dict, "ColorModel"))
    if colormodel == ""
        # Bail out to ImageMagick
        @warn "QuartzImageIO found empty colormodel string"
        CFRelease(imgsrc)
        return nothing
    end
    imtype = CFStringGetCString(CGImageSourceGetType(imgsrc))
    alphacode, storagedepth = alpha_and_depth(imgsrc)

    # Get image description string
    # This is unused after the update to Images 0.6, but it would be nice to have again.
    imagedescription = ""
    if imtype == "public.tiff"
        tiffdict = CFDictionaryGetValue(dict, "{TIFF}")
        imagedescription = tiffdict != C_NULL ?
            CFStringGetCString(CFDictionaryGetValue(tiffdict,"ImageDescription")) :
            nothing
    end
    CFRelease(dict)

    # Allocate the buffer and get the pixel data
    sz = imframes > 1 ?
        (convert(Int, imwidth), convert(Int, imheight), convert(Int, imframes)) :
        (convert(Int, imwidth), convert(Int, imheight))
    ufixedtype = Dict(10 => N6f10, 12 => N4f12, 14 => N2f14, 16 => N0f16)
    T = pixeldepth <= 8 ? N0f8 : ufixedtype[pixeldepth]
    if colormodel == "Gray" && alphacode == 0 && storagedepth == 1
        buf = Array{Gray{T}}(undef, sz)
        fillgray!(reshape(reinterpret(T, buf), sz), imgsrc)
    elseif colormodel == "Gray" && alphacode ∈ [1, 3]
        buf = Array{GrayA{T}}(undef, sz)
        fillgrayalpha!(reshape(reinterpret(T, buf), (2, sz...)), imgsrc)
    elseif colormodel == "Gray" && alphacode ∈ [2, 4]
        buf = Array{AGray{T}}(undef, sz)
        fillgrayalpha!(reshape(reinterpret(T, buf), (2, sz...)), imgsrc)
    elseif colormodel == "RGB" && alphacode ∈ [1, 3]
        buf = Array{RGBA{T}}(undef, sz)
        fillcolor!(reshape(reinterpret(T, buf), (4, sz...)), imgsrc, storagedepth)
    elseif colormodel == "RGB" && alphacode ∈ [2, 4]
        buf = Array{ARGB{T}}(undef, sz)
        fillcolor!(reshape(reinterpret(T, buf), (4, sz...)), imgsrc, storagedepth)
    elseif colormodel == "RGB" && alphacode == 0
        buf = Array{RGB{T}}(undef, sz)
        fillcolor!(reshape(reinterpret(T, buf), (3, sz...)), imgsrc, storagedepth)
    elseif colormodel == "RGB" && alphacode ∈ [5, 6]
        buf = alphacode == 5 ? Array{RGB4{T}}(undef, sz) : Array{RGB1{T}}(undef, sz)
        fillcolor!(reshape(reinterpret(T, buf), (4, sz...)), imgsrc, storagedepth)
    else
        @warn "Unknown colormodel ($colormodel) and alphacode ($alphacode) found by QuartzImageIO"
        CFRelease(imgsrc)
        return nothing
    end
    CFRelease(imgsrc)
    # TODO: Override this flip with the Exif information, if available.  See:
    # https://github.com/JuliaIO/ImageMagick.jl/blob/f5fd22dbe5564e57e710de54254425e34fd6571c/src/ImageMagick.jl#L125
    permutedims(buf, [2; 1; 3:ndims(buf)])
end

function alpha_and_depth(imgsrc)
    CGimg = CGImageSourceCreateImageAtIndex(imgsrc, 0)  # Check only first frame
    alphacode = CGImageGetAlphaInfo(CGimg)
    bitspercomponent = CGImageGetBitsPerComponent(CGimg)
    bitsperpixel = CGImageGetBitsPerPixel(CGimg)
    CGImageRelease(CGimg)
    # Alpha codes documented here:
    # https://developer.apple.com/library/mac/documentation/graphicsimaging/reference/CGImage/Reference/reference.html#//apple_ref/doc/uid/TP30000956-CH3g-459700
    # Dividing bits per pixel by bits per component tells us how many
    # color + alpha components we have in each pixel
    alphacode, convert(Int, div(bitsperpixel, bitspercomponent))
end

function fillgray!(buffer::AbstractArray{T, 2}, imgsrc) where T
    imwidth, imheight = size(buffer, 1), size(buffer, 2)
    CGimg = CGImageSourceCreateImageAtIndex(imgsrc, 0)
    imagepixels = CopyImagePixels(CGimg)
    pixelptr = CFDataGetBytePtr(imagepixels, eltype(buffer))
    imbuffer = unsafe_wrap(Array, pixelptr, (imwidth, imheight), own=false)
    buffer[:, :] .= imbuffer
    CFRelease(imagepixels)
    CGImageRelease(CGimg)
end

# Image stack
function fillgray!(buffer::AbstractArray{T, 3}, imgsrc) where T
    imwidth, imheight, nimages = size(buffer, 1), size(buffer, 2), size(buffer, 3)
    for i in 1:nimages
        CGimg = CGImageSourceCreateImageAtIndex(imgsrc, i - 1)
        imagepixels = CopyImagePixels(CGimg)
        pixelptr = CFDataGetBytePtr(imagepixels, T)
        imbuffer = unsafe_wrap(Array, pixelptr, (imwidth, imheight), own=false)
        buffer[:, :, i] .= imbuffer
        CFRelease(imagepixels)
        CGImageRelease(CGimg)
    end
end

function fillgrayalpha!(buffer::AbstractArray{UInt8, 3}, imgsrc)
    imwidth, imheight = size(buffer, 2), size(buffer, 3)
    CGimg = CGImageSourceCreateImageAtIndex(imgsrc, 0)
    imagepixels = CopyImagePixels(CGimg)
    pixelptr = CFDataGetBytePtr(imagepixels, UInt16)
    imbuffer = unsafe_wrap(Array, pixelptr, (imwidth, imheight), own=false)
    buffer[1, :, :] .= imbuffer .& 0xff
    buffer[2, :, :] .= div.(imbuffer .& 0xff00, 256)
    CFRelease(imagepixels)
    CGImageRelease(CGimg)
end
fillgrayalpha!(buffer::AbstractArray{N0f8, 3}, imgsrc) =
    fillgrayalpha!(reinterpret(UInt8, buffer), imgsrc)

function fillcolor!(buffer::AbstractArray{T, 3}, imgsrc, nc) where T
    imwidth, imheight = size(buffer, 2), size(buffer, 3)
    CGimg = CGImageSourceCreateImageAtIndex(imgsrc, 0)
    imagepixels = CopyImagePixels(CGimg)
    pixelptr = CFDataGetBytePtr(imagepixels, T)
    imbuffer = unsafe_wrap(Array, pixelptr, (nc, imwidth, imheight), own=false)
    buffer[:, :, :] .= imbuffer
    CFRelease(imagepixels)
    CGImageRelease(CGimg)
end

function fillcolor!(buffer::AbstractArray{T, 4}, imgsrc, nc) where T
    imwidth, imheight, nimages = size(buffer, 2), size(buffer, 3), size(buffer, 4)
    for i in 1:nimages
        CGimg = CGImageSourceCreateImageAtIndex(imgsrc, i - 1)
        imagepixels = CopyImagePixels(CGimg)
        pixelptr = CFDataGetBytePtr(imagepixels, T)
        imbuffer = unsafe_wrap(Array, pixelptr, (nc, imwidth, imheight), own=false)
        buffer[:, :, :, i] .= imbuffer
        CFRelease(imagepixels)
        CGImageRelease(CGimg)
    end
end

## Saving Images ###############################################################

# For supported pixel formats, see Table 2-1 in:
# https://developer.apple.com/library/content/documentation/GraphicsImaging/Conceptual/drawingwithquartz2d/dq_context/dq_context.html#//apple_ref/doc/uid/TP30001066-CH203-BCIBHHBB

""" `save_(f, img, image_type)`

- f is the file to save to, of type `FileIO.File{DataFormat}`
- image_type should be one of Apple's image types (eg. "public.jpeg")
- permute_horizontal, if true, will transpose the image (flip x and y)
- mapi is the mapping to apply to the data before saving. Defaults to `identity`.
  A useful alternative value is `clamp01nan`.
"""
function save_(f::File{R}, img::AbstractArray;
               permute_horizontal=true, mapi=identity) where R <: DataFormat
    # Setup buffer
    local imgm
    try
        imgm = map(x -> mapCG(mapi(x)), img)
    catch
        @warn """QuartzImageIO: Mapping to the storage type failed.
                Perhaps your data had out-of-range values?
                Try `map(clamp01nan, img)` to clamp values to a valid range."""
        rethrow()
    end
    permute_horizontal && (imgm = permutedims_horizontal(imgm))
    ndims(imgm) > 3 && error("QuartzImageIO: At most 3 dimensions are supported for saving.")
    contig = to_contiguous(imgm)
    buf = to_explicit(contig)
    # Color type and order
    T = eltype(img)
    bitmap_info = zero(UInt32)
    if T <: Union{Gray, GrayA}
        # Note: GrayA is not an Apple supported combination, so we map to Gray
        bitmap_info |= kCGImageAlphaNone
        colorspace = CGColorSpaceCreateWithName("kCGColorSpaceGenericGray")
        components = 1
    elseif T <: Union{RGB, RGB4, HSV, HSL, Lab}
        bitmap_info |= kCGImageAlphaNoneSkipLast
        colorspace = CGColorSpaceCreateWithName("kCGColorSpaceSRGB")
        components = 4
    elseif T <: RGBA
        bitmap_info |= kCGImageAlphaPremultipliedLast
        colorspace = CGColorSpaceCreateWithName("kCGColorSpaceSRGB")
        components = 4
    elseif T <: ARGB
        bitmap_info |= kCGImageAlphaPremultipliedFirst
        colorspace = CGColorSpaceCreateWithName("kCGColorSpaceSRGB")
        components = 4
    elseif T <: Union{BGRA, ABGR}
        error("QuartzImageIO can only handle RGB byte orders")
    else
        error("QuartzImageIO: tried to save unknown buffer type $T")
    end
    # Image bit depth
    S = eltype(buf)
    if S <: Union{Int8, UInt8}
        bits_per_component = 8
        bitmap_info |= kCGBitmapByteOrderDefault
    elseif S <: Union{Int16, UInt16}
        bits_per_component = 16
        bitmap_info |= kCGBitmapByteOrder16Little
    elseif S <: Union{Int32, UInt32}
        # does this even exist?
        bits_per_component = 32
        bitmap_info |= kCGBitmapByteOrder32Little
    elseif S <: Float32
        bits_per_component = 32
        bitmap_info |= kCGBitmapFloatComponents
        bitmap_info |= kCGBitmapByteOrder32Little
    end
    # Image size
    width, height = size(contig)
    nframes = size(contig, 3)
    bytes_per_row = width*components*bits_per_component ÷ 8
    # Output type
    apple_format_names = Dict(format"BMP" => "com.microsoft.bmp",
                              format"GIF" => "com.compuserve.gif",
                              format"JPEG" => "public.jpeg",
                              format"PNG" => "public.png",
                              format"TIFF" => "public.tiff",
                              format"TGA" => "com.truevision.tga-image")
    image_type = apple_format_names[R]
    # Ready to save
    out_url = CFURLCreateWithFileSystemPath(filename(f))
    out_dest = CGImageDestinationCreateWithURL(out_url, image_type, nframes)
    if ndims(imgm) == 2
        bmp_context = CGBitmapContextCreate(buf, width, height, bits_per_component,
                                            bytes_per_row, colorspace, bitmap_info)
        out_image = CGBitmapContextCreateImage(bmp_context)
        CFRelease(bmp_context)
        CGImageDestinationAddImage(out_dest, out_image)
        CGImageRelease(out_image)
    else
        for i in 1:nframes
            bmp_context = CGBitmapContextCreate(buf[:,:,i], width, height,
                                                bits_per_component, bytes_per_row,
                                                colorspace, bitmap_info)
            out_image = CGBitmapContextCreateImage(bmp_context)
            CFRelease(bmp_context)
            CGImageDestinationAddImage(out_dest, out_image)
            CGImageRelease(out_image)
        end
    end
    CGImageDestinationFinalize(out_dest)
    CFRelease(out_dest)
    CFRelease(colorspace)
    CFRelease(out_url)
    nothing
end

function save_(io::Stream, img::AbstractArray; permute_horizontal=true, mapi = clamp01nan)
    write(io, getblob(img, permute_horizontal, mapi))
end

function getblob(img::AbstractArray, permute_horizontal, mapi)
    # In theory we could save the image directly to a buffer via
    # CGImageDestinationCreateWithData - TODO. But I couldn't figure out how
    # to get the length of the CFMutableData object. So I take the inefficient
    # route of saving the image to a temporary file for now.
    temp_file = joinpath(tempdir(), "QuartzImageIO_temp.png")
    save(File(format"PNG", temp_file), img,
         permute_horizontal=permute_horizontal, mapi=mapi)
    read(open(temp_file))
end

# Element-mapping function. Converts to RGB/RGBA and uses
# N0f8 "inner" element type.
const Color1{T} = Color{T, 1}
const Color2{T, C<:Color1} = TransparentColor{C, T, 2}
const Color3{T} = Color{T, 3}
const Color4{T, C<:Color3} = TransparentColor{C, T, 4}

mapCG(c::Color1) = mapCG(convert(Gray, c))
mapCG(c::Gray{T}) where T = convert(Gray{N0f8}, c)
mapCG(c::Gray{T}) where T <: Normed = c

mapCG(c::Color2) = mapCG(convert(Gray, c))
mapCG(c::GrayA{T}) where T = convert(Gray{N0f8}, c)
mapCG(c::GrayA{T}) where T <: Normed = convert(Gray, c)

# Note: macOS does not handle 3 channel buffers, only 4,
# but we can tell it to use or skip that 4th (alpha) channel
mapCG(c::Color3) = mapCG(convert(RGBA, c))
mapCG(c::RGB{T}) where T = convert(RGBA{N0f8}, c)
mapCG(c::RGB{T}) where T <: Normed = convert(RGBA{T}, c)
mapCG(c::RGB4{T}) where T <: Real = convert(RGBA{T}, c)

mapCG(c::Color4) = mapCG(convert(RGBA, c))
mapCG(c::RGBA{T}) where T = convert(RGBA{N0f8}, c)
mapCG(c::RGBA{T}) where T <: Normed = c

mapCG(x::UInt8) = reinterpret(N0f8, x)
mapCG(x::Bool) = convert(N0f8, x)
mapCG(x::AbstractFloat) = convert(N0f8, x)
mapCG(x::Normed) = x

# Make the data contiguous in memory, because writers don't handle stride.
to_contiguous(A::Array) = A
to_contiguous(A::AbstractArray) = collect(A)
to_contiguous(A::BitArray) = convert(Array{N0f8}, A)
to_contiguous(A::ColorView) = to_contiguous(channelview(A))

to_explicit(A::Array{C}) where C <: Colorant = to_explicit(channelview(A))
to_explicit(A::Base.ReinterpretArray{T}) where T = to_explicit(copyto!(Array{T}(undef, size(A)), A))
to_explicit(A::Array{T}) where T <: Normed = rawview(A)
to_explicit(A::Array{Float32}) = A
to_explicit(A::Array{T}) where T <: AbstractFloat = to_explicit(convert(Array{N0f8}, A))

permutedims_horizontal(img::AbstractVector) = img
function permutedims_horizontal(img)
    # Vertical-major is hard-coded here
    p = [2; 1; 3:ndims(img)]
    permutedims(img, p)
end

## OSX Framework Wrappers ######################################################

# Commented out functions remain here because they might be useful for future
# debugging.

const foundation = Libdl.find_library(["/System/Library/Frameworks/Foundation.framework/Resources/BridgeSupport/Foundation"])
const imageio = Libdl.find_library(["/System/Library/Frameworks/ImageIO.framework/ImageIO"])

const kCFNumberSInt8Type = 1
const kCFNumberSInt16Type = 2
const kCFNumberSInt32Type = 3
const kCFNumberSInt64Type = 4
const kCFNumberFloat32Type = 5
const kCFNumberFloat64Type = 6
const kCFNumberCharType = 7
const kCFNumberShortType = 8
const kCFNumberIntType = 9
const kCFNumberLongType = 10
const kCFNumberLongLongType = 11
const kCFNumberFloatType = 12
const kCFNumberDoubleType = 13
const kCFNumberCFIndexType = 14
const kCFNumberNSIntegerType = 15
const kCFNumberCGFloatType = 16
const kCFNumberMaxType = 16

# enum defined at https://developer.apple.com/library/mac/documentation/GraphicsImaging/Reference/CGImage/index.html#//apple_ref/c/tdef/CGImageAlphaInfo
include("CG_const.jl")

# Objective-C and NS wrappers
oms(id, uid) = ccall(:objc_msgSend, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}), id, selector(uid))

ogc(id) = ccall((:objc_getClass, "Cocoa.framework/Cocoa"), Ptr{Cvoid}, (Ptr{UInt8}, ), id)

selector(sel::String) = ccall(:sel_getUid, Ptr{Cvoid}, (Ptr{UInt8}, ), sel)

NSString(init::String) = ccall(:objc_msgSend, Ptr{Cvoid},
                               (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{UInt8}, UInt64),
                               oms(ogc("NSString"), "alloc"),
                               selector("initWithCString:encoding:"), init, 4)

# NSLog(str::String, obj) = ccall((:NSLog, foundation), Ptr{Cvoid},
#                                 (Ptr{Cvoid}, Ptr{Cvoid}), NSString(str), obj)

# NSLog(str::String) = ccall((:NSLog, foundation), Ptr{Cvoid},
#                            (Ptr{Cvoid}, ), NSString(str))

# NSLog(obj::Ptr) = ccall((:NSLog, foundation), Ptr{Cvoid}, (Ptr{Cvoid}, ), obj)

# Core Foundation
# General

# CFRetain(CFTypeRef::Ptr{Cvoid}) = CFTypeRef != C_NULL &&
#     ccall(:CFRetain, Void, (Ptr{Cvoid}, ), CFTypeRef)

CFRelease(CFTypeRef::Ptr{Cvoid}) = CFTypeRef != C_NULL &&
    ccall(:CFRelease, Nothing, (Ptr{Cvoid}, ), CFTypeRef)

# function CFGetRetainCount(CFTypeRef::Ptr{Cvoid})
#     CFTypeRef == C_NULL && return 0
#     ccall(:CFGetRetainCount, Clonglong, (Ptr{Cvoid}, ), CFTypeRef)
# end

CFShow(CFTypeRef::Ptr{Cvoid}) = CFTypeRef != C_NULL &&
    ccall(:CFShow, Nothing, (Ptr{Cvoid}, ), CFTypeRef)

# function CFCopyDescription(CFTypeRef::Ptr{Cvoid})
#     CFTypeRef == C_NULL && return C_NULL
#     ccall(:CFCopyDescription, Ptr{Cvoid}, (Ptr{Cvoid}, ), CFTypeRef)
# end

# CFCopyTypeIDDescription(CFTypeID::Cint) = CFTypeRef != C_NULL &&
#     ccall(:CFCopyTypeIDDescription, Ptr{Cvoid}, (Cint, ), CFTypeID)

# function CFGetTypeID(CFTypeRef::Ptr{Cvoid})
#     CFTypeRef == C_NULL && return nothing
#     ccall(:CFGetTypeID, Culonglong, (Ptr{Cvoid}, ), CFTypeRef)
# end

# CFURLCreateWithString(filename) =
#     ccall(:CFURLCreateWithString, Ptr{Cvoid},
#           (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}), C_NULL, NSString(filename), C_NULL)

CFURLCreateWithFileSystemPath(filename::String) =
    ccall(:CFURLCreateWithFileSystemPath, Ptr{Cvoid},
          (Ptr{Cvoid}, Ptr{Cvoid}, Cint, Bool), C_NULL, NSString(filename), 0, false)

# CFDictionary

# CFDictionaryGetKeysAndValues(CFDictionaryRef::Ptr{Cvoid}, keys, values) =
#     CFDictionaryRef != C_NULL &&
#     ccall(:CFDictionaryGetKeysAndValues, Void,
#           (Ptr{Cvoid}, Ptr{Ptr{Cvoid}}, Ptr{Ptr{Cvoid}}), CFDictionaryRef, keys, values)

function CFDictionaryGetValue(CFDictionaryRef::Ptr{Cvoid}, key)
    CFDictionaryRef == C_NULL && return C_NULL
    ccall(:CFDictionaryGetValue, Ptr{Cvoid},
          (Ptr{Cvoid}, Ptr{Cvoid}), CFDictionaryRef, key)
end

CFDictionaryGetValue(CFDictionaryRef::Ptr{Cvoid}, key::String) =
    CFDictionaryGetValue(CFDictionaryRef::Ptr{Cvoid}, NSString(key))

# CFNumber
function CFNumberGetValue(CFNum::Ptr{Cvoid}, numtype)
    CFNum == C_NULL && return nothing
    out = Cint[0]
    ccall(:CFNumberGetValue, Bool, (Ptr{Cvoid}, Cint, Ptr{Cint}), CFNum, numtype, out)
    out[1]
end

CFNumberGetValue(CFNum::Ptr{Cvoid}, ::Type{Int8}) =
    CFNumberGetValue(CFNum, kCFNumberSInt8Type)

CFNumberGetValue(CFNum::Ptr{Cvoid}, ::Type{Int16}) =
    CFNumberGetValue(CFNum, kCFNumberSInt16Type)

CFNumberGetValue(CFNum::Ptr{Cvoid}, ::Type{Int32}) =
    CFNumberGetValue(CFNum, kCFNumberSInt32Type)

CFNumberGetValue(CFNum::Ptr{Cvoid}, ::Type{Int64}) =
    CFNumberGetValue(CFNum, kCFNumberSInt64Type)

CFNumberGetValue(CFNum::Ptr{Cvoid}, ::Type{Float32}) =
    CFNumberGetValue(CFNum, kCFNumberFloat32Type)

CFNumberGetValue(CFNum::Ptr{Cvoid}, ::Type{Float64}) =
    CFNumberGetValue(CFNum, kCFNumberFloat64Type)

CFNumberGetValue(CFNum::Ptr{Cvoid}, ::Type{UInt8}) =
    CFNumberGetValue(CFNum, kCFNumberCharType)

#CFBoolean
CFBooleanGetValue(CFBoolean::Ptr{Cvoid}) =
    CFBoolean != C_NULL &&
    ccall(:CFBooleanGetValue, Bool, (Ptr{Cvoid}, ), CFBoolean)

# CFString
function CFStringGetCString(CFStringRef::Ptr{Cvoid})
    CFStringRef == C_NULL && return ""
    buffer = Array{UInt8}(undef, 1024)  # does this need to be bigger for Open Microscopy TIFFs?
    res = ccall(:CFStringGetCString, Bool, (Ptr{Cvoid}, Ptr{UInt8}, UInt, UInt16),
                CFStringRef, buffer, length(buffer), 0x0600)
    res == C_NULL && return ""
    return unsafe_string(pointer(buffer))
end

# These were unsafe, can return null pointers at random times.
# See Apple Developer Docs
# CFStringGetCStringPtr(CFStringRef::Ptr{Cvoid}) =
#     ccall(:CFStringGetCStringPtr, Ptr{UInt8}, (Ptr{Cvoid}, UInt16), CFStringRef, 0x0600)
#
# getCFString(CFStr::Ptr{Cvoid}) = CFStringGetCStringPtr(CFStr) != C_NULL ?
#     bytestring(CFStringGetCStringPtr(CFStr)) : ""

# Core Graphics
# CGImageSource

CGImageSourceCreateWithURL(myURL::Ptr{Cvoid}) =
    ccall((:CGImageSourceCreateWithURL, imageio), Ptr{Cvoid},
          (Ptr{Cvoid}, Ptr{Cvoid}), myURL, C_NULL)

CGImageSourceCreateWithData(data::Ptr{Cvoid}) =
    ccall((:CGImageSourceCreateWithData, imageio), Ptr{Cvoid},
          (Ptr{Cvoid}, Ptr{Cvoid}), data, C_NULL)

CGImageSourceGetType(CGImageSourceRef::Ptr{Cvoid}) =
    ccall(:CGImageSourceGetType, Ptr{Cvoid}, (Ptr{Cvoid}, ), CGImageSourceRef)

CGImageSourceGetStatus(CGImageSourceRef::Ptr{Cvoid}) =
    ccall(:CGImageSourceGetStatus, UInt32, (Ptr{Cvoid}, ), CGImageSourceRef)

CGImageSourceGetStatusAtIndex(CGImageSourceRef::Ptr{Cvoid}, n) =
    ccall(:CGImageSourceGetStatusAtIndex, Int32,
          (Ptr{Cvoid}, Csize_t), CGImageSourceRef, n) #Int32?

# CGImageSourceCopyProperties(CGImageSourceRef::Ptr{Cvoid}) =
#     ccall(:CGImageSourceCopyProperties, Ptr{Cvoid},
#           (Ptr{Cvoid}, Ptr{Cvoid}), CGImageSourceRef, C_NULL)

CGImageSourceCopyPropertiesAtIndex(CGImageSourceRef::Ptr{Cvoid}, n) =
    ccall(:CGImageSourceCopyPropertiesAtIndex, Ptr{Cvoid},
          (Ptr{Cvoid}, Csize_t, Ptr{Cvoid}), CGImageSourceRef, n, C_NULL)

CGImageSourceGetCount(CGImageSourceRef::Ptr{Cvoid}) =
    ccall(:CGImageSourceGetCount, Csize_t, (Ptr{Cvoid}, ), CGImageSourceRef)

CGImageSourceCreateImageAtIndex(CGImageSourceRef::Ptr{Cvoid}, i) =
    ccall(:CGImageSourceCreateImageAtIndex, Ptr{Cvoid},
          (Ptr{Cvoid}, UInt64, Ptr{Cvoid}), CGImageSourceRef, i, C_NULL)

# CGImageGet

CGImageGetAlphaInfo(CGImageRef::Ptr{Cvoid}) =
    ccall(:CGImageGetAlphaInfo, UInt32, (Ptr{Cvoid}, ), CGImageRef)

# Use this to detect if image contains floating point values
# CGImageGetBitmapInfo(CGImageRef::Ptr{Cvoid}) =
#     ccall(:CGImageGetBitmapInfo, UInt32, (Ptr{Cvoid}, ), CGImageRef)

CGImageGetBitsPerComponent(CGImageRef::Ptr{Cvoid}) =
    ccall(:CGImageGetBitsPerComponent, Csize_t, (Ptr{Cvoid}, ), CGImageRef)

CGImageGetBitsPerPixel(CGImageRef::Ptr{Cvoid}) =
    ccall(:CGImageGetBitsPerPixel, Csize_t, (Ptr{Cvoid}, ), CGImageRef)

# CGImageGetBytesPerRow(CGImageRef::Ptr{Cvoid}) =
#     ccall(:CGImageGetBytesPerRow, Csize_t, (Ptr{Cvoid}, ), CGImageRef)

CGImageGetColorSpace(CGImageRef::Ptr{Cvoid}) =
    ccall(:CGImageGetColorSpace, UInt32, (Ptr{Cvoid}, ), CGImageRef)

# CGImageGetDecode(CGImageRef::Ptr{Cvoid}) =
#     ccall(:CGImageGetDecode, Ptr{Float64}, (Ptr{Cvoid}, ), CGImageRef)

# CGImageGetHeight(CGImageRef::Ptr{Cvoid}) =
#     ccall(:CGImageGetHeight, Csize_t, (Ptr{Cvoid}, ), CGImageRef)

# CGImageGetRenderingIntent(CGImageRef::Ptr{Cvoid}) =
#     ccall(:CGImageGetRenderingIntent, UInt32, (Ptr{Cvoid}, ), CGImageRef)

# CGImageGetShouldInterpolate(CGImageRef::Ptr{Cvoid}) =
#     ccall(:CGImageGetShouldInterpolate, Bool, (Ptr{Cvoid}, ), CGImageRef)

# CGImageGetTypeID() =
#     ccall(:CGImageGetTypeID, Culonglong, (),)

# CGImageGetWidth(CGImageRef::Ptr{Cvoid}) =
#     ccall(:CGImageGetWidth, Csize_t, (Ptr{Cvoid}, ), CGImageRef)

CGImageRelease(CGImageRef::Ptr{Cvoid}) =
    ccall(:CGImageRelease, Nothing, (Ptr{Cvoid}, ), CGImageRef)

# Get pixel data
# See: https://developer.apple.com/library/mac/qa/qa1509/_index.html
CGImageGetDataProvider(CGImageRef::Ptr{Cvoid}) =
    ccall(:CGImageGetDataProvider, Ptr{Cvoid}, (Ptr{Cvoid}, ), CGImageRef)

CGDataProviderCopyData(CGDataProviderRef::Ptr{Cvoid}) =
    ccall(:CGDataProviderCopyData, Ptr{Cvoid}, (Ptr{Cvoid}, ), CGDataProviderRef)

CopyImagePixels(inImage::Ptr{Cvoid}) =
    CGDataProviderCopyData(CGImageGetDataProvider(inImage))

CFDataGetBytePtr(CFDataRef::Ptr{Cvoid}, ::Type{T}) where T =
    ccall(:CFDataGetBytePtr, Ptr{T}, (Ptr{Cvoid}, ), CFDataRef)

# CFDataGetLength(CFDataRef::Ptr{Cvoid}) =
#     ccall(:CFDataGetLength, Ptr{Int64}, (Ptr{Cvoid}, ), CFDataRef)

CFDataCreate(bytes::Array{UInt8,1}) =
    ccall(:CFDataCreate, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{UInt8}, Csize_t),
          C_NULL, bytes, length(bytes))

### For output #################################################################

""" `check_null(x)`

Triggers an error if `x` is `NULL`, else returns `x` """
function check_null(x)
    # Poor-man's error handling. TODO: raise more specific exceptions.
    if x == C_NULL
        error("C call returned NULL")
    else
        x
    end
end

function CGImageDestinationCreateWithURL(url::CFURLRef,
                                         filetype::String,
                                         count::Integer,
                                         options::CFDictionaryRef=C_NULL)
    check_null(ccall((:CGImageDestinationCreateWithURL, imageio),
                     CGImageDestinationRef,
                     (CFURLRef, CFStringRef, Csize_t, CFDictionaryRef),
                     url, NSString(filetype), count, options))
end

function CGImageDestinationAddImage(dest::CGImageDestinationRef,
                                    image::CGImageRef,
                                    properties::CFDictionaryRef=C_NULL)
    # Returns NULL.
    # From the Apple docs: "The function logs an error if you add more images
    # than what you specified when you created the image destination. "
    # Maybe we should catch that somehow?
    ccall((:CGImageDestinationAddImage, imageio),
          Ptr{Cvoid},
          (CGImageDestinationRef, CGImageRef, CFDictionaryRef),
          dest, image, properties)
end

mutable struct WritingImageFailed <: Exception end

function CGImageDestinationFinalize(dest::CGImageDestinationRef)
    rval = ccall((:CGImageDestinationFinalize, imageio),
                 Bool, (CGImageDestinationRef,), dest)
    # See https://developer.apple.com/library/mac/documentation/GraphicsImaging/Reference/CGImageDestination/index.html#//apple_ref/c/func/CGImageDestinationFinalize
    if !rval throw(WritingImageFailed()) end
end

CGColorSpaceCreateDeviceRGB() =
    check_null(ccall((:CGColorSpaceCreateDeviceRGB, imageio), CGColorSpaceRef, ()))

# Valid colorspace names are (excluding deprecated values):
# kCGColorSpaceGenericRGBLinear, kCGColorSpaceGenericCMYK, kCGColorSpaceGenericXYZ,
# kCGColorSpaceGenericGrayGamma2_2, kCGColorSpaceExtendedGray, kCGColorSpaceLinearGray,
# kCGColorSpaceExtendedLinearGray, kCGColorSpaceSRGB, kCGColorSpaceLinearSRGB,
# kCGColorSpaceExtendedLinearSRGB, kCGColorSpaceDCIP3, kCGColorSpaceDisplayP3,
# kCGColorSpaceAdobeRGB1998, kCGColorSpaceACESCGLinear, kCGColorSpaceITUR_709,
# kCGColorSpaceITUR_2020, kCGColorSpaceROMMRGB
CGColorSpaceCreateWithName(name::String) =
    check_null(ccall((:CGColorSpaceCreateWithName, imageio), CGColorSpaceRef,
                     (CFStringRef,), NSString(name)))

function CGBitmapContextCreate(data, # void*
                               width, height, # size_t
                               bitsPerComponent, bytesPerRow, # size_t
                               space::CGColorSpaceRef,
                               bitmapInfo) # uint32_t
    check_null(ccall((:CGBitmapContextCreate, imageio),
                     CGContextRef,
                     (Ptr{Cvoid}, Csize_t, Csize_t, Csize_t, Csize_t,
                      CGColorSpaceRef, UInt32),
                     data, width, height, bitsPerComponent, bytesPerRow,
                     space, bitmapInfo))
end

function CGBitmapContextCreateImage(context_ref::CGContextRef)
    check_null(ccall((:CGBitmapContextCreateImage, imageio),
                     CGImageRef,
                     (CGContextRef,), context_ref))
end

function __init__()
    if !Sys.isapple()
        # change this to an error in future
        @warn("""
              QuartzImageIO.jl can only be used on Apple macOS. Suggested usage is
                  @static if Sys.isapple()
                      using QuartzImageIO
                      # QuartzImageIO specific code goes here
                  end
              """)
    end
end

end # Module


text:
#-------------------------------------------------------------------------------
#"""
#Containing all relevant parameters:
#"""
@everywhere mutable struct Params
    structure_net_vf ::Array{Int,1};
    structure_net_qf ::Array{Int,1};
    structure_net_qf2::Array{Int,1};
    structure_net_pl ::Array{Int,1};
    init_range       ::Float64;
    lr_pl            ::Float64;
    lr_vf            ::Float64;
    lr_qf            ::Float64;
    σ                ::Array{Float64,1};
    dt               ::Float64;
    t_horizon        ::Float64;
    batch_size       ::Int;
    l_traj           ::Int;
    n_traj           ::Int;
    n_iter           ::Int;
    γ_disc           ::Float64;
    u_range          ::Array{Float64,1};
    s_init           ::Array{Float64,1};
    b_t              ::Array{Float64,1};
    l_s_t            ::Int;
    l_out            ::Int;
    λ_critic         ::Float64;
    λ_actor          ::Float64;
    eps              ::Float64;
    δ                ::Float64;
    cg_damping       ::Float64;
    τ                ::Float64;
    buffer_size      ::Int;
    σ_off_policy     ::Float64;
    clip_σ           ::Float64;
    d_update         ::Int;
    Params() = new(); #to instanziate easily!
end
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
#"""
#Containing the local saving of one gradient
#"""
@everywhere mutable struct Net
    grad  ::Array{Any,1}
    Net() = new()
end

"""
Initialize the Network
"""
function init_net(θ::Any)
    n_net = Net()
    n_net.grad = copy(θ) .* 0
    return n_net
end

"""
Copy the Network
"""
function copy_net(o_net::Net)
    c_net = Net()
    c_net.grad = copy(o_net.grad)
    return c_net
end
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
#"""
#Containing all data of one trajectory
#"""
@everywhere mutable struct Trajectory
    #Struct of the Data collected in one Episode
    s_t    ::Array{Float64,2};   #states
    s_t2   ::Array{Float64,2};   #new_states
    a_t    ::Array{Float64,2};   #actions
    p_a_t  ::Array{Float64,2};   #Probability of actions
    r_t    ::Array{Float64,1};   #rewards
    val_t  ::Array{Float64,1};   #used value function results
    out_t  ::Array{Float64,2};   #NN-output
    A_t    ::Array{Float64,1};   #Local Advantage
    R_t    ::Array{Float64,1};   #Local discounted reward
    len    ::Int;                #Actual length of the Trajectory

    Trajectory() = new(); #for easy instanziation
end

#"""
#Create new Instance of trajectory
#"""
@everywhere function init_traj(params::Params,l_traj::Int)
    #Get from param struct
    n_traj = params.n_traj;
    l_s_t  = params.l_s_t;
    l_out  = params.l_out;

    #Preallocate required data
    new_traj          = Trajectory()
    new_traj.s_t      = zeros(l_traj,l_s_t)
    new_traj.s_t2     = zeros(l_traj,l_s_t)
    new_traj.a_t      = zeros(l_traj,l_out)
    new_traj.p_a_t    = zeros(l_traj,l_out)
    new_traj.r_t      = zeros(l_traj)
    new_traj.val_t    = zeros(l_traj)
    new_traj.out_t    = zeros(l_traj,l_out)
    new_traj.A_t      = zeros(l_traj)
    new_traj.R_t      = zeros(l_traj)
    new_traj.len      = 0

    return new_traj
end

#"""
#Copy the trajectory
#"""
@everywhere function copy_traj(o_traj::Trajectory,params::Params)
    #Get from param struct
    l_traj = params.l_traj;

    #start coping:
    c_traj          = Trajectory()
    c_traj.s_t      = copy(o_traj.s_t)
    c_traj.s_t2     = copy(o_traj.s_t2)
    c_traj.a_t      = copy(o_traj.a_t)
    c_traj.p_a_t    = copy(o_traj.p_a_t)
    c_traj.r_t      = copy(o_traj.r_t)
    c_traj.out_t    = copy(o_traj.out_t)
    c_traj.A_t      = copy(o_traj.A_t)
    c_traj.R_t      = copy(o_traj.R_t)
    c_traj.val_t    = copy(o_traj.val_t)
    c_traj.len      = copy(o_traj.len)

    return c_traj
end
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
#"""
#Containing all data of all trajectories
#"""
@everywhere mutable struct Traj_Data
    traj      ::Array{Trajectory,1} #All data of Trajectories
    timesteps ::Int
    n_traj    ::Int

    Traj_Data() = new(); #for easy instanziation
end

"""
Initialize the trajectory data:
"""
function init_all_traj(params::Params)
    n_traj = params.n_traj;
    l_traj = params.l_traj;
    l_s_t  = params.l_s_t;

    #all_traj = Traj_Data(Array{Trajectory,1}(undef,n_traj),0,n_traj);
    all_traj = Traj_Data()
    all_traj.traj = Array{Trajectory,1}(undef,n_traj)
    all_traj.timesteps = 0
    all_traj.n_traj = n_traj

    #Create instance of trajectory:
    traj = init_traj(params,l_traj)

    #Fill the big struct with the small one:
    for i_traj = 1:n_traj
        all_traj.traj[i_traj] = copy_traj(traj,params)
    end

    return all_traj;
end

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

#Struct Running Mean_Square
#-------------------------------------------------------------------------------
@everywhere mutable struct RunMeanStd
    mean  ::Array{Float64,1}
    var   ::Array{Float64,1}
    std   ::Array{Float64,1}
    count ::Int
end

#Init the thing:
function init_rms(s_t::Array{Float64,1})
    len = length(s_t)
    rms = RunMeanStd(zeros(len) , zeros(len) .+1 , zeros(len) .+1 , 2)
    return rms
end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
export @P_str

struct PotID
    path::Tuple{Vararg{Symbol}}
end

"""
    P"[/]your/pot/path"

The path can be either relative or absolute path. If a relative path is provided, it will be resolved to an absolute path based on the current context.

!!! note
    We don't validate the path for you during construction. A [`PotNotRegisteredError`](@ref) will be thrown when you try to send messages to an unregistered path.
"""
macro P_str(s)
    PotID(s)
end

function Base.show(io::IO, p::PotID)
    if isempty(getfield(p, :path))
        print(io, "/")
    else
        for x in getfield(p, :path)
            print(io, '/')
            print(io, x)
        end
    end
end

function PotID(s::String)
    if length(s) > 0
        if s[1] == '/'
            PotID(Tuple(Symbol(x) for x in split(s, '/';keepempty=false)))
        else
            self() / s
        end
    else
        PotID(())
    end
end

function Base.:(/)(p::PotID, s::String)
    PotID((getfield(p, :path)..., Symbol(s)))
end

const ROOT = P"/"
const LOGGER = P"/log"
const SCHEDULER = P"/scheduler"
const USER = P"/user"


text:
using SpectroscopicUnits
using Unitful
using Test

@testset "SpectroscopicUnits.jl" begin
    @testset "convert Hz to length" begin
        @test ustrip(convert_unit(m, 1e15Hz)) ≈ 2.9979e-7 atol = 1e-11
        @test ustrip(convert_unit(cm, 1e15Hz)) ≈ 2.9979e-5 atol = 1e-9
        @test ustrip(convert_unit(mm, 1e15Hz)) ≈ 0.00029979 atol = 1e-8
        @test ustrip(convert_unit(μm, 1e15Hz)) ≈ 0.29979 atol = 1e-5
        @test ustrip(convert_unit(nm, 1e15Hz)) ≈ 299.79246 atol = 1e-5
    end

    @testset "convert radHz to length" begin
        @test ustrip(convert_unit(m, 2.0*pi*1e15radHz)) ≈ 2.9979e-7 atol = 1e-11
        @test ustrip(convert_unit(cm, 2.0*pi*1e15radHz)) ≈ 2.9979e-5 atol = 1e-9
        @test ustrip(convert_unit(mm, 2.0*pi*1e15radHz)) ≈ 0.00029979 atol = 1e-8
        @test ustrip(convert_unit(μm, 2.0*pi*1e15radHz)) ≈ 0.29979 atol = 1e-5
        @test ustrip(convert_unit(nm, 2.0*pi*1e15radHz)) ≈ 299.79246 atol = 1e-5
    end

    @testset "convert eV to length" begin
        @test ustrip(convert_unit(nm, 3.2eV)) ≈ 387.451 atol = 1e-3
        @test ustrip(convert_unit(μm, 3.2eV)) ≈ 0.38745 atol = 1e-5
    end
end


text:
import ODLCommonTools.eqn_deepcopy
"""
  SimpleODEMod.eqn_deepcopy

  This function performs a proper deepcopy (unlike julia's builtin deepcopy)
    on an Euler equation object.
  It preserves reference topology (i.e. q & q_vec pointing to same array in DG schemes).

    Inputs:
      eqn
      mesh
      sbp
      opts

    Outputs:
      eqn_copy

    One reason for doing this is this case:
      a = rand(2,2)
      b = a
      a[3] = 8
      b[3] == 8
      this is because 'a[3] =' is actually setindex!

"""
function eqn_deepcopy(mesh::AbstractMesh{Tmsh}, sbp, eqn::SimpleODEData_{Tsol, Tres}, opts::Dict) where {Tmsh, Tsol, Tres}

  # The eqn object has over 100 fields, so it is necessary to write a better approach for
  #   copying than explicitly copying every named field
  # This is the second write of eqn_deepcopy; the first version explicitly copied each field.
  #   Was done for SimpleODE and Advection before Euler made it obvious that that was intractable.

  # 1: call constructor on eqn_copy

  param_tuple = getAllTypeParams(mesh, eqn, opts)

  eqn_copy = SimpleODEData_{param_tuple...}(mesh, sbp, opts)
  copy!(eqn_copy, eqn)

  return eqn_copy

end


text:
# firstpassagedistributions.jl

using LinearAlgebra
using Distributions
using Roots  # needed for numerically finding quantiles
using ExponentialAction  # More numerically stable than built-in exp(::Matrix)

"""
    fpdistribution(T, A, p0)

Creates a new first-passage time distribution from the transient transition rate matrix T,
the absorbing transition rate matrix A, and the vector of initial probabilities of the
transient states p0.

Checks to make sure that the inputs meet the requirements for being a first passage time
distribution for a continuous-time, discrete-state Markov process.
"""
struct fpdistribution{T1, T2, T3} <: ContinuousUnivariateDistribution
    T::T1  # transient matrix
    A::T2  # absorbing matrix
    p0::T3  # initial condition

    # Internal constructor function
    fpdistribution(T::T1, A::T2, p0::T3) where {T1, T2, T3} = begin
        # At least one column of T sum to be less than zero
        @assert any(sum(T, dims=1) .< 0) "Transient T matrix incorrect:\n$T\n"
        # Checks to make sure that the total exit rate from T is equal to the total rate of
        # entry into A
        @assert isapprox(vec(-ones(1, size(T, 1)) * T), vec(A'*ones(size(A, 1)))) "Absorbing matrix A incorrect"
        # Size of p0 should correspond to the number of transient states
        @assert size(p0, 1) == size(T, 1) == size(A, 2) "Dimension mismatch with T, A, and/or p0"
        # p0 should be a probability distribution
        @assert isapprox(sum(p0), 1.0) "Initial condition p0 not a probability distribution"
        new{T1, T2, T3}(T, A, p0)
    end
end

# Helper functions
Distributions.@distr_support fpdistribution 0 +Inf
Distributions.minimum(d::fpdistribution) = 0.0
Distributions.maximum(d::fpdistribution) = +Inf
Distributions.mean(d::fpdistribution) = -sum(inv(d.T) * d.p0)
Distributions.var(d::fpdistribution) = 2*sum(d.T^(-2) * d.p0) - mean(d)^2

"""
    pdf(d::fpdistribution, t::Real)

Return the probability density function of d evaluated at the value t.
"""
Distributions.pdf(d::fpdistribution, t::Real) = begin
    #insupport(d, t) ? sum(d.A * exp(t * d.T) * d.p0) : zero(t)
    insupport(d, t) ? sum(d.A * expv(t, d.T, d.p0)) : zero(t)
end

"""
    pdf(d::fpdistribution, t, dims)

Return the conditional probability density at `t` for the absorbing dimensions given in
`dims`.
"""
Distributions.pdf(d::fpdistribution, t::Real, dims) = begin
    #insupport(d, t) ? getindex(d.A * exp(t * d.T) * d.p0 ./ splittingprobabilities(d), dims) : zero(t)
    insupport(d, t) ? getindex(d.A * expv(t, d.T, d.p0) ./ splittingprobabilities(d), dims) : zero(t)
end

"""
    logpdf(d::fpdistribution, t::Real)

Returns the log probability density of `d` evaluated at `t`.
"""
Distributions.logpdf(d::fpdistribution, t::Real) = begin
    log(pdf(d, t))
end

"""
    logpdf(d::fpdistribution, t::Real)

Returns the sum of the log probability densities of `d` evaluated at each value in `t`.
"""
Distributions.logpdf(d::fpdistribution, t::AbstractVector) = begin
    sum(log.(pdf(d, x) for x in t))
end

"""
    logpdf(d::fpdistribution, t::Real)

Returns the log probability density of `d` evaluated at `t` for the absorbing dimensions in
`dims`.
"""
Distributions.logpdf(d::fpdistribution, t::Real, dims) = begin
    log.(pdf(d, t, dims))
end

"""
    cdf(d::fpdistribution, t::Real)

Returns the cumulative distribution function of `d` evaluated at `t`.
"""
Distributions.cdf(d::fpdistribution, t::Real) = begin
    insupport(d, t) ? 1 - sum(expv(t, d.T, d.p0)) : zero(t)
end

Distributions.cdf(d::fpdistribution, t::Real, dims) = begin
    error("cdf not implemented for multiple absorbing states.")
end

"""
    quantile(d::fpdistribution, p)

Return the ``p``-th quantile for the first-passage time distribution `d`.

The quantile function for these distributions has no closed-form solution, so this method
uses the `find_zero` function from the `Roots` package to find the quantile numerically.
"""
Distributions.quantile(d::fpdistribution, p) = begin
    if p <= zero(p)
        return zero(p)
    elseif p >= one(p)
        return Inf
    else
        # Use the mean of the distribution as the initial guess
        find_zero((x -> cdf(d, x) - p), mean(d))
    end
end

Distributions.quantile(d::fpdistribution, p, dims) = begin
    error("Quantile function not implemented for multiple absorbing states.")
end

Distributions.rand(d::fpdistribution, rng::AbstractVector{<:Real}) = begin
    quantile(d, rng)
end

"""
    splittingprobabilities(T, A, p0)

Compute the splitting probabilities, i.e., the probabilities of being absorbed into each of
the absorbing states.
"""
function splittingprobabilities(d::fpdistribution)
    return -d.A * inv(d.T) * d.p0
end

"""
    quasistationary(d::fpdistribution)

Return the quasi-stationary distribution of the distribution `d`, i.e., the probability
distribution over transient states conditional on not yet being absorbed.
"""
function quasistationary(d::fpdistribution)
    vals, vecs = eigen(d.T)
    idx = argmax(vals)
    vec = vecs[:,idx]
    return vec / sum(vec)
end


text:
module AD

using ..NiLang
using NiLangCore
using Reexport
@reexport using NiLangCore.ADCore
using MLStyle

import ..NiLang: ⊕, ⊖, NEG, CONJ, ROT, IROT, SWAP

NiLangCore.ADCore.GVar(x::Dup) = Dup(GVar(x.x), GVar(x.twin))
(invg::Type{Inv{GVar}})(x::Dup) = Dup(invg(x.x), invg(x.twin))

include("instructs_basic.jl")
include("instructs_ext.jl")
include("simple_hessian.jl")
include("taylor.jl")

for op in [:>, :<, :>=, :<=, :isless]
    @eval Base.$op(a::Bundle, b::Bundle) = $op(value(a), value(b))
    @eval Base.$op(a::Bundle, b) = $op(value(a), b)
    @eval Base.$op(a, b::Bundle) = $op(a, value(b))
end
#Base.conj(x::GVar) = GVar(x.x', x.g')
end


text:
# Used for objectives and solvers where the gradient is available/exists
mutable struct OnceDifferentiable{TF, TDF, TX} <: AbstractObjective
    f # objective
    df # (partial) derivative of objective
    fdf # objective and (partial) derivative of objective
    F::TF # cache for f output
    DF::TDF # cache for df output
    x_f::TX # x used to evaluate f (stored in F)
    x_df::TX # x used to evaluate df (stored in DF)
    f_calls::Vector{Int}
    df_calls::Vector{Int}
end

### Only the objective
# Ambiguity
OnceDifferentiable(f, x::AbstractArray,
                   F::Real = real(zero(eltype(x))),
                   DF::AbstractArray = alloc_DF(x, F); inplace = true, autodiff = :finite,
                   chunk::ForwardDiff.Chunk = ForwardDiff.Chunk(x)) =
    OnceDifferentiable(f, x, F, DF, autodiff, chunk)
#OnceDifferentiable(f, x::AbstractArray, F::AbstractArray; autodiff = :finite) =
#    OnceDifferentiable(f, x::AbstractArray, F::AbstractArray, alloc_DF(x, F))
function OnceDifferentiable(f, x::AbstractArray,
                   F::AbstractArray, DF::AbstractArray = alloc_DF(x, F);
                   inplace = true, autodiff = :finite)
    f! = f!_from_f(f, F, inplace)

    OnceDifferentiable(f!, x::AbstractArray, F::AbstractArray, DF, autodiff)
end

function OnceDifferentiable(f, x_seed::AbstractArray{T},
                            F::Real,
                            DF::AbstractArray,
                            autodiff, chunk) where T
    # When here, at the constructor with positional autodiff, it should already
    # be the case, that f is inplace.
    if  typeof(f) <: Union{InplaceObjective, NotInplaceObjective}

        fF = make_f(f, x_seed, F)
        dfF = make_df(f, x_seed, F)
        fdfF = make_fdf(f, x_seed, F)

        return OnceDifferentiable(fF, dfF, fdfF, x_seed, F, DF)
    else
        if is_finitediff(autodiff)

            # Figure out which Val-type to use for FiniteDiff based on our
            # symbol interface.
            fdtype = finitediff_fdtype(autodiff)
            df_array_spec = DF
            x_array_spec = x_seed
            return_spec = typeof(F)
            gcache = FiniteDiff.GradientCache(df_array_spec, x_array_spec, fdtype, return_spec)

            function g!(storage, x)
                FiniteDiff.finite_difference_gradient!(storage, f, x, gcache)
                return
            end
            function fg!(storage, x)
                g!(storage, x)
                return f(x)
            end
        elseif is_forwarddiff(autodiff)
            gcfg = ForwardDiff.GradientConfig(f, x_seed, chunk)
            g! = (out, x) -> ForwardDiff.gradient!(out, f, x, gcfg)

            fg! = (out, x) -> begin
                gr_res = DiffResults.DiffResult(zero(T), out)
                ForwardDiff.gradient!(gr_res, f, x, gcfg)
                DiffResults.value(gr_res)
            end
        else
            error("The autodiff value $autodiff is not support. Use :finite or :forward.")
        end

        return OnceDifferentiable(f, g!, fg!, x_seed, F, DF)
    end
end

has_not_dep_symbol_in_ad = Ref{Bool}(true)
OnceDifferentiable(f, x::AbstractArray, F::AbstractArray, autodiff::Symbol, chunk::ForwardDiff.Chunk = ForwardDiff.Chunk(x)) =
OnceDifferentiable(f, x, F, alloc_DF(x, F), autodiff, chunk)
function OnceDifferentiable(f, x::AbstractArray, F::AbstractArray,
                            autodiff::Bool, chunk::ForwardDiff.Chunk = ForwardDiff.Chunk(x))
    if autodiff == false
        throw(ErrorException("It is not possible to set the `autodiff` keyword to `false` when constructing a OnceDifferentiable instance from only one function. Pass in the (partial) derivative or specify a valid `autodiff` symbol."))
    elseif has_not_dep_symbol_in_ad[]
        @warn("Setting the `autodiff` keyword to `true` is deprecated. Please use a valid symbol instead.")
        has_not_dep_symbol_in_ad[] = false
    end
    OnceDifferentiable(f, x, F, alloc_DF(x, F), :forward, chunk)
end
function OnceDifferentiable(f, x_seed::AbstractArray, F::AbstractArray, DF::AbstractArray,
    autodiff::Symbol , chunk::ForwardDiff.Chunk = ForwardDiff.Chunk(x_seed))
    if  typeof(f) <: Union{InplaceObjective, NotInplaceObjective}
        fF = make_f(f, x_seed, F)
        dfF = make_df(f, x_seed, F)
        fdfF = make_fdf(f, x_seed, F)
        return OnceDifferentiable(fF, dfF, fdfF, x_seed, F, DF)
    else
        if is_finitediff(autodiff)
            # Figure out which Val-type to use for FiniteDiff based on our
            # symbol interface.
            fdtype = finitediff_fdtype(autodiff)
            # Apparently only the third input is aliased.
            j_finitediff_cache = FiniteDiff.JacobianCache(copy(x_seed), copy(F), copy(F), fdtype)
            if autodiff == :finiteforward
                # These copies can be done away with if we add a keyword for
                # reusing arrays instead for overwriting them.
                Fx = copy(F)
                DF = copy(DF)

                x_f, x_df = x_of_nans(x_seed), x_of_nans(x_seed)
                f_calls, j_calls = [0,], [0,]
                function j_finiteforward!(J, x)
                    # Exploit the possibility that it might be that x_f == x
                    # then we don't have to call f again.

                    # if at least one element of x_f is different from x, update
                    if any(x_f .!= x)
                        f(Fx, x)
                        f_calls .+= 1
                    end

                    FiniteDiff.finite_difference_jacobian!(J, f, x, j_finitediff_cache, Fx)
                end
                function fj_finiteforward!(F, J, x)
                    f(F, x)
                    FiniteDiff.finite_difference_jacobian!(J, f, x, j_finitediff_cache, F)
                end

                return OnceDifferentiable(f, j_finiteforward!, fj_finiteforward!, Fx, DF, x_f, x_df, f_calls, j_calls)
            end

            function fj_finitediff!(F, J, x)
                f(F, x)
                FiniteDiff.finite_difference_jacobian!(J, f, x, j_finitediff_cache)
                F
            end
            function j_finitediff!(J, x)
                F_cache = copy(F)
                fj_finitediff!(F_cache, J, x)
            end

            return OnceDifferentiable(f, j_finitediff!, fj_finitediff!, x_seed, F, DF)

        elseif is_forwarddiff(autodiff)

            jac_cfg = ForwardDiff.JacobianConfig(f, F, x_seed, chunk)
            ForwardDiff.checktag(jac_cfg, f, x_seed)

            F2 = copy(F)
            function j_forwarddiff!(J, x)
                ForwardDiff.jacobian!(J, f, F2, x, jac_cfg, Val{false}())
            end
            function fj_forwarddiff!(F, J, x)
                jac_res = DiffResults.DiffResult(F, J)
                ForwardDiff.jacobian!(jac_res, f, F2, x, jac_cfg, Val{false}())
                DiffResults.value(jac_res)
            end

            return OnceDifferentiable(f, j_forwarddiff!, fj_forwarddiff!, x_seed, F, DF)
        else
            error("The autodiff value $(autodiff) is not supported. Use :finite or :forward.")
        end
    end
end

### Objective and derivative
function OnceDifferentiable(f, df,
                   x::AbstractArray,
                   F::Real = real(zero(eltype(x))),
                   DF::AbstractArray = alloc_DF(x, F);
                   inplace = true)

    df! = df!_from_df(df, F, inplace)

    fdf! = make_fdf(x, F, f, df!)

    OnceDifferentiable(f, df!, fdf!, x, F, DF)
end

function OnceDifferentiable(f, j,
                   x::AbstractArray,
                   F::AbstractArray,
                   J::AbstractArray = alloc_DF(x, F);
                   inplace = true)

    f! = f!_from_f(f, F, inplace)
    j! = df!_from_df(j, F, inplace)
    fj! = make_fdf(x, F, f!, j!)

    OnceDifferentiable(f!, j!, fj!, x, F, J)
end

### Objective, derivative and combination
function OnceDifferentiable(f, df, fdf,
    x::AbstractArray,
    F::Real = real(zero(eltype(x))),
    DF::AbstractArray = alloc_DF(x, F);
    inplace = true)

    # f is never "inplace" since F is scalar
    df! = df!_from_df(df, F, inplace)
    fdf! = fdf!_from_fdf(fdf, F, inplace)

    x_f, x_df = x_of_nans(x), x_of_nans(x)

    OnceDifferentiable{typeof(F),typeof(DF),typeof(x)}(f, df!, fdf!,
    copy(F), copy(DF),
    x_f, x_df,
    [0,], [0,])
end

function OnceDifferentiable(f, df, fdf,
                            x::AbstractArray,
                            F::AbstractArray,
                            DF::AbstractArray = alloc_DF(x, F);
                            inplace = true)

    f = f!_from_f(f, F, inplace)
    df! = df!_from_df(df, F, inplace)
    fdf! = fdf!_from_fdf(fdf, F, inplace)

    x_f, x_df = x_of_nans(x), x_of_nans(x)

    OnceDifferentiable(f, df!, fdf!, copy(F), copy(DF), x_f, x_df, [0,], [0,])
end


text:
include("vector.jl")
include("points/point.jl")
include("points/sort.jl")
include("axisalignedbox.jl")
include("polytopes/polytope.jl")
include("polytopes/interpolate.jl")
include("polytopes/jacobian.jl")
include("polytopes/faces.jl")
include("polytopes/edges.jl")
include("polytopes/measure.jl")
include("polytopes/triangulate.jl")

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
#####
##### Velocity gradients
#####

# Diagonal
@inline ∂x_u(i, j, k, grid, u) = ∂xᶜᵃᵃ(i, j, k, grid, u)
@inline ∂y_v(i, j, k, grid, v) = ∂yᵃᶜᵃ(i, j, k, grid, v)
@inline ∂z_w(i, j, k, grid, w) = ∂zᵃᵃᶜ(i, j, k, grid, w)

# Off-diagonal
@inline ∂x_v(i, j, k, grid, v) = ∂xᶠᵃᵃ(i, j, k, grid, v)
@inline ∂x_w(i, j, k, grid, w) = ∂xᶠᵃᵃ(i, j, k, grid, w)

@inline ∂y_u(i, j, k, grid, u) = ∂yᵃᶠᵃ(i, j, k, grid, u)
@inline ∂y_w(i, j, k, grid, w) = ∂yᵃᶠᵃ(i, j, k, grid, w)

@inline ∂z_u(i, j, k, grid, u) = ∂zᵃᵃᶠ(i, j, k, grid, u)
@inline ∂z_v(i, j, k, grid, v) = ∂zᵃᵃᶠ(i, j, k, grid, v)

#####
##### Strain components
#####

# ccc strain components
@inline Σ₁₁(i, j, k, grid, u) = ∂xᶜᵃᵃ(i, j, k, grid, u)
@inline Σ₂₂(i, j, k, grid, v) = ∂yᵃᶜᵃ(i, j, k, grid, v)
@inline Σ₃₃(i, j, k, grid, w) = ∂zᵃᵃᶜ(i, j, k, grid, w)

@inline tr_Σ(i, j, k, grid, u, v, w) =
    Σ₁₁(i, j, k, grid, u) + Σ₂₂(i, j, k, grid, v) + Σ₃₃(i, j, k, grid, w)

# ffc
@inline Σ₁₂(i, j, k, grid::AbstractGrid{FT}, u, v) where FT =
    FT(0.5) * (∂yᵃᶠᵃ(i, j, k, grid, u) + ∂xᶠᵃᵃ(i, j, k, grid, v))

# fcf
@inline Σ₁₃(i, j, k, grid::AbstractGrid{FT}, u, w) where FT =
    FT(0.5) * (∂zᵃᵃᶠ(i, j, k, grid, u) + ∂xᶠᵃᵃ(i, j, k, grid, w))

# cff
@inline Σ₂₃(i, j, k, grid::AbstractGrid{FT}, v, w) where FT =
    FT(0.5) * (∂zᵃᵃᶠ(i, j, k, grid, v) + ∂yᵃᶠᵃ(i, j, k, grid, w))

@inline Σ₁₂²(i, j, k, grid, u, v) = Σ₁₂(i, j, k, grid, u, v)^2
@inline Σ₁₃²(i, j, k, grid, u, w) = Σ₁₃(i, j, k, grid, u, w)^2
@inline Σ₂₃²(i, j, k, grid, v, w) = Σ₂₃(i, j, k, grid, v, w)^2

#####
##### Renamed functions for consistent function signatures
#####

@inline ∂x_u(i, j, k, grid, u, v, w) = ∂x_u(i, j, k, grid, u)
@inline ∂x_v(i, j, k, grid, u, v, w) = ∂x_v(i, j, k, grid, v)
@inline ∂x_w(i, j, k, grid, u, v, w) = ∂x_w(i, j, k, grid, w)

@inline ∂y_u(i, j, k, grid, u, v, w) = ∂y_u(i, j, k, grid, u)
@inline ∂y_v(i, j, k, grid, u, v, w) = ∂y_v(i, j, k, grid, v)
@inline ∂y_w(i, j, k, grid, u, v, w) = ∂y_w(i, j, k, grid, w)

@inline ∂z_u(i, j, k, grid, u, v, w) = ∂z_u(i, j, k, grid, u)
@inline ∂z_v(i, j, k, grid, u, v, w) = ∂z_v(i, j, k, grid, v)
@inline ∂z_w(i, j, k, grid, u, v, w) = ∂z_w(i, j, k, grid, w)

@inline Σ₁₁(i, j, k, grid, u, v, w) = Σ₁₁(i, j, k, grid, u)
@inline Σ₂₂(i, j, k, grid, u, v, w) = Σ₂₂(i, j, k, grid, v)
@inline Σ₃₃(i, j, k, grid, u, v, w) = Σ₃₃(i, j, k, grid, w)

@inline Σ₁₂(i, j, k, grid, u, v, w) = Σ₁₂(i, j, k, grid, u, v)
@inline Σ₁₃(i, j, k, grid, u, v, w) = Σ₁₃(i, j, k, grid, u, w)
@inline Σ₂₃(i, j, k, grid, u, v, w) = Σ₂₃(i, j, k, grid, v, w)

# Symmetry relations
const Σ₂₁ = Σ₁₂
const Σ₃₁ = Σ₁₃
const Σ₃₂ = Σ₂₃

# Trace and squared strains
@inline tr_Σ²(ijk...) = Σ₁₁(ijk...)^2 +  Σ₂₂(ijk...)^2 +  Σ₃₃(ijk...)^2

@inline Σ₁₂²(i, j, k, grid, u, v, w) = Σ₁₂²(i, j, k, grid, u, v)
@inline Σ₁₃²(i, j, k, grid, u, v, w) = Σ₁₃²(i, j, k, grid, u, w)
@inline Σ₂₃²(i, j, k, grid, u, v, w) = Σ₂₃²(i, j, k, grid, v, w)

#####
##### Same-location velocity products
#####

# ccc
@inline ∂x_u²(ijk...) = ∂x_u(ijk...)^2
@inline ∂y_v²(ijk...) = ∂y_v(ijk...)^2
@inline ∂z_w²(ijk...) = ∂z_w(ijk...)^2

# ffc
@inline ∂x_v²(ijk...) = ∂x_v(ijk...)^2
@inline ∂y_u²(ijk...) = ∂y_u(ijk...)^2

@inline ∂x_v_Σ₁₂(ijk...) = ∂x_v(ijk...) * Σ₁₂(ijk...)
@inline ∂y_u_Σ₁₂(ijk...) = ∂y_u(ijk...) * Σ₁₂(ijk...)

# fcf
@inline ∂z_u²(ijk...) = ∂z_u(ijk...)^2
@inline ∂x_w²(ijk...) = ∂x_w(ijk...)^2

@inline ∂x_w_Σ₁₃(ijk...) = ∂x_w(ijk...) * Σ₁₃(ijk...)
@inline ∂z_u_Σ₁₃(ijk...) = ∂z_u(ijk...) * Σ₁₃(ijk...)

# cff
@inline ∂z_v²(ijk...) = ∂z_v(ijk...)^2
@inline ∂y_w²(ijk...) = ∂y_w(ijk...)^2
@inline ∂z_v_Σ₂₃(ijk...) = ∂z_v(ijk...) * Σ₂₃(ijk...)
@inline ∂y_w_Σ₂₃(ijk...) = ∂y_w(ijk...) * Σ₂₃(ijk...)

#####
##### Tracer gradients
#####

@inline ∂x_c²(ijk...) = ∂xᶠᵃᵃ(ijk...)^2
@inline ∂y_c²(ijk...) = ∂yᵃᶠᵃ(ijk...)^2
@inline ∂z_c²(ijk...) = ∂zᵃᵃᶠ(ijk...)^2

#####
##### Normalized gradients
#####

# ccc
const norm_∂x_u = ∂x_u
const norm_∂y_v = ∂y_v
const norm_∂z_w = ∂z_w

# ffc
@inline norm_∂x_v(i, j, k, grid, v) =
    Δᶠx_ffc(i, j, k, grid) / Δᶠy_ffc(i, j, k, grid) * ∂xᶠᵃᵃ(i, j, k, grid, v)

@inline norm_∂y_u(i, j, k, grid, u) =
    Δᶠy_ffc(i, j, k, grid) / Δᶠx_ffc(i, j, k, grid) * ∂yᵃᶠᵃ(i, j, k, grid, u)

# fcf
@inline norm_∂x_w(i, j, k, grid, w) =
    Δᶠx_fcf(i, j, k, grid) / Δᶠz_fcf(i, j, k, grid) * ∂xᶠᵃᵃ(i, j, k, grid, w)

@inline norm_∂z_u(i, j, k, grid, u) =
    Δᶠz_fcf(i, j, k, grid) / Δᶠx_fcf(i, j, k, grid) * ∂zᵃᵃᶠ(i, j, k, grid, u)

# cff
@inline norm_∂y_w(i, j, k, grid, w) =
    Δᶠy_cff(i, j, k, grid) / Δᶠz_cff(i, j, k, grid) * ∂yᵃᶠᵃ(i, j, k, grid, w)

@inline norm_∂z_v(i, j, k, grid, v) =
    Δᶠz_cff(i, j, k, grid) / Δᶠy_cff(i, j, k, grid) * ∂zᵃᵃᶠ(i, j, k, grid, v)

# tracers
@inline norm_∂x_c(i, j, k, grid, c) = Δᶠx_fcc(i, j, k, grid) * ∂xᶠᵃᵃ(i, j, k, grid, c)
@inline norm_∂y_c(i, j, k, grid, c) = Δᶠy_cfc(i, j, k, grid) * ∂yᵃᶠᵃ(i, j, k, grid, c)
@inline norm_∂z_c(i, j, k, grid, c) = Δᶠz_ccf(i, j, k, grid) * ∂zᵃᵃᶠ(i, j, k, grid, c)

#####
##### Strain operators
#####

# ccc
@inline norm_Σ₁₁(i, j, k, grid, u) = norm_∂x_u(i, j, k, grid, u)
@inline norm_Σ₂₂(i, j, k, grid, v) = norm_∂y_v(i, j, k, grid, v)
@inline norm_Σ₃₃(i, j, k, grid, w) = norm_∂z_w(i, j, k, grid, w)

@inline norm_tr_Σ(i, j, k, grid, u, v, w) =
    norm_Σ₁₁(i, j, k, grid, u) + norm_Σ₂₂(i, j, k, grid, v) + norm_Σ₃₃(i, j, k, grid, w)

# ffc
@inline norm_Σ₁₂(i, j, k, grid::AbstractGrid{T}, u, v) where T =
    T(0.5) * (norm_∂y_u(i, j, k, grid, u) + norm_∂x_v(i, j, k, grid, v))

# fcf
@inline norm_Σ₁₃(i, j, k, grid::AbstractGrid{T}, u, w) where T =
    T(0.5) * (norm_∂z_u(i, j, k, grid, u) + norm_∂x_w(i, j, k, grid, w))

# cff
@inline norm_Σ₂₃(i, j, k, grid::AbstractGrid{T}, v, w) where T =
    T(0.5) * (norm_∂z_v(i, j, k, grid, v) + norm_∂y_w(i, j, k, grid, w))

@inline norm_Σ₁₂²(i, j, k, grid, u, v) = norm_Σ₁₂(i, j, k, grid, u, v)^2
@inline norm_Σ₁₃²(i, j, k, grid, u, w) = norm_Σ₁₃(i, j, k, grid, u, w)^2
@inline norm_Σ₂₃²(i, j, k, grid, v, w) = norm_Σ₂₃(i, j, k, grid, v, w)^2

# Consistent function signatures for convenience:
@inline norm_∂x_v(i, j, k, grid, u, v, w) = norm_∂x_v(i, j, k, grid, v)
@inline norm_∂x_w(i, j, k, grid, u, v, w) = norm_∂x_w(i, j, k, grid, w)

@inline norm_∂y_u(i, j, k, grid, u, v, w) = norm_∂y_u(i, j, k, grid, u)
@inline norm_∂y_w(i, j, k, grid, u, v, w) = norm_∂y_w(i, j, k, grid, w)

@inline norm_∂z_u(i, j, k, grid, u, v, w) = norm_∂z_u(i, j, k, grid, u)
@inline norm_∂z_v(i, j, k, grid, u, v, w) = norm_∂z_v(i, j, k, grid, v)

@inline norm_Σ₁₁(i, j, k, grid, u, v, w) = norm_Σ₁₁(i, j, k, grid, u)
@inline norm_Σ₂₂(i, j, k, grid, u, v, w) = norm_Σ₂₂(i, j, k, grid, v)
@inline norm_Σ₃₃(i, j, k, grid, u, v, w) = norm_Σ₃₃(i, j, k, grid, w)

@inline norm_Σ₁₂(i, j, k, grid, u, v, w) = norm_Σ₁₂(i, j, k, grid, u, v)
@inline norm_Σ₁₃(i, j, k, grid, u, v, w) = norm_Σ₁₃(i, j, k, grid, u, w)
@inline norm_Σ₂₃(i, j, k, grid, u, v, w) = norm_Σ₂₃(i, j, k, grid, v, w)

# Symmetry relations
const norm_Σ₂₁ = norm_Σ₁₂
const norm_Σ₃₁ = norm_Σ₁₃
const norm_Σ₃₂ = norm_Σ₂₃

# Trace and squared strains
@inline norm_tr_Σ²(ijk...) = norm_Σ₁₁(ijk...)^2 +  norm_Σ₂₂(ijk...)^2 +  norm_Σ₃₃(ijk...)^2

@inline norm_Σ₁₂²(i, j, k, grid, u, v, w) = norm_Σ₁₂²(i, j, k, grid, u, v)
@inline norm_Σ₁₃²(i, j, k, grid, u, v, w) = norm_Σ₁₃²(i, j, k, grid, u, w)
@inline norm_Σ₂₃²(i, j, k, grid, u, v, w) = norm_Σ₂₃²(i, j, k, grid, v, w)

#####
##### Same-location velocity products
#####

# ccc
@inline norm_∂x_u²(ijk...) = norm_∂x_u(ijk...)^2
@inline norm_∂y_v²(ijk...) = norm_∂y_v(ijk...)^2
@inline norm_∂z_w²(ijk...) = norm_∂z_w(ijk...)^2

# ffc
@inline norm_∂x_v²(ijk...) = norm_∂x_v(ijk...)^2
@inline norm_∂y_u²(ijk...) = norm_∂y_u(ijk...)^2

@inline norm_∂x_v_Σ₁₂(ijk...) = norm_∂x_v(ijk...) * norm_Σ₁₂(ijk...)
@inline norm_∂y_u_Σ₁₂(ijk...) = norm_∂y_u(ijk...) * norm_Σ₁₂(ijk...)

# fcf
@inline norm_∂z_u²(ijk...) = norm_∂z_u(ijk...)^2
@inline norm_∂x_w²(ijk...) = norm_∂x_w(ijk...)^2

@inline norm_∂x_w_Σ₁₃(ijk...) = norm_∂x_w(ijk...) * norm_Σ₁₃(ijk...)
@inline norm_∂z_u_Σ₁₃(ijk...) = norm_∂z_u(ijk...) * norm_Σ₁₃(ijk...)

# cff
@inline norm_∂z_v²(ijk...) = norm_∂z_v(ijk...)^2
@inline norm_∂y_w²(ijk...) = norm_∂y_w(ijk...)^2

@inline norm_∂z_v_Σ₂₃(ijk...) = norm_∂z_v(ijk...) * norm_Σ₂₃(ijk...)
@inline norm_∂y_w_Σ₂₃(ijk...) = norm_∂y_w(ijk...) * norm_Σ₂₃(ijk...)

#####
##### Tracer gradients squared
#####

@inline norm_∂x_c²(ijk...) = norm_∂x_c(ijk...)^2
@inline norm_∂y_c²(ijk...) = norm_∂y_c(ijk...)^2
@inline norm_∂z_c²(ijk...) = norm_∂z_c(ijk...)^2


text:
using Pkg;

Pkg.add("JUDI")
Pkg.add("SegyIO")
Pkg.add(PackageSpec(url="https://github.com/tierra-colada/ColadaSeismic.jl.git"))

max_stars_repo_name: tierra-colada/ColadaSeismic

max_issues_repo_name: tierra-colada/ColadaSeismic

max_forks_repo_name: tierra-colada/ColadaSeismic


text:
using TintiNet
using Test

@testset "TintiNet.jl" begin
    # Write your tests here.
end


text:
# Raw memory management
# Copied from CUDAdrv.jl

# NOTE: Mostly implemented for compatibility and ease-of-use, however we're
# currently not nearly as featureful as the original, which should change
# where possible.

export Mem

module Mem

using ..HSARuntime
using ..HSARuntime.HSA
import HSARuntime: check, get_region, AGENTS

## buffer type

struct Buffer
    ptr::Ptr{Cvoid}
    bytesize::Int
    agent::HSAAgent
    coherent::Bool
end

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

function view(buf::Buffer, bytes::Int)
    bytes > buf.bytesize && throw(BoundsError(buf, bytes))
    return Buffer(buf.ptr+bytes, buf.bytesize-bytes, buf.agent)
end

## refcounting

const refcounts = Dict{Buffer, Int}()

function refcount(buf::Buffer)
    get(refcounts, Base.unsafe_convert(Ptr{Cvoid}, buf), 0)
end

"""
    retain(buf)

Increase the refcount of a buffer.
"""
function retain(buf::Buffer)
    refcount = get!(refcounts, buf, 0)
    refcounts[buf] = refcount + 1
    return
end

"""
    release(buf)

Decrease the refcount of a buffer. Returns `true` if the refcount has dropped to 0, and
some action needs to be taken.
"""
function release(buf::Buffer)
    haskey(refcounts, buf) || error("Release of unmanaged $buf")
    refcount = refcounts[buf]
    @assert refcount > 0 "Release of dead $buf"
    refcounts[buf] = refcount - 1
    return refcount==1
end

## memory info

"""
    info()

Returns a tuple of two integers, indicating respectively the free and total amount of memory
(in bytes) available for allocation on the agent.
"""
function info()
    free_ref = Ref{Csize_t}()
    total_ref = Ref{Csize_t}()
    # FIXME: I'm not sure HSA has an API for this...
    return convert(Int, free_ref[]), convert(Int, total_ref[])
end

"""
    free()

Returns the free amount of memory (in bytes), available for allocation on the agent.
"""
free() = info()[1]

"""
    total()

Returns the total amount of memory (in bytes), available for allocation on the agent.
"""
total() = info()[2]

"""
    used()

Returns the used amount of memory (in bytes), allocated on the agent.
"""
used() = total()-free()

"""
    pointerinfo(ptr::Ptr)
    pointerinfo(buf::Buffer)
    pointerinfo(a::Array)

Retrieve information about the allocation referenced by the given pointer.
"""
function pointerinfo(ptr::Ptr)
    ptrinfo = Ref{HSA.PointerInfo}()

    HSA.amd_pointer_info(Ptr{Cvoid}(ptr), ptrinfo, C_NULL, Ptr{UInt32}(C_NULL), C_NULL) |> check

    return ptrinfo[]
end
pointerinfo(buf::Buffer) = pointerinfo(buf.ptr)
pointerinfo(a::Array) = pointerinfo(pointer(a))

## Page-locking

"""
    lock(ptr::Ptr, bytesize::Integer, agent::HSAAgent)
    lock(ptr, bytesize)
    lock(a::Array, agent)
    lock(a)

Page-lock a host pointer allocated by the OS allocator and return a new pointer from
the given `agent`. For more information, see `hsa_amd_memory_lock()`.

See also: [`unlock`](@ref)
"""
function lock(ptr::Ptr, bytesize::Integer, agent::HSAAgent)
    plocked = Ref{Ptr{Cvoid}}()
    HSA.amd_memory_lock(Ptr{Cvoid}(ptr), bytesize, Ref(agent.agent), 1, plocked) |> check
    return plocked[]
end
lock(ptr, bytesize) = lock(ptr, bytesize, get_default_agent())
lock(a::Array, agent::HSAAgent) = lock(pointer(a), sizeof(a), agent)
lock(a::Array) = lock(pointer(a), sizeof(a), get_default_agent())

"""
    unlock(ptr::Ptr)
    unlock(a::Array)

Unlock the host pointer previously page-locked with [`lock`](@ref).
NB: `ptr` should be the original locked host pointer and not the pointer returned by `lock`!
"""
function unlock(ptr::Ptr)
    HSA.amd_memory_unlock(Ptr{Cvoid}(ptr)) |> check
end
unlock(a::Array) = unlock(pointer(a))

## generic interface (for documentation purposes)

"""
Allocate linear memory on the device and return a buffer to the allocated memory. The
allocated memory is suitably aligned for any kind of variable. The memory will not be freed
automatically, use [`free(::Buffer)`](@ref) for that.
"""
function alloc end

"""
Free device memory.
"""
function free end

"""
Initialize device memory with a repeating value.
"""
function set! end

"""
Upload memory from host to device.
Executed asynchronously on `queue` if `async` is true.
"""
function upload end
@doc (@doc upload) upload!

"""
Download memory from device to host.
Executed asynchronously on `queue` if `async` is true.
"""
function download end
@doc (@doc download) download!

"""
Transfer memory from device to device.
Executed asynchronously on `queue` if `async` is true.
"""
function transfer end
@doc (@doc transfer) transfer!

## pointer-based

"""
    alloc(bytes::Integer; coherent=false)
    alloc(agent::HSAAgent, bytes::Integer; coherent=false)

Allocate `bytesize` bytes of HSA-managed memory.
Allocations are not coherent by default, meaning that the allocated buffer is
only accessible from the given agent (the default agent if not specified).

If `coherent` is set to `true`, the allocated buffer will be accessible from
all HSA agents, including the host CPU.
This, even though convenient, can sometimes be slower than explicit memory transfers.
"""
function alloc(agent::HSAAgent, bytesize::Integer; coherent=false)
    bytesize == 0 && return Buffer(C_NULL, 0, agent, coherent)

    ptr_ref = Ref{Ptr{Cvoid}}()

    region_kind = if coherent
        :finegrained
    else
        :coarsegrained
    end

    region = get_region(agent, region_kind)
    check(HSA.memory_allocate(region[], bytesize, ptr_ref))
    return Buffer(ptr_ref[], bytesize, agent, coherent)
end
alloc(bytesize; kwargs...) = alloc(get_default_agent(), bytesize; kwargs...)

function free(buf::Buffer)
    if buf.ptr != C_NULL
        check(HSA.memory_free(buf.ptr))
    end
    return
end

"""
    set!(buf::Buffer, value::UInt32, len::Integer)

Write `len` copies of the 32-bit `value` at the start of `buf`.
"""
function set!(buf::Buffer, value::UInt32, len::Integer)
    HSA.amd_memory_fill(buf.ptr, value, len) |> check
end

"""
    upload!(dst::Buffer, src, nbytes::Integer)

Upload `nbytes` memory from `src` at the host to `dst` on the device.
"""
function upload!(dst::Buffer, src::Ptr{T}, nbytes::Integer) where T
    # Page-locking coherent memory results in a ReadOnlyMemoryError.
    plocked = if dst.coherent
        src
    else
        lock(src, nbytes, dst.agent)
    end

    HSA.memory_copy(Ptr{T}(dst.ptr), Ptr{T}(plocked), nbytes) |> check

    if !dst.coherent
        unlock(src)
    end
end

"""
    download!(dst::Ref, src::Buffer, nbytes::Integer)

Download `nbytes` memory from `src` on the device to `dst` on the host.
"""
function download!(dst::Ptr{T}, src::Buffer, nbytes::Integer) where T
    plocked = if src.coherent
        dst
    else
        lock(dst, nbytes, src.agent)
    end

    HSA.memory_copy(Ptr{T}(plocked), Ptr{T}(src.ptr), nbytes) |> check

    if !src.coherent
        unlock(dst)
    end
end

"""
    transfer!(dst::Buffer, src::Buffer, nbytes::Integer)

Transfer `nbytes` of device memory from `src` to `dst`.
"""
function transfer!(dst::Buffer, src::Buffer, nbytes::Integer)
    HSA.memory_copy(dst.ptr, src.ptr, nbytes) |> check
end

## array based

"""
    alloc(src::AbstractArray; alloc_kwargs...)

Allocate space to store the contents of `src`.
"""
function alloc(src::AbstractArray; alloc_kwargs...)
    return alloc(sizeof(src); alloc_kwargs...)
end

"""
    upload!(dst::Buffer, src::AbstractArray)

Upload the contents of an array `src` to `dst`.
"""
function upload!(dst::Buffer, src::AbstractArray)
    GC.@preserve src upload!(dst, pointer(src), sizeof(src))
end

"""
    upload(src::AbstractArray; alloc_kwargs...)::Buffer

Allocates space for and uploads the contents of an array `src`, returning a Buffer.
For the allocation keywoard arguments see [`alloc`](@ref).
"""
function upload(src::AbstractArray; alloc_kwargs...)
    dst = alloc(src; alloc_kwargs...)
    upload!(dst, src)
    return dst
end

"""
    download!(dst::AbstractArray, src::Buffer)

Downloads memory from `src` to the array at `dst`. The amount of memory downloaded is
determined by calling `sizeof` on the array, so it needs to be properly preallocated.
"""
function download!(dst::AbstractArray, src::Buffer)
    GC.@preserve dst download!(pointer(dst), src, sizeof(dst))
    return
end

## type based

function check_type(::Type{Buffer}, T)
    if isa(T, UnionAll) || T.abstract || !isconcretetype(T)
        throw(ArgumentError("cannot represent abstract or non-leaf object"))
    end
    Base.datatype_pointerfree(T) || throw(ArgumentError("cannot handle non-ptrfree objects"))
    sizeof(T) == 0 && throw(ArgumentError("cannot represent singleton objects"))
end

"""
    alloc(T::Type, [count::Integer=1]; alloc_kwargs...)

Allocate space for `count` objects of type `T`.
"""
function alloc(::Type{T}, count::Integer=1; alloc_kwargs...) where {T}
    check_type(Buffer, T)

    return alloc(sizeof(T)*count; alloc_kwargs...)
end

"""
    download(::Type{T}, src::Buffer, [count::Integer=1])::Vector{T}

Download `count` objects of type `T` from the device at `src`, returning a vector.
"""
function download(::Type{T}, src::Buffer, count::Integer=1) where {T}
    dst = Vector{T}(undef, count)
    download!(dst, src)
    return dst
end

# Pretty-printing
function Base.show(io::IO, ptrinfo::HSA.PointerInfo)
    println(io, "Pointer type: $(ptrinfo.type)")
    println(io, "Owner: $(AGENTS[ptrinfo.agentOwner.handle])")
    println(io, "Agent base address: $(ptrinfo.agentBaseAddress)")
    println(io, "Host base address: $(ptrinfo.hostBaseAddress)")
    print(io, "Size (bytes): $(ptrinfo.sizeInBytes)")
end

end # module Mem


text:
using SFML

window = RenderWindow(VideoMode(800, 600), "Square", window_defaultstyle)
set_framerate_limit(window, 60)

event = Event()

square = RectangleShape()
set_size(square, Vector2f(40, 40))
set_position(square, Vector2f(400, 300))
set_origin(square, Vector2f(20, 20))
set_fillcolor(square, SFML.magenta)
set_outlinecolor(square, Color(4, 200, 127))
set_outline_thickness(square, 2)

while isopen(window)
    while pollevent(window, event)
        if get_type(event) == EventType.CLOSED
            close(window)
        end
    end
    clear(window, SFML.white)
    draw(window, square)
    display(window)
end


text:
## Author: Victor H. Aguiar

## questions: slack

## Import packages, make sure to install them.
## Julia version: v"1.5.2"
import JuMP
import Ipopt

## define a model
robinson=JuMP.Model(Ipopt.Optimizer)
#Knitro
# define variables
#coconuts
@JuMP.variable(robinson,c>=0)
#labor
@JuMP.variable(robinson,l>=0)
#wages
@JuMP.variable(robinson,w>=0)

# define parameters
# you must change this values to fit the dataset given to you
α1=.4
α2=(1-α1)
β = .1
A=50
Lbar=100.0

# Solve model
# each term of the summation is an equation of the system of equations
tax=0
@JuMP.NLobjective(robinson,Min,(α1/(α2)*(Lbar-l)*w-c)^2+(A*l^β*(1-tax)-c)^2+(A*β*l^(β-1)*(1-tax)-w)^2)
JuMP.optimize!(robinson)

JuMP.value(c)
JuMP.value(l)
JuMP.value(w)

taxvec=collect(0:.01:.99)
revenuevec=zeros(length(taxvec))
for t in 1:length(taxvec)
    tax=taxvec[t]
    @JuMP.NLobjective(robinson,Min,(α1/(α2)*(Lbar-l)*w-c)^2+(A*l^β*(1-tax)-c)^2+(A*β*l^(β-1)*(1-tax)-w)^2)
    JuMP.optimize!(robinson)
    revenuevec[t]=A*JuMP.value(l)^β*tax
end

sol=findmax(revenuevec)
taxvec[sol[2]]


text:
using Flux
using Transformers
using Transformers.Basic

struct TransformersModel
    embed::Union{WordPositionEmbed,CompositeEmbedding}
    encoder::Union{Encoder,Bert}
    decoder::Decoder
    generator::Generator
end

Flux.@functor TransformersModel
Flux.trainable(transformers::TransformersModel) = (transformers.embed, transformers.encoder, transformers.decoder, transformers.generator)

function (transformers::TransformersModel)(
    vocabulary::Vocabulary{T},
    inputs::Vector{T},
    outputs::Vector{T}
)::AbstractMatrix{<:AbstractFloat} where T
    # Encode inputs.
    input_indices = vocabulary(inputs) |> gpu
    if typeof(transformers.embed) <: CompositeEmbedding
        input_indices = (
            tok = input_indices,
            segment = fill(1, length(input_indices))
        )
    end
    input_embedding = transformers.embed(input_indices)
    encoded_embedding = transformers.encoder(input_embedding)

    # Decode outputs.
    output_indices = vocabulary(outputs) |> gpu
    if typeof(transformers.embed) <: CompositeEmbedding
        output_indices = (
            tok = output_indices,
            segment = fill(1, length(output_indices))
        )
    end
    output_embedding = transformers.embed(output_indices)
    decoded_embedding = transformers.decoder(output_embedding, encoded_embedding)

    # Calculate probabilities for next token.
    return transformers.generator(decoded_embedding)
end

function TransformersModel(
    embed::Union{WordPositionEmbed,CompositeEmbedding},
    size::Integer,
    head::Integer,
    hs::Integer,
    ps::Integer,
    layer::Integer,
    vocab_size::Integer;
    act=relu,
    pdrop=0.1,
    length_normalization::AbstractFloat=0.0
)
    return TransformersModel(
        embed,
        Encoder(size, head, hs, ps, layer; act=act,pdrop=pdrop),
        Decoder(size, head, hs, ps, layer; act=act,pdrop=pdrop),
        Generator(size, vocab_size)
    )
end


text:
function inex04()
#
#  test polynomial suggested by Goedecker
#
    inex(4)
end


text:
using ITF1788tests
using Test

ITF1788tests.regenerate_test

@testset "ITF1788tests.jl" begin

end