text:
const GEO_FILE = "cities1000"
const GEO_SOURCE = "http://download.geonames.org/export/dump"
const DATA_DIR = joinpath(dirname(dirname(pathof(ReverseGeocode))),"data")

"""
Column names in the geonames dumpfile (from http://download.geonames.org/export/dump):
"""
const COLUMNS = [:geonameid, :name, :asciiname, :alternatenames, :latitude, :longitude,
                :feature_class, :feature_code, :country_code, :cc2, :admin1_code, :admin2_code,
                :admin3_code, :admin4_code, :population, :elevation, :dem, :timezone, :modification_date]

"""
    Geocoder(;data_dir="./data", geo_file="cities1000"))

Geocoder structure that holds the reference points and their labels (city name and country code).
"""
struct Geocoder
    tree::NNTree
    info::Array{NamedTuple{(:city, :country_code),Tuple{String, String}}}
    country_codes::Dict{String, String}

    function Geocoder(;data_dir::String=DATA_DIR, geo_file::String=GEO_FILE)
        if ! isfile(joinpath(data_dir,"$geo_file.csv"))
            download_data(;data_dir=data_dir, geo_file=geo_file)
        end
        points, info = read_data(;data_dir=data_dir, geo_file=geo_file)

        tree = KDTree(points)
        country_codes = Dict(CSV.File(joinpath(data_dir,"country_codes.csv"); delim="\t", header=false))

        new(tree, info, country_codes)
    end
end

"""
    read_data(;data_dir="./data", geo_file="cities1000")

Load coordinates, country codes and city names from the `.csv` saved export of the geonames file.
Make sure to call `download_data()` before `read_data()`.
"""
function read_data(;data_dir::String=DATA_DIR, geo_file::String=GEO_FILE)
    data = CSV.File(joinpath(data_dir,"$geo_file.csv"); delim="\t", header=true, types=[String, Float64, Float64, String])

    n = length(data)
    points = Array{Float64}(undef, 2, n)
    info = Array{NamedTuple{(:city, :country_code),Tuple{String, String}}}(undef, n)
    for (i,row) in enumerate(data)
        points[:,i] .= row.latitude, row.longitude
        info[i] = (city = row.name, country_code = row.country_code)
    end
    points, info
end

"""
    download_data(;data_dir="./data", geo_file="cities1000", header=COLUMNS)

Download dump from [geonames.org](http://download.geonames.org/export/dump/). This function
fetches a file of cities with a population > 1000 (and seats of administrations of ceratain country subdivisions,
other options are population 500, 5000, 15000, see geonames.org for details).
The dump is unpacked and city name, coordinates and country code are saved
in a `.csv` file for use in the Geocoder.
"""
function download_data(;data_dir::String=DATA_DIR, geo_file::String=GEO_FILE, header=COLUMNS)
    # Download the source file
    download("$GEO_SOURCE/$geo_file.zip", joinpath(data_dir,"$geo_file.zip"))
    # extract the csv and drop unnecessary columns
    r = ZipFile.Reader(joinpath(data_dir,"$geo_file.zip"))
    data = CSV.File(read(r.files[1]); delim="\t", header=header, select=[:name,:latitude,:longitude,:country_code])
    close(r)
    # save needed data as csv
    CSV.write(joinpath(data_dir,"$geo_file.csv"), data; delim="\t")
    # clean up
    rm(joinpath(data_dir,"$geo_file.zip"))

    @info "Reference dataset sucessfuly saved in $data_dir."
end

"""
Download and resave the country codes csv from geonames.
Country codes are part of the package so this function does not usually need to run during install.
"""
function download_country_codes(;data_dir::String=DATA_DIR)
    download("http://download.geonames.org/export/dump/countryInfo.txt", "$data_dir/countryInfo.txt")
    country_info = CSV.File(joinpath(data_dir,"countryInfo.txt"); delim="\t", header=false, select=[1,5], datarow=51)
    country_codes = Dict([(c.Column1, c.Column5) for c in country_info])
    CSV.write(joinpath(data_dir,"country_codes.csv"), country_codes, delim="\t", header=false)
end

"""
    decode(gc::Geocoder, points) => Array{NamedTuple{(:country, :country_code, :city), Tuple{String, String, String}}}

Return country name and city for collection of points. Points should be either an array of staticaly sized arrays
(e.g. StaticArrays) or a Matrix (see NearestNeighbors.jl documentation for details).

The country and city is determined by the nearest neighbor search in the labelled list of city locations from geonames.org.
As such, the results may not be exactly accurate (e.g searches for points close to borders or in the middle of nowhere).

Nearest neighbor search uses the euclidian metric in the space of lat/lon coordinates.

```
# Example
julia> gc = Geocoder();
julia> ReverseGeocode.decode(gc, [SA[49.5863897, 17.2627342], SA[63.3342550, 12.0280064]])
2-element Array{Tuple{String,String},1}:
 (country="Czechia", country_code="CZ", city="Olomouc")
 (country="Norway", country_code="NO", city="Meråker")
```
"""
function decode(gc::Geocoder, points::Union{AbstractVector{<:AbstractVector{<:Real}}, AbstractMatrix{<:Real}})::Vector{NamedTuple{(:country, :country_code, :city), Tuple{String, String, String}}}
    idxs, dist = nn(gc.tree, points)
    infos = [gc.info[idx] for idx in idxs]
    tags = [(country=gc.country_codes[i.country_code], country_code=i.country_code, city=i.city) for i in infos]
end

"""
    decode(gc::Geocoder, point) => NamedTuple{(:country, :country_code, :city), Tuple{String, String, String}}

Decode for single point. If processing many points, preferably use `decode(gc, points)` instead of using this method in a loop.
"""
function decode(gc::Geocoder, point::AbstractVector{<:Real})
    decode(gc, [point])[1]
end


text:
using COPInstances
using Test

@testset "COPInstances.jl" begin
    # Write your tests here.
    @test 1 == 1
end

ext: jl


text:
# RUN: julia -e "import Brutus; Brutus.lit(:emit_translated)" --startup-file=no %s 2>&1 | FileCheck %s

@noinline add(x, y) = x + y
function caller(x, y)
    add(x, y)
end
emit(caller, Int64, Int64)

# CHECK:   func nested @"Tuple{typeof(Main.caller), Int64, Int64}"(%arg0: !jlir<"typeof(Main.caller)">, %arg1: !jlir.Int64, %arg2: !jlir.Int64) -> !jlir.Int64 attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     "jlir.goto"()[^bb1] : () -> ()
# CHECK-NEXT:   ^bb1:  // pred: ^bb0
# CHECK-NEXT:     %0 = "jlir.constant"() {value = #jlir.Main.add} : () -> !jlir<"typeof(Main.add)">
# CHECK-NEXT:     %1 = "jlir.invoke"(%0, %arg1, %arg2) {methodInstance = #jlir<"add(Int64, Int64) from add(Any, Any)">} : (!jlir<"typeof(Main.add)">, !jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     "jlir.return"(%1) : (!jlir.Int64) -> ()
# CHECK-NEXT:   }
# CHECK-NEXT: }

max_stars_repo_name: femtomc/brutus

max_issues_repo_name: femtomc/brutus

max_forks_repo_name: femtomc/brutus


text:
using Documenter, ModelingToolkit

makedocs(
    sitename="ModelingToolkit.jl",
    authors="Chris Rackauckas",
    modules=[ModelingToolkit],
    clean=true,doctest=false,
    format = Documenter.HTML(#analytics = "UA-90474609-3",
                             assets = ["assets/favicon.ico"],
                             canonical="https://mtk.sciml.ai/stable/"),
    pages=[
        "Home" => "index.md",
        "Tutorials" => Any[
            "tutorials/symbolic_functions.md",
            "tutorials/ode_modeling.md",
            "tutorials/higher_order.md",
            "tutorials/nonlinear.md",
            "tutorials/modelingtoolkitize.md",
            "tutorials/auto_parallel.md",
            "tutorials/converting_to_C.md"
        ],
        "Systems" => Any[
            "systems/AbstractSystem.md",
            "systems/ODESystem.md",
            "systems/SDESystem.md",
            "systems/JumpSystem.md",
            "systems/NonlinearSystem.md",
            "systems/OptimizationSystem.md",
            "systems/ControlSystem.md",
            "systems/ReactionSystem.md",
            "systems/PDESystem.md",
            "systems/DependencyGraphs.md"
        ],
        "Comparison Against SymPy" => "comparison.md",
        "highlevel.md",
        "build_function.md",
        "IR.md"
    ]
)

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


text:
function reflection_coefficient(ω::T, amps::AverageWave{T}, medium::Medium, specie::Specie;
        θin::T = 0.0, scheme::Symbol = :trapezoidal) where T <: AbstractFloat

    k = ω/medium.c
    M = amps.hankel_order
    σ = k .* integration_scheme(amps.x; scheme=scheme) # multiple by k to be the same as the non-dimensional version

    R = T(2)*specie.num_density/(cos(θin)*k^2)*sum(
        im^T(m)*exp(im*k*amps.x[j]*cos(θin) - im*θin*m)*amps.amplitudes[j,m+M+1,1]*σ[j]
    for m=-M:M, j in eachindex(amps.x))

    return R
end


text:
# support for single-process Jus apps, where the program is a client of itself

@kwdef struct ClientVar
    config
    id
    value
    metadata
end

"""
    local_server(namespace)

Start a server running in a new task.
Returns (config::Config, task::Task, ws::FakeWS)
"""
function local_server(namespace::AbstractString)
    c = Config()
    ws = FakeWS()
    other_ws = other(ws)
    t = @task try
        serve(c, ws)
    catch err
        println("@@@@@@@@\n@@@@@@@@ CAUGHT SERVER ERROR")
        println(typeof(err))
        println(!isopen(other_ws))
        if err isa HTTP.WebSockets.WebSocketError && !isopen(other_ws)
            println("@@@@@@@@\n@@@@@@@@ OK CLOSE")
        else
            @error join(["$(err)", stacktrace(catch_backtrace())...], "\n")
        end
    end
    bind(ws.input, t)
    bind(ws.output, t)
    schedule(t)
    output(other_ws, (namespace, secret = ""))
    c, t, other_ws
end


text:
"""
General Interface for Pricing

		Price=pricer(mcProcess,rfCurve,mcBaseData,payoff_)

Where:\n
		mcProcess          = Process to be simulated.
		rfCurve  = Zero Rate Data.
		mcBaseData = Basic properties of MonteCarlo simulation
		payoff_ = Payoff(s) to be priced

		Price     = Price of the derivative

"""
function distribution(mcProcess::BaseProcess, rfCurve::AbstractZeroRateCurve, mcConfig::MonteCarloConfiguration, T::num_) where {num_ <: Number}
    set_seed(mcConfig)
    S = simulate(mcProcess, rfCurve, mcConfig, T)

    return S[:, end]
end


text:
type DiscretizedSwap{ST <: SwapType} <: DiscretizedAsset
  nominal::Float64
  swapT::ST
  fixedResetTimes::Vector{Float64}
  fixedPayTimes::Vector{Float64}
  floatingResetTimes::Vector{Float64}
  floatingPayTimes::Vector{Float64}
  args::VanillaSwapArgs
  common::DiscretizedAssetCommon
end

function DiscretizedSwap{DC <: DayCount, ST <: SwapType}(nominal::Float64, swapT::ST, referenceDate::Date, dc::DC, fixedPayDates::Vector{Date}, fixedResetDates::Vector{Date},
                        floatingPayDates::Vector{Date}, floatingResetDates::Vector{Date}, args::VanillaSwapArgs)
  fixed_n = length(fixedPayDates)
  float_n = length(floatingPayDates)

  fixedResetTimes = zeros(fixed_n)
  fixedPayTimes = zeros(fixed_n)
  floatingResetTimes = zeros(float_n)
  floatingPayTimes = zeros(float_n)

  for i = 1:fixed_n
    fixedResetTimes[i] = year_fraction(dc, referenceDate, fixedResetDates[i])
    fixedPayTimes[i] = year_fraction(dc, referenceDate, fixedPayDates[i])
  end

  for i = 1:float_n
    floatingResetTimes[i] = year_fraction(dc, referenceDate, floatingResetDates[i])
    floatingPayTimes[i] = year_fraction(dc, referenceDate, floatingPayDates[i])
  end

  DiscretizedSwap(nominal, swapT, fixedResetTimes, fixedPayTimes, floatingResetTimes, floatingPayTimes, args, DiscretizedAssetCommon())
end

# methods #
function mandatory_times(discretizedSwap::DiscretizedSwap)
  # get times
  times = vcat(discretizedSwap.fixedResetTimes[discretizedSwap.fixedResetTimes .>= 0.0], discretizedSwap.fixedPayTimes[discretizedSwap.fixedPayTimes .>= 0.0],
               discretizedSwap.floatingResetTimes[discretizedSwap.floatingResetTimes .>= 0.0], discretizedSwap.floatingPayTimes[discretizedSwap.floatingPayTimes .>= 0.0])

  return times
end

function reset!(dSwap::DiscretizedSwap, sz::Int)
  dSwap.common.values = zeros(sz)
  adjust_values!(dSwap)

  return dSwap
end

