text:
# Execution control
# modeled after: CUDAdrv/src/execution.jl

export ROCDim, ROCModule, ROCFunction, roccall

## kernel allocations

struct KernelMetadata
    kern::UInt64
    buf::Mem.Buffer
end

## module and function

const MAX_EXCEPTIONS = 256
const EXE_TO_MODULE_MAP = IdDict{Any,WeakRef}()
mutable struct ROCModule{E}
    exe::RuntimeExecutable{E}
    metadata::Vector{KernelMetadata}
    exceptions::Mem.Buffer
end
function ROCModule(exe)
    metadata = KernelMetadata[]
    exceptions = Mem.alloc(sizeof(ExceptionEntry)*MAX_EXCEPTIONS; coherent=true)
    mod = ROCModule(exe, metadata, exceptions)
    _exe = exe.exe
    EXE_TO_MODULE_MAP[_exe] = WeakRef(mod)
    hsaref!()
    finalizer(mod) do x
        # FIXME: Free all metadata
        Mem.free(mod.exceptions)
        delete!(EXE_TO_MODULE_MAP, _exe)
        hsaunref!()
    end
end
mutable struct ROCFunction
    mod::ROCModule
    entry::String
end

"""
    ROCDim3(x)
    ROCDim3((x,))
    ROCDim3((x, y))
    ROCDim3((x, y, x))

A type used to specify dimensions, consisting of 3 integers for the `x`, `y`,
and `z` dimension, respectively. Unspecified dimensions default to `1`.

Often accepted as argument through the `ROCDim` type alias, eg. in the case of
[`roccall`](@ref) or [`launch`](@ref), allowing to pass dimensions as a plain
integer or a tuple without having to construct an explicit `ROCDim3` object.
"""
struct ROCDim3
    x::Cuint
    y::Cuint
    z::Cuint
end

ROCDim3(dims::Integer)             = ROCDim3(dims,    Cuint(1), Cuint(1))
ROCDim3(dims::NTuple{1,<:Integer}) = ROCDim3(dims[1], Cuint(1), Cuint(1))
ROCDim3(dims::NTuple{2,<:Integer}) = ROCDim3(dims[1], dims[2],  Cuint(1))
ROCDim3(dims::NTuple{3,<:Integer}) = ROCDim3(dims[1], dims[2],  dims[3])

# Type alias for conveniently specifying the dimensions
# (e.g. `(len, 2)` instead of `ROCDim3((len, 2))`)
const ROCDim = Union{Integer,
                    Tuple{Integer},
                    Tuple{Integer, Integer},
                    Tuple{Integer, Integer, Integer}}

function Base.getindex(dims::ROCDim3, idx::Int)
    return idx == 1 ? dims.x :
           idx == 2 ? dims.y :
           idx == 3 ? dims.z :
           error("Invalid dimension: $idx")
end

"""
    roccall(f::ROCFunction, types, values...;
            queue::RuntimeQueue, signal::RuntimeEvent,
            groupsize::ROCDim, gridsize::ROCDim)

`ccall`-like interface for launching a ROC function `f` on a GPU.

For example:

    vadd = ROCFunction(md, "vadd")
    a = rand(Float32, 10)
    b = rand(Float32, 10)
    ad = Mem.upload(a)
    bd = Mem.upload(b)
    c = zeros(Float32, 10)
    cd = Mem.alloc(c)

    roccall(vadd, (Ptr{Cfloat},Ptr{Cfloat},Ptr{Cfloat}), ad, bd, cd;
             gridsize=10)
    Mem.download!(c, cd)

The `groupsize` and `gridsize` arguments control the launch configuration, and should both
consist of either an integer, or a tuple of 1 to 3 integers (omitted dimensions default to
1). The `types` argument can contain both a tuple of types, and a tuple type, the latter
being slightly faster.
"""
roccall

# we need a generated function to get a tuple of converted arguments (using unsafe_convert),
# without having to inspect the types at runtime
@generated function roccall(f::ROCFunction, tt::Type, args::Vararg{Any,N}; queue::RuntimeQueue,
                             signal::RuntimeEvent, groupsize::ROCDim=1, gridsize::ROCDim=groupsize) where N

    # the type of `tt` is Type{Tuple{<:DataType...}}
    types = tt.parameters[1].parameters

    ex = Expr(:block)
    push!(ex.args, :(Base.@_inline_meta))

    # convert the argument values to match the kernel's signature (specified by the user)
    # (this mimics `lower-ccall` in julia-syntax.scm)
    converted_args = Vector{Symbol}(undef, length(args))
    arg_ptrs = Vector{Symbol}(undef, length(args))
    for i in 1:length(args)
        converted_args[i] = gensym()
        arg_ptrs[i] = gensym()
        push!(ex.args, :($(converted_args[i]) = Base.cconvert($(types[i]), args[$i])))
        push!(ex.args, :($(arg_ptrs[i]) = Base.unsafe_convert($(types[i]), $(converted_args[i]))))
    end

    append!(ex.args, (quote
        #GC.@preserve $(converted_args...) begin
            launch(queue, signal, f, groupsize, gridsize, ($(arg_ptrs...),))
        #end
    end).args)

    return ex
end

"""
    launch(queue::RuntimeQueue, signal::RuntimeEvent, f::ROCFunction,
           groupsize::ROCDim, gridsize::ROCDim, args...)

Low-level call to launch a ROC function `f` on the GPU, using `groupsize` and
`gridsize` as the grid and block configuration, respectively. The kernel is
launched on `queue` and is waited on by `signal`.

Arguments to a kernel should either be bitstype, in which case they will be
copied to the internal kernel parameter buffer, or a pointer to device memory.

This is a low-level call, preferably use [`roccall`](@ref) instead.
"""
@inline function launch(queue::RuntimeQueue, signal::RuntimeEvent, f::ROCFunction,
                        groupsize::ROCDim, gridsize::ROCDim, args...)
    groupsize = ROCDim3(groupsize)
    gridsize = ROCDim3(gridsize)
    (groupsize.x>0 && groupsize.y>0 && groupsize.z>0) ||
        throw(ArgumentError("Group dimensions must be non-zero"))
    (gridsize.x>0 && gridsize.y>0 && gridsize.z>0) ||
        throw(ArgumentError("Grid dimensions must be non-zero"))
    (groupsize.x<=_max_group_size.x+1 &&
     groupsize.y<=_max_group_size.y+1 &&
     groupsize.z<=_max_group_size.z+1) ||
        throw(ArgumentError("Group dimensions too large"))

    _launch(queue, signal, f, groupsize, gridsize, args...)
end

## HSA object preservation while a kernel is active

const SIGNAL_PRESERVED = IdDict{HSASignal, Vector{Any}}()

function preserve!(sig::HSASignal, @nospecialize(x))
    set = get!(()->Any[], SIGNAL_PRESERVED, sig)
    push!(set, x)
end
preserve!(sig::HSAStatusSignal, @nospecialize(x)) = preserve!(sig.signal, x)
preserve!(ev::RuntimeEvent, @nospecialize(x)) = preserve!(ev.event, x)

unpreserve!(sig::HSASignal) = delete!(SIGNAL_PRESERVED, sig)
unpreserve!(sig::HSAStatusSignal) = unpreserve!(sig.signal)
unpreserve!(ev::RuntimeEvent) = unpreserve!(ev.event)

# we need a generated function to get an args array,
# without having to inspect the types at runtime
@generated function _launch(queue::RuntimeQueue, signal::RuntimeEvent, f::ROCFunction,
                            groupsize::ROCDim3, gridsize::ROCDim3,
                            args::NTuple{N,Any}) where N

    all(isbitstype, args.parameters) ||
        throw(ArgumentError("Arguments to kernel should be bitstype."))

    ex = Expr(:block)
    push!(ex.args, :(Base.@_inline_meta))

    # put arguments in Ref boxes so that we can get a pointers to them
    arg_refs = Vector{Symbol}(undef, N)
    for i in 1:N
        arg_refs[i] = gensym()
        push!(ex.args, :($(arg_refs[i]) = Base.RefValue(args[$i])))
    end

    append!(ex.args, (quote
        GC.@preserve $(arg_refs...) begin
            # create kernel instance
            kern = create_kernel(get_device(queue), f.mod.exe, f.entry, args)

            # launch kernel
            lock($RT_LOCK)
            try
                push!($_active_kernels[queue.queue], signal.event)
            finally
                unlock($RT_LOCK)
            end
            launch_kernel(queue, kern, signal;
                          groupsize=groupsize, gridsize=gridsize)

            # preserve kernel and arguments
            $preserve!(signal, kern)
            for arg in args
                $preserve!(signal, arg)
            end
        end
    end).args)

    return ex
end


text:
##
##!------------------------------------------------------------------------------------------------!
##  common.jl
##
##  Tests for common.jl
##
##  Creation date: 18/08/15
##  Author: Arindam Biswas <arindam dot b at eml dot cc>
##!------------------------------------------------------------------------------------------------!
##

ext: jl


text:
function weno3!(dU::Vector{Float64},k::Float64,A::Vector{Float64},Q::Vector{Float64},
    beta::Float64,rho::Float64,μ::Float64,diff::Float64,ts::Float64,xs::Float64,
    Qs::Float64,As::Float64,ISA::Vector{Float64},ISQ::Vector{Float64},
    PA::Vector{Function},PQ::Vector{Function},wA::Vector{Float64},wQ::Vector{Float64},
    c::Vector{Float64},a::Vector{Float64},p::Vector{Float64},x::Vector{Float64},
    Al::Vector{Float64},Ar::Vector{Float64},Ql::Vector{Float64},Qr::Vector{Float64},
    Hl::Vector{Float64},Hr::Vector{Float64},Fl::Vector{Float64},Fr::Vector{Float64},
    Pi::Vector{Float64},σ::Vector{Float64},eigs::Vector{Float64},abso::Vector{Float64},
    J::Matrix{Float64},Aint::Vector{Float64},Qint::Vector{Float64},Pret::Vector{Float64},
    s::Vector{Float64})

    # parameters
    ϵ = k^2;

    # allocation
    Aint[:] .= A[2:end-1];
    Qint[:] .= Q[2:end-1];

    # polynomial reconstruction loop (interior points only)
    # println("Reconstruction functions:")
    for j = 1:length(Aint)
        if j == 1
            c[1] = ϵ;
            c[2] = 0.25;
            c[3] = 0.75-ϵ;
        elseif j == length(A)
            c[1] = ϵ;
            c[2] = 0.25;
            c[3] = 0.75-ϵ
        else
            c[1] = 0.25;
            c[2] = 0.5;
            c[3] = 0.25;
        end
        smoothinds!(ISA,Aint,c,k,j,σ)
        smoothinds!(ISQ,Qint,c,k,j,σ)
        weights!(wA,c,ϵ,ISA,s)
        weights!(wQ,c,ϵ,ISQ,s)
        PA[j] = reconstruct(wA,c,Aint,j,k,σ)
        PQ[j] = reconstruct(wQ,c,Qint,j,k,σ)
    end

    # rhs loop
    # println("RHS loop:")
    for j = 1:length(Aint)
        if j == 1 # leftmost interior point
            # cell boundaries
            x[1] = 0.;
            x[2] = k;
            # reconstructed values at cell boundaries
            Al[1] = A[1];
            PA[j](Pret,x[1],Pi);
            Al[2] = Pret[1];
            PA[j](Pret,x[2],Pi);
            Ar[1] = Pret[1];
            PA[j+1](Pret,x[2],Pi);
            Ar[2] = Pret[1];
            Ql[1] = Q[1];
            PQ[j](Pret,x[1],Pi);
            Ql[2] = Pret[1];
            PQ[j](Pret,x[2],Pi);
            Qr[1] = Pret[1];
            PQ[j+1](Pret,x[2],Pi);
            Qr[2] = Pret[1];
            # numerical viscosity coefficient based on spectral radius of flux Jacobian
            abseigs1d(abso,eigs,J,Ql[1],Al[1],beta,rho,ts,xs,Qs,As);
            p[1] = maximum(abso);
            abseigs1d(abso,eigs,J,Ql[2],Al[2],beta,rho,ts,xs,Qs,As);
            p[2] = maximum(abso);
            a[1] = maximum(p);
            abseigs1d(abso,eigs,J,Qr[1],Ar[1],beta,rho,ts,xs,Qs,As);
            p[1] = maximum(abso);
            abseigs1d(abso,eigs,J,Qr[2],Ar[2],beta,rho,ts,xs,Qs,As);
            p[2] = maximum(abso);
            a[2] = maximum(p);
        elseif j == length(Aint) # rightmost interior point
            x[1] = k*(length(A)-3.);
            x[2] = k*(length(A)-2.);

            PA[j-1](Pret,x[1],Pi);
            Al[1] = Pret[1];
            PA[j](Pret,x[1],Pi);
            Al[2] = Pret[1];
            PA[j](Pret,x[2],Pi);
            Ar[1] = Pret[1];
            Ar[2] = A[end];
            PQ[j-1](Pret,x[1],Pi);
            Ql[1] = Pret[1];
            PQ[j](Pret,x[1],Pi);
            Ql[2] = Pret[1];
            PQ[j](Pret,x[2],Pi);
            Qr[1] = Pret[1];
            Qr[2] = Q[end];

            abseigs1d(abso,eigs,J,Ql[1],Al[1],beta,rho,ts,xs,Qs,As);
            p[1] = maximum(abso);
            abseigs1d(abso,eigs,J,Ql[2],Al[2],beta,rho,ts,xs,Qs,As);
            p[2] = maximum(abso);
            a[1] = maximum(p);
            abseigs1d(abso,eigs,J,Qr[1],Ar[1],beta,rho,ts,xs,Qs,As);
            p[1] = maximum(abso);
            abseigs1d(abso,eigs,J,Qr[2],Ar[2],beta,rho,ts,xs,Qs,As);
            p[2] = maximum(abso);
            a[2] = maximum(p);
        else # other interior points
            x[1] = k*(j-1);
            x[2] = k*j;

            PA[j-1](Pret,x[1],Pi);
            Al[1] = Pret[1];
            PA[j](Pret,x[1],Pi);
            Al[2] = Pret[1];
            PA[j](Pret,x[2],Pi);
            Ar[1] = Pret[1];
            PA[j+1](Pret,x[2],Pi);
            Ar[2] = Pret[1];
            PQ[j-1](Pret,x[1],Pi);
            Ql[1] = Pret[1];
            PQ[j](Pret,x[1],Pi);
            Ql[2] = Pret[1];
            PQ[j](Pret,x[2],Pi);
            Qr[1] = Pret[1];
            PQ[j+1](Pret,x[2],Pi);
            Qr[2] = Pret[1];

            abseigs1d(abso,eigs,J,Ql[1],Al[1],beta,rho,ts,xs,Qs,As);
            p[1] = maximum(abso);
            abseigs1d(abso,eigs,J,Ql[2],Al[2],beta,rho,ts,xs,Qs,As);
            p[2] = maximum(abso);
            a[1] = maximum(p);
            abseigs1d(abso,eigs,J,Qr[1],Ar[1],beta,rho,ts,xs,Qs,As);
            p[1] = maximum(abso);
            abseigs1d(abso,eigs,J,Qr[2],Ar[2],beta,rho,ts,xs,Qs,As);
            p[2] = maximum(abso);
            a[2] = maximum(p);
        end
        # numerical fluxes (local Lax-Friedrichs)
        F1d(Fl,Ql[1],Al[1],beta,rho,ts,xs,Qs,As)
        F1d(Fr,Ql[2],Al[2],beta,rho,ts,xs,Qs,As)
        Hl[1] = 0.5*(Fl[1] + Fr[1] - a[1]*(Al[2] - Al[1]))
        Hl[2] = 0.5*(Fl[2] + Fr[2] - a[1]*(Ql[2] - Ql[1]))
        F1d(Fr,Qr[2],Ar[2],beta,rho,ts,xs,Qs,As)
        F1d(Fl,Qr[1],Ar[1],beta,rho,ts,xs,Qs,As)
        Hr[1] = 0.5*(Fl[1] + Fr[1] - a[2]*(Ar[2] - Ar[1]))
        Hr[2] = 0.5*(Fl[2] + Fr[2] - a[2]*(Qr[2] - Qr[1]))
        # 3rd-order WENO spatial derivatives (plus source terms)
        dU[j] = -1./k*(Hr[1] - Hl[1]);
        dU[length(Aint)+j] = -1./k*(Hr[2] - Hl[2]) - ts/As*diff*π*μ/rho*Q[j+1]/A[j+1];
    end
    # println("Calculated dU: $dU")
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
pushfirst!(LOAD_PATH, joinpath(@__DIR__, ".."))
using Documenter
using LaTeXStrings
using RRTMGP