function pre_adjust_values_impl!(dSwap::DiscretizedSwap)
  # Floating payments
  @simd for i = 1:length(dSwap.floatingResetTimes)
    @inbounds t = dSwap.floatingResetTimes[i]
    if t >= 0.0 && is_on_time(dSwap, t)
      bond = DiscretizedDiscountBond()
      initialize!(bond, dSwap.common.method, dSwap.floatingPayTimes[i])
      rollback!(bond, dSwap.common.time)

      nominal = dSwap.nominal
      T = dSwap.args.floatingAccrualTimes[i]
      spread = dSwap.args.floatingSpreads[i]

      accruedSpread = nominal * T * spread
      for j = 1:length(dSwap.common.values)
        @inbounds coup = nominal * (1.0 - bond.common.values[j]) + accruedSpread * bond.common.values[j]

        if isa(dSwap.swapT, Payer)
          @inbounds dSwap.common.values[j] += coup
        else
          @inbounds dSwap.common.values[j] -= coup
        end
      end
    end
  end

  # Fixed Payments
  @simd for i = 1:length(dSwap.fixedResetTimes)
    t = dSwap.fixedResetTimes[i]
    if t >= 0.0 && is_on_time(dSwap, t)
      bond = DiscretizedDiscountBond()
      initialize!(bond, dSwap.common.method, dSwap.fixedPayTimes[i])
      rollback!(bond, dSwap.common.time)

      @inbounds fixedCoup = dSwap.args.fixedCoupons[i]

      for j = 1:length(dSwap.common.values)
        @inbounds coup = fixedCoup * bond.common.values[j]
        if isa(dSwap.swapT, Payer)
          @inbounds dSwap.common.values[j] -= coup
        else
          @inbounds dSwap.common.values[j] += coup
        end
      end
    end
  end

  return dSwap
end

function post_adjust_values_impl!(dSwap::DiscretizedSwap)
  # fixed coupons whose reset time is in the past won't be managed in pre_adjust_values
  @simd for i = 1:length(dSwap.fixedPayTimes)
    @inbounds t = dSwap.fixedPayTimes[i]
    @inbounds _reset = dSwap.fixedResetTimes[i]
    if t >= 0.0 && is_on_time(dSwap, t) && _reset < 0.0
      fixedCoup = dSwap.args.fixedCoupons[i]
      if isa(dSwap.swapT, Payer)
        dSwap.common.values -= fixedCoup
      else
        dSwap.common.values += fixedCoup
      end
    end
  end

  # the same applies to floating payments whose rate is already fixed
  @simd for i = 1:length(dSwap.floatingPayTimes)
    @inbounds t = dSwap.floatingPayTimes[i]
    @inbounds _reset = dSwap.floatingResetTimes[i]
    if t >= 0.0 && is_on_time(dSwap, t) && _reset < 0.0
      @inbounds currentFloatingCoup = dSwap.args.floatingCoupons[i]

      if isa(dSwap.swapT, Payer)
        dSwap.common.values += currentFloatingCoup
      else
        dSwap.common.values -= currentFloatingCoup
      end
    end
  end

  return dSwap
end


text:
@muladd function perform_step!(integrator,cache::PCEulerConstantCache)
  @unpack t,dt,uprev,u,W,p,f,g = integrator
  @unpack theta,eta,ggprime = integrator.alg
  dW = W.dW
  if !is_diagonal_noise(integrator.sol.prob) || typeof(W.dW) <: Number
    gdW_n = g(uprev,p,t)*dW
  else
    gdW_n = g(uprev,p,t).*dW
  end

  f_n = f(uprev,p,t)
  ubar = uprev + dt*f_n + gdW_n

  fbar_n = f_n - eta*ggprime(uprev, p, t)

  tnp1 = t + dt
  fbar_np1 = f(ubar, p, tnp1) - eta*ggprime(ubar, p, tnp1)

  if !is_diagonal_noise(integrator.sol.prob) || typeof(dW) <: Number
    gdW_np1 = g(ubar,p,tnp1)*dW
  else
    gdW_np1 = g(ubar,p,tnp1).*dW
  end

  u = uprev + (theta*dt)*fbar_np1 + ((1-theta)*dt)*fbar_n + eta*gdW_np1 + (1-eta)*gdW_n
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::PCEulerCache)
  @unpack utmp,ftmp,gtmp,gdWtmp, bbprimetmp = cache
  @unpack t,dt,uprev,u,W,p,f,g = integrator
  @unpack theta,eta,ggprime = integrator.alg
  dW = W.dW

  f(ftmp,uprev,p,t)
  g(gtmp,uprev,p,t)
  ggprime(bbprimetmp,uprev,p,t)

  if is_diagonal_noise(integrator.sol.prob)
    @.. gtmp *= dW
  else
    mul!(gdWtmp,gtmp,dW)
  end

  @.. utmp = uprev + ftmp*dt + gdWtmp
  @.. u = uprev + (1-theta)*(ftmp- eta*bbprimetmp)*dt + (1-eta)*gdWtmp

  tnp1 = t + dt
  f(ftmp,utmp,p,tnp1)
  g(gtmp,utmp,p,tnp1)
  ggprime(bbprimetmp,utmp,p,tnp1)

  if is_diagonal_noise(integrator.sol.prob)
    @.. gtmp *= dW
  else
    mul!(gdWtmp,gtmp,dW)
  end

  @.. u += theta*(ftmp - eta*bbprimetmp)*dt + eta*gdWtmp
end

max_stars_repo_name: JuliaDiffEq/StochasticDiffEq


text:
"""
	ARCHModel <: StatisticalModel
"""
abstract type ARCHModel <: StatisticalModel end

# this makes predict.(am, :variance, 1:3) work
Base.Broadcast.broadcastable(am::ARCHModel) = Ref(am)

"""
	VolatilitySpec{T}
Abstract supertype of UnivariateVolatilitySpec{T} and MultivariateVolatilitySpec{T} .
"""
abstract type VolatilitySpec{T} end

"""
    MeanSpec{T}
Abstract supertype that mean specifications inherit from.
"""
abstract type MeanSpec{T} end

struct NumParamError <: Exception
    expected::Int
    got::Int
end

function showerror(io::IO, e::NumParamError)
    print(io, "incorrect number of parameters: expected $(e.expected), got $(e.got).")
end

nobs(am::ARCHModel) = length(am.data)
islinear(am::ARCHModel) = false
isfitted(am::ARCHModel) = am.fitted

function confint(am::ARCHModel, level::Real=0.95)
    hcat(coef(am), coef(am)) .+ stderror(am)*quantile(Normal(),(1. -level)/2.)*[1. -1.]
end

score(am::ARCHModel) = sum(scores(am), dims=1)

function vcov(am::ARCHModel)
	S = scores(am)
	V = S'S/nobs(am)
    J = informationmatrix(am; expected=false) #Note: B&W use expected information.
    Ji = try
        inv(J)
    catch e
        if e isa LinearAlgebra.SingularException
            @warn "Fisher information is singular; vcov matrix is inaccurate."
            pinv(J)
        else
            rethrow(e)
        end
    end
    v = Ji*V*Ji/nobs(am) #Huber sandwich
    all(diag(v).>0) || @warn "non-positive variance encountered; vcov matrix is inaccurate."
    v
end

function show(io::IO, spec::VolatilitySpec)
    println(io, modname(typeof(spec)), " specification.\n\n", length(spec.coefs) > 0 ? CoefTable(spec.coefs, coefnames(typeof(spec)), ["Parameters:"]) : "No estimable parameters.")
end

stderror(am::ARCHModel) = sqrt.(abs.(diag(vcov(am))))

"""
    fit!(am::ARCHModel; algorithm=BFGS(), autodiff=:forward, kwargs...)

Fit the uni- or multivariate ARCHModel specified by `am`, modifying `am` in place.
Keyword arguments are passed on to the optimizer.
"""
function fit!(am::ARCHModel; kwargs...) end

"""
    fit(am::ARCHModel; algorithm=BFGS(), autodiff=:forward, kwargs...)

Fit the uni- or multivariate ARCHModel specified by `am` and return the result in a new instance of
`ARCHModel`. Keyword arguments are passed on to the optimizer.
"""
function fit(am::ARCHModel; kwargs...) end

"""
    simulate!(am::ARCHModel; warmup=100)
Simulate an ARCHModel, modifying `am` in place.
"""
function simulate! end

"""
    simulate(am::ARCHModel; warmup=100)
	simulate(am::ARCHModel, T; warmup=100)
    simulate(spec::UnivariateVolatilitySpec, T; warmup=100, dist=StdNormal(), meanspec=NoIntercept())
Simulate a length-T time series from a UnivariateARCHModel.
	simulate(spec::MultivariateVolatilitySpec, T; warmup=100, dist=MultivariateStdNormal(), meanspec=[NoIntercept() for i = 1:d])
Simulate a length-T time series from a MultivariateARCHModel.
"""
function simulate end

function simulate!(am::ARCHModel; warmup=100)
	am.fitted = false
    _simulate!(am.data, am.spec; warmup=warmup, dist=am.dist, meanspec=am.meanspec)
    am
end

function simulate(am::ARCHModel, nobs; warmup=100)
	am2 = deepcopy(am)
	simulate(am2.spec, nobs; warmup=warmup, dist=am2.dist, meanspec=am2.meanspec)
end

simulate(am::ARCHModel; warmup=100) = simulate(am, size(am.data)[1]; warmup=warmup)


text:
module VkExt

using VulkanCore
include("VkExt.VkPhysicalDeviceFeatures.jl")

const VK_TRUE = vk.VkBool32(vk.VK_TRUE)
const VK_FALSE = vk.VkBool32(vk.VK_FALSE)

# instance
function createInstance(info::vk.VkInstanceCreateInfo)::vk.VkInstance
    outInstance = Ref{vk.VkInstance}()
    err = vk.vkCreateInstance(Ref(info), C_NULL, outInstance)
    if err != vk.VK_SUCCESS
        error(err, "Failed to create instance!")
    end
    return outInstance[]
end

function enumerateInstanceExtensionProperties()::Tuple{Vector{vk.VkExtensionProperties}, UInt32}
    extensionCount = Ref{UInt32}(0)
    vk.vkEnumerateInstanceExtensionProperties(C_NULL, extensionCount, C_NULL)
    availableExtensions = Vector{vk.VkExtensionProperties}(undef, extensionCount[])
    vk.vkEnumerateInstanceExtensionProperties(C_NULL, extensionCount, availableExtensions)
    return availableExtensions, extensionCount[]
end

function enumerateInstanceLayerProperties()::Tuple{Vector{vk.VkLayerProperties}, UInt32}
    layerCount = Ref{UInt32}(0)
    vk.vkEnumerateInstanceLayerProperties(layerCount, C_NULL);
    avaliableLayers = Vector{vk.VkLayerProperties}(undef, layerCount[])
    vk.vkEnumerateInstanceLayerProperties(layerCount, avaliableLayers)
    return avaliableLayers, layerCount[]
end

function vkCreateDebugReportCallbackEXT(instance, callbackInfoRef, allocatorRef, callbackRef)
    fnptr = vk.vkGetInstanceProcAddr(instance, "vkCreateDebugReportCallbackEXT") |> vk.PFN_vkCreateDebugReportCallbackEXT
    ccall(fnptr, vk.VkResult, (vk.VkInstance, Ptr{vk.VkDebugReportCallbackCreateInfoEXT}, Ptr{vk.VkAllocationCallbacks},
                               Ptr{vk.VkDebugReportCallbackEXT}), instance, callbackInfoRef, allocatorRef, callbackRef)
end

function createDebugReportCallbackEXT(instance::vk.VkInstance, createInfo::vk.VkDebugReportCallbackCreateInfoEXT)::vk.VkDebugReportCallbackEXT
    callback = Ref{vk.VkDebugReportCallbackEXT}()
    err = vkCreateDebugReportCallbackEXT(instance, Ref(createInfo), C_NULL, callback)
    if err != vk.VK_SUCCESS
        error(err, "Create debug report callback ext failed!")
    end
    return callback[]
end

function destroyDebugReportCallbackEXT(instance::vk.VkInstance, callback::vk.VkDebugReportCallbackEXT)
    fnptr = vk.vkGetInstanceProcAddr(instance, "vkDestroyDebugReportCallbackEXT") |> vk.PFN_vkDestroyDebugReportCallbackEXT
    err = ccall(fnptr, vk.VkResult, (vk.VkInstance, vk.VkDebugReportCallbackEXT, Ptr{vk.VkAllocationCallbacks}),
                                     instance, callback, C_NULL)
    if err != vk.VK_SUCCESS
       error(err, "Destroy debug report callback ext failed!")
    end
end
# physical device
function enumeratePhysicalDevices(vkInstance::vk.VkInstance)
    physicalDeviceCount = Ref{Cuint}(0)
    vk.vkEnumeratePhysicalDevices(vkInstance, physicalDeviceCount, C_NULL)
    if (physicalDeviceCount[] == 0)
        error("Failed to find GPUs with Vulkan support!")
    end
    physicalDevices = Vector{vk.VkPhysicalDevice}(undef, physicalDeviceCount[])
    vk.vkEnumeratePhysicalDevices(vkInstance, physicalDeviceCount, physicalDevices)
    return physicalDevices
end

function enumeratePhysicalDeviceGroups(vkInstance::vk.VkInstance)
    groupCount = Ref{Cuint}(0)
    vk.vkEnumeratePhysicalDeviceGroups(vkInstance, groupCount, C_NULL)
    if (groupCount[] == 0)
        error("Failed to find GPU groups with Vulkan support!")
    end
    groupProps = Vector{vk.VkPhysicalDeviceGroupProperties}(undef, groupCount[])
    vk.vkEnumeratePhysicalDeviceGroups(vkInstance, groupCount, groupProps)
    return groupProps
end

function getProperties(phyDevice::vk.VkPhysicalDevice)::vk.VkPhysicalDeviceProperties
    deviceProperties = Ref{vk.VkPhysicalDeviceProperties}()
    vk.vkGetPhysicalDeviceProperties(phyDevice, deviceProperties)
    return deviceProperties[]
end

function enumerateDeviceExtensionProperties(phyDevice::vk.VkPhysicalDevice)
    extensionCount = Ref{UInt32}(0)
    vk.vkEnumerateDeviceExtensionProperties(phyDevice, C_NULL, extensionCount, C_NULL)

    availableExtensions = Vector{vk.VkExtensionProperties}(undef, extensionCount[])
    vk.vkEnumerateDeviceExtensionProperties(phyDevice, C_NULL, extensionCount, availableExtensions)
    return availableExtensions
end

function getQueueFamilyProperties(phyDevice::vk.VkPhysicalDevice)::Vector{vk.VkQueueFamilyProperties}
    queueFamilyCount = Ref{Cuint}(0)
    vk.vkGetPhysicalDeviceQueueFamilyProperties(phyDevice, queueFamilyCount, C_NULL)

    queueFamilies = Vector{vk.VkQueueFamilyProperties}(undef, queueFamilyCount[])
    vk.vkGetPhysicalDeviceQueueFamilyProperties(phyDevice, queueFamilyCount, queueFamilies)
    return queueFamilies
end

function createDevice(phyDevice::vk.VkPhysicalDevice, createInfo::Ref{vk.VkDeviceCreateInfo})::vk.VkDevice
    logicalDevice = Ref{vk.VkDevice}()
    err = vk.vkCreateDevice(phyDevice, createInfo, C_NULL, logicalDevice)
    if err != vk.VK_SUCCESS
         error(err, "failed to create logical device!")
    end
    return logicalDevice[]
end

function getSurfaceSupportKHR(phyDevice::vk.VkPhysicalDevice, queueFamilyIndex::UInt32, surface::vk.VkSurfaceKHR)
    presentSupport = Ref{vk.VkBool32}(false)
    err = vk.vkGetPhysicalDeviceSurfaceSupportKHR(phyDevice, queueFamilyIndex, surface, presentSupport)
    if err != vk.VK_SUCCESS
        error(err, "failed to get physical device surface support!")
    end
    return presentSupport[]
end

function getSurfaceFormatsKHR(phyDevice::vk.VkPhysicalDevice, surface::vk.VkSurfaceKHR)::Vector{vk.VkSurfaceFormatKHR}
    surfaceFormats = Vector{vk.VkSurfaceFormatKHR}()
    formatCount = Ref{UInt32}()
    vk.vkGetPhysicalDeviceSurfaceFormatsKHR(phyDevice, surface, formatCount, C_NULL)
    if (formatCount[] != 0)
        resize!(surfaceFormats, formatCount[])
        vk.vkGetPhysicalDeviceSurfaceFormatsKHR(phyDevice, surface, formatCount, surfaceFormats)
    end
    return surfaceFormats
end

function getSurfacePresentModesKHR(phyDevice::vk.VkPhysicalDevice, surface::vk.VkSurfaceKHR)::Vector{vk.VkPresentModeKHR}
    presentModes = Vector{vk.VkPresentModeKHR}(undef, 1)
    presentModeCount = Ref{UInt32}()
    vk.vkGetPhysicalDeviceSurfacePresentModesKHR(phyDevice, surface, presentModeCount, C_NULL)
    if (presentModeCount != 0)
        resize!(presentModes, presentModeCount[])
        vk.vkGetPhysicalDeviceSurfacePresentModesKHR(phyDevice, surface, presentModeCount, presentModes)
    end
    return presentModes
end

function getMemoryProperties(phyDevice::vk.VkPhysicalDevice)::vk.VkPhysicalDeviceMemoryProperties
    memoryProperties = Ref{vk.VkPhysicalDeviceMemoryProperties}()
    vk.vkGetPhysicalDeviceMemoryProperties(phyDevice, memoryProperties)
    return memoryProperties[]
end

# logical device
function createShaderModule(logicalDevice::vk.VkDevice, code::Ptr{UInt8}, codeSize::Int64)::vk.VkShaderModule
    createInfo = Ref(vk.VkShaderModuleCreateInfo(
        vk.VK_STRUCTURE_TYPE_SHADER_MODULE_CREATE_INFO, #sType::VkStructureType
        C_NULL, #pNext::Ptr{Cvoid}
        0, #flags::VkShaderModuleCreateFlags
        codeSize, #codeSize::Csize_t
        code, #pCode::Ptr{UInt32}
    ))

    shaderModule = Ref{vk.VkShaderModule}()
    if (vk.vkCreateShaderModule(logicalDevice, createInfo, C_NULL, shaderModule) != vk.VK_SUCCESS)
        error("Failed to create shader module!")
    end
    return shaderModule[]
end

function createCommandPool(logicalDevice::vk.VkDevice, createInfo::vk.VkCommandPoolCreateInfo)::vk.VkCommandPool
    commandPool = Ref{vk.VkCommandPool}()
    if (vk.vkCreateCommandPool(logicalDevice, Ref(createInfo), C_NULL, commandPool) != vk.VK_SUCCESS)
        error("Failed to create command pool!")
    end
    return commandPool[]
end

function createSemaphore(logicalDevice::vk.VkDevice,
                            createInfo::vk.VkSemaphoreCreateInfo = vk.VkSemaphoreCreateInfo(
                                                                        vk.VK_STRUCTURE_TYPE_SEMAPHORE_CREATE_INFO, #sType::VkStructureType
                                                                        C_NULL, #pNext::Ptr{Cvoid}
                                                                        0 #flags::VkSemaphoreCreateFlags
                                                                    ))
    semaphore = Ref{vk.VkSemaphore}()
    if (vk.vkCreateSemaphore(logicalDevice, Ref(createInfo), C_NULL, semaphore) != vk.VK_SUCCESS)
        error("Failed to create semaphore!")
    end
    return semaphore[]
end

function getSwapchainImagesKHR(logicalDevice::vk.VkDevice, swapchain::vk.VkSwapchainKHR)::Vector{vk.VkImage}
    imageCount = Ref{UInt32}()
    swapChainImages = Vector{vk.VkImage}()
    vk.vkGetSwapchainImagesKHR(logicalDevice, swapchain, imageCount, C_NULL)
    resize!(swapChainImages, imageCount[])
    vk.vkGetSwapchainImagesKHR(logicalDevice, swapchain, imageCount, swapChainImages)
    return swapChainImages
end

function allocateMemory(logicalDevice::vk.VkDevice, allocateInfo::Ref{vk.VkMemoryAllocateInfo})::vk.VkDeviceMemory
    memory = Ref{vk.VkDeviceMemory}()
    vk.vkAllocateMemory(logicalDevice, allocateInfo, C_NULL, memory)
    return memory[]
end

function getImageMemoryRequirements(logicalDevice::vk.VkDevice, image::vk.VkImage)::vk.VkMemoryRequirements
    memoryRequirements = Ref{vk.VkMemoryRequirements}()
    vk.vkGetImageMemoryRequirements(logicalDevice, image, memoryRequirements)
    return memoryRequirements[]
end

function getBufferMemoryRequirements(logicalDevice::vk.VkDevice, buffer::vk.VkBuffer)::vk.VkMemoryRequirements
    memoryRequirements = Ref{vk.VkMemoryRequirements}()
    vk.vkGetBufferMemoryRequirements(logicalDevice, buffer, memoryRequirements)
    return memoryRequirements[]
end

function createBuffer(logicalDevice::vk.VkDevice, createInfo::Ref{vk.VkBufferCreateInfo})::vk.VkBuffer
    buffer = Ref{vk.VkBuffer}()
    if vk.vkCreateBuffer(logicalDevice, createInfo, C_NULL, buffer) != vk.VK_SUCCESS
        error("Failed to create buffer!")
    end
    return buffer[]
end

function mapMemory(logicalDevice::vk.VkDevice, memory::vk.VkDeviceMemory, offset::vk.VkDeviceSize, size::vk.VkDeviceSize, flags::vk.VkMemoryMapFlags)::Ptr{Cvoid}
    pData = Ref{Ptr{Cvoid}}()
    if vk.vkMapMemory(logicalDevice, memory, offset, size, flags, pData) != vk.VK_SUCCESS
        error("Failed to map memory!")
    end
    return pData[]
end

function waitForFences(logicalDevice::vk.VkDevice, fences::Vector{vk.VkFence}, waitAll::vk.VkBool32, timeout::UInt64)
    if vk.vkWaitForFences(logicalDevice, length(fences), pointer(fences), waitAll, timeout) != vk.VK_SUCCESS
        error("Failed to wait for fences!")
    end
end

function createFence(logicalDevice::vk.VkDevice)::vk.VkFence
    fenceInfo = Ref(vk.VkFenceCreateInfo(
        vk.VK_STRUCTURE_TYPE_FENCE_CREATE_INFO, #sType::VkStructureType
        C_NULL, #pNext::Ptr{Cvoid}
        0 #flags::VkFenceCreateFlags
    ))
    outFence = Ref{vk.VkFence}()
    err = vk.vkCreateFence(logicalDevice, fenceInfo, C_NULL, outFence)
    if err == vk.VK_ERROR_OUT_OF_HOST_MEMORY
        error("Queue::findFreeFence(): Ran out of memory trying to allocate a fence for CommandBuffer submission. \n" *
              "Use Queue::catchUp to prevent oversized backlogs of submissions.")
    elseif err != vk.VK_SUCCESS
        error("Failed to create fence!")
    end
    return outFence[]
end

# common
function ClearValue(color::vk.VkClearColorValue)
    return vk.VkClearValue(color)
end

function ClearValue(depthStencil::vk.VkClearDepthStencilValue)
    vk.VkClearValue(vk.VkClearColorValue((depthStencil.depth, depthStencil.stencil, 0, 0)))
end

end #module


text:
function findRequiredIdxs(A::LagrangeInterpolation, t)
  idxs = sortperm(A.t,by=x->abs(t-x))
  idxs[1:(A.n+1)]
end

function spline_coefficients(n, d, k, u::Number)
  N = zeros(n)
  if u == k[1]
    N[1] = one(u)
  elseif u == k[end]
    N[end] = one(u)
  else
    i = findfirst(x->x>u,k) - 1
    N[i] = one(u)
    for deg = 1:d
      N[i-deg] = (k[i+1]-u)/(k[i+1]-k[i-deg+1]) * N[i-deg+1]
      for j = (i-deg+1):(i-1)
        N[j] = (u-k[j])/(k[j+deg]-k[j]) * N[j] + (k[j+deg+1]-u)/(k[j+deg+1]-k[j+1]) * N[j+1]
      end
      N[i] = (u-k[i])/(k[i+deg]-k[i]) * N[i]
    end
  end
  N
end

function spline_coefficients(n, d, k, u::AbstractVector)
  N = zeros(eltype(u),n,n)
  for i = 1:n
    N[i,:] .= spline_coefficients(n,d,k,u[i])
  end
  N
end

# helper function for data manipulation
function munge_data(u::AbstractVector{<:Real}, t::AbstractVector{<:Real})
  return u, t
end

function munge_data(u::AbstractVector, t::AbstractVector)
  Tu = Base.nonmissingtype(eltype(u))
  Tt = Base.nonmissingtype(eltype(t))
  @assert length(t) == length(u)
  non_missing_indices = collect(i for i in 1:length(t) if !ismissing(u[i]) && !ismissing(t[i]))
  newu = Tu.([u[i] for i in non_missing_indices])
  newt = Tt.([t[i] for i in non_missing_indices])

  return newu, newt
end

function munge_data(U::StridedMatrix, t::AbstractVector)
  TU = Base.nonmissingtype(eltype(U))
  Tt = Base.nonmissingtype(eltype(t))
  @assert length(t) == size(U,2)
  non_missing_indices = collect(i for i in 1:length(t) if !any(ismissing,U[:,i]) && !ismissing(t[i]))
  newUs = [TU.(U[:,i]) for i in non_missing_indices]
  newt= Tt.([t[i] for i in non_missing_indices])

  return hcat(newUs...), newt
end


text:
DATA_SIZE = 360
index_timetype = Date(2000, 1,1) + Day.(0:(DATA_SIZE - 1))
vec1 = randn(DATA_SIZE)
vec2 = randn(DATA_SIZE)
vec3 = randn(DATA_SIZE)
df = DataFrame(Index = index_timetype, vec1 = vec1,vec2 = vec2, vec3 = vec3)
ts = TSx.TS(df, 1)

functions = [mean, median, sum, minimum, maximum, std]
cols = [2, 3, :vec1, :vec3]          # col 1 is index
windowsize = [1, 5, 100, DATA_SIZE]

@test typeof(TSx.rollapply(mean, ts, 2, 5)) == TSx.TS

for fun in functions
    res = @test TSx.rollapply(fun, ts, 1, 5).coredata[!, 2] == RollingFunctions.rolling(fun, df[!, 2], 5)
end

for fun in functions
    colname = "vec1_rolling_$(fun)"
    res = @test TSx.rollapply(fun, ts, :vec1, 5).coredata[!, colname] == RollingFunctions.rolling(fun, df[!, :vec1], 5)

end

for fun in functions
    colname = "vec3_rolling_$(fun)"
    res = @test TSx.rollapply(fun, ts, :vec3, 100).coredata[!, colname] == RollingFunctions.rolling(fun, df[!, :vec3], 100)
end

for fun in functions
    colname = "vec3_rolling_$(fun)"
    res = @test TSx.rollapply(fun, ts, :vec3, DATA_SIZE).coredata[!, colname] == RollingFunctions.rolling(fun, df[!, :vec3], DATA_SIZE)
end

for fun in functions
    res = @test_throws ErrorException TSx.rollapply(fun, ts, 2, 10000)
end


text:
using JuliaWebAPI
using Test

include("clnt.jl")