makedocs(
    sitename = "RRTMGP.jl",
    doctest = false,
    strict = false,
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        mathengine = MathJax(Dict(:TeX => Dict(:equationNumbers => Dict(:autoNumber => "AMS"), :Macros => Dict()))),
    ),
    clean = true,
    modules = [RRTMGP],
    pages = Any[
        "Home" => "index.md",
        "Mathematical Formulation" => "MathFormulation.md",
        "Optics" => Any[
            "Angular Discretization" => "optics/AngularDiscretizations.md"
            "Atmospheric State" => "optics/AtmosphericStates.md"
            "Boundary Conditions" => "optics/BCs.md"
            "Fluxes" => "optics/Fluxes.md"
            "Utilities for gray radiation simulation" => "optics/GrayUtils.md"
            "Lookup Tables" => "optics/LookUpTables.md"
            "Optics" => "optics/Optics.md"
            "Solver struct" => "optics/RTE.md"
            "Source Functions" => "optics/Sources.md"
            "Volume Mixing Ratios" => "optics/Vmrs.md"
        ],
        "RTE Solver" => "RTESolver.md",
        "References" => "References.md",
    ],
)

deploydocs(
    repo = "github.com/CliMA/RRTMGP.jl.git",
    target = "build",
    push_preview = true,
    devbranch = "main",
    forcepush = true,
)

max_stars_repo_name: climate-machine/JRRTMGP

max_issues_repo_name: climate-machine/JRRTMGP

max_forks_repo_name: climate-machine/JRRTMGP


text:
# ---
# title: 857. Minimum Cost to Hire K Workers
# id: problem857
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Heap
# link: <https://leetcode.com/problems/minimum-cost-to-hire-k-workers/description/>
# hidden: true
# ---
#
# There are `N` workers.  The `i`-th worker has a `quality[i]` and a minimum
# wage expectation `wage[i]`.
#
# Now we want to hire exactly `K` workers to form a _paid group_.  When hiring a
# group of K workers, we must pay them according to the following rules:
#
#   1. Every worker in the paid group should be paid in the ratio of their quality compared to other workers in the paid group.
#   2. Every worker in the paid group must be paid at least their minimum wage expectation.
#
# Return the least amount of money needed to form a paid group satisfying the
# above conditions.
#
#
#
# **Example 1:**
#
#
#
#     Input: quality = [10,20,5], wage = [70,50,30], K = 2
#     Output: 105.00000
#     **Explanation** : We pay 70 to 0-th worker and 35 to 2-th worker.
#
#
# **Example 2:**
#
#
#
#     Input: quality = [3,1,10,10,1], wage = [4,8,2,2,7], K = 3
#     Output: 30.66667
#     **Explanation** : We pay 4 to 0-th worker, 13.33333 to 2-th and 3-th workers seperately.
#
#
#
#
# **Note:**
#
#   1. `1 <= K <= N <= 10000`, where `N = quality.length = wage.length`
#   2. `1 <= quality[i] <= 10000`
#   3. `1 <= wage[i] <= 10000`
#   4. Answers within `10^-5` of the correct answer will be considered correct.
#
#
## @lc code=start
using LeetCode

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