const opts = Base.JLOptions()
const inline_flag = opts.can_inline == 1 ? `` : `--inline=no`
const cov_flag = (opts.code_coverage == 1) ? `--code-coverage=user` :
                 (opts.code_coverage == 2) ? `--code-coverage=all` :
                 ``

function spawn_srvr()
    srvrscript = joinpath(dirname(@__FILE__), "srvr.jl")
    srvrcmd = `$(joinpath(Sys.BINDIR, "julia")) $cov_flag $inline_flag $srvrscript --runsrvr`
    println("spawining $srvrcmd")
    srvrproc = @static (VERSION < v"0.7.0-") ? spawn(srvrcmd) : run(srvrcmd, wait=false)
end

function kill_spawned_srvr(srvrproc)
    kill(srvrproc)
end

function test_remotesrvr()
    srvrproc = spawn_srvr()
    println("started server process, running client")
    run_clnt(JuliaWebAPI.JSONMsgFormat(), JuliaWebAPI.ZMQTransport("127.0.0.1", 9999, ZMQ.REQ, false))
    println("stopping server process")
    kill_spawned_srvr(srvrproc)
end

# run tests if invoked with run flag
!isempty(ARGS) && (ARGS[1] == "--runclntsrvr") && test_remotesrvr()


text:
module ConservationLawsParticles

using RecursiveArrayTools
using DocStringExtensions

include("utils.jl")
include("densities.jl")
include("models.jl")
include("interactions.jl")
include("diffusion.jl")
include("velocities.jl")
include("plotting.jl")
include("examples.jl")

end # module


text:
import Pkg
Pkg.activate("..")

using CrimsonSkyline
using Distributions: Normal, LogNormal, Poisson
using StatsBase: mean

function log_rate_model(t::Trace, ::Tuple)
    loc = sample(t, "loc", Normal())
    scale = sample(t, "scale", LogNormal())
    sample(t, "log_rate", Normal(loc, scale))
end

function count_model(t::Trace, log_rate::CompiledSF)
    log_rate_value = sample(t, "log_rate", log_rate)
    sample(t, "count", Poisson(exp(log_rate_value)))
end

function main()
    compiled_log_rate_model = compile(
        TypedSF(log_rate_model, Tuple, Float64),
        (),
        ["log_rate"]
    )
    conditioned_count_model = condition(
        count_model,
        Dict("count" => 30)
    )
    results = inference(
        conditioned_count_model, Metropolis();
        params = (compiled_log_rate_model,),
        inference_params = Dict("num_iterations" => 55000, "burn" => 5000, "thin" => 500)
    )
    @info "Posterior rate: $(exp(mean(results, "log_rate")))"

end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
# ---
# title: 1131. Maximum of Absolute Value Expression
# id: problem1131
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math, Bit Manipulation
# link: <https://leetcode.com/problems/maximum-of-absolute-value-expression/description/>
# hidden: true
# ---
#
# Given two arrays of integers with equal lengths, return the maximum value of:
#
# `|arr1[i] - arr1[j]| + |arr2[i] - arr2[j]| + |i - j|`
#
# where the maximum is taken over all `0 <= i, j < arr1.length`.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr1 = [1,2,3,4], arr2 = [-1,4,5,6]
#     Output: 13
#
#
# **Example 2:**
#
#
#
#     Input: arr1 = [1,-2,-5,0,10], arr2 = [0,-2,-1,-7,-4]
#     Output: 20
#
#
#
#
# **Constraints:**
#
#   * `2 <= arr1.length == arr2.length <= 40000`
#   * `-10^6 <= arr1[i], arr2[i] <= 10^6`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using AllocationOpt
using Test
using GraphQLClient

@testset "AllocationOpt.jl" begin
    include("../src/exceptions.jl")
    include("../src/domainmodel.jl")
    include("../src/query.jl")
    include("../src/service.jl")
    include("../src/ActionQueue.jl")

    gateway_url = "https://gateway.thegraph.com/network"

    # Tests
    include("domainmodel.jl")
    include("query.jl")
    include("service.jl")
    include("actionqueue.jl")

    @testset "network_state" begin
        client = Client(gateway_url)
        id = query(client, "indexers"; query_args=Dict("where" => Dict("stakedTokens_gte" => "100000000000000000000000")), output_fields="id").data["indexers"][1]["id"]

        @test_throws AllocationOpt.UnknownIndexerError network_state(
            "0x6ac85b9d834b51b14a7b0ed849bb5199e",
            String[],
            String[],
            String[],
            String[],
            gateway_url,
        )
        @test_throws AllocationOpt.BadSubgraphIpfsHashError network_state(
            id, String[""], String[], String[], String[], gateway_url
        )
    end

    @testset "optimize_indexer" begin
        client = Client(gateway_url)
        id = query(client, "indexers"; query_args=Dict("where" => Dict("stakedTokens_gte" => "100000000000000000000000")), output_fields="id").data["indexers"][1]["id"]
        all_hashes = [
            deployment["ipfsHash"] for deployment in
            query(client, "subgraphDeployments"; query_args=Dict("first" => 1000, "where" => Dict("signalledTokens_gte" => "100000000000000000000000")), output_fields="ipfsHash").data["subgraphDeployments"]
        ]
        ipfshash = all_hashes[1]
        another_ipfshash = all_hashes[2]

        # Indexer stake should be equal to the sum of allocations
        # query indexer's stake
        indexer = query(client, "indexers"; query_args=Dict("where" => Dict("id" => id)), output_fields=["delegatedTokens", "stakedTokens"]).data["indexers"][1]
        stake = togrt(indexer["delegatedTokens"]) + togrt(indexer["stakedTokens"])

        # run optimize_indexer
        repo, optindexer = network_state(
            id, String[ipfshash], String[], String[], String[], gateway_url
        )
        allocs = optimize_indexer(optindexer, repo, 0.0, 10000)
        # Sum allocation amounts
        ω = sum(values(allocs))
        @test isapprox(ω, stake; atol=1e-6)
        # Length of allocations is 1
        @test length(allocs) == 1
        # run optimize_indexer
        repo, optindexer = network_state(
            id,
            String[ipfshash, another_ipfshash],
            String[],
            String[],
            String[],
            gateway_url,
        )
        allocs = optimize_indexer(optindexer, repo, 0.0, 10000)
        # Sum allocation amounts
        ω = sum(values(allocs))
        @test isapprox(ω, stake; atol=1e-6)
        # Length of allocations ≤ length of whitelist
        @test length(allocs) ≤ 2

        # Should handle CSV input
        client = Client(gateway_url)
        id = query(client, "indexers"; query_args=Dict("where" => Dict("stakedTokens_gte" => "100000000000000000000000")), output_fields="id").data["indexers"][1]["id"]
        indexer = query(client, "indexers"; query_args=Dict("where" => Dict("id" => id)), output_fields=["delegatedTokens", "stakedTokens"]).data["indexers"][1]
        stake = togrt(indexer["delegatedTokens"]) + togrt(indexer["stakedTokens"])
        cols = read_filterlists("example.csv")
        repo, optindexer = network_state(id, cols..., gateway_url)
        allocs = optimize_indexer(optindexer, repo, 0.0, 10000)
        ω = sum(values(allocs))
        @test isapprox(ω, stake; atol=1e-6)
    end
end


text:
include("fixtures.jl")

@testset "Window Functions" begin
    @testset "Rectangular Window" begin
        b = rect_window() # width = .5
        @test_throws Exception b(prevfloat(0.))
        @test b(0) == 1.
        @test b(prevfloat(0.5)) == 1.
        @test b(0.5) == 1.
        @test b(nextfloat(0.5)) == 0.
        @test b(100) == 0.
    end

    @testset "Triangular Window" begin
        b = tri_window() # width = 1.
        @test_throws Exception b(prevfloat(0.))
        @test b(0) == 1.
        @test b(.2) == .8
        @test b(.5) == .5
        @test b(1) == 0
    end

    @testset "Lanczos Window" begin
        w = lanczos_window() # lobes=3, width=1
        @test w(0) == 1
        @test w(.5) > .5
        @test w(1) ≈ 0 atol=1e-15
        @test w(1.5) < -.1
        @test w(2) ≈ 0 atol=1e-15
        @test w(2.5) > .01
        @test w(3) ≈ 0 atol=1e-15

        w = lanczos_window(3, width = 2) # lobes=3, width=2
        @test w(0) == 1
        @test w(1) > .5
        @test w(2) ≈ 0 atol=1e-15
        @test w(3) < -.1
        @test w(4) ≈ 0 atol=1e-15
        @test w(5) > .01
        @test w(6) ≈ 0 atol=1e-15
    end

    @testset "Invalid width" begin
        @test_throws AssertionError rect_window(0)
        @test_throws AssertionError rect_window(0.)
    end

    @testset "API Example for Predefined Window Functions" begin
        nuresample(xin, yin, [3.3, 4.1], rect_window(.7),
            required_points_per_slice=1, upsampling_function=rect_window(1.))
    end
end


text:
import StaticArrays, Printf
import LinearAlgebra: dot, cross

export DipoleParameters,
       doppler_temperature,
       dipole_temperature,
       SOLSYSDIR_ECL_THETA,
       SOLSYSDIR_ECL_PHI,
       SOLSYSSPEED_M_S,
       SPEED_OF_LIGHT_M_S,
       T_CMB,
       SPEED_OF_LIGHT_M_S,
       PLANCK_H_MKS,
       BOLTZMANN_K_MKS

# These data have been taken from Planck 2015 I, Table 1 (10.1051/0004-6361/201527101)
const SOLSYSDIR_ECL_THETA = 1.7656051330336222
const SOLSYSDIR_ECL_PHI = 2.9958842149922833
const SOLSYSSPEED_M_S = 370082.2332
const T_CMB = 2.72548

const SPEED_OF_LIGHT_M_S = 2.99792458e8
const PLANCK_H_MKS = 6.62606896e-34
const BOLTZMANN_K_MKS = 1.3806504e-23

@doc raw"""
This structure encodes the information needed to estimate the contribution of the
temperature of the CMB dipole caused by the motion of the Solar System with respect
to the CMB rest frame.

# Fields

- `solsysdir_theta_rad`: colatitude (in radians) of the dipole axis
- `solsysdir_phi_rad`: longitude (in radians) of the dipole axis
- `solsys_speed_m_s`: speed (in m/s) of the reference frame with respect to the
  CMB
- `solsys_velocity_m_s`: velocity 3-vector (in m/s) of the reference frame with
  respect to the CMB
- `tcmb_k`: thermodynamic temperature (in Kelvin) of the CMB

# Object creation

The simplest way to create a `DipoleParameters` object is to call the
constructor without arguments. In this case, the dipole axis will be provided in
Ecliptic coordinates, using the estimate by Planck 2015, and the best COBE
estimate for the CMB monopole temperature will be used. Otherwise, you can pass
any of the following keywords to set the parameters:

- `theta_rad` defaults to `SOLSYSDIR_ECL_THETA`
- `phi_rad` defaults to `SOLSYSDIR_ECL_PHI`
- `speed_m_s` defaults to `SOLSYSSPEED_M_S`
- `t_k` defaults to `T_CMB`

Here is an example where we produce a dipole that is 10% stronger than Planck's:

````julia
dip = DipoleParameters(speed_m_s = SOLSYS_SPEED_VEC_M_S * 1.10)
````
"""
struct DipoleParameters
    solsysdir_theta_rad::Float64
    solsysdir_phi_rad::Float64
    solsys_speed_m_s::Float64
    solsys_velocity_m_s::StaticArrays.StaticVector
    tcmb_k::Float64

    DipoleParameters(;
        theta_rad = SOLSYSDIR_ECL_THETA,
        phi_rad = SOLSYSDIR_ECL_PHI,
        speed_m_s = SOLSYSSPEED_M_S,
        t_k = T_CMB) = new(theta_rad,
        phi_rad,
        speed_m_s,
        StaticArrays.SVector{3}(speed_m_s * [
            sin(theta_rad) * cos(phi_rad),
            sin(theta_rad) * sin(phi_rad),
            cos(theta_rad),
        ]),
        t_k)
end

function Base.show(io::IO, params::DipoleParameters)
    if get(io, :compact, false)
        Printf.@printf(io, "DipoleParameters([%.4f, %.4f, %.4f], %.4f)",
            params.solsys_velocity_m_s[1],
            params.solsys_velocity_m_s[2],
            params.solsys_velocity_m_s[3],
            params.tcmb_k,
        )
    else
        Printf.@printf(io, """DipoleParameters(
    theta_rad=%.6e,
    phi_rad=%.6e,
    speed_m_s=%.6e,
    t_k=%.6e,
)""",
            params.solsysdir_theta_rad,
            params.solsysdir_phi_rad,
            params.solsys_speed_m_s,
            params.tcmb_k,
        )
    end
end

function doppler_temperature(velocity_m_s, dir, tcmb_k)
    betavec = velocity_m_s / SPEED_OF_LIGHT_M_S
    gamma = 1 / sqrt(1 - dot(betavec, betavec))

    tcmb_k * (1 / (gamma * (1 - dot(betavec, dir))) - 1)
end

function doppler_temperature(velocity_m_s,
    dir,
    tcmb_k, freq_hz,
)
    fact = PLANCK_H_MKS * freq_hz / (BOLTZMANN_K_MKS * tcmb_k)
    expfact = exp(fact)
    q = (fact / 2) * (expfact + 1) / (expfact - 1)

    betavec = velocity_m_s / SPEED_OF_LIGHT_M_S
    dotprod = dot(betavec, dir)
    tcmb_k * (dotprod + q * dotprod^2)
end