text:
# Autogenerated wrapper script for Mineos_jll for armv7l-linux-musleabihf-libgfortran5
export cucss2sac, eigcon, eigen2asc, endi, green, minos_bran, simpledit, syndat

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Mineos")
JLLWrappers.@declare_executable_product(cucss2sac)
JLLWrappers.@declare_executable_product(eigcon)
JLLWrappers.@declare_executable_product(eigen2asc)
JLLWrappers.@declare_executable_product(endi)
JLLWrappers.@declare_executable_product(green)
JLLWrappers.@declare_executable_product(minos_bran)
JLLWrappers.@declare_executable_product(simpledit)
JLLWrappers.@declare_executable_product(syndat)
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_executable_product(
        cucss2sac,
        "bin/cucss2sac",
    )

    JLLWrappers.@init_executable_product(
        eigcon,
        "bin/eigcon",
    )

    JLLWrappers.@init_executable_product(
        eigen2asc,
        "bin/eigen2asc",
    )

    JLLWrappers.@init_executable_product(
        endi,
        "bin/endi",
    )

    JLLWrappers.@init_executable_product(
        green,
        "bin/green",
    )

    JLLWrappers.@init_executable_product(
        minos_bran,
        "bin/minos_bran",
    )

    JLLWrappers.@init_executable_product(
        simpledit,
        "bin/simpledit",
    )

    JLLWrappers.@init_executable_product(
        syndat,
        "bin/syndat",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Test
using StationXML

@testset "Types" begin
    # Types with restrictions on what values they can take
    @testset "Restricted types" begin
        # Types with optional unit but which can take only one value
        @testset "Units" begin
            for (T, unit) in (
                    (StationXML.Second, "SECONDS"),
                    (StationXML.Voltage, "VOLTS"),
                    (StationXML.Angle, "DEGREES"),
                    (StationXML.Azimuth, "DEGREES"),
                    (StationXML.Dip, "DEGREES"),
                    (StationXML.ClockDrift, "SECONDS/SAMPLE"),
                    (StationXML.Frequency, "HERTZ"),
                    (StationXML.SampleRate, "SAMPLES/S"),
                    (StationXML.Latitude, "DEGREES"),
                    (StationXML.Longitude, "DEGREES"))
                obj = T(0)
                @test obj.unit === missing
                obj.unit = unit
                @test obj.unit == unit
                @test_throws ArgumentError obj.unit = "invalid unit string"
                @test_throws ArgumentError T(; value=0, unit="invalid unit string")
            end
        end

        # Types which can take only a range of values
        @testset "Value ranges" begin
            # Float- and FloatNoUnit-derived types
            for (T, vals) in (
                (StationXML.Angle, (-361, 361)),
                (StationXML.Azimuth, (-1, 361)),
                (StationXML.Dip, (-91, 91)),
                (StationXML.ClockDrift, (-1,)),
                (StationXML.Latitude, (-91, 91)),
                (StationXML.Longitude, (-361, 361)),
                    )
                for val in vals
                    @test_throws ArgumentError T(val)
                    t = T(0)
                    @test_throws ArgumentError t.value = val
                end
            end
            # Enumerated types
            @test_throws ArgumentError StationXML.Nominal("WEIRD_VALUE")
            @test_throws ArgumentError StationXML.Email("not_an_email_address")
            @test_throws ArgumentError StationXML.PhoneNumber(area_code=1,
                phone_number="not a phone number")
            @test_throws ArgumentError StationXML.RestrictedStatus("WEIRD_VALUE")
            # Network
            @test_throws ArgumentError StationXML.Network(code="AN", total_number_stations=-1)
            @test_throws ArgumentError StationXML.Network(code="AN", selected_number_stations=-1)
            @test_throws ArgumentError StationXML.Coefficient(value=1, number=-1)
        end
    end

    # Types which only permit certain combinations of fields to be set
    @testset "Field combinations" begin
        # All or none of the frequency_* fields for Sensitivity
        let input_units = StationXML.Units("M/S"), output_units = StationXML.Units("V")
            @test_throws ArgumentError StationXML.Sensitivity(value=1, frequency=0,
                input_units=input_units, output_units=output_units,
                frequency_start=1)
            @test_throws ArgumentError StationXML.Sensitivity(value=1, frequency=0,
                input_units=input_units, output_units=output_units,
                frequency_start=1, frequency_end=2)
            @test_throws ArgumentError StationXML.Sensitivity(value=1, frequency=0,
                input_units=input_units, output_units=output_units,
                frequency_end=1, frequency_db_variation=2)
        end
        # ResponseStage stage number and only one type of function
        let iu = StationXML.Units("M/S"), ou = StationXML.Units("V"),
                    pz = StationXML.PolesZeros(input_units=iu,
                                               output_units=ou,
                                               pz_transfer_function_type="LAPLACE (HERTZ)",
                                               normalization_frequency=1),
                coefs = StationXML.Coefficients(input_units=iu, output_units=ou,
                                                cf_transfer_function_type="ANALOG (HERTZ)"),
                resplist = StationXML.ResponseList(input_units=iu, output_units=ou),
                fir = StationXML.FIR(input_units=iu, output_units=ou, symmetry="NONE"),
                poly = StationXML.Polynomial(input_units=iu, output_units=ou,
                                             frequency_lower_bound=StationXML.Frequency(0),
                                             frequency_upper_bound=StationXML.Frequency(1),
                                             approximation_lower_bound=0,
                                             approximation_upper_bound=1,
                                             maximum_error=1)
            # Only one type allowed
            all_args = (:poles_zeros=>pz, :coefficients=>coefs,
                    :response_list=>resplist, :fir=>fir, :polynomial=>poly)
            for args in (all_args[1:2], all_args[1:3], all_args[1:4], all_args)
                @test_throws ArgumentError StationXML.ResponseStage(; number=1, args...)
            end
            for arg in all_args
                @test StationXML.ResponseStage(; (:number=>1, arg)...) isa StationXML.ResponseStage
            end
            # Stage number >= 1
            @test_throws ArgumentError StationXML.ResponseStage(poles_zeros=pz, number=-1)
            # No stage gain with polynomial
            @test_logs (:warn, "stage_gain cannot be specified for polynomial stages in StationXML v1.1"
                ) StationXML.ResponseStage(polynomial=poly, number=1,
                    stage_gain=StationXML.Gain(1, 1))
        end
    end
end


text:
@defcomp PAGE09SeaLevelRise begin

  # Parameters

    rt_g_globaltemperature = Parameter(index=[time], unit="degreeC")
    sltemp_SLtemprise = Parameter(unit="m-degreeC", default=1.7333333333333334)
    sla_SLbaselinerise = Parameter(unit="m", default=1.00)
    sltau_SLresponsetime = Parameter(unit="years", default=1000.)
    s0_initialSL = Parameter(unit="m", default=0.15)
    y_year = Parameter(index=[time], unit="year")
    y_year_0 = Parameter(unit="year")

  # Variables

    es_equilibriumSL = Variable(index=[time], unit="m")
    s_sealevel = Variable(index=[time], unit="m")
    expfs_exponential = Variable(index=[time], unit="unitless")
    yp_timestep = Variable(index=[time], unit="years")

    function run_timestep(p, v, d, t)
        s0_initialSL_2009 = quantile(TriangularDist(0.1, 0.2, 0.15), cdf(TriangularDist(0.17, 0.21, 0.19), p.s0_initialSL))
        sltau_SLresponsetime_2009 = quantile(TriangularDist(500, 1500, 1000), cdf(Gamma(16.0833333333333333, 24.), p.sltau_SLresponsetime))

        if is_first(t)
            v.yp_timestep[t] = p.y_year[1] - p.y_year_0
            v.es_equilibriumSL[t] = p.sltemp_SLtemprise * p.rt_g_globaltemperature[t] + p.sla_SLbaselinerise
            v.expfs_exponential[t] = exp(-v.yp_timestep[t] / sltau_SLresponsetime_2009)
            v.s_sealevel[t] = s0_initialSL_2009 + (v.es_equilibriumSL[t] - s0_initialSL_2009) * (1 - v.expfs_exponential[t])

        else
            v.yp_timestep[t] = p.y_year[t] - p.y_year[t - 1]
            v.es_equilibriumSL[t] = p.sltemp_SLtemprise * p.rt_g_globaltemperature[t] + p.sla_SLbaselinerise
            v.expfs_exponential[t] = exp(-v.yp_timestep[t] / sltau_SLresponsetime_2009)
            v.s_sealevel[t] = v.s_sealevel[t - 1] + (v.es_equilibriumSL[t] - v.s_sealevel[t - 1]) * (1 - v.expfs_exponential[t])

        end
    end
end


text:
include("../various_mle_things/incremental.jl")

function unitvec(n, i)
    v = zeros(n, 1)
    v[i] = 1
    return v
end

"""
compute determinant of submatrix of Σθ formed by deleting ith row and column efficiently
given determinant of full matrix
"""
#function det_loocv(holeskyΣθ::IncrementalCholeskyc, detcholeskyΣθ::Float64, i::Int64)
function det_loocv(choleskyΣθ, detcholeskyΣθ::Float64, i::Int64)
    n = size(choleskyΣθ)[1]
    #println("i: ", i)
    @assert i >= 1 && n >= i
    #R = get_chol(choleskyΣθ)
    R = choleskyΣθ.U
    return detcholeskyΣθ * norm(R'\ unitvec(n, i))^2
end

"""
compute determinant of submatrix of Σθ formed by deleting ith row and column efficiently
given determinant of full matrix
"""
#function det_loocv(holeskyΣθ::IncrementalCholeskyc, detcholeskyΣθ::Float64, i::Int64)
function det_loocv_IC(choleskyΣθ::IncrementalCholesky, detcholeskyΣθ::Float64, i::Int64)
    n = size(choleskyΣθ)[1]
    #println("i: ", i)
    @assert i >= 1 && n >= i
    #R = get_chol(choleskyΣθ)
    L = get_chol(choleskyΣθ).L
    return detcholeskyΣθ * norm(L\ unitvec(n, i))^2
end

"""
stably compute log determinant of submatrix of Σθ formed by deleting ith row and column efficiently
given determinant of full matrix
"""
#function det_loocv(holeskyΣθ::IncrementalCholeskyc, detcholeskyΣθ::Float64, i::Int64)
function logdet_loocv_IC(choleskyΣθ::IncrementalCholesky, logdetcholeskyΣθ::Float64, i::Int64)
    n = size(choleskyΣθ)[1]
    #println("i: ", i)
    @assert i >= 1 && n >= i
    #R = get_chol(choleskyΣθ)
    L = get_chol(choleskyΣθ).L
    return logdetcholeskyΣθ + 2*log(norm(L\ unitvec(n, i)))
end

function logdet_loocv(choleskyΣθ::Cholesky{Float64,Array{Float64,2}}, logdetcholeskyΣθ::Float64, i::Int64)
    n = size(choleskyΣθ)[1]
    #println("i: ", i)
    @assert i >= 1 && n >= i
    #R = get_chol(choleskyΣθ)
    R = choleskyΣθ.U
    return logdetcholeskyΣθ + 2*log(norm(R'\ unitvec(n, i)))
end

"""
Compute determinant of submatrix of - [Sigmainv X; X' 0]
"""
function det_XΣX_loocv(X, choleskyΣθ, detcholeskyΣθ::Float64, i)
    n = size(choleskyΣθ, 1)
    p = size(X, 2)
    @assert n == size(X, 1)
    @assert i >=1 && n>=i
    ei = unitvec(n, i)
    #R11 = choleskyΣθ.U #upper triangular factor
    R12 = choleskyΣθ.L\X
    cholR22 = cholesky(R12'*R12)
    a = choleskyΣθ.L\ei
    b = - (cholR22.L \ (R12'*a))
    eWe =  norm(a)^2 - norm(b)^2
    detW = (-1)^n * detcholeskyΣθ * det(cholR22)
    det_Σθ_minus_i = det_loocv(choleskyΣθ, detcholeskyΣθ, i)
    return eWe * detW / det_Σθ_minus_i
end

"""
Compute determinant of submatrix of - [Sigmainv X; X' 0]. IC stands for IncrementalCholesky
"""
function det_XΣX_loocv_IC(X, choleskyΣθ::IncrementalCholesky, detcholeskyΣθ::Float64, i)
    n = size(choleskyΣθ, 1)
    p = size(X, 2)
    @assert n == size(X, 1)
    @assert i >=1 && n>=i
    chol = get_chol(choleskyΣθ)
    ei = unitvec(n, i)
    #R11 = choleskyΣθ.U #upper triangular factor
    R12 = get_chol(choleskyΣθ).L\X
    cholR22 = cholesky(R12'*R12)
    a = get_chol(choleskyΣθ).L\ei
    b = - (cholR22.L \ (R12'*a))
    eWe =  norm(a)^2 - norm(b)^2
    detW = (-1)^n * detcholeskyΣθ * det(cholR22)
    det_Σθ_minus_i = det_loocv_IC(choleskyΣθ, detcholeskyΣθ, i)
    return eWe * detW / det_Σθ_minus_i
end

"""
Stably compute log of absolute value of determinant of submatrix of - [Sigmainv X; X' 0], where
the determinant is computed up to a sign, and the log of the absolute value of
the determinant is taken
"""
function logdet_XΣX_loocv(X, choleskyΣθ, logdetcholeskyΣθ::Float64, i)
    n = size(choleskyΣθ, 1)
    p = size(X, 2)
    @assert n == size(X, 1)
    @assert i >=1 && n>=i
    ei = unitvec(n, i)
    #R11 = choleskyΣθ.U #upper triangular factor
    R12 = choleskyΣθ.L\X
    cholR22 = cholesky(R12'*R12)
    a = choleskyΣθ.L\ei
    b = - (cholR22.L \ (R12'*a))
    #println("norm a squared: ", norm(a)^2)
    #println("norm b squared: ", norm(b)^2)
    logeWe =  log(abs(norm(a)^2 - norm(b)^2))
    #println("logeWe: ", logeWe)
    #logdetW = (-1)^n * detcholeskyΣθ * det(cholR22)
    #println("logdet(cholR22): ", logdet(cholR22))
    logdetW = logdetcholeskyΣθ + logdet(cholR22)
    #println("logdetW: ", logdetW)
    logdet_Σθ_minus_i = logdet_loocv(choleskyΣθ, logdetcholeskyΣθ, i) #submatrix of PD matrix is PD
    return logeWe + logdetW - logdet_Σθ_minus_i
end

function logdet_XΣX_loocv_IC(X, choleskyΣθ::IncrementalCholesky, logdetcholeskyΣθ::Float64, i)
    n = size(choleskyΣθ)[1]
    p = size(X, 2)
    @assert n == size(X, 1)
    @assert i >=1 && n>=i
    ei = unitvec(n, i)
    #R11 = choleskyΣθ.U #upper triangular factor
    R12 = get_chol(choleskyΣθ).L\X
    cholR22 = cholesky(R12'*R12)
    a = get_chol(choleskyΣθ).L\ei
    b = - (cholR22.L \ (R12'*a))
    #println("norm a squared: ", norm(a)^2)
    #println("norm b squared: ", norm(b)^2)
    logeWe =  log(abs(norm(a)^2 - norm(b)^2))
    #println("logeWe: ", logeWe)
    #logdetW = (-1)^n * detcholeskyΣθ * det(cholR22)
    #println("logdet(cholR22): ", logdet(cholR22))
    logdetW = logdetcholeskyΣθ + logdet(cholR22)
    #println("logdetW: ", logdetW)
    logdet_Σθ_minus_i = logdet_loocv_IC(choleskyΣθ, logdetcholeskyΣθ, i) #submatrix of PD matrix is PD
    return logeWe + logdetW - logdet_Σθ_minus_i
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using DiffEqBase, DiffEqJump, DiffEqBiological
using BenchmarkTools

doBenchmark = true

# number of save times to look at
num_save_vec = 25000

# diffusion model
N = 32
println("N = ", N)
diffnetwork = "@reaction_network difftype begin\n"
for i in 1:(N-1)
    diffnetwork *= "\t K, X$(i) --> X$(i+1)\n"
    diffnetwork *= "\t K, X$(i+1) --> X$(i)\n"
end
diffnetwork *= "end K"
rs = eval( parse(diffnetwork) )
tf = 10.
params = (1.,)
prob   = DiscreteProblem(10*ones(Int64,N), (0.0, tf), params)

methods = (Direct(),)

if doBenchmark
    # methods with different numbers of points
    println("Saving at fixed times:\n")
    for num_save in num_save_vec
        dtsave = tf / num_save

        for method in methods
            println("Solving with method: ", typeof(method), ", using SSAStepper and ", num_save, " points")
            jump_prob = JumpProblem(prob, method, rs, save_positions=(false,false))
            sol = solve(jump_prob, SSAStepper(), saveat=dtsave);
            @btime solve($jump_prob, SSAStepper(), saveat=$dtsave)
        end

        println()
    end
end

max_stars_repo_name: isaacsas/DiffEqJumpExtensions

max_issues_repo_name: isaacsas/DiffEqJumpExtensions

max_forks_repo_name: isaacsas/DiffEqJumpExtensions


text:
using GridOperations; const GO = GridOperations
using Base.Test
using ColorTypes

@testset "Colorize" begin
    arr = [3.0, 4.0, 5.0, NaN]
    colored = GO.colorize(arr)
    @test colored == [
        RGB{Float64}(0.267004, 0.004874, 0.329415),
        RGB{Float64}(0.1281485, 0.565107, 0.5508925),
        RGB{Float64}(0.993248, 0.906157, 0.143936),
        RGB{Float64}(1.0, 1.0, 1.0)
    ]
end

@testset "Progressive morphological filter" begin
    srand(1)
    zmin = rand(3.0f0:0.1f0:5.0f0, 30, 40)
    max_window_radius = 9.0
    slope = 0.151234
    dh_max = 2.0
    dh_min = 0.4
    nodata = -9999.0
    boundarymask = trues(zmin)
    zmax_pmf, flags = GO.pmf_filter(zmin, boundarymask, max_window_radius, slope, dh_max, dh_min, 1.0)
    @test count(flags) === 754
    @test minimum(zmax_pmf) ≈ 3.4
    @test maximum(zmax_pmf) ≈ 3.702468  # large Float32 inaccuracy, should be 3.7
    @test hash(zmax_pmf) === 0x520f02f7bbe47b0e
end


text:
include("setup.jl")

@testset "Downloads.jl" begin
    @testset "libcurl configuration" begin
        julia = "$(VERSION.major).$(VERSION.minor)"
        @test Curl.USER_AGENT == "curl/$(Curl.CURL_VERSION) julia/$julia"
        if VERSION > v"1.6-"
            @test Curl.SYSTEM_SSL == Sys.iswindows() | Sys.isapple()
        end
    end

    @testset "API coverage" begin
        value = "Julia is great!"
        base64 = "SnVsaWEgaXMgZ3JlYXQh"
        url = "$server/base64/$base64"
        headers = ["Foo" => "Bar"]
        # test with one argument
        path = download(url)
        @test isfile(path)
        @test value == read(path, String)
        rm(path)
        # with headers
        path = download(url, headers=headers)
        @test isfile(path)
        @test value == read(path, String)
        rm(path)
        # test with two arguments
        arg_writers() do path, output
            @arg_test output begin
                @test output == download(url, output)
            end
            @test isfile(path)
            @test value == read(path, String)
            rm(path)
            # with headers
            @arg_test output begin
                @test output == download(url, output, headers=headers)
            end
            @test isfile(path)
            @test value == read(path, String)
            rm(path)
        end

        # not an API test, but a convenient place to test this
        @testset "follow redirects" begin
            redirect = "$server/redirect-to?url=$(url_escape(url))"
            path = download(redirect)
            @test isfile(path)
            @test value == read(path, String)
            rm(path)
        end
    end

    @testset "get request" begin
        url = "$server/get"
        json = download_json(url)
        @test json["url"] == url
        resp = request(url)
        @test resp isa Response
        @test resp.proto == "https"
        @test resp.status == 200
    end

    # https://github.com/JuliaLang/Downloads.jl/issues/131
    @testset "head request" begin
        url = server * "/image/jpeg"
        output = IOBuffer()
        resp = request(url; method="HEAD", output=output)
        @test resp isa Response
        @test resp.proto == "https"
        @test resp.status == 200
        @test isempty(take!(output)) # no output from a `HEAD`
        len = parse(Int, Dict(resp.headers)["content-length"])

        # when we make a `GET` instead of a `HEAD`, we get a body with the content-length
        # returned from the `HEAD` request.
        resp = request(url; method="GET", output=output)
        bytes = take!(output)
        @test length(bytes) == len
    end

    @testset "put request" begin
        url = "$server/put"
        data = "Hello, world!"
        resp, json = request_json(url, input=IOBuffer(data))
        @test json["url"] == url
        @test json["data"] == data
    end

    @testset "post request" begin
        url = "$server/post"
        data = "Hello, world!"
        resp, json = request_json(url, input=IOBuffer(data), method="POST")
        @test json["url"] == url
        @test json["data"] == data
    end

    @testset "put from file" begin
        url = "$server/put"
        file = tempname()
        write(file, "Hello, world!")
        resp, json = request_json(url, input=file)
        @test json["url"] == url
        @test json["data"] == read(file, String)
        rm(file)
    end

    @testset "redirected get" begin
        url = "$server/get"
        redirect = "$server/redirect-to?url=$(url_escape(url))"
        json = download_json(url)
        @test json["url"] == url
    end

    @testset "redirected put" begin
        url = "$server/put"
        redirect = "$server/redirect-to?url=$(url_escape(url))"
        data = "Hello, world!"
        resp, json = request_json(redirect, input=IOBuffer(data))
        @test json["url"] == url
        @test json["data"] == data
    end

    @testset "redirected post" begin
        url = "$server/post"
        redirect = "$server/redirect-to?url=$(url_escape(url))"
        data = "Hello, world!"
        resp, json = request_json(redirect, input=IOBuffer(data), method="POST")
        @test json["url"] == url
        @test json["data"] == data
    end

    @testset "redirected put from file" begin
        url = "$server/put"
        redirect = "$server/redirect-to?url=$(url_escape(url))"
        file = tempname()
        write(file, "Hello, world!")
        resp, json = request_json(redirect, input=file)
        @test json["url"] == url
        @test json["data"] == read(file, String)
    end

    @testset "headers" begin
        url = "$server/headers"

        @testset "set headers" begin
            headers = ["Foo" => "123", "Header" => "VaLuE", "Empty" => ""]
            json = download_json(url, headers = headers)
            for (key, value) in headers
                @test header(json["headers"], key) == value
            end
            @test header(json["headers"], "Accept") == "*/*"
            @test header(json["headers"], "User-Agent") == Curl.USER_AGENT
        end

        @testset "override default header" begin
            headers = [
                "Accept"     => "application/tar"
                "User-Agent" => "MyUserAgent/1.0"
            ]
            json = download_json(url, headers = headers)
            @test header(json["headers"], "Accept") == "application/tar"
            @test header(json["headers"], "User-Agent") == "MyUserAgent/1.0"
        end

        @testset "override default header with empty value" begin
            headers = [
                "Accept"     => ""
                "User-Agent" => ""
            ]
            json = download_json(url, headers = headers)
            @test header(json["headers"], "Accept") == ""
            @test header(json["headers"], "User-Agent") == ""
        end

        @testset "delete default header" begin
        headers = [
            "Accept"     => nothing
            "User-Agent" => nothing
        ]
        json = download_json(url, headers = headers)
            @test !("Accept" in keys(json["headers"]))
            @test !("User-Agent" in keys(json["headers"]))
        end

        @testset "HTTP/2 user agent bug" begin
            json = download_json(url)
            @test header(json["headers"], "User-Agent") == Curl.USER_AGENT
            @sync for _ = 1:2
                @async begin
                    json = download_json(url)
                    @test header(json["headers"], "User-Agent") == Curl.USER_AGENT
                end
            end
        end
    end

    @testset "session support" begin
        downloader = Downloader()

        # This url will redirect to /cookies, which echoes the set cookies as json
        set_cookie_url = "$server/cookies/set?k1=v1&k2=v2"
        cookies = download_json(set_cookie_url, downloader=downloader)
        @test get(cookies, "k1", "") == "v1"
        @test get(cookies, "k2", "") == "v2"

        # As the handle is destroyed, subsequent requests have no cookies
        cookie_url = "$server/cookies"
        cookies = download_json(cookie_url, downloader=downloader)
        @test isempty(cookies)
    end

    @testset "netrc support" begin
        user = "gVvkQiHN62"
        passwd = "dlctfSMTno8n"
        auth_url = "$server/basic-auth/$user/$passwd"
        resp = request(auth_url)
        @test resp isa Response
        @test resp.status == 401  # no succesful authentication

        # Setup .netrc
        hostname = match(r"^\w+://([^/]+)"i, server).captures[1]
        netrc = tempname()
        open(netrc, "w") do io
            write(io, "machine $hostname login $user password $passwd\n")
        end

        # Setup config to point to custom .netrc (normally in ~/.netrc)
        downloader = Downloads.Downloader()
        downloader.easy_hook = (easy, info) ->
            Curl.setopt(easy, Curl.CURLOPT_NETRC_FILE, netrc)

        resp = request(auth_url, throw=false, downloader=downloader)
        @test resp isa Response
        @test resp.status == 200  # succesful authentication

        # Cleanup
        rm(netrc)
    end

    @testset "file protocol" begin
        @testset "success" begin
            path = tempname()
            data = rand(UInt8, 256)
            write(path, data)
            temp = download("file://$path")
            @test data == read(temp)
            output = IOBuffer()
            resp = request("file://$path", output = output)
            @test resp isa Response
            @test resp.proto == "file"
            @test resp.status == 0
            @test take!(output) == data
            rm(path)
        end
        @testset "failure" begin
            path = tempname()
            @test_throws RequestError download("file://$path")
            @test_throws RequestError request("file://$path")
            output = IOBuffer()
            resp = request("file://$path", output = output, throw = false)
            @test resp isa RequestError
        end
    end

    @testset "errors" begin
        @test_throws ArgumentError download("ba\0d")
        @test_throws ArgumentError download("good", "ba\0d")

        err = @exception download("xyz://domain.invalid")
        @test err isa RequestError
        @test err.code != 0
        @test startswith(err.message, "Protocol \"xyz\" not supported")
        @test err.response.proto === nothing

        err = @exception request("xyz://domain.invalid", input = IOBuffer("Hi"))
        @test err isa RequestError
        @test err.code != 0
        @test startswith(err.message, "Protocol \"xyz\" not supported")
        @test err.response.proto === nothing

        err = @exception download("https://domain.invalid")
        @test err isa RequestError
        @test err.code != 0
        @test startswith(err.message, "Could not resolve host")
        @test err.response.proto === nothing

        err = @exception request("https://domain.invalid", input = IOBuffer("Hi"))
        @test err isa RequestError
        @test err.code != 0
        @test startswith(err.message, "Could not resolve host")
        @test err.response.proto === nothing

        err = @exception download("$server/status/404")
        @test err isa RequestError
        @test err.code == 0 && isempty(err.message)
        @test err.response.status == 404
        @test contains(err.response.message, r"^HTTP/\d+(?:\.\d+)?\s+404\b")
        @test err.response.proto === "https"

        resp = request("$server/get", input = IOBuffer("Hi"))
        @test resp isa Response
        @test resp.status == 405
        @test contains(resp.message, r"^HTTP/\d+(?:\.\d+)?\s+405\b")
        @test err.response.proto === "https"

        path = tempname()
        @test_throws RequestError download("$server/status/404", path)
        @test !ispath(path)
    end

    @testset "concurrent requests" begin
        mine = Downloader()
        for downloader in (nothing, mine)
            have_lsof = Sys.which("lsof") !== nothing
            count_tcp() = Base.count(x->contains("TCP",x), split(read(`lsof -p $(getpid())`, String), '\n'))
            if have_lsof
                n_tcp = count_tcp()
            end
            delay = 2
            count = 100
            url = "$server/delay/$delay"
            t = @elapsed @sync for id = 1:count
                @async begin
                    json = download_json("$url?id=$id", downloader = downloader)
                    @test get(json["args"], "id", nothing) == ["$id"]
                end
            end
            @test t < 0.9*count*delay
            if have_lsof
                @test n_tcp == count_tcp()
            end
        end
    end

    @testset "request API" begin
        @testset "basic request usage" begin
            for status in (200, 300, 400)
                url = "$server/status/$status"
                resp, body = request_body(url)
                @test resp.url == url
                @test resp.status == status
                test_response_string(resp.message, status)
                @test all(hdr isa Pair{String,String} for hdr in resp.headers)
                headers = Dict(resp.headers)
                @test "content-length" in keys(headers)
            end
        end

        @testset "custom headers" begin
            url = "$server/response-headers?FooBar=VaLuE"
            resp, body = request_body(url)
            @test resp.url == url
            @test resp.status == 200
            test_response_string(resp.message, 200)
            headers = Dict(resp.headers)
            @test headers["foobar"] == "VaLuE"
        end

        @testset "url for redirect" begin
            url = "$server/get"
            redirect = "$server/redirect-to?url=$(url_escape(url))"
            resp, json = request_json(redirect)
            @test resp.url == url
            @test resp.status == 200
            test_response_string(resp.message, 200)
            @test json["url"] == url
        end

        @testset "timeouts" begin
            url = "$server/delay/2"
            @testset "download" begin
                @test_throws ArgumentError download(url, devnull, timeout = -1)
                @test_throws ArgumentError download(url, devnull, timeout = 0)
                @test_throws RequestError download(url, devnull, timeout = 1)
                @test_throws RequestError download(url, devnull, timeout = 0.5)
                @test download(url, devnull, timeout = 100) == devnull
                @test download(url, devnull, timeout = Inf) == devnull
            end
            @testset "request(throw = true)" begin
                @test_throws ArgumentError request(url, timeout = -1)
                @test_throws ArgumentError request(url, timeout = 0)
                @test_throws RequestError request(url, timeout = 1)
                @test_throws RequestError request(url, timeout = 0.5)
                @test request(url, timeout = 100) isa Response
                @test request(url, timeout = Inf) isa Response
            end
            @testset "request(throw = false)" begin
                @test request(url, throw = false, timeout = 1) isa RequestError
                @test request(url, throw = false, timeout = 0.5) isa RequestError
                @test request(url, throw = false, timeout = 100) isa Response
                @test request(url, throw = false, timeout = Inf) isa Response
            end
        end

        @testset "progress" begin
            url = "https://httpbingo.org/drip"
            progress = []
            dl_funcs = [
                download,
                (url; progress) ->
                    request(url, output=devnull, progress=progress)
            ]
            p_funcs = [
                (prog...) -> push!(progress, prog),
                (total, now) -> push!(progress, (total, now)),
                (total, now, _, _) -> push!(progress, (total, now)),
            ]
            for f in dl_funcs, p in p_funcs
                @testset "request" begin
                    empty!(progress)
                    f(url; progress = p)
                    @test progress[1][1] == 0
                    @test progress[1][2] == 0
                    @test progress[end][1] == 10
                    @test progress[end][2] == 10
                    @test issorted(p[1] for p in progress)
                    @test issorted(p[2] for p in progress)
                end
            end
        end
    end

    save_env = get(ENV, "JULIA_SSL_NO_VERIFY_HOSTS", nothing)
    delete!(ENV, "JULIA_SSL_NO_VERIFY_HOSTS")

    @testset "bad TLS" begin
        urls = [
            "https://wrong.host.badssl.com"
            "https://untrusted-root.badssl.com"
        ]
        @testset "bad TLS is rejected" for url in urls
            resp = request(url, throw=false)
            @test resp isa RequestError
            @test resp.code == Curl.CURLE_PEER_FAILED_VERIFICATION
        end
        @testset "easy hook work-around" begin
            local url
            easy_hook = (easy, info) -> begin
                Curl.set_ssl_verify(easy, false)
                @test info.url == url
            end
            # downloader-specific easy hook
            downloader = Downloader()
            downloader.easy_hook = easy_hook
            for outer url in urls
                resp = request(url, throw=false, downloader=downloader)
                @test resp isa Response
                @test resp.status == 200
            end
            # default easy hook
            Downloads.EASY_HOOK[] = easy_hook
            Downloads.DOWNLOADER[] = nothing
            for outer url in urls
                resp = request(url, throw=false)
                @test resp isa Response
                @test resp.status == 200
            end
            Downloads.DOWNLOADER[] = nothing
            Downloads.EASY_HOOK[] = nothing
        end
        ENV["JULIA_SSL_NO_VERIFY_HOSTS"] = "**.badssl.com"
        @testset "SSL no verify override" for url in urls
            resp = request(url, throw=false)
            @test resp isa Response
            @test resp.status == 200
        end
        delete!(ENV, "JULIA_SSL_NO_VERIFY_HOSTS")
    end

    @testset "SNI required" begin
        url = "https://juliahub.com" # anything served by CloudFront
        # secure verified host request
        resp = request(url, throw=false, downloader=Downloader())
        @test resp isa Response
        @test resp.status == 200
        # insecure unverified host request
        ENV["JULIA_SSL_NO_VERIFY_HOSTS"] = "**"
        resp = request(url, throw=false, downloader=Downloader())
        @test resp isa Response
        @test resp.status == 200
    end

    if save_env !== nothing
        ENV["JULIA_SSL_NO_VERIFY_HOSTS"] = save_env
    else
        delete!(ENV, "JULIA_SSL_NO_VERIFY_HOSTS")
    end

    @__MODULE__() == Main && @testset "ftp download" begin
        url = "ftp://xmlsoft.org/libxslt/libxslt-1.1.33.tar.gz"
        file = Downloads.download(url)
        @test isfile(file)
        @test filesize(file) == 3444093
        head = String(read!(open(file), Vector{UInt8}(undef, 16)))
        @test head == "\x1f\x8b\b\0\xa5T.\\\x02\x03\xec]{s۶"
    end

    @testset "grace cleanup" begin
        dl = Downloader(grace=1)
        Downloads.download("https://httpbingo.org/drip"; downloader=dl)
        Downloads.download("https://httpbingo.org/drip"; downloader=dl)
    end
end

Downloads.DOWNLOADER[] = nothing
GC.gc(true)


text:
include("../common/PWGrid_v02.jl")
include("../common/ortho_gram_schmidt.jl")
include("../common/wrappers_fft.jl")

include("EnergiesT.jl")
include("PotentialsT.jl")
include("gen_dr.jl")
include("op_K.jl")
include("op_V_loc.jl")
include("op_H.jl")
include("calc_rho.jl")
include("calc_grad.jl")
include("calc_Energies.jl")
include("KS_solve_Emin_sd.jl")
include("KS_solve_Emin_cg.jl")
include("Poisson_solve.jl")
include("LDA_VWN.jl")
include("Kprec.jl")

include("../common/gen_dr_center.jl")
include("../common/calc_strfact_v2.jl")
include("../common/calc_ewald_v2.jl")

include("diag_lobpcg.jl")
include("KS_solve_SCF.jl")

include("andersonmix.jl")
include("KS_solve_SCF_andersonmix.jl")

include("get_ub_lb_lanczos.jl")
include("KS_solve_ChebySCF.jl")
include("chebyfilt.jl")
include("norm_matrix_induced.jl")

function test_main( Ns, xx; method="SCF" )

    const LatVecs = 16.0*diagm( ones(3) )

    pw = PWGrid( Ns, LatVecs )

    const Ω  = pw.Ω
    const r  = pw.r
    const G  = pw.gvec.G
    const G2 = pw.gvec.G2
    const Npoints = prod(Ns)
    const Ngwx = pw.gvecw.Ngwx

    Xpos = zeros( 3, 2 )
    Xpos[:,1] = [0.0, 0.0, 0.0]
    Xpos[:,2] = [xx, 0.0, 0.0]

    println(Xpos')

    Nspecies = 1
    atmsymb = ["H"] # unique list of atomic symbols
    atm2species = [1, 1]    # mapping from atom to species
    Zv = [1.0]    # only valence ?

    Sf = calc_strfact( Xpos, Nspecies, atm2species, pw.gvec.G )

    E_nn = calc_ewald( pw, Sf, Xpos, Nspecies, atm2species, Zv )
    #@printf("E_nn = %18.10f\n", E_nn)

    Vg = zeros(ComplexF64,Npoints)
    prefactor = -4*pi/Ω
    for ig=2:Npoints
        Vg[ig] = prefactor/G2[ig]
    end
    V_ionic = real( G_to_R(Ns, Vg .* Sf) ) * Npoints

    # Need to sum up over Nspecies for more than one species type
    # We simply need reshape because we only have one species type here.
    V_ionic = reshape( V_ionic, (Npoints) )

    #println("sum(Sf) = ", sum(Sf))
    #println("size(V_ionic) = ", size(V_ionic))

    const Nstates = 1
    Focc = [2.0]

    if method == "CG"
        psi, Energies, Potentials = KS_solve_Emin_cg(
            pw, V_ionic, Focc, Nstates, NiterMax=1000, E_NN=E_nn )
        #
        Y = ortho_gram_schmidt(psi)
        mu = Y' * op_H( pw, Potentials, Y )
        evals, evecs = eig(mu)
        psi = Y*evecs
    elseif method == "SCF"
        Energies, Potentials, psi, evals = KS_solve_SCF( pw, V_ionic, Focc, Nstates, E_NN=E_nn )
    elseif method == "SCF_andersonmix"
        Energies, Potentials, psi, evals = KS_solve_SCF_andersonmix( pw, V_ionic, Focc, Nstates, E_NN=E_nn )
    elseif method == "ChebySCF"
        Energies, Potentials, psi, evals = KS_solve_ChebySCF( pw, V_ionic, Focc, Nstates, β=0.8, E_NN=E_nn )
    end

    for st = 1:Nstates
        @printf("State # %d, Energy = %f\n", st, real(evals[st]))
    end

    @printf("E_nn = %18.10f\n", E_nn)
    @printf("E total = %18.10f\n", Energies.Total)
    print_Energies(Energies)

end

@time test_main( [64,64,64], 1.5, method="SCF" )
@time test_main( [64,64,64], 1.5, method="SCF_andersonmix" )

#@time test_main( [64,64,64], 1.5, method="ChebySCF" )
#@time test_main( [64,64,64], 1.5, method="CG" )

#for xx in [0.5, 1.0, 1.25, 1.50, 1.75, 2.0, 4.0, 6.0]
#    test_main( [64,64,64], xx )
#end


text:
"""
    SegmentedSum{V <: AbstractVector{<:Number}} <: AbstractAggregation

[`AbstractAggregation`](@ref) implementing segmented sum aggregation:

``
f(\\{x_1, \\ldots, x_k\\}) = \\sum_{i = 1}^{k} x_i
``

Stores a vector of parameters `ψ` that are filled into the resulting matrix in case an empty bag is encountered.

See also: [`AbstractAggregation`](@ref), [`AggregationStack`](@ref), [`sum_aggregation`](@ref),
    [`SegmentedMax`](@ref), [`SegmentedMean`](@ref), [`SegmentedPNorm`](@ref), [`SegmentedLSE`](@ref).
"""
struct SegmentedSum{V <: AbstractVector{<:Number}} <: AbstractAggregation
    ψ::V
end

Flux.@functor SegmentedSum

SegmentedSum(T::Type, d::Int) = SegmentedSum(zeros(T, d))
SegmentedSum(d::Int) = SegmentedSum(Float32, d)

Flux.@forward SegmentedSum.ψ Base.getindex, Base.length, Base.size, Base.firstindex, Base.lastindex,
        Base.first, Base.last, Base.iterate, Base.eltype

Base.vcat(as::SegmentedSum...) = reduce(vcat, as |> collect)
function Base.reduce(::typeof(vcat), as::Vector{<:SegmentedSum})
    SegmentedSum(reduce(vcat, [a.ψ for a in as]))
end

function (a::SegmentedSum)(x::Maybe{AbstractMatrix{T}}, bags::AbstractBags,
                              w::Optional{AbstractVecOrMat{T}}=nothing) where T
    _check_agg(a, x)
    segmented_sum_forw(x, a.ψ, bags, w)
end

segmented_sum_forw(::Missing, ψ::AbstractVector, bags::AbstractBags, w) = repeat(ψ, 1, length(bags))
function segmented_sum_forw(x::AbstractMatrix, ψ::AbstractVector, bags::AbstractBags, w::Optional{AbstractVecOrMat})
    t = promote_type(eltype(x), eltype(ψ))
    y = zeros(t, size(x, 1), length(bags))
    @inbounds for (bi, b) in enumerate(bags)
        if isempty(b)
            for i in eachindex(ψ)
                y[i, bi] = ψ[i]
            end
        else
            for j in b
                for i in 1:size(x, 1)
                    y[i, bi] += _weight(w, i, j, t) * x[i, j]
                end
            end
        end
    end
    y
end

function segmented_sum_back(Δ, y, x, ψ, bags, w)
    dx = similar(x)
    dψ = zero(ψ)
    dw = isnothing(w) ? ZeroTangent() : zero(w)
    @inbounds for (bi, b) in enumerate(bags)
        if isempty(b)
            for i in eachindex(ψ)
                dψ[i] += Δ[i, bi]
            end
        else
            for j in b
                for i in 1:size(x, 1)
                    dx[i, j] = _weight(w, i, j, eltype(x)) * Δ[i, bi]
                    ∇dw_segmented_sum!(dw, Δ, x, y, w, i, j, bi)
                end
            end
        end
    end
    dx, dψ, NoTangent(), dw
end

function segmented_sum_back(Δ, y, x::Missing, ψ, bags, w::Nothing)
    dψ = zero(ψ)
    @inbounds for (bi, b) in enumerate(bags)
        for i in eachindex(ψ)
            dψ[i] += Δ[i, bi]
        end
    end
    ZeroTangent(), dψ, NoTangent(), ZeroTangent()
end

∇dw_segmented_sum!(dw::ZeroTangent, Δ, x, y, w::Nothing, i, j, bi) = nothing
function ∇dw_segmented_sum!(dw::AbstractVector, Δ, x, y, w::AbstractVector, i, j, bi)
    dw[j] += Δ[i, bi] * (x[i, j])
end
function ∇dw_segmented_sum!(dw::AbstractMatrix, Δ, x, y, w::AbstractMatrix, i, j, bi)
    dw[i, j] += Δ[i, bi] * (x[i, j])
end

function ChainRulesCore.rrule(::typeof(segmented_sum_forw), args...)
    y = segmented_sum_forw(args...)
    grad = Δ -> (NoTangent(), segmented_sum_back(Δ, y, args...)...)
    y, grad
end


text:
abstract type NPlayerNavigationCost <: PlayerCost end

"Returns an the input cost object, e.g. ::QuadCost"
function inputcost end
"Returns an iterable of `::SoftConstr` for the inputs."
function inputconstr end
"Returns  a `::QuadCost`."
function statecost end
"Returns an interable of `::SoftConstr` for the states."
function stateconstr end
"The `::ProximityCost` for a player with any other player."
function proximitycost end
"The terminal `::GoalCost` to encourage reaching some goal"
function goalcost end

"--------------------------------- Implementation ---------------------------------"

function quadraticize!(qcache::QuadCache, pc::NPlayerNavigationCost,
                       g::GeneralGame, x::SVector, u::SVector, t::AbstractFloat)
    nx = n_states(g)
    nu = n_controls(g)

    xi = xindex(g)[player_id(pc)]
    ui = uindex(g)[player_id(pc)]

    # reset the cache
    reset!(qcache)

    # quadratic input cost
    quad!(qcache.R, qcache.r, inputcost(pc), u[ui], ui)
    # input constraints
    for constrᵢ in inputconstr(pc)
        quad!(qcache.R, qcache.r, constrᵢ, u, ui)
    end

    # quadratic state cost
    quad!(qcache.Q, qcache.l, statecost(pc), x[xi], xi)
    # state constraints
    for constrᵢ in stateconstr(pc)
        quad!(qcache.Q, qcache.l, constrᵢ, x, xi)
    end

    # pairwise proximity cost
    xyi_ego = xyindex(g)[player_id(pc)]
    for (j, xyi_other) in enumerate(xyindex(g))
        j != player_id(pc) || continue
        quad!(qcache.Q, qcache.l, proximitycost(pc), x, xyi_ego, xyi_other, j)
    end

    # the goal cost
    quad!(qcache.Q, qcache.l, goalcost(pc), x[xi], xi, t)

    return QuadraticPlayerCost(SVector(qcache.l), SMatrix(qcache.Q),
                               SVector(qcache.r), SMatrix(qcache.R))
end

function (pc::NPlayerNavigationCost)(g::AbstractGame, x::SVector, u::SVector, t::Float64)
    # extract the states and inputs for this player
    xi = xindex(g)[player_id(pc)]
    ui = uindex(g)[player_id(pc)]
    xᵢ = x[xi]
    uᵢ = u[ui]
    # setup the cost: each player wan't to:
    cost = 0.0

    # control cost: only cares about own control
    cost += inputcost(pc)(uᵢ)
    for constrᵢ in inputconstr(pc)
        cost += constrᵢ(u, ui)
    end

    # running cost for states (e.g. large steering)
    cost += statecost(pc)(xᵢ)
    for constrᵢ in stateconstr(pc)
        cost += constrᵢ(x, xi)
    end

    # proximity constraint
    # ego positions
    xys = xyindex(g)
    xp_ego = x[xys[player_id(pc)]]
    for (j, xy_other) in enumerate(xys)
        j != player_id(pc) || continue
        # other positions
        xp_other = x[xys[j]]
        cost += proximitycost(pc)(xp_ego, xp_other, j)
    end

    # goal state cost cost:
    cost += goalcost(pc)(xᵢ, t)

    return cost
end

function generate_nplayer_navigation_game(DynType::Type, CostModelType::Type,
                                          T_horizon::Float64, ΔT::Float64,
                                          goals::Vararg{SVector}; kwargs...)
    # the number of players
    np = length(goals)
    # the time at which the goal cost is activated
    t_final = T_horizon - 1.5*ΔT
    h = Int(T_horizon/ΔT)
    # setting up the dynamics
    dyn = ProductSystem(Tuple(DynType{ΔT}() for i in 1:np))

    xids = xindex(dyn)
    uids = uindex(dyn)
    costs = SVector{np}([CostModelType(i, goals[i], t_final, np; kwargs...) for i in 1:np])

    return GeneralGame(h, uids, dyn, costs)
end


text:
#!/usr/bin/env julia

using Winston

target = "getMarketsInfo"
(timing_data, labels) = readcsv("data/timing-" * target * ".csv", header=true)

num_markets = timing_data[:,1];
time_elapsed = timing_data[:,2];

P = plot(num_markets, time_elapsed, "b-");
xlabel(labels[1]);
ylabel(labels[2]);
title(target)

savefig(P, "data/timing-" * target * ".png")


text:
for N in [Float64, Float32, Rational{Int}]
    # BallInf approximation of a 3D unit ball in the 1-norm centered at [1,2,0]
    b = Ball1(N[1, 2, 0], N(1))
    bi = ballinf_approximation(b)
    biexp = BallInf(N[1, 2, 0], N(1))
    @test bi.center ≈ biexp.center
    @test bi.radius ≈ biexp.radius

    # BallInf approximation of a 2D polygon whose vertices are
    # (0,1), (1,1), (2,-1), (-1,0)
    p = HPolygon{N}()
    addconstraint!(p, LinearConstraint(N[0, 1], N(1)))
    addconstraint!(p, LinearConstraint(N[-1, 1], N(1)))
    addconstraint!(p, LinearConstraint(N[-1, -3], N(1)))
    addconstraint!(p, LinearConstraint(N[2, 1], N(3)))
    bi = ballinf_approximation(p)
    biexp = BallInf(N[0.5, 0], N(1.5))
    @test bi.center ≈ biexp.center
    @test bi.radius ≈ biexp.radius

    # empty set
    E = EmptySet{N}(2)
    @test ballinf_approximation(E) == E

    if N == Float64
        # robustness (see issue #2532): the set has two contradicting constraints
        # => the set is empty, but requires to use high precision
        # x >= 1.0000000000000002
        # x <= 1
        s1 = HalfSpace([-1.0], -1.0000000000000002)
        s2 = HalfSpace([1.0], 1.0)
        P = s1 ∩ s2
        @test ballinf_approximation(P) == BallInf([1.0], 0.0)

        s1big = HalfSpace([-big(1.0)], -big(1.0000000000000002))
        s2big = HalfSpace([big(1.0)], big(1.0))
        Pbig = HPolytope([s1big, s2big])
        @test ballinf_approximation(Pbig) == EmptySet{BigFloat}(1)
    end
end


text:
# This code generates the 2l+l dimensional real irreducible representations of O(3).
#
# For a function defined on the sphere (f: S² → R) and represented in the basis of the
# real spherical harmonics, these matrix representations rotate/reflect the function
# about an axis in 3D space.
#
# The routines in this file are based upon:
#
#   Rotation Matrices for Real Spherical Harmonics, Direct Determination by Recursion
#   Joseph Ivanic and Klaus Ruedenberg
#   J. Phys. Chem., Vol. 100, No. 15, 1996
#
# The article contains numerous errors, and the "Additions and Corrections"
# — J. Phys. Chem. A, Vol. 102, No. 45, 1998 — also contains at least one error,
# for the matrices generated using Table 2 do not constitute an O(3) group
# homomorphism. A few empirical tests reveal that the term "sqrt(1 - δ(m,-1))" should
# probably be changed to "sqrt(1 + δ(m,-1))" in Table 2 for "V_mm'" and "m < 0".
#
# This file contains a reference implementation of the (corrected) methods in the
# article. The code in this file is not intended for actual use, however, as an
# optimized implementation exists in irreps_O3.jl.

function R(R₀, l::Int, m::Int, n::Int)
    #println("R: l = $l, m = $m, n = $n")
    if abs(m) > l || abs(n) > l
        0
    elseif l == 1
        R₀(m, n)
    elseif l > 1
        u(l, m, n) * U(R₀, l, m, n) +
        v(l, m, n) * V(R₀, l, m, n) +
        w(l, m, n) * W(R₀, l, m, n)
    end
end

δ(i::Int, j::Int) = i == j ? 1 : 0

function u(l::Int, m::Int, n::Int)
    if abs(n) < l
        sqrt(((l + m) * (l - m)) / ((l + n) * (l - n)))
    elseif abs(n) == l
        sqrt(((l + m) * (l - m)) / ((2l) * (2l - 1)))
    end
end

function v(l::Int, m::Int, n::Int)
    if abs(n) < l
        (1/2) *
        sqrt(((1 + δ(m, 0)) * (l + abs(m) - 1) * (l + abs(m))) / ((l + n) * (l - n))) *
        (1 - 2δ(m, 0))
    elseif abs(n) == l
        (1/2) *
        sqrt(((1 + δ(m, 0)) * (l + abs(m) - 1) * (l + abs(m))) / ((2l) * (2l - 1))) *
        (1 - 2δ(m, 0))
    end
end

function w(l::Int, m::Int, n::Int)
    if abs(n) < l
        (-1/2) *
        sqrt(((l - abs(m) - 1) * (l - abs(m))) / ((l + n) * (l - n))) *
        (1 - δ(m, 0))
    elseif abs(n) == l
        (-1/2) *
        sqrt(((l - abs(m) - 1) * (l - abs(m))) / ((2l) * (2l - 1))) *
        (1 - δ(m, 0))
    end
end

function P(R₀, l::Int, i::Int, m::Int, n::Int)
    #println("P: l = $l, i = $i, m = $m, n = $n")
    if abs(n) < l
        R₀(i, 0) * R(R₀, l - 1, m, n)
    elseif n == l
        R₀(i, 1) * R(R₀, l - 1, m,  l - 1) - R₀(i, -1) * R(R₀, l - 1, m, -l + 1)
    elseif n == -l
        R₀(i, 1) * R(R₀, l - 1, m, -l + 1) + R₀(i, -1) * R(R₀, l - 1, m,  l - 1)
    end
end

function U(R₀, l::Int, m::Int, n::Int)
    #println("U: l = $l, m = $m, n = $n")
    if m == 0
        P(R₀, l, 0, 0, n)
    elseif m > 0
        P(R₀, l, 0, m, n)
    elseif m < 0
        P(R₀, l, 0, m, n)
    end
end

function V(R₀, l::Int, m::Int, n::Int)
    #println("V: l = $l, m = $m, n = $n")
    if m == 0
        P(R₀, l, 1, 1, n) + P(R₀, l, -1, -1, n)
    elseif m > 0
        P(R₀, l, 1, m - 1, n) * sqrt(1 + δ(m, 1)) - P(R₀, l, -1, -m + 1, n) * (1 - δ(m, 1))
    elseif m < 0
        P(R₀, l, 1, m + 1, n) * (1 - δ(m, -1)) + P(R₀, l, -1, -m - 1, n) * sqrt(1 + δ(m, -1))
    end
end

function W(R₀, l::Int, m::Int, n::Int)
    #println("W: l = $l, m = $m, n = $n")
    if m == 0
        0
    elseif m > 0
        P(R₀, l, 1, m + 1, n) + P(R₀, l, -1, -m - 1, n)
    elseif m < 0
        P(R₀, l, 1, m - 1, n) - P(R₀, l, -1, -m + 1, n)
    end
end


text:
include("precedence.jl")
include("extensional.jl")


text:
using StableDQMC
using Test, Random, LinearAlgebra

@testset "StableDQMC.jl" begin

    @testset "Opt-in JacobiSVD/GenericSVD" begin
        @test_throws UndefVarError gesvj
        @test_throws UndefVarError gesvj!
        @test_throws UndefVarError genericsvd
        @test_throws UndefVarError genericsvd!
        StableDQMC.addJacobiSVD()
        StableDQMC.addGenericSVD()
    end

    @testset "UDT Type" begin
        Random.seed!(1234)
        x = rand(10,10)

        @test typeof(udt(x)) == UDT{Float64,Float64,Array{Float64,2}}
        @test typeof(udt!(copy(x))) == UDT{Float64,Float64,Array{Float64,2}}

        F = udt(x)
        u,d,t = F

        @test u * u' ≈ I # unitarity of u
        @test isreal(d)
        @static if VERSION < v"1.1.0-DEV.792"
            eachcol(A::AbstractMatrix) = (view(A, :, i) for i in axes(A, 2))
        end
        @test sort!(count.(iszero,eachcol(t))) == 0:9 # triangularity of t (istriu(t) == false due to pivoting)
        @test Matrix(F) ≈ x
        @test u * Diagonal(d) * t ≈ x

        @test size(F) == size(x)
        @test typeof(similar(F)) == UDT{Float64,Float64,Array{Float64,2}}

        N = udt(3)
        @test N.U == reshape([1.0], 1,1)
        @test N.D == [3.0]
        @test N.T == reshape([1.0], 1,1)

        # operations
        @test inv(F) ≈ inv(x)
        G = fact_mult(F, F)
        @test Matrix(G) ≈ x*x
    end

    @testset "SVD decomposition variants" begin
        for dtype in (Float64, ComplexF64)
            for decomp in (gesdd, gesvd, gesvj, genericsvd)
                x = rand(dtype, 5,5)

                F = decomp(x)
                @test typeof(F) <: SVD
                @test Matrix(F) ≈ x
                @test F.U * F.U' ≈ I
                @test isreal(F.S)
                @test F.Vt * F.Vt' ≈ I
            end
        end

        a, b = rand(5,5), rand(5,5)
        A, B = svd(a), svd(b)
        @test Matrix(fact_mult(A,B)) ≈ a*b
    end

    @testset "QR / UDT operations" begin
        for dtype in (Float64, ComplexF64)
            x = rand(dtype, 5,5)
            X = udt(x)
            y = rand(dtype, 5,5)
            Y = udt(y)
            res = similar(x)

            inv_I_plus_x = inv(I + x)
            inv_I_plus_xydagger = inv(I + x * y')
            inv_sum_xy = inv(x + y)

            @test inv_one_plus(X) ≈ inv_I_plus_x
            @test inv_one_plus!(res, X) ≈ inv_I_plus_x
            @test res ≈ inv_I_plus_x

            @test typeof(udt_inv_one_plus(X, Y)) <: UDT
            @test Matrix(udt_inv_one_plus(X, Y)) ≈ inv_I_plus_xydagger
            @test inv_one_plus(X,Y) ≈ inv_I_plus_xydagger
            @test inv_one_plus!(res, X,Y) ≈ inv_I_plus_xydagger
            @test res ≈ inv(I + x * y')

            @test typeof(udt_inv_sum(X,Y)) <: UDT
            @test Matrix(udt_inv_sum(X,Y)) ≈ inv_sum_xy
            @test inv_sum(X,Y) ≈ inv_sum_xy
            @test inv_sum!(res, X,Y) ≈ inv_sum_xy
            @test res ≈ inv_sum_xy

            # Loh et al variants
            @test typeof(udt_inv_one_plus_loh(X)) <: UDT
            @test Matrix(udt_inv_one_plus_loh(X)) ≈ inv_I_plus_x
            @test inv_one_plus_loh(X) ≈ inv_I_plus_x
            @test inv_one_plus_loh!(res,X) ≈ inv_I_plus_x
            @test res ≈ inv_I_plus_x

            @test typeof(udt_inv_sum_loh(X,Y)) <: UDT
            @test Matrix(udt_inv_sum_loh(X,Y)) ≈ inv_sum_xy
            @test inv_sum_loh(X,Y) ≈ inv_sum_xy
            @test inv_sum_loh!(res,X,Y) ≈ inv_sum_xy
            @test res ≈ inv_sum_xy
        end
    end

    @testset "SVD / USVt operations" begin
        for dtype in (Float64, ComplexF64)
            x = rand(dtype, 5,5)
            X = svd(x)
            y = rand(dtype, 5,5)
            Y = svd(y)
            res = similar(x)

            inv_I_plus_x = inv(I + x)
            inv_sum_xy = inv(x + y)

            @test inv_one_plus(X) ≈ inv_I_plus_x
            @test Matrix(svd_inv_one_plus_loh(X)) ≈ inv_I_plus_x
            @test inv_one_plus_loh(X) ≈ inv_I_plus_x
            @test inv_one_plus_loh!(res, X) ≈ inv_I_plus_x
            @test res ≈ inv_I_plus_x

            @test Matrix(svd_inv_sum(X, Y)) ≈ inv_sum_xy
            @test inv_sum(X, Y) ≈ inv_sum_xy
            @test inv_sum!(res, X, Y) ≈ inv_sum_xy
            @test res ≈ inv_sum_xy

            @test Matrix(svd_inv_sum_loh(X, Y)) ≈ inv_sum_xy
            @test inv_sum_loh(X, Y) ≈ inv_sum_xy
            @test inv_sum_loh!(res, X, Y) ≈ inv_sum_xy
            @test res ≈ inv_sum_xy
        end
    end

    @testset "Helpers" begin
        x = [1 2; 3 4]
        @test cond(svdvals(x)) ≈ cond(x)

        R, svs = calc_Bchain(x, 5)
        @test R ≈ x^5
        for i in 1:5
            @test svs[i,:] ≈ log.(svdvals(x^i))
        end

        # SVD / UDV
        R, svs = calc_Bchain_svd(x, 5)
        xsvd = svd(x^5)
        @test R.U ≈ xsvd.U
        @test R.S ≈ xsvd.S
        @test R.Vt ≈ xsvd.Vt
        for i in 1:5
            @test svs[i,:] ≈ log.(svdvals(x^i))
        end

        # QR / UDT
        R, svs = calc_Bchain_qr(x, 5)
        xudt = udt(x^5)
        @test R.U ≈ xudt.U
        @test R.D ≈ xudt.D
        @test R.T ≈ xudt.T
        for i in 1:5
            @test svs[i,:] ≈ log.(udt(x^i).D)
        end

        g = calc_tdgf_qr(x,2,3)
        @test g ≈ inv(x^-2 + x^3) # [B^-N1 + B^N2]^-1
        g = calc_tdgf_svd(x,2,3)
        @test g ≈ inv(x^-2 + x^3) # [B^-N1 + B^N2]^-1
    end

end


text:
__precompile__()
module CloudWatchLogs

using AWSCore: AWSConfig, AWSException
using AWSCore.Services: logs
using Dates
using MbedTLS: MbedException
using Memento
using Mocking
using TimeZones
using UUIDs

export CloudWatchLogStream, LogEvent, submit_log, submit_logs
export create_group, delete_group, create_stream, delete_stream
export CloudWatchLogHandler
export StreamNotFoundException, LogSubmissionException

const LOGGER = getlogger(@__MODULE__)

# Info on limits:
# https://docs.aws.amazon.com/AmazonCloudWatch/latest/logs/cloudwatch_limits_cwl.html

# 1 MB (maximum). This limit cannot be changed.
const MAX_BATCH_SIZE = 1048576

# 256 KB (maximum). This limit cannot be changed.
const MAX_EVENT_SIZE = 262144

# https://docs.aws.amazon.com/AmazonCloudWatchLogs/latest/APIReference/API_PutLogEvents.html#CWL-PutLogEvents-request-logEvents
const MAX_BATCH_LENGTH = 10000

# 5 requests per second per log stream. This limit cannot be changed.
const PUTLOGEVENTS_RATE_LIMIT = 0.2
const PUTLOGEVENTS_DELAYS =
    ExponentialBackOff(n=10, first_delay=PUTLOGEVENTS_RATE_LIMIT, factor=1.1)

const GENERIC_AWS_DELAYS = ExponentialBackOff(n=10, first_delay=0.2, factor=2, jitter=0.2)

__init__() = Memento.register(LOGGER)

include("exceptions.jl")
include("event.jl")
include("stream.jl")
include("handler.jl")

end


text:
using BinaryBuilder

name = "DecFP"
version = v"2.0.3" # --> This is a lie to build for experimental platforms, upstream version is 2.0 Update 2

# Collection of sources required to build DecFP
sources = [
    ArchiveSource("https://www.netlib.org/misc/intel/IntelRDFPMathLib20U2.tar.gz",
                  "93c0c78e0989df88f8540bf38d6743734804cef1e40706fd8fe5c6a03f79e173"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/IntelRDFPMathLib20U2
patch --binary -p1 < $WORKSPACE/srcdir/patches/decfp.patch
cd LIBRARY
if [[ $nbits == 64 ]]; then
     _HOST_ARCH=x86_64
else
     _HOST_ARCH=x86
fi
if [[ $target == *-w64-* ]]; then
    _HOST_OS="Windows_NT"
    objext="obj"
elif [[ $target == *-freebsd* ]]; then
    _HOST_OS="FreeBSD"
    objext="o"
else
    # Linux and macOS/Darwin
    _HOST_OS="Linux"
    objext="o"
fi
CFLAGS_OPT="-fPIC -fsigned-char -D__ENABLE_BINARY80__=0"
if [[ $target == *"-musl"* ]]; then
    CFLAGS_OPT+=" -D__QNX__"
elif [[ $target == *"freebsd"* ]]; then
    CFLAGS_OPT+=" -D__linux"
fi
make CC_NAME=cc CFLAGS_OPT="$CFLAGS_OPT" CFLAGS="$CFLAGS_OPT" _HOST_OS="$_HOST_OS" AR_CMD="ar rv" _HOST_ARCH=$_HOST_ARCH CALL_BY_REF=0 GLOBAL_RND=0 GLOBAL_FLAGS=0 UNCHANGED_BINARY_FLAGS=0
mkdir -p "${libdir}"
$CC $LDFLAGS -shared -o "${libdir}/libbid.${dlext}" *.${objext}
install_license ../eula.txt
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libbid", :libbid)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
kernelM(l::Int,x, y) = x * diffsbesselj(l,x) * sbesselj(l,y) - y * sbesselj(l,x) * diffsbesselj(l,y)

function scattering_field(wavemode::EffectiveRegularWaveMode{T,3,Acoustic{T,3},WithoutSymmetry{3}}) where T
    L = wavemode.basis_order
    L1 = wavemode.basis_field_order

    i2s = 1:size(wavemode.eigenvectors,2)
    i3s = 1:size(wavemode.eigenvectors,3)

    nn1_to_n = [
            basisorder_to_basislength(Acoustic{T,3},l-1) + m + l + 1
    for l = 0:L for m = -l:l for l1 = 0:L1 for m1 = -l1:l1];

    nn1_to_n1 = [
            basisorder_to_basislength(Acoustic{T,3},l1-1) + m1 + l1 + 1
    for l = 0:L for m = -l:l for l1 = 0:L1 for m1 = -l1:l1];

    indices_same_n = collect(groupby(nn1 -> nn1_to_n[nn1] , 1:length(nn1_to_n)))

    v = regular_basis_function(wavemode.wavenumber, wavemode.medium)

    function scat_field(x::AbstractVector{T})
        vs =  v(L1, x)

        vecs = [
            sum(wavemode.eigenvectors[ind,i2,i3] .* vs[nn1_to_n1[ind]]) # note here:  vs[nn1_to_n1[ind]]] == vs
        for ind in indices_same_n, i2 in i2s, i3 in i3s];

        # Sum over different eigvectors
        vecs = sum(vecs, dims = 3)[:,:];

        # for only 1 species the code below gives the same
        # nn1_to_n1 = [
        #         basisorder_to_basislength(Acoustic{T,3},l1-1) + m1 + l1 + 1
        # for l = 0:L for m = -l:l for l1 = 0:L1 for m1 = -l1:l1];
        #
        # vecs2 = [
        #     wavemode.eigenvectors[i[1],i[2],i[3]] * vs[nn1_to_n1[i[1]]]
        # for i in CartesianIndices(wavemode.eigenvectors)];
        #
        # vecs2 = sum(vecs2, dims = 3);
        # vecs2 = [sum(vecs2[inds]) for inds in indices_same_n]
        # vecs ≈ vecs2
        return vecs
    end

end

function scattering_field(wavemode::EffectiveRegularWaveMode{T,3,Acoustic{T,3},Sym}) where {T, Sym<:AbstractAzimuthalSymmetry{3}}
    L = wavemode.basis_order
    L1 = wavemode.basis_field_order

    i2s = 1:size(wavemode.eigenvectors,2)
    i3s = 1:size(wavemode.eigenvectors,3)

    nl1_to_n = [
            basisorder_to_basislength(Acoustic{T,3},l-1) + m + l + 1
    for l = 0:L for m = -l:l for l1 = abs(m):L1];

    nl1_to_n1 = [
            basisorder_to_basislength(Acoustic{T,3},l1-1) - m + l1 + 1
    for l = 0:L for m = -l:l for l1 = abs(m):L1];

    indices_same_n = collect(groupby(nl1 -> nl1_to_n[nl1] , 1:length(nl1_to_n)))

    v = regular_basis_function(wavemode.wavenumber, wavemode.medium)

    function scat_field(x::AbstractVector{T})
        vs =  v(L1, x)

        vecs = [
            sum(wavemode.eigenvectors[ind,i2,i3] .* vs[nl1_to_n1[ind]])
        for ind in indices_same_n, i2 in i2s, i3 in i3s];

        # Sum over different eigvectors
        vecs = sum(vecs, dims = 3)[:,:];

        return vecs
    end

end

"""
    scattering_field(wavemode::EffectiveRegularWaveMode{T,3,Acoustic{T,3},RadialSymmetry{3}})

In general ``F_{nn'} = \\delta_{m',-m} \\delta_{\\ell,\\ell'}(-1)^m F_{(\\ell,0),(\\ell,0)}``, so the scattering field is given by

``\\langle f_n \\rangle (x) = \\sum_{n'} F_{nn'} \\mathrm v_{n'}(k_* x) =  F_{(\\ell,0),(\\ell,0)} (-1)^m \\mathrm v_{(\\ell,-m)}(k_* x)``
"""
function scattering_field(wavemode::EffectiveRegularWaveMode{T,3,Acoustic{T,3},RadialSymmetry{3}}) where T
    L = wavemode.basis_order
    len = basisorder_to_basislength(Acoustic{T,3},L)

    ls, ms = spherical_harmonics_indices(L)
    inds = lm_to_spherical_harmonic_index.(ls,-ms)

    v = regular_basis_function(wavemode.wavenumber, wavemode.medium)

    function scat_field(x::AbstractVector{T})
        vs =  v(L, x)
        vs = vs[inds] .* (-one(T)) .^ ms

        vecs = [
            sum(wavemode.eigenvectors[ls[n]+1,i2,:]) * vs[n]
        for n = 1:len, i2 = 1:size(wavemode.eigenvectors,2)];

        return vecs
    end
end

function material_scattering_coefficients(wavemode::EffectiveRegularWaveMode{T,3,Acoustic{T,3},WithoutSymmetry{3}}) where T

    # Unpacking parameters
    k = wavemode.ω / wavemode.medium.c

    species = wavemode.material.species
    S = length(species)
    rs = outer_radius.(species)

    R = outer_radius(wavemode.material.shape)
    k_eff = wavemode.wavenumber

    L = wavemode.basis_order
    L1 = wavemode.basis_field_order
    dL = L + L1

    # Pre-calculating
    Ml1s = [
        (R - rs[j]) * kernelM(l1, k*(R - rs[j]), k_eff*(R - rs[j])) * number_density(species[j])
    for l1 = 0:L1, j in eachindex(species)] ./ (k_eff^T(2) - k^T(2))

    nn1_to_l1 = [l1 for l = 0:L for m = -l:l for l1 = 0:L1 for m1 = -l1:l1];

    # Sum over species and different eigvectors
    vecs = [
        wavemode.eigenvectors[i] * Ml1s[nn1_to_l1[i[1]]+1,i[2]]
    for i in CartesianIndices(wavemode.eigenvectors)]

    vec = reshape(sum(vecs, dims = (2,3)),:)

    # Sum over wave basis components
    gaunt = [
            gaunt_coefficient(dl,dm,l,m,l1,m1)
        for l = 0:L for m = -l:l
        for l1 = 0:L1 for m1 = -l1:l1
    for dl = 0:dL for dm = -dl:dl]

    # gaunt[(dn),(n,l1)] = gaunt_coefficient(dl,dm,l,m,l1,m1)
    gaunt = reshape(gaunt, (:,length(vec)))

    return gaunt * vec
end

function material_scattering_coefficients(wavemode::EffectiveRegularWaveMode{T,3,Acoustic{T,3},Sym}) where {T,Sym<:AbstractAzimuthalSymmetry{3}}

    # Unpacking parameters
    k = wavemode.ω / wavemode.medium.c

    species = wavemode.material.species
    S = length(species)
    rs = outer_radius.(species)

    R = outer_radius(wavemode.material.shape)
    k_eff = wavemode.wavenumber

    L = wavemode.basis_order
    L1 = wavemode.basis_field_order
    dL = L + L1

    # Pre-calculating
    Ml1s = [
        (R - rs[j]) * kernelM(l1, k*(R - rs[j]), k_eff*(R - rs[j])) * number_density(species[j])
    for l1 = 0:L1, j in eachindex(species)] ./ (k_eff^T(2) - k^T(2))

    nl1_to_l1 = [l1 for l = 0:L for m = -l:l for l1 = abs(m):L1];

    # Sum over species and different eigvectors
    vecs = [
        wavemode.eigenvectors[i] * Ml1s[nl1_to_l1[i[1]]+1,i[2]]
    for i in CartesianIndices(wavemode.eigenvectors)]

    vec = reshape(sum(vecs, dims = (2,3)),:)

    # Sum over wave basis components
    gaunt = [
            gaunt_coefficient(dl,dm,l,m,l1,-m)
        for l = 0:L for m = -l:l
        for l1 = abs(m):L1
    for dl = 0:dL for dm = -dl:dl]

    gaunt = reshape(gaunt, (:,length(vec)))

    return gaunt * vec
end

function material_scattering_coefficients(wavemode::EffectiveRegularWaveMode{T,3,Acoustic{T,3},RadialSymmetry{3}}) where T

    # Unpacking parameters
    k = wavemode.ω / wavemode.medium.c
    k_eff = wavemode.wavenumber
    R = outer_radius(wavemode.material.shape)

    species = wavemode.material.species
    S = length(species)
    rs = outer_radius.(species)

    L = wavemode.basis_order

    Fscat0 = sum(
        sqrt(T(4pi)) * wavemode.eigenvectors[i] * (2*(i[1]-1) + 1) * T(-1)^(i[1]-1) *
        (R - rs[i[2]]) * kernelM(i[1]-1, k*(R - rs[i[2]]), k_eff*(R - rs[i[2]])) * number_density(species[i[2]])
    for i in CartesianIndices(wavemode.eigenvectors)) / (k_eff^T(2) - k^T(2))

    return [Fscat0]
end

max_stars_repo_name: arturgower/EffectiveWaves


text:
module TestHamerly

using ParallelKMeans
using ParallelKMeans: chunk_initialize, double_argmax
using Test
using Random

@testset "initialize" begin
    X = permutedims([1.0 2; 2 1; 4 5; 6 6])
    centroids = permutedims([1.0 2; 4 5; 6 6])
    nrow, ncol = size(X)
    containers = ParallelKMeans.create_containers(Hamerly(), 3, nrow, ncol, 1)

    ParallelKMeans.chunk_initialize(Hamerly(), containers, centroids, X, 1:ncol, 1)
    @test containers.lb == [18.0, 20.0, 5.0, 5.0]
    @test containers.ub == [0.0, 2.0, 0.0, 0.0]
end

@testset "double argmax" begin
    @test double_argmax([0.5, 0, 0]) == (1, 2, 0.5, 0.0)
end

@testset "singlethread linear separation" begin
    # with the same amount of iterations answer should be the same as in Lloyd case
    Random.seed!(2020)

    X = rand(3, 100)
    res = kmeans(Hamerly(), X, 3; n_threads = 1, tol = 1e-10, max_iters = 10, verbose = false)

    @test res.totalcost ≈ 14.16198704459199
    @test !res.converged
    @test res.iterations == 10

    Random.seed!(2020)
    X = rand(3, 100)
    res = kmeans(Hamerly(), X, 3; n_threads = 1, tol = 1e-10, max_iters = 1000, verbose = false)

    @test res.totalcost ≈ 14.161987044591992
    @test res.converged
    @test res.iterations == 11
end

@testset "multithread linear separation quasi two threads" begin
    Random.seed!(2020)

    X = rand(3, 100)
    res = kmeans(Hamerly(), X, 3; n_threads = 2, tol = 1e-10, max_iters = 10, verbose = false)

    @test res.totalcost ≈ 14.16198704459199
    @test !res.converged
    @test res.iterations == 10

    Random.seed!(2020)
    X = rand(3, 100)
    res = kmeans(Hamerly(), X, 3; n_threads = 2, tol = 1e-10, max_iters = 1000, verbose = false)

    @test res.totalcost ≈ 14.161987044591992
    @test res.converged
    @test res.iterations == 11
end

end # module


text:
input = read("2020/data/day_8.txt", String)

test_input = """
nop +0
acc +1
jmp +4
acc +3
jmp -3
acc -99
acc +1
jmp -4
acc +6
"""

mutable struct State
    total::Int
    index::Int
    used_codes::Array{Int,1}
end

function apply_rule!(state::State, rule::AbstractString)
    push!(state.used_codes, state.index)

    key, str_val = split(rule)
    val = parse(Int, str_val)

    if key == "nop"
        state.index += 1
    elseif key == "jmp"
        state.index += val
    elseif key == "acc"
        state.total += val
        state.index += 1
    end
    # @info state
end

function run_program(rules)
    state = State(0, 1, [])
    for unused in 1:length(rules)  # theoretical max iterations
        if state.index ∈ state.used_codes
            return (state, "ERROR: Infinite Loop")
        elseif state.index == length(rules) + 1
            return state, "SUCCESS"
        else
            rule = rules[state.index]
            apply_rule!(state, rule)
        end
    end
end

function part_1(input)
    rules = split(strip(input), "\n")

    state, exit_code = run_program(rules)
    @assert exit_code == "ERROR: Infinite Loop"
    state
end
@assert part_1(test_input).total == 5
@info part_1(input)

function part_2(input)
    state = State(0, 1, [])
    rules = split(strip(input), "\n")

    for i in 1:length(rules)
        if occursin("nop", rules[i]) | occursin("jmp", rules[i])
            new_rules = copy(rules)
            if occursin("nop", rules[i])
                new_rules[i] = replace(rules[i], "nop" => "jmp")
            else
                new_rules[i] = replace(rules[i], "jmp" => "nop")
            end
            state, exit_code = run_program(new_rules)
            if exit_code == "SUCCESS"
                return state, i
            end
        end
    end
end
@info part_2(input)

max_stars_repo_name: RensDimmendaal/advent-of-code-julia

max_issues_repo_name: RensDimmendaal/advent-of-code-julia

max_forks_repo_name: RensDimmendaal/advent-of-code-julia


text:
# make this keyworded?
init(::NEqProblem, ::LineSearch, x) = (z=copy(x), d=copy(x), Fx=copy(x), Jx=x*x')
# the bang is just potentially inplace x and state. nonbang copies these
function solve(problem::NEqProblem, x, method::LineSearch=LineSearch(Newton(), Static(1)), options=NEqOptions(), state=init(problem, method, x))
    t0 = time()

    # Unpack
    scheme, linesearch = modelscheme(method), algorithm(method)
    # Unpack important objectives
    F = problem.R.F
    FJ = problem.R.FJ
    # Unpack state
    z, d, Fx, Jx = state
    T = eltype(Fx)

    # Set up MeritObjective. This defines the least squares
    # objective for the line search.
    merit = MeritObjective(problem, F, FJ, Fx, Jx, d)
    meritproblem = OptimizationProblem(merit, nothing, Euclidean(0), nothing, mstyle(problem), nothing)
    # Evaluate the residual and Jacobian
    Fx, Jx = FJ(Fx, Jx, x)
    ρF0, ρ2F0 = norm(Fx, Inf),  norm(Fx, 2)

    stoptol = T(options.f_reltol)*ρF0 + T(options.f_abstol)
    if ρF0 < stoptol
        return ConvergenceInfo(method, (solution=x, best_residual=Fx, ρF0=ρF0, ρ2F0=ρ2F0, ρs=T(NaN), iter=0, time=time()-t0), options)
    end

    # Create variable for norms but keep the first ones for printing purposes.
    ρs, ρ2F = ρF0, ρ2F0

    iter = 1
    while iter ≤ options.maxiter
        # Shift z into x
        if mstyle isa InPlace
            x .= z
        else
            x = copy(z)
        end
        # Update the search direction
        if mstyle isa InPlace
            d = scheme.linsolve(d, Jx, -Fx)
        else
            d = scheme.linsolve(Jx, -Fx)
        end
        # Need to restrict to static and backtracking here because we don't allow
        # for methods that calculate the gradient of the line objective.
        #
        # For non-linear systems of equations we choose the sum-of-
        # squares merit function. Some useful things to remember is:
        #
        # f(y) = 1/2*|| F(y) ||^2 =>
        # ∇_df = -d'*J(x)'*F(x)
        #
        # where we remember the notation x means the current iterate and y is any
        # proposal. This means that if we step in the Newton direction such that d
        # is defined by
        #
        # J(x)*d = -F(x) => -d'*J(x)' = F(x)' =>
        # ∇_df = -F(x)'*F(x) = -f(x)*2
        #
        # φ = LineObjective!(F, ∇fz, z, x, d, fx, dot(∇fx, d))
        φ = LineObjective(meritproblem, z, z, x, d, (ρ2F^2)/2, -ρ2F^2)

        # Perform line search along d
        α, f_α, ls_success = find_steplength(mstyle, linesearch, φ, T(1))

        # # Calculate final step vector and update the state

        # Step in the direction α*d
        z = retract(problem, z, x, d, α)

        # Update residual and jacobian
        Fx, Jx = FJ(Fx, Jx, z)

        # Update 2-norm for line search conditions: ϕ(0) and ϕ'(0)
        ρ2F = norm(Fx, 2)

        # Update the largest successive change in the iterate
        ρs = mapreduce(x->abs(x[1]-x[2]), max, zip(x,z)) # norm(x.-z, Inf)

        if ρ2F < stoptol || ρs <= 1e-12
            break
        end
        iter += 1
    end
    return ConvergenceInfo(method, (solution=z, best_residual=Fx, ρF0=ρF0, ρ2F0=ρ2F0, ρs=ρs, iter=iter, time=time()-t0), options)
end


text:
"""
    store perf info of each variable
"""
struct SplitInfoGPU{V,M,B}
    gains::V
    gainsL::V
    gainsR::V
    ∑Ls::M
    ∑Rs::M
    bins::B
end

"""
    Carries training information for a given tree node
"""
struct TrainNodeGPU{T<:AbstractFloat, S, I<:AbstractVector, V<:AbstractVector}
    parent::S
    depth::S
    ∑::V
    gain::T
    𝑖::I
    𝑗::I
end

struct TreeNodeGPU{T<:AbstractFloat, S, B<:Bool}
    left::S
    right::S
    feat::S
    cond::T
    gain::T
    pred::Vector{T}
    split::B
end

TreeNodeGPU(left::S, right::S, feat::S, cond::T, gain::T, K) where {T<:AbstractFloat, S} = TreeNodeGPU(left, right, feat, cond, gain, zeros(T,K), true)
TreeNodeGPU(pred::Vector{T}) where {T} = TreeNodeGPU(UInt32(0), UInt32(0), UInt32(0), zero(T), zero(T), pred, false)

"""
    single tree is made of a a vector of nodes
"""
struct TreeGPU{T<:AbstractFloat, S}
    nodes::Vector{TreeNodeGPU{T,S,Bool}}
end

# gradient-boosted tree is formed by a vector of trees
struct GBTreeGPU{T<:AbstractFloat, S}
    trees::Vector{TreeGPU{T,S}}
    params::EvoTypes
    metric::Metric
    K::S
    levels
end


text:
getreward(rewardcall::Type{pre_partial}, rewardfunc, mReward, tStateVectors,
	vChoiceOne, vChoiceTwo,
	cnt, ix, j, jprime, l, lprime) =
	rewardfunc(mReward[cnt.i_state],
					   getindex.(tStateVectors, (j, l, ix.I...)),
					   (vChoiceOne[jprime], vChoiceTwo[lprime]))

getreward(rewardcall::Type{jit}, rewardfunc, mReward, tStateVectors,
	vChoiceOne, vChoiceTwo,
	cnt, ix, j, jprime, l, lprime) =
	reward = rewardfunc(getindex.(tStateVectors, (j, l, ix.I...)),
						(vChoiceOne[jprime], vChoiceTwo[lprime]))

getreward(rewardcall::Type{pre}, rewardfunc, mReward, tStateVectors,
	vChoiceOne, vChoiceTwo,
	cnt, ix, j, jprime, l, lprime) =  mReward[cnt.i_choice, cnt.i_state]

_solve(p::DDP,
	mTransition::Array{Float64,2}, mReward::Union{Array{Float64}, Nothing},
	opts::SolverOptions) =
	_solve(p, mTransition, mReward, opts, p.tStateVectors,
		getchoicevars(p.tStateVectors, p.tChoiceVectors),
		getnonchoicevars(p.tStateVectors, p.tChoiceVectors),)

function _solve(p::DDP,
				mTransition::Array{Float64,2}, mReward::Union{Array{Float64}, Nothing},
				opts::SolverOptions,
				tStateVectors,
				tChoiceVectors,
				tOtherStateVectors)

	@unpack β, rewardfunc = p
	@unpack disp, disp_each_iter, max_iter, epsilon, rewardcall,
	monotonicity, concavity, intdim = opts
	monotonicity1, monotonicity2 = monotonicity
	concavity1, concavity2 = concavity

    (nChoiceOne, nChoiceTwo) = length.(tChoiceVectors)
    (vChoiceOne, vChoiceTwo) = tChoiceVectors

	nChoices = nChoiceOne * nChoiceTwo
    nStates = prod(length.(tStateVectors))
    nOtherStates = prod(length.(tOtherStateVectors))

	mValFun    = zeros((nChoices, nOtherStates)) # matrix so that can multiply with mTransition
	mValFunNew = zeros(nStates) # vector so that can index easily

    mPolFunInd1 = zeros(Int16, nStates)
    mPolFunInd2 = zeros(Int16, nStates)

    mβEV = zeros(nChoices, size(mTransition, 1)) # depends on intdim

    # VFI initialization
    maxDifference::Float64 = Inf # to initialize
    iteration::Int64 = 0 # initialize counter
	if β == 0.0
        tolerance = Inf
    elseif β == 1.0
        throw(ArgumentError("method invalid for beta = 1"))
    else
        tolerance = epsilon * (1-β) / (2*β)
    end

    # initialize for less memory allocation
	mValFunDiff = zeros(nStates)
    mValFunDiffAbs = zeros(nStates)

    liquidationvalue::Float64 = -Inf
    inactionvalue::Float64 = -Inf
    reward::Float64 = -Inf
	βEV::Float64 = -Inf

	valueHighSoFarOne::Float64 = -Inf
	valueProvisionalOne::Float64 = -Inf

	valueHighSoFarTwo::Float64 = -Inf
	valueProvisionalTwo::Float64 = -Inf

    iChoice1::Int16 = 0
    iChoice2::Int16 = 0

    vChoice2Start = ones(Int64, nChoiceOne)
    vChoice1Start = 1
	iChoice1inner = 0

	i::Int64 = 0 # outer loop for other state vars
	cnt = initialize_counter()

	# will need beta times transpose of transition matrix
	mTransition_βT = β * transpose(mTransition)

    # VFI
    while maxDifference > tolerance

        mul!(mβEV, mValFun, mTransition_βT)

		i = 0
		reset!(cnt)

		# @inbounds
		for ix in CartesianIndices(length.(tOtherStateVectors)) # other states
			i = i + 1
			cnt.i_exogstate += 1

            if monotonicity2
                vChoice2Start[:] .= 1
            end

            for l = 1:nChoiceTwo

				if monotonicity1
	                vChoice1Start = 1
	            end

                for j = 1:nChoiceOne # first state

                    valueHighSoFarTwo = -Inf
                    valueProvisionalTwo = -Inf
                    iChoice1 = 0
                    iChoice2 = 0

					cnt.i_state += 1
					cnt.i_choice = 0

					cnt.i_statechoice += (vChoice2Start[j]-1)*nChoiceOne # compensate if start later
					cnt.i_choice += (vChoice2Start[j]-1)*nChoiceOne # compensate if start later

                    for lprime = vChoice2Start[j]:nChoiceTwo

						# get optimal second choice variable conditional on first
						valueHighSoFarOne = -Inf
						valueProvisionalOne = -Inf
						iChoice1inner  = 0

						# find highest value for second state var
						cnt.i_statechoice += vChoice1Start-1 # compensate if start later
						cnt.i_choice += vChoice1Start-1
						for jprime = vChoice1Start:nChoiceOne

							cnt.i_statechoice += 1
							cnt.i_choice += 1

							reward = getreward(rewardcall, rewardfunc,
								mReward, tStateVectors,
								vChoiceOne, vChoiceTwo,
								cnt, ix, j, jprime, l, lprime)

							valueProvisionalOne = reward + mβEV[cnt.i_choice,
								getcounter(cnt, intdim)]

							if valueProvisionalOne >= valueHighSoFarOne
    							   valueHighSoFarOne = valueProvisionalOne
    							   iChoice1inner = jprime
	                        elseif concavity1
								adj = nChoiceOne - jprime # adjust counter if finish early
								cnt.i_statechoice += adj
								cnt.i_choice += adj
    							break
							end

						end # jprime

                        valueProvisionalTwo = valueHighSoFarOne

                        if valueProvisionalTwo > valueHighSoFarTwo
                            valueHighSoFarTwo = valueProvisionalTwo
                            iChoice1 = iChoice1inner
                            iChoice2 = lprime
                        elseif concavity2
							adj = (nChoiceTwo - lprime)*nChoiceOne # adjust counter if finish early
							cnt.i_statechoice += adj
							cnt.i_choice += adj
                            break
                        end

                    end #lprime

                    mValFunNew[cnt.i_state] = valueHighSoFarTwo
                    mPolFunInd1[cnt.i_state] = iChoice1
                    mPolFunInd2[cnt.i_state] = iChoice2

					if monotonicity1
                        vChoice1Start = iChoice1inner
                    end
					if monotonicity2
					  vChoice2Start[j] = iChoice2
					end

                end # j (first choice)

            end # l (second choice)

        end #i (stochastic variables)

        mValFunDiff .= mValFunNew .- @view(mValFun[:])
        mValFunDiffAbs .= abs.(mValFunDiff)
        maxDifference  = maximum(mValFunDiffAbs)

        # mcqueen!(mValFunNew, mValFun, β) #--> does not even converge somehow

        copyto!(mValFun, mValFunNew)

        # howards improvement?

        if disp
            display_iter(iteration, disp_each_iter, maxDifference)
        end

		iteration += 1

        if iteration > max_iter
            println("WARNING: maximum iterations exceeded")
            # println("Parameters used:")
			# println(p.params)
            break
        end

    end #while

    # println("Iteration ", string(iteration), ": Sup change is = ", string(maxDifference))

    mPolFun1 = vChoiceOne[mPolFunInd1]
    mPolFun2 = vChoiceTwo[mPolFunInd2]

    nNodes = length.(tStateVectors)
    meshPolFun1 = reshape(mPolFun1, tuple(nNodes...))
    meshPolFun2 = reshape(mPolFun2, tuple(nNodes...))

    meshValFun = reshape(mValFun, tuple(nNodes...))
    # meshExit   = reshape(mExit, tuple(nNodes...))

	return meshValFun, (meshPolFun1, meshPolFun2)
end # solve