@doc raw"""
    doppler_temperature(velocity_m_s, dir, tcmb_k)
    doppler_temperature(velocity_m_s, dir, tcmb_k, freq_hz)

Compute the temperature caused by the motion by `velocity_m_s` (a 3-vector, in
m/s) through the Doppler effect. If `freq_hz` is specified, the computation
includes quadrupolar relativistic corrections.

# See also

If you need to compute the solar dipole caused by the motion of the Solar System
with respect to the CMB rest frame, it is easier to use
[`dipole_temperature`](@ref).

# Example

````jldoctest
julia>

julia>
````
"""
doppler_temperature

################################################################################

function dipole_temperature(dir;
    params::DipoleParameters = DipoleParameters(),)
    doppler_temperature(params.solsys_velocity_m_s, dir, params.tcmb_k)
end

function dipole_temperature(dir,
    freq_hz::Number;
    params::DipoleParameters = DipoleParameters(),)
    doppler_temperature(params.solsys_velocity_m_s, dir, params.tcmb_k, freq_hz)
end

@doc raw"""
    dipole_temperature(dir; params::DipoleParameters = DipoleParameters())
    dipole_temperature(dir, freq_hz; params::DipoleParameters = DipoleParameters())

Compute the temperature of the dipole along the direction `dir`, which should be
a 3 - element array containing the XYZ components of the pointing vector.The
vector must be expressed in the same coordinates as the vector used to specify
the dipole in `params`, which is an object of type `DipoleParameters`.(Hint:if
you created this object calling the constructor without arguments, as in
`DipoleParameters()`, Ecliptic coordinates will be used.)

If `freq_hz` is specified, a relativistic frequency - dependent correction is
applied.

# See also

If you need to compute the temperature caused by any other kinetic component
(e.g., the orbital dipole caused by the motion of the spacecraft), use
[`doppler_temperature`](@ref).

# Example

````jldoctest
julia> dipole_temperature([0, 0, 1])
-0.0006532169921239991

julia> dipole_temperature([0, 0, 1], params=DipoleParameters(speed_m_s=371e3))
-0.0006548416782232279

julia> dipole_temperature([0, 0, 1], 30e9, params=DipoleParameters(speed_m_s=371e3))
-0.0006527515338213527
````
"""
dipole_temperature


text:
@info("Simulation: Chained Equations vs Chained ALM")

let
    Random.seed!(65536)
    ROW = 500
    COL = 500
    input_rank = 50
    input_sample = 80
    truth_matrix        = rand([(FixedRankMatrix(Distributions.Gaussian(0, 1), rank = input_rank), ROW, COL)])
    sample_model        = provide(Sampler{BernoulliModel}(), rate = input_sample / 100)
    input_matrix        = sample_model.draw(truth_matrix)
    manual_type_matrix  = Array{Symbol}(undef, ROW, COL)
    manual_type_matrix .= :Gaussian
    complete(MICE(), A = input_matrix, type_assignment = manual_type_matrix)
end


text:
module TextSearch

using NBInclude

nbinclude("text_search.ipynb")

end # module


text:
// Test scope

int main() {
    int a = 5;
    return f();
}

int f() {
    return a;
}

max_stars_repo_name: fredefox/javalette

max_issues_repo_name: fredefox/javalette

max_forks_repo_name: fredefox/javalette


text:
"""
    Streams:
        Notes:
            The whereclause and expressions in the sink types allow the use of 2 variables, RANDOM and ROW.
            These can be used to do random sampling or split up the table.
            The DataMovement service allows for a maximum of 1000000 rows to be written to an output table at one time. Larger results
            can be written through multiple calls to the stream function using ROW filters in the whereclause, the sqlwhereclause or the
            rechunk_datalake_csv function
        Functions:
            - stream(source,sink) returns StreamResponse
            - sinktosource(sink,streamresponse)
            - getdatatypes()
        Types:
            Sources:
            - SqlSourceParameters
            - CsvNotebookFileShareSourceParameters
            - CsvDataLakeSourceParameters
            - ExternalSqlSourceParameters
            Sinks:
            - SqlSinkParameters
            - CsvNotebookFileShareSinkParameters
            - CsvDataLakeSinkParameters
    Sql commands:
        Functions:
            - sqltodf(datastorename::String,sqlquery::SqlQuery)
            - sqltocsv(datastorename::String,sqlquery::SqlQuery,filename::String,folderpathfromroot::Union{String,Void}=nothing)
            - sqltransformation(datastorename::String,sqlquery::SqlQuery,sinktablename::String)
            - sqldeleterows!(datastorename::String,tablename::String;whereclause=nothing)
        Types:
            - SqlQuery
            - SqlJoin
    Datalake commands:
        Functions:
            - delete_datalake_file!(datastorename::String,tablename::String,filename_including_partition::String)
            - list_datalake_table_files_with_partitions(datastorename::String,tablename::String)
            - get_lines_in_datalake_csv(datastorename::String,tablename::String,filename_including_partition::AbstractString)
            - rechunk_datalake_csv(datastorename::String,fromtablename::String,filename_including_partition::AbstractString,totablename::String)
    Schema manager commands:
        Functions:
            - create_destination_table(;storename::String=nothing,tabledefinition::DestinationTableDefinition=nothing)
            - get_table_definition(;storename::String=nothing,tablename::String=nothing)
            - add_destination_table_indexes(;storename::String=nothing,tablename::String=nothing,tableindexes::Array{DestinationTableDefinitionIndex,1}=nothing)
            - save_table_definition(;tabledefinition::DestinationTableDefinition=nothing,filename::String=nothing)
            - load_table_definition(;filename::String=nothing)
            - create_table_mapping(;storename=nothing,tablename=nothing,mappingname=nothing,notmapped=Array{String,1}(),source_dest_name_pairs=Array{Tuple{String,String}}(0))
            - create_processed_table!(datastorename::String,tablename::String,columnnames::Array{String,1},columntypes::Array{String,1};partitionpath::Union{String,Void}=nothing)
            - delete_processed_table!(datastorename::String,tablename::String)
        Types:
            - DestinationTableDefinitionColumn
            - DestinationTableDefinition
            - DestinationTableDefinitionIndex
"""

module NeuroData
    using NeuroJulia
    using DataFrames
    using JSON
    using CSV

    include(Pkg.dir() * "/NeuroJulia/src/NeuroData/NeuroDataTypes.jl")
    include(Pkg.dir() * "/NeuroJulia/src/NeuroData/NeuroDataSqlQuery.jl")
    include(Pkg.dir() * "/NeuroJulia/src/NeuroData/NeuroDataSourceSink.jl")
    include(Pkg.dir() * "/NeuroJulia/src/NeuroData/NeuroDataStream.jl")
    include(Pkg.dir() * "/NeuroJulia/src/NeuroData/NeuroDataSqlCommands.jl")
    include(Pkg.dir() * "/NeuroJulia/src/NeuroData/NeuroDataDatalakeCommands.jl")
    include(Pkg.dir() * "/NeuroJulia/src/NeuroData/NeuroDataSchemaManager.jl")
end

max_stars_repo_name: SnowdenNeuroverse/NeuroJulia

max_issues_repo_name: SnowdenNeuroverse/NeuroJulia

max_forks_repo_name: SnowdenNeuroverse/NeuroJulia


text:
#
# Draft ACEModel Interface:
#
# * A model incorporates both the parameterisation and the parameters.
# * E.g. a LinearACEModel would know about the basis and the coefficients
# * We can then perform the following operations:
#    - evaluate the model at a given configuration, with given parameters
#    - take gradient w.r.t. the configuration
#    - take gradient w.r.t. the parameters
#
# An advantage will be that we can do the parameter reduction trick and then
# use a fast evaluator to obtain the gradient w.r.t. configuration which is
# very expensive if we naively take derivatives w.r.t. the basis and then
# apply parameters. (this is not implemented below and will be the next step)
# Another advantage is that model, and full model construction are all stored
# together for future inspection.

struct LinearACEModel{TB, TP, TEV} <: AbstractACEModel
   basis::TB
   c::Vector{TP}
   evaluator::TEV
end

struct NaiveEvaluator end

function LinearACEModel(basis::SymmetricBasis, c = zeros(length(basis));
               evaluator = :standard)
   if evaluator == :naive
      ev = NaiveEvaluator()
   elseif evaluator == :standard
      ev = PIEvaluator(basis, c)
   elseif evaluator == :recursive
      error("Recursive evaluator not yet implemented")
   else
      error("unknown evaluator")
   end
   return LinearACEModel(basis, c, ev)
end

params(m::LinearACEModel) = copy(m.c)

function set_params!(m::LinearACEModel, c)
   m.c[:] .= c
   set_params!(m.evaluator, m.basis, c)
   return m
end

set_params!(::NaiveEvaluator, args...) = nothing

# ------------------- dispatching on the evaluators

alloc_temp(m::LinearACEModel) = alloc_temp(m.evaluator, m)

evaluate(m::LinearACEModel, X::AbstractConfiguration) =
      evaluate!(alloc_temp(m), m, X)

evaluate!(tmp, m::LinearACEModel, X::AbstractConfiguration) =
      evaluate!(tmp, m::LinearACEModel, m.evaluator, X::AbstractConfiguration)

alloc_temp_d(m::LinearACEModel, X::AbstractConfiguration) =
      alloc_temp_d(m, length(X))

alloc_temp_d(m::LinearACEModel, N::Integer) =
      alloc_temp_d(m.evaluator, N, m)

# this one seems generic and doesn't need to be dispatched?
alloc_grad_config(m::LinearACEModel, X::AbstractConfiguration) =
      Vector{gradtype(m.basis)}(undef, length(X))

grad_config(m::LinearACEModel, X::AbstractConfiguration) =
      grad_config!(alloc_grad_config(m, X), alloc_temp_d(m, X), m, X)

grad_config!(g, tmpd, m::LinearACEModel, X::AbstractConfiguration) =
      grad_config!(g, tmpd, m, m.evaluator, X)

grad_params(m::LinearACEModel, X::AbstractConfiguration) =
      grad_params!(alloc_B(m.basis), alloc_temp(m.basis), m, X)

function grad_params!(g, tmp, m::LinearACEModel, X::AbstractConfiguration)
   evaluate!(g, tmp, m.basis, X)
   return g
end

function grad_params_config(m::LinearACEModel, X::AbstractConfiguration)
   tmpd = alloc_temp_d(m.basis, length(X))
   dB = alloc_dB(m.basis, length(X))
   return grad_params_config!(dB, tmpd, m, X)
end

grad_params_config!(dB, tmpd, m::LinearACEModel, X::AbstractConfiguration)  =
      evaluate_d!(dB, tmpd, m.basis, X)

# ------------------- implementation of naive evaluator

alloc_temp(::NaiveEvaluator, m::LinearACEModel) =
   ( tmpbasis = alloc_temp(m.basis),
     B = alloc_B(m.basis)
   )

function evaluate!(tmp, m::LinearACEModel, ::NaiveEvaluator,
                    X::AbstractConfiguration)
   evaluate!(tmp.B, tmp.tmpbasis, m.basis, X)
   return sum(prod, zip(m.c, tmp.B))
end

alloc_temp_d(::NaiveEvaluator, N::Integer, m::LinearACEModel) =
      ( tmpdbasis = alloc_temp_d(m.basis, N),
        B = alloc_B(m.basis),
        dB = alloc_dB(m.basis, N)
      )

function grad_config!(g, tmpd, m::LinearACEModel, ::NaiveEvaluator,
                     X::AbstractConfiguration)
   evaluate_d!(tmpd.dB, tmpd.tmpdbasis, m.basis, X)
   fill!(g, zero(eltype(g)))
   for ix = 1:length(X), ib = 1:length(m.basis)
      g[ix] += m.c[ib] * tmpd.dB[ib, ix]
   end
   return g
end


text:
using Test
using CSVFiles
using DataFrames
using Distributions

regenerate = false # do a large MC run, to regenerate information needed for std. errors
samplesize = 1000 # normal MC sample size (takes ~5 seconds)
confidence = 2.576 # 99% CI by default; use 1.96 to apply a 95% CI, but expect more spurious errors

# Monte Carlo distribution information
# Filled in from a run with regenerate = true
information = Dict(
    :td => Dict(:transform => x -> log(x), :mu => 19.268649852193303, :sigma => 1.0189429437682878),
    :tpc => Dict(:transform => x -> log(-x), :mu => 16.723536381185546, :sigma => 0.9100375935397476),
    :tac => Dict(:transform => x -> log(x), :mu => 17.35110886191308, :sigma => 0.3586397403122257),
    :te => Dict(:transform => x -> log(x), :mu => 19.30193679432187, :sigma => 1.048364875171592),
    :c_co2concentration => Dict(:transform => x -> x, :mu => 913073.0003103315, :sigma => 71674.18878428967),
    :ft => Dict(:transform => x -> x, :mu => 8.637901336461365, :sigma => 0.4315112767385383),
    :rt_g => Dict(:transform => x -> x, :mu => 5.888175260413073, :sigma => 1.7673873280904946),
    :sealevel => Dict(:transform => x -> x, :mu => 1.6056649448621665, :sigma => 0.6124270241884567)
)

compare = DataFrame(load(joinpath(@__DIR__, "validationdata/PAGE09montecarloquantiles.csv")))

if regenerate
    println("Regenerating MC distribution information")

    # Perform a large MC run and extract statistics
    do_monte_carlo_runs(100_000)
    df = DataFrame(load(joinpath(@__DIR__, "../output/mimipagemontecarlooutput.csv")))

     for ii in 1:nrow(compare)
        name = Symbol(compare[ii, :Variable_Name])
        if kurtosis(df[name]) > 2.9 # exponential distribution
            if name == :tpc # negative across all quantiles
                print("    :$name => Dict(:transform => x -> log(-x), :mu => $(mean(log(-df[df[name] .< 0, name]))), :sigma => $(std(log(-df[df[name] .< 0, name]))))")
            else
                print("    :$name => Dict(:transform => x -> log(x), :mu => $(mean(log(df[df[name] .> 0, name]))), :sigma => $(std(log(df[df[name] .> 0, name]))))")
            end
        else
            print("    :$name => Dict(:transform => x -> x, :mu => $(mean(df[name])), :sigma => $(std(df[name])))")
        end
        if ii != nrow(compare)
            println(",")
        end
    end
else
    println("Performing MC sample")
    # Perform a small MC run
    do_monte_carlo_runs(samplesize)
    df = DataFrame(load(joinpath(@__DIR__, "../output/mimipagemontecarlooutput.csv")))
end

# Compare all known quantiles
for ii in 1:nrow(compare)
    name = Symbol(compare[ii, :Variable_Name])
    transform = information[name][:transform]
    distribution = Normal(information[name][:mu], information[name][:sigma])
    for qval in [.05, .10, .25, .50, .75, .90, .95]
        estimated = transform(quantile(collect(Missings.skipmissing(df[!, name])), qval)) # perform transform *after* quantile, so captures effect of all values
        stderr = sqrt(qval * (1 - qval) / (samplesize * pdf(distribution, estimated)^2))

        expected = transform(compare[ii, Symbol("perc_$(trunc(Int, qval * 100))")])

        #println("$name x $qval: $estimated ≈ $expected rtol=$(ceil(confidence * stderr, -trunc(Int, log10(stderr))))")
        @test estimated ≈ expected rtol=ceil(confidence * stderr; digits = -trunc(Int, log10(stderr)))
    end
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
module RocketTestActorTest

using Test
using Rocket

import Rocket: test_actor, check_isvalid, isreceived, iscompleted, isfailed, check_data_equals, check_error_equals

@testset "TestActor" begin

    println("Testing: actor TestActor")

    @testset begin
        actor = test_actor(Int)

        @test check_isvalid(actor)
    end

    @testset begin
        actor1 = test_actor(Any)
        actor2 = test_actor(Any)

        push!(Rocket.data(actor1), Rocket.ErrorTestEvent("err"))
        push!(Rocket.data(actor2), Rocket.CompleteTestEvent())

        @test_throws Rocket.DataEventIncorrectTypeException check_isvalid(actor1)
        @test_throws Rocket.DataEventIncorrectTypeException check_isvalid(actor2)
    end

    @testset begin
        actor1 = test_actor(Any)
        actor2 = test_actor(Any)

        push!(Rocket.errors(actor1), Rocket.DataTestEvent("data"))
        push!(Rocket.errors(actor2), Rocket.CompleteTestEvent())

        @test_throws Rocket.ErrorEventIncorrectTypeException check_isvalid(actor1)
        @test_throws Rocket.ErrorEventIncorrectTypeException check_isvalid(actor2)
    end

    @testset begin
        actor1 = test_actor(Any)
        actor2 = test_actor(Any)

        push!(Rocket.completes(actor1), Rocket.DataTestEvent("data"))
        push!(Rocket.completes(actor2), Rocket.ErrorTestEvent("err"))

        @test_throws Rocket.CompleteEventIncorrectTypeException check_isvalid(actor1)
        @test_throws Rocket.CompleteEventIncorrectTypeException check_isvalid(actor2)
    end

    @testset begin
        actor = test_actor(Any)

        d1 = Rocket.DataTestEvent(0)
        d2 = Rocket.DataTestEvent(0)

        push!(Rocket.data(actor), d2)
        push!(Rocket.data(actor), d1)

        @test_throws Rocket.DataEventIncorrectTimestampsOrderException check_isvalid(actor)
    end

    @testset begin
        actor = test_actor(Any)

        error!(actor, "error1")
        error!(actor, "error2")

        @test_throws Rocket.MultipleErrorEventsException check_isvalid(actor)
    end

    @testset begin
        actor = test_actor(Any)

        complete!(actor)
        complete!(actor)

        @test_throws Rocket.MultipleCompleteEventsException check_isvalid(actor)
    end

    @testset begin
        actor = test_actor(Any)

        complete!(actor)
        error!(actor, "error")

        @test_throws Rocket.ErrorAfterCompleteEventException check_isvalid(actor)
    end

    @testset begin
        actor = test_actor(Any)

        error!(actor, "error")
        complete!(actor)

        @test_throws Rocket.CompleteAfterErrorEventException check_isvalid(actor)
    end

    @testset begin
        actor = test_actor(Any)

        error!(actor, "error")
        next!(actor, 0)

        @test_throws Rocket.NextAfterErrorEventException check_isvalid(actor)
    end

    @testset begin
        actor = test_actor(Any)

        complete!(actor)
        next!(actor, 0)

        @test_throws Rocket.NextAfterCompleteEventException check_isvalid(actor)
    end

    @testset begin
        actor = test_actor(Int)

        next!(actor, "string")

        @test_throws Rocket.UnacceptableNextEventDataTypeException check_isvalid(actor)
    end

    @testset begin
        actor1 = test_actor(Int)

        next!(actor1, 0)
        next!(actor1, 1)
        next!(actor1, 2)
        complete!(actor1)

        @test check_isvalid(actor1)
        @test isreceived(actor1)
        @test iscompleted(actor1)
        @test !isfailed(actor1)
        @test check_data_equals(actor1, [ 0, 1, 2 ])

        actor2 = test_actor(String)

        next!(actor2, "0")
        next!(actor2, "1")
        next!(actor2, "2")
        error!(actor2, "error")

        @test check_isvalid(actor2)
        @test isreceived(actor2)
        @test !iscompleted(actor2)
        @test isfailed(actor2)
        @test check_data_equals(actor2, [ "0", "1", "2" ])
        @test check_error_equals(actor2, "error")

        actor3 = test_actor(Union{Int, String})

        next!(actor3, 0)
        next!(actor3, "1")
        next!(actor3, "2")
        next!(actor3, 3)

        @test check_isvalid(actor3)
        @test isreceived(actor3)
        @test !iscompleted(actor3)
        @test !isfailed(actor3)
        @test check_data_equals(actor3, [ 0, "1", "2", 3 ])

        actor4 = test_actor(Any)

        @test check_isvalid(actor4)
        @test !isreceived(actor4)
        @test !iscompleted(actor4)
        @test !isfailed(actor4)
        @test check_data_equals(actor4, [ ])

        actor5 = test_actor(Any)

        next!(actor5, rand())

        @test_throws Rocket.DataEventsEqualityFailedException check_data_equals(actor5, [ rand() ])

        actor6 = test_actor(Any)

        error!(actor6, rand())

        @test_throws Rocket.ErrorEventEqualityFailedException check_error_equals(actor6, rand())
    end

    @testset begin
        @test test_actor(Int) isa Rocket.TestActor
        @test test_actor()    isa Rocket.TestActorFactory
    end
end

end


text:
module OrthoPoly
using Polynomials

macro recur_poly_gen(name, rule, p₁, extra_args...)
    generator = gensym()
    @eval begin
        if isempty($extra_args)
            ($name)(n::Integer, x = Poly([0, 1])) = ($generator)(n::Integer, x)[1]
        else
            ($name)(n::Integer, $(extra_args...), x = Poly([0, 1])) = ($generator)(n::Integer, x, $(extra_args...))[1]
        end
        function ($generator)(n::Integer, x::T, $(extra_args...)) where T
            if n == 0
                one(x), zero(x)
            elseif n == 1
                $p₁, one(x)
            else
                p₋₁, p₋₂ = ($generator)(n-1, x, $(extra_args...))
                $rule, p₋₁
            end
        end
    end
end

# Jacobi polynomials are orthogonal on the interval [-1,1] with respect
# to the weight (1-x)^α(1+x)^β. Legendre and Chebyshev polynomials are
# special cases.
# TODO: General Jacobi, Gegenbauer and associated Legendre polynomials
# (the latter are not polynomials for odd values of m).
@recur_poly_gen(chebyshevT, 2x*p₋₁ - p₋₂, x) # Chebyshev polynomials, first kind
@recur_poly_gen(chebyshevU, 2x*p₋₁ - p₋₂, 2x) # Chebyshev polynomials, second kind
@recur_poly_gen(legendre, ((2n - 1)x*p₋₁ - (n - 1)p₋₂)/n, x) # Legendre polynomials

# Hermite polynomials are orthogonal on the interval (-∞,∞) with
# weight exp(-x²/a) where a = 2 for probabilists and a = 1 for
# physicists.
@recur_poly_gen(hermite_prob, x*p₋₁ - (n - 1)p₋₂, x) # Probabilists' Hermite polynomials
@recur_poly_gen(hermite, 2x*p₋₁ - 2(n - 1)p₋₂, 2x) # Physicists' Hermite polynomials

# Associated Laguerre polynomials are orthogonal on the interval [0,∞)
# with respect to the weight x^α exp(-x), ∀ α ∈ ℝ.
@recur_poly_gen(laguerre, ((2n - 1 - x)p₋₁ - (n - 1)p₋₂)/n, 1 - x) # Laguerre polynomials
@recur_poly_gen(laguerre_assoc, ((2n - 1 + α - x)p₋₁ - (n - 1 + α)p₋₂)/n, 1 + α - x, α)

export chebyshevT, chebyshevU, legendre, laguerre, laguerre_assoc, hermite, hermite_prob

end # module


text:
# ---
# title: 1019. Next Greater Node In Linked List
# id: problem1019
# author: Indigo
# date: 2022-03-19
# difficulty: Medium
# categories: Linked List, Stack
# link: <https://leetcode.com/problems/next-greater-node-in-linked-list/description/>
# hidden: true
# ---
#
# We are given a linked list with `head` as the first node.  Let's number the
# nodes in the list: `node_1, node_2, node_3, ...` etc.
#
# Each node may have a _next larger_ **value** : for `node_i`,
# `next_larger(node_i)` is the `node_j.val` such that `j > i`, `node_j.val >
# node_i.val`, and `j` is the smallest possible choice.  If such a `j` does not
# exist, the next larger value is `0`.
#
# Return an array of integers `answer`, where `answer[i] =
# next_larger(node_{i+1})`.
#
# Note that in the example **inputs**  (not outputs) below, arrays such as
# `[2,1,5]` represent the serialization of a linked list with a head node value
# of 2, second node value of 1, and third node value of 5.
#
#
#
# **Example 1:**
#
#
#
#     Input: [2,1,5]
#     Output: [5,5,0]
#
#
# **Example 2:**
#
#
#
#     Input: [2,7,4,3,5]
#     Output: [7,0,5,5,0]
#
#
# **Example 3:**
#
#
#
#     Input: [1,7,5,1,9,2,5,1]
#     Output: [7,9,9,9,0,5,0,0]
#
#
#
#
# **Note:**
#
#   1. `1 <= node.val <= 10^9` for each node in the linked list.
#   2. The given list has length in the range `[0, 10000]`.
#
#
## @lc code=start
using LeetCode

function next_larger_nodes(head::ListNode{Int})
    head = reverse_list(head)
    stk = Int[]
    res = Int[]
    while !isnothing(head)
        while !isempty(stk) && stk[end] <= val(head)
            pop!(stk)
        end
        pushfirst!(res, isempty(stk) ? 0 : stk[end])
        push!(stk, val(head))
        head = next(head)
    end
    return res
end
## @lc code=end


text:
argtypes(x::Combined{T, A}) where {T, A} = A
argtypes(x) = Any

function data_limits(x)
    error("No datalimits for $(typeof(x)) and $(argtypes(x))")
end

function data_limits(x::Atomic)
    isempty(x.plots) ? atomic_limits(x) : data_limits(x.plots)
end

"""
Data limits calculate a minimal boundingbox from the data points in a plot.
This doesn't include any transformations, markers etc.
"""
atomic_limits(x::Atomic{<: Tuple{Arg1}}) where Arg1 = FRect3D(to_value(x[1]))

# TODO makes this generically work
atomic_limits(x::Atomic{<: Tuple{<: AbstractVector{<: NTuple{N, <: Number}}}}) where N = FRect3D(Point{N, Float32}.(to_value(x[1])))
function atomic_limits(x::Atomic{<: Tuple{<: AbstractVector{<: NTuple{2, T}}}}) where T <: VecTypes
    FRect3D(reinterpret(T, to_value(x[1])))
end

function atomic_limits(x::Atomic{<: Tuple{X, Y, Z}}) where {X, Y, Z}
    xyz_boundingbox(to_value.(x[1:3])...)
end

function atomic_limits(x::Atomic{<: Tuple{X, Y}}) where {X, Y}
    xyz_boundingbox(to_value.(x[1:2])...)
end

_isfinite(x) = isfinite(x)
_isfinite(x::VecTypes) = all(isfinite, x)
scalarmax(x::AbstractArray, y::AbstractArray) = max.(x, y)
scalarmax(x, y) = max(x, y)
scalarmin(x::AbstractArray, y::AbstractArray) = min.(x, y)
scalarmin(x, y) = min(x, y)

extrema_nan(itr::Pair) = (itr[1], itr[2])
extrema_nan(itr::ClosedInterval) = (minimum(itr), maximum(itr))

function extrema_nan(itr)
    vs = iterate(itr)
    vs === nothing && return (NaN, NaN)
    v, s = vs
    vmin = vmax = v
    # find first finite value
    while vs !== nothing && !_isfinite(v)
        v, s = vs
        vmin = vmax = v
        vs = iterate(itr, s)
    end
    while vs !== nothing
        x, s = vs
        vs = iterate(itr, s)
        _isfinite(x) || continue
        vmax = scalarmax(x, vmax)
        vmin = scalarmin(x, vmin)
    end
    return (vmin, vmax)
end

function xyz_boundingbox(x, y, z = (0 => 0))
    minmax = extrema_nan.((x, y, z))
    mini, maxi = first.(minmax), last.(minmax)
    FRect3D(mini, maxi .- mini)
end

const ImageLike{Arg} = Union{Heatmap{Arg}, Image{Arg}}
function data_limits(x::ImageLike{<: Tuple{X, Y, Z}}) where {X, Y, Z}
    xyz_boundingbox(to_value.((x[1], x[2]))...)
end

function data_limits(x::Volume)
    xyz_boundingbox(to_value.((x[1], x[2], x[3]))...)
end

function text_limits(x::VecTypes)
    p = to_ndim(Vec3f0, x, 0.0)
    FRect3D(p, p)
end
function text_limits(x::AbstractVector)
    FRect3D(x)
end
function atomic_limits(x::Text{<: Tuple{Arg1}}) where Arg1
    boundingbox(x)
end

function data_limits(x::Annotations)
    # data limits is supposed to not include any transformation.
    # for the annotation, we use the model matrix directly, so we need to
    # to inverse that transformation for the correct limits
    bb = data_limits(x.plots[1])
    inv(modelmatrix(x)) * bb
end

Base.isfinite(x::Rect) = all(isfinite.(minimum(x))) &&  all(isfinite.(maximum(x)))

function data_limits(plots::Vector)
    isempty(plots) && return
    bb = FRect3D()
    plot_idx = iterate(plots)
    while plot_idx !== nothing
        plot, idx = plot_idx
        plot_idx = iterate(plots, idx)
        # axis shouldn't be part of the data limit
        isaxis(plot) && continue
        isa(plot, Legend) && continue
        bb2 = data_limits(plot)
        isfinite(bb) || (bb = bb2)
        isfinite(bb2) || continue
        bb = union(bb, bb2)
    end
    bb
end

data_limits(s::Scene) = data_limits(plots_from_camera(s))
data_limits(plot::Combined) = data_limits(plot.plots)


text:
"""
```
  lsq_constrsparsereg_admm(
    X             :: AbstractMatrix{T},
    y             :: AbstractVector{T},
    ρ             :: Number = zero(T);
    proj          :: Function = x -> x,
    obswt         :: Vector{T} = ones(T, length(y)),
    penwt         :: Vector{T} = ones(T, size(X, 2)),
    β0            :: Vector{T} = zeros(T, size(X, 2)),
    admmmaxite    :: Int = 10000,
    admmabstol    :: Float64 = 1e-4,
    admmreltol    :: Float64 = 1e-4,
    admmscale     :: Float64 = 1 / length(y),
    admmvaryscale :: Bool = false
    )
```
Fit constrained lasso at a fixed tuning parameter value by applying the
alternating direction method of multipliers (ADMM) algorithm.

### Arguments
- `X`       : predictor matrix.
- `y`       : response vector.
- `ρ`       : tuning parameter. Default 0.

### Optional arguments
- `proj`         : projection onto the constraint set. Default is identity (no constraint).
- `obswt`        : observation weights. Default is `[1 1 1 ... 1]`.
- `penwt`        : predictor penalty weights. Default is `[1 1 1 ... 1]`.
- `β0`           : starting point.
- `admmmaxite`   : maximum number of iterations for ADMM. Default is `10000`.
- `admmabstol`   : absolute tolerance for ADMM.
- `admmreltol`   : relative tolerance for ADMM.
- `admmscale`    : ADMM scale parameter. Default is `1/n`.
- `admmvaryscale`: dynamically chance the ADMM scale parameter. Default is false.

### Returns
- `β`       : estimated coefficents.
"""
function lsq_constrsparsereg_admm(
    X::AbstractMatrix{T},
    y::AbstractVector{T},
    ρ::Number        = zero(T);
    proj::Function   = x -> x,
    obswt::Vector{T} = ones(T, length(y)),
    penwt::Vector{T} = ones(T, size(X, 2)),
    β0::Vector{T}    = zeros(T, size(X, 2)),
    admmmaxite::Int       = 10000,
    admmabstol::Float64   = 1e-4,
    admmreltol::Float64   = 1e-4,
    admmscale::Float64    = 1 / length(y),
    admmvaryscale::Bool   = false
    ) where T <: Union{Float64}

    n, p = size(X)
    β = copy(β0)
    z = proj(β)
    u = β - z
    v = similar(β)
    zold = similar(z)
    primalres = similar(z)

    # allocate working arrays
    Xaug = vcat(X, Matrix{Float64}(I, p, p))
    yaug = vcat(y, fill(1, p))
    obswtaug = vcat(obswt, fill(1, p))
    λ = [ρ / (n + p)]

    for iter in 1:admmmaxite

        # update working arrays
        for j in 1:p
          admmscaleinv = 1 / √admmscale
          if iter == 1 || admmvaryscale
            Xaug[n+j, j] = admmscaleinv
          end
          yaug[n+j] = (z[j] - u[j]) * admmscaleinv
        end
        # update β - lasso
        path = glmnet(Xaug, yaug;
                weights = obswtaug, lambda = λ, penalty_factor = penwt,
                standardize = false, intercept = false)
        copyto!(β, path.betas)
        # update z - projection to constraint set
        v .= β .+ u
        copyto!(zold, z)
        z = proj(v) ## could not do in-place; z = proj!(v) not working
        # update scaled dual variables u
        dualresnorm = norm((z - zold) / admmscale)
        # dualresnorm = 0.0
        # for j in 1:p
        #   dualresnorm += ((z[j] - zold[j]) / admmscale)^2
        # end
        # dualresnorm ^= 1/2
        primalres .= β .- z
        primalresnorm = norm(primalres)
        u .+= primalres
        # check convergence criterion
        if (primalresnorm <= √p * admmabstol
                + admmreltol * max(norm(β), norm(z))) &&
                (dualresnorm <= √n * admmabstol
                + admmreltol * norm(u / admmscale))
            break
        end
        # update ADMM scale parameter if requested
        if admmvaryscale
            if primalresnorm / dualresnorm > 10
                admmscale /= 2
                u /= 2
            elseif primalresnorm / dualresnorm < 0.1
                admmscale *= 2
                u *= 2
            end
        end
    end

    return β

end # end of function

"""
```
  lsq_constrsparsereg_admm(
    X             :: AbstractMatrix{T},
    y             :: AbstractVector{T},
    ρlist         :: Vector;
    proj          :: Function = x -> x,
    obswt         :: Vector{T} = ones(T, length(y)),
    penwt         :: Vector{T} = ones(T, size(X, 2)),
    admmmaxite    :: Int = 10000,
    admmabstol    :: Float64 = 1e-4,
    admmreltol    :: Float64 = 1e-4,
    admmscale     :: Float64 = 1 / length(y),
    admmvaryscale :: Bool = false
    )
```
Fit constrained lasso at fixed tuning parameter values by applying the
alternating direction method of multipliers (ADMM) algorithm.

### Arguments
- `X`       : predictor matrix.
- `y`       : response vector.
- `ρlist`   : a vector of tuning parameter values.

### Optional arguments
- `proj`         : projection onto the constraint set. Default is identity (no constraint).
- `obswt`        : observation weights. Default is `[1 1 1 ... 1]`.
- `penwt`        : predictor penalty weights. Default is `[1 1 1 ... 1]`.
- `admmmaxite`   : maximum number of iterations for ADMM. Default is `10000`.
- `admmabstol`   : absolute tolerance for ADMM.
- `admmreltol`   : relative tolerance for ADMM.
- `admmscale`    : ADMM scale parameter. Default is `1/n`.
- `admmvaryscale`: dynamically chance the ADMM scale parameter. Default is false.

### Returns
- `βpath`       : estimated coefficents along the grid of `ρ` values.
"""
function lsq_constrsparsereg_admm(
    X::AbstractMatrix{T},
    y::AbstractVector{T},
    ρlist;
    proj::Function   = x -> x,
    obswt::Vector = ones(T, length(y)),
    penwt::Vector = ones(T, size(X, 2)),
    admmmaxite::Int       = 10000,
    admmabstol::Float64   = 1e-4,
    admmreltol::Float64   = 1e-4,
    admmscale::Float64    = 1 / length(y),
    admmvaryscale::Bool   = false
    ) where T

    βpath = zeros(T, size(X, 2), length(ρlist))
    for i in eachindex(ρlist)
      βpath[:, i] = lsq_constrsparsereg_admm(X, y, ρlist[i]; proj = proj,
            obswt = obswt, penwt = penwt, admmmaxite = admmmaxite,
            admmabstol = admmabstol, admmreltol = admmreltol,
            admmscale = admmscale, admmvaryscale = admmvaryscale)
    end
    return βpath
end


text:
using TensorToolbox
using LinearAlgebra
using Arpack
using TensorDecompositions

function lranmf_mu(X, r ; max_iter=200, tol = 1.0E-3, verbose = false)
    # lranmf_mu is effective when r << min(n,m)
    # input is a onnegative matrix
    # proposed by Guoxu Zhou in 2012
    # https://ieeexplore.ieee.org/document/6166354

    n, m = size(X)
    epsilon = 0.0001

    # See step1 in section II-B in the paper
    if r == min(n, m)
        svd_X = svd(X)
    else
        svd_X = svds(X; nsv=r, ritzvec=true)[1]
    end
    Achil = svd_X.U * diagm(svd_X.S)
    Bchil = svd_X.V

    # See step2 in section II-B in the paper
    A = rand(n, r)
    B = rand(m, r)
    cost_at_init = norm(X - A*B')
    previous_cost = cost_at_init
    for iter = 1:max_iter
        B .= B .* ( max.( Bchil*(Achil' * A), epsilon ) ) ./ ( B*(A'*A) )
        A .= A .* ( max.( Achil*(Bchil' * B), epsilon ) ) ./ ( A*(B'*B) )

        if tol > 0 && iter % 10 == 0
            cost = norm(X - A*B')
            if verbose
                println("iter: $iter cost: $cost")
            end
            if (previous_cost - cost) / cost_at_init < tol
                break
            end
            previous_cost = cost
        end
    end

    # A * B' is rank-r matrix
    return A, B'
end

max_stars_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_issues_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_forks_repo_name: gkazunii/Legendre-tucker-rank-reduction


text:
struct Device <: LavaAbstraction
    handle::Vk.Device
    extensions::Vector{String}
    features::Vk.PhysicalDeviceFeatures2
    queues::QueueDispatch
    pipeline_ht::HashTable{Pipeline}
    pipeline_layout_ht::HashTable{PipelineLayout}
    pipeline_layouts::Dictionary{Vk.PipelineLayout,PipelineLayout}
    pending_pipelines::Vector{Vk.GraphicsPipelineCreateInfo}
    shader_cache::ShaderCache
    transfer_ops::Vector{Vk.SemaphoreSubmitInfoKHR}
    command_pools::CommandPools
end

vk_handle_type(::Type{Device}) = Vk.Device

function Device(physical_device::Vk.PhysicalDevice, extensions, queue_config, features::Vk.PhysicalDeviceFeatures2; surface = nothing, next = C_NULL)
    infos = queue_infos(QueueDispatch, physical_device, queue_config)
    info = Vk.DeviceCreateInfo(
        infos,
        [],
        extensions;
        next,
    )

    handle = unwrap(create(Device, physical_device, info))
    queues = QueueDispatch(handle, infos; surface)
    Device(handle, extensions, features, queues, HashTable{Pipeline}(), HashTable{PipelineLayout}(), Dictionary(), [], ShaderCache(handle), [], CommandPools(handle))
end

function request_command_buffer(device::Device, usage::Vk.QueueFlag)
    index = get_queue_family(device.queues, usage)
    pool = request_pool!(device.command_pools, index)
    handle = first(unwrap(Vk.allocate_command_buffers(device, Vk.CommandBufferAllocateInfo(pool, Vk.COMMAND_BUFFER_LEVEL_PRIMARY, 1))))
    cb = SimpleCommandBuffer(handle, index)
    start_recording(cb)
    cb
end

queue_family_indices(device::Device) = queue_family_indices(device.queues)

submit(device, args...; kwargs...) = submit(device.queues, args...; kwargs...)

function create_pipelines!(device::Device)
    batch_create!(device.pipeline_ht, device.pending_pipelines) do infos
        (handles, _) = unwrap(Vk.create_graphics_pipelines(device, infos))
        map(zip(handles, infos)) do (handle, info)
            Pipeline(handle, PipelineType(Vk.PIPELINE_BIND_POINT_GRAPHICS), pipeline_layout(device, info.layout))
        end
    end
end

function pipeline_layout(device::Device, resources)
    info = Vk.PipelineLayoutCreateInfo(
        [resources.gset.set.layout],
        [Vk.PushConstantRange(Vk.SHADER_STAGE_ALL, 0, sizeof(PushConstantData))],
    )
    get!(device.pipeline_layout_ht, info) do info
        handle = unwrap(Vk.create_pipeline_layout(device, info))
        layout = PipelineLayout(handle, info.set_layouts, info.push_constant_ranges)
        insert!(device.pipeline_layouts, handle, layout)
        layout
    end
end

pipeline_layout(device::Device, handle::Vk.PipelineLayout) = device.pipeline_layouts[handle]

function Base.show(io::IO, device::Device)
    print(io, Device, "($(device.handle))")
end


text:
# ---
# title: 1039. Minimum Score Triangulation of Polygon
# id: problem1039
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/minimum-score-triangulation-of-polygon/description/>
# hidden: true
# ---
#
# Given `N`, consider a convex `N`-sided polygon with vertices labelled `A[0],
# A[i], ..., A[N-1]` in clockwise order.
#
# Suppose you triangulate the polygon into `N-2` triangles.  For each triangle,
# the value of that triangle is the **product**  of the labels of the vertices,
# and the _total score_ of the triangulation is the sum of these values over all
# `N-2` triangles in the triangulation.
#
# Return the smallest possible total score that you can achieve with some
# triangulation of the polygon.
#
#
#
# **Example 1:**
#
#
#
#     Input: [1,2,3]
#     Output: 6
#     Explanation: The polygon is already triangulated, and the score of the only triangle is 6.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2019/05/01/minimum-score-
# triangulation-of-polygon-1.png)
#
#
#
#     Input: [3,7,4,5]
#     Output: 144
#     Explanation: There are two triangulations, with possible scores: 3*7*5 + 4*5*7 = 245, or 3*4*5 + 3*4*7 = 144.  The minimum score is 144.
#
#
# **Example 3:**
#
#
#
#     Input: [1,3,1,4,1,5]
#     Output: 13
#     Explanation: The minimum score triangulation has score 1*1*3 + 1*1*4 + 1*1*5 + 1*1*1 = 13.
#
#
#
#
# **Note:**
#
#   1. `3 <= A.length <= 50`
#   2. `1 <= A[i] <= 100`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using CSV
using DataFrames
using DelimitedFiles
using YAML
using Query

using SLiDE

# ******************************************************************************************
#   READ BLUENOTE OUTPUT RESULTS TO CHECK.
# ******************************************************************************************
y = read_file(joinpath("dev","buildstream","check_share.yml"));
bluenote = Dict(k => sort(read_file(joinpath(y["SharePath"]..., ensurearray(v["name"])...)))
    for (k,v) in y["ShareInput"])
bluenote = Dict(Symbol(k) => edit_with(v, Rename.(names(v), Symbol.(y["ShareInput"][k]["col"])))
    for (k,v) in bluenote)

# ******************************************************************************************
#   READ SETS AND SLiDE SUPPLY/USE DATA.
# ******************************************************************************************
y = read_file(joinpath("dev","buildstream","setlist.yml"));
set = Dict(Symbol(k) => sort(read_file(joinpath(y["SetPath"]..., ensurearray(v)...)))[:,1]
    for (k,v) in y["SetInput"])

# Define edits to perform on each input DataFrame.
x = Dict()
x[:gsp] = [
    Map(joinpath("crosswalk","gsp.csv"), [:gsp_code], [:windc_code], [:si], [:s], :left);
    Map(joinpath("..","coresets","yr.csv"), [:yr], [:yr], [:yr], [:yr], :left);
    Map(joinpath("..","coresets","r","state.csv"), [:r], [:r], [:r], [:r], :left);
    Drop.([:yr,:r,:s], missing, "==");
    Order([:yr,:r,:gdpcat,:s,:si,:value], [Int, String, String, String, Int, Float64])
];
x[:pce] = [
    Map(joinpath("crosswalk","pce.csv"), [:pg], [:g], [:pg], [:g], :left);
    Map(joinpath("..","coresets","yr.csv"), [:yr], [:yr], [:yr], [:yr], :left);
    Map(joinpath("..","coresets","r","state.csv"), [:r], [:r], [:r], [:r], :left);
    Drop.([:yr,:r,:g], missing, "==");
    Order([:yr,:r,:g,:value], [Int, String, String, Float64])
];
x[:utd] = [
    Add(:share, 0.0);
    Map(joinpath("crosswalk","naics.csv"), [:naics_code], [:windc_code], [:n], [:s], :left);
    Map(joinpath("..","coresets","yr.csv"), [:yr], [:yr], [:yr], [:yr], :left);
    Map(joinpath("..","coresets","r","state.csv"), [:r], [:r], [:r], [:r], :left);
    Drop.([:yr,:r], missing, "==");
    Order([:yr,:r,:s,:t,:value,:n], [Int, String, String, String, Float64, Int])
];
x[:cfs] = [
    Map(joinpath("crosswalk","sctg.csv"), [:sctg_code], [:windc_code], [:sg], [:g], :left);
    Drop(:g, missing, "==")
    # Order([:orig_state,:dest_state,:n,:g,:units,:value],[String,String,Int,String,String,Float64])
]

# Read share info.
shr = Dict()
shr[:cfs] = read_file(joinpath("data", "output", "cfs_state.csv"))
shr[:gsp] = read_file(joinpath("data", "output", "gsp_state.csv"))
shr[:pce] = read_file(joinpath("data", "output", "pce.csv"))
shr[:sgf] = read_file(joinpath("data", "output", "sgf.csv"))
shr[:utd] = read_file(joinpath("data", "output", "utd.csv"))

[global shr[k] = edit_with(shr[k], x[k]) for k in keys(shr) if k in keys(x)]

function filter_with(df::DataFrame, set::Dict)
    cols = find_oftype(df, Not(AbstractFloat))

    cols_sets = intersect(cols, collect(keys(set)));
    vals_sets = [set[k] for k in cols_sets]
    list_sets = NamedTuple{Tuple(cols_sets,)}(vals_sets,)

    df = fill_zero(list_sets, df);
    return df
end

# # ******************************************************************************************
# """
# ## PCE Shares
# """
# shr[:pce][!,:share] .= shr[:pce][:,:value] ./ sum_over(shr[:pce], :r; keepkeys = true)
# # PCE Share. Regional shares of final consumption
# #!  pce_map(r,g,yr) = sum(map(g,pg), pce_raw_units(yr,r,pg,"millions of us dollars (USD)"));
# #!  pce_shr(yr,r,g) = pce_map(r,g,yr) / sum(r.local, pce_map(r,g,yr));

# # ******************************************************************************************
# # # Share of total trade by region.
# # #!  usatrd_(yr,r,s,t) = sum(map(n,s), usatrd(r,n,yr,t));
# # #!  usatrd_shr(yr,r,s,t)$(NOT notinc(s) and sum(r.local, usatrd_(yr,r,s,t)))
# # #!      = usatrd_(yr,r,s,t) / sum(r.local, usatrd_(yr,r,s,t));
# # #!  usatrd_shr(yr,r,s,t)$(NOT notinc(s) AND NOT sum(r.local, usatrd_(yr,r,s,t)))
# # #!      = sum(yr.local, usatrd_(yr,r,s,t)) / sum((r.local,yr.local), usatrd_(yr,r,s,t));
# # shr[:utd] = fill_zero(sum_over(shr[:utd], :n; values_only = false))

# # shr[:utd][!,:share] .= shr[:utd][:,:value] ./ sum_over(shr[:utd], :r; keepkeys = true)
# # shr[:utd][isnan.(shr[:utd][:,:share]),:share] .= (sum_over(shr[:utd], :yr; keepkeys = true) ./
# #     sum_over(shr[:utd], [:yr,:r]; keepkeys = true))[isnan.(shr[:utd][:,:share])]

# # # Sectors not included in USA Trade Data.
# # #!  notinc(s) = yes$(NOT sum(n, map(n,s)));
# # shr[:notinc] = DataFrame(s = setdiff(set[:i], shr[:utd][:,:s]))

# # """
# # ## GSP Shares
# # """
# # df = fill_zero(copy(shr[:sgf]))

# # Labor share of value added from national dataset
# #!  PARAMETER lshr_0(yr,s) "Labor share of value added from national dataset";
# #!  lshr_0(yr,s)$(va_0(yr,'compen',s) + va_0(yr,'surplus',s)) =
# #!      va_0(yr,'compen',s) / (va_0(yr,'compen',s) + va_0(yr,'surplus',s));
# # !!!! I would move to partitionbea.jl
# io[:va0] = edit_with(io[:va0], Rename(:j,:s))
# io[:lshr0] = io[:va0] |> @filter(_.va == set[:va][1])  |> DataFrame
# io[:lshr0] = edit_with(io[:lshr0], Drop(:va, "all", "=="))

# io[:lshr0][!,:value] = io[:va0][io[:va0][:,:va] .== "compen",:value] ./
#     (io[:va0][io[:va0][:,:va] .== "compen",:value] + io[:va0][io[:va0][:,:va] .== "surplus",:value])

# # Calculated gross state product.
# shr[:gsp] = fill_zero(sum_over(shr[:gsp], :si; values_only = false))
# shr[:gsp][!,:calc] .= shr[:gsp][:,:value]

# shr[:gsp][shr[:gsp][:,:gdpcat] .== "gdp",:calc] .=
#     sum_over(shr[:gsp][findall(in(["cmp","gos","taxsbd"]), shr[:gsp][:,:gdpcat]),:], :gdpcat)
# shr[:gsp][!,:diff] = shr[:gsp][:,:calc] - shr[:gsp][:,:value]

# #!  region_shr(yr,r,s)$(sum(r.local, gsp0(yr,r,s,'Reported'))) =
# #!      gsp0(yr,r,s,'Reported') / sum(r.local,  gsp0(yr,r,s,'Reported'));
# #!  region_shr(yr,r,'use')$sum((r.local,s), region_shr(yr,r,s)) =
# #!      sum(s, region_shr(yr,r,s)) / sum((r.local,s), region_shr(yr,r,s));
# #!  region_shr(yr,r,'oth')$sum((r.local,s), region_shr(yr,r,s)) =
# #!      sum(s, region_shr(yr,r,s)) / sum((r.local,s), region_shr(yr,r,s));

# """
# Regional share of value added.
# """
# shr[:region] = shr[:gsp][shr[:gsp][:,:gdpcat] .== "gdp", [:yr,:r,:s,:value]]
# shr[:region][!,:share] .= shr[:region][:,:value] ./ sum_over(shr[:region], :r; keepkeys = true)

# # !!!! Not matching
# for SECT in ["use", "oth"]
#     df_temp = fill_zero((yr = set[:yr], r = set[:r], s = SECT))
#     df_temp[!,:share] .= sum_over(shr[:region], :s) ./
#         sum_over(sum_over(shr[:region], :r; keepkeys = true, values_only = false), :s)

#     global shr[:region] = [shr[:region]; df_temp]
# end

# # ******************************************************************************************
# # # # CFS
# # # df_cfs = copy(shr[:cfs])
# # # shr[:cfs] = sort(edit_with(shr[:cfs], x[:cfs]))

# # # # Local supply-demand. Trade that remains within the same region.
# # # #!  PARAMETER d0(r,g) "Local supply-demand (CFS)";
# # # #!  d0_(r,n,sg) = cfs2012_units(r,r,n,sg,"millions of us dollars (USD)");
# # # #!  d0(r,g) = sum(map(sg,g), sum(n, d0_(r,n,sg)));
# # # shr[:d0] = shr[:cfs][shr[:cfs][:,:orig_state] .== shr[:cfs][:,:dest_state],:]
# # # shr[:d0] = edit_with(shr[:d0], Rename(:orig_state,:r))[:,[:r,:n,:g,:units,:value]]
# # # shr[:d0] = sum_over(shr[:d0], :n; values_only = false)

# # # # Interstate trade (CFS)
# # # #!  PARAMETER mrt0(r,r,g) "Interstate trade (CFS)";
# # # #!  mrt0_(r,rr,n,sg)$(NOT SAMEAS(r,rr)) = cfs2012_units(r,rr,n,sg,"millions of us dollars (USD)");
# # # #!  mrt0(r,rr,g) = sum(map(sg,g), sum(n, mrt0_(r,rr,n,sg)));
# # # shr[:mrt0] = shr[:cfs][shr[:cfs][:,:orig_state] .!= shr[:cfs][:,:dest_state],:]
# # # shr[:mrt0] = sum_over(shr[:mrt0], :n; values_only = false)

# # # # National exports (CFS)
# # # #!  PARAMETER xn0(r,g) "National exports (CFS)";
# # # #!  xn0(r,g) = sum(rr, mrt0(r,rr,g));
# # # shr[:x0] = sum_over(shr[:mrt0], :dest_state; values_only = false)
# # # shr[:x0] = edit_with(shr[:x0], Rename(:orig_state, :r))

# # # # National demand (CFS)
# # # #   PARAMETER mn0(r,g) "National demand (CFS)";
# # # #   mn0(r,g) = sum(rr, mrt0(rr,r,g));
# # # shr[:mn0] = sum_over(shr[:mrt0], :orig_state; values_only = false)
# # # shr[:mn0] = edit_with(shr[:mn0], Rename(:dest_state, :r))

# # # # Regional purchase coefficient
# # # #   PARAMETER rpc(*,g) "Regional purchase coefficient";
# # # #   rpc(r,g)$(d0(r,g) + mn0(r,g)) = d0(r,g) / (d0(r,g) + mn0(r,g));
# # # shr[:mn0], shr[:d0] = fill_zero(shr[:mn0], shr[:d0]);
# # # shr[:rpc] = shr[:d0]
# # # shr[:rpc][!,:value] .= shr[:d0][:,:value] ./ (shr[:mn0][:,:value] + shr[:d0][:,:value])

# # # UTI = 0.9
# # # shr[:rpc] = [shr[:rpc];
# # #     DataFrame(permute((
# # #         r = set[:r], g = :uti, units = "billions of us dollard (USD)", value = UTI)))]

# # # # # # ii_zero = shr[:mn0] + shr[:d0] .== 0.0

# # # # # # df = copy(shr[:cfs])

# # # # # # first(df,3)

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
@testset "cauchyendgame" begin
     PolyImpl.@polyvar x
     F = PolySystem([(x - 2.0)^4])
     G, solutions = totaldegree(F)

     H = StraightLineHomotopy(G,F)

     results = solve(H, solutions)
     for result in results
          @test norm([2.0] - result.solution) ≈ 0.0 atol=1e-3
     end

     # check that show doesn't throw
     @test all(x -> length(string(x)) > 0, results)

     @test all(x -> !isnull(x.convergent_cluster), results)
end