text:
using FiniteDifferences
using BSON: @save, @load
import Base:zero

extract(X, n, b) = reduce(vcat, (xt[n,b] for xt in X))
unzip(a) = (getfield.(a, x) for x in fieldnames(eltype(a)))
count_params(M) = sum(length(p) for p in Flux.params(M))

function fd_test_grads(loss, params; order=10)
    zygote_grads = gradient(loss, params)

    for (id, p) in enumerate(params)
        function fd_loss!(ptild)
            old_p = copy(p)
            p .= ptild  # This changes the param in model
            J = loss()
            p .= old_p  # return to original
            return J
        end

        fd_grad = grad(central_fdm(order, 1), fd_loss!, p)[1];
        error = norm(fd_grad - zygote_grads[p])
        println("Param ", id, " \tsize ", size(p), "\t gradient err = ", error)
    end
end

# Why is inv not implemented for CUDA arrays ....
using CUDA, LinearAlgebra
function cuinv(m::AbstractArray{T}) where T
    if size(m, 1) != size(m, 2) throw(ArgumentError("Matrix not square.")) end
    A = CuArray(m)
    B = CuArray(Matrix{T}(I(size(A, 1))))
    A, ipiv = CUDA.CUSOLVER.getrf!(A)
    Matrix{T}(CUDA.CUSOLVER.getrs!('N', A, ipiv, B))
end

function estimate_lipschitz_lower(model; batches=1, seq_len=3000, plot_res=false, maxIter=300, step_size=1E-2, clip_at=1E-2, init_var=1E-3)

    # solver = PeacemanRachford(tol=1E-8, verbose=true)
    # (model::typeof(model))(x0, u) = model(x0, u, solver)

    x0 = init_state(model, batches)

    # u1 = collect(randn(model.nu, batches) for t in 1:seq_len)
    # u2 = collect(randn(model.nu, batches) for t in 1:seq_len)
    u1 = init_var * randn(seq_len, model.nu, batches)
    u2 = u1 .+ 1E-4 * randn(seq_len, model.nu, batches)

    unpack(ui) = [ui[t,:,:] for t in 1:seq_len]

    θ = Flux.Params([u1, u2, x0])

    opt = Flux.Optimiser(ClipValue(clip_at),
                        Flux.Optimise.ADAM(step_size),
                        Flux.Optimise.ExpDecay(1.0, 0.1, 200, 0.001)
                        )

    Lip() = norm(model(x0, unpack(u1))[1] - model(x0, unpack(u2))[1]) / (norm(u1 - u2))

    # Lip() = sum((model(x0, unpack(u1))[1] - model(x0, unpack(u2))[1]).^2) / sum((u1 - u2).^2)

    lips = []
    for iter in 1:maxIter
        L, back = Zygote.pullback(() -> -Lip(), θ)
        ∇L = back(one(L))
        update!(opt, θ, ∇L)
        append!(lips, L)
        println("Iter: ", iter, "\t L: ", -L, "\t η: ", opt[3].eta)
    end

    if plot_res == true

        du = unpack(u2) - unpack(u1)
        # Δu = 0.05 * du / norm(du)
        Δu = unpack(u2) - unpack(u1)
        yest1 = model(x0, unpack(u1))[1]
        yest2 = model(x0, unpack(u2))[1]
        # yest2 = model(x0, unpack(u1) + Δu)[1]
        Δy = yest1 - yest2

        data = (dy = Δy, du = Δu, Lip = norm(Δy) / norm(Δu),  u1 = unpack(u1), u2 = unpack(u2), yest1 = yest1, yest2 = yest2, model = model)
        @save "./results/sys_id/f16_adversarial_perturbations/ffren_10.bson" data

        p1 = plot(extract(Δy, 2, 1))
        plot!(extract(Δu, 1, 1))
        plot!(extract(Δu, 2, 1))
        plot(p1, p2, layout=(2, 1))
    end
    return maximum(-lips)
end

# return vectorized parameter list
function vec(param_list)
    return vcat((p[:] for p in param_list)...)
end

# Inverse operation of vec. Stores vectorized parameter into parameter list.
function store_vec!(param_list, θ)
    counter = 1
    for p in param_list
        θᵢ = θ[counter:counter + length(p) - 1]
        p .= reshape(θᵢ, size(p))
        counter += length(θᵢ)
    end
end

max_stars_repo_name: imanchester/REN

max_issues_repo_name: imanchester/REN

max_forks_repo_name: imanchester/REN


text:
using JuliaCommunityStatistics
using Test

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

ext: jl


text:
using Parameters

@component struct Test1
    p::Int
end

Test1() = Test1(0)

@component @with_kw struct Test2
    p::Int = 1
end

@pooled_component struct Test3
    p::Int
end
Test3() = Test3(1)

@pooled_component @with_kw struct Test4
    p::Int = 1
end

const c1 = Overseer.component_type(Test1){Test1}()
const c2 = Overseer.component_type(Test2){Test2}()
const c3 = Overseer.component_type(Test3){Test3}()
const c4 = Overseer.component_type(Test4){Test4}()

const entities1 = [Entity(i) for i in 2:2:10]
const entities2 = [Entity(i) for i in 10:3:20]
const entities3 = [Entity(i) for i in 3:10]
const entities4 = [Entity(1)]

@component struct ParametricComp{T}
    x::T
end

@component @with_kw struct ParametricCompKw{T}
    x::T = 1.0
end

@component struct ParametricCompFunc{T}
    x::T
    function ParametricCompFunc()
        new{Float64}(1.0)
    end
    function ParametricCompFunc{T}() where {T}
        new{T}(T(1.0))
    end
end

@testset "Basic Component definitions" begin
    @test Overseer.component_type(Test1) == Component
    @test Overseer.component_type(Test2) == Component
    @test Overseer.component_type(Test3) == PooledComponent
    @test Overseer.component_type(Test4) == PooledComponent

    for (c, es) in zip((c1, c2, c3, c4), (entities1, entities2, entities3, entities4))
        for e in es
            c[e] = eltype(c)()
        end
    end

    for (c, es) in zip((c1, c2, c3, c4), (entities1, entities2, entities3, entities4))
        for e in es
            @test in(e, c)
        end
    end

    @test ParametricComp <: Overseer.ComponentData

    @test ParametricCompKw <: Overseer.ComponentData
    @test ParametricCompKw().x == 1.0

    @test ParametricCompFunc <: Overseer.ComponentData
    @test ParametricCompFunc().x == 1.0
    @test ParametricCompFunc{Int}().x === 1

end

@testset "Component iteration" begin
    t = 0
    for e in @entities_in(((c1 && c3) || c4) && !c2)
        t += e.id
    end
    @test t == 4+6+8+1

    t = 0
    for e in @entities_in((c1 || c3) && !c2)
        if e in c1
            t += e.id
        end
        if e in c3
            t += c3[e].p
        end
    end
    @test t == 27

    t = 0
    for e in @entities_in(c1)
        t += e.id
    end
    @test t == sum(2:2:10)
end

@testset "Component and entity manipulation" begin
    @test pop!(c1, Entity(10)) == Test1()

    @test length(c1) == length(entities1) - 1

    @test pop!(c2, Entity(10)) == Test2()

    @test length(c2) == length(entities2) - 1
    @test c1[1] == Test1()
    @test c3[1] == Test3()

    c2[Entity(13)] = Test2(50)
    @test c2[Entity(13)] == Test2(50)

    c3[Entity(13)] = Test3(50)
    @test c3[Entity(13)] == Test3(50)

    pop!(c3, Entity(13))
    @test !in(Entity(13), c3)

    empty!(c1)
    @test isempty(c1)

    empty!(c3)
    @test isempty(c3)

    # swap_ordering
    c2[Entity(12)] = Test2()

    @test_throws BoundsError swap_order!(c2, Entity(14), Entity(15))
    @test_throws BoundsError swap_order!(c2, Entity(13), Entity(14))

    orig1 = c2[Entity(12)]
    orig2 = c2[Entity(13)]

    orig_id1 = c2.indices[12]
    orig_id2 = c2.indices[13]

    swap_order!(c2, Entity(12), Entity(13))
    @test c2[Entity(12)] == orig1
    @test c2[Entity(13)] == orig2

    @test c2.indices[13] == orig_id1
    @test c2.indices[12] == orig_id2

    c3[Entity(12)] = Test3()
    c3[Entity(13)] = Test3(50)

    @test_throws BoundsError swap_order!(c3, Entity(14), Entity(15))
    @test_throws BoundsError swap_order!(c3, Entity(13), Entity(14))

    orig1 = c3[Entity(12)]
    orig2 = c3[Entity(13)]

    orig_id1 = c3.indices[12]
    orig_id2 = c3.indices[13]

    swap_order!(c3, Entity(12), Entity(13))
    @test c3[Entity(12)] == orig1
    @test c3[Entity(13)] == orig2

    @test c3.indices[13] == orig_id1
    @test c3.indices[12] == orig_id2
end

@component @with_kw struct Spatial_t
    position::NTuple{3, Float64} = (1.0,1.0,1.0)
end
@component struct Spring_t
    center::NTuple{3, Float64}
end

@component struct Spatial1_t
    position::NTuple{3, Float64}
end

@component mutable struct Spatial2_t
    position::NTuple{3, Float64}
end

@testset "iteration syntax" begin
    m = Ledger()
    for i = 1:5
        Entity(m, Spatial_t(), Spring_t((i, 1.0, 2.0)), Spatial1_t((3.0, 1.0, 2.0)), Spatial2_t((1.0, 1.0, 1.0)))
    end
    for e in @entities_in(m, Spatial_t && Spring_t)
        e[Spatial_t] = Spatial_t(e.position .+ e.center)
    end
    @test m[Spatial_t][Entity(1)].position == (2.0, 2.0, 3.0)
    for e in @entities_in(m, Spatial_t && Spatial1_t)
        e[Spatial_t] = Spatial_t(e[Spatial_t].position .+ e[Spatial1_t].position)
    end
    @test m[Spatial_t][Entity(1)].position ==  (5.0, 3.0, 5.0)
    @test_throws ErrorException begin
        for e in @entities_in(m, Spatial_t && Spatial1_t)
            e[Spatial_t] = Spatial_t(e.position .+ e[Spatial1_t].position)
        end
    end
    for e in @entities_in(m, Spatial2_t )
        e.position = 2 .* e.position
    end
    @test m[Spatial2_t][Entity(1)].position == (2.0, 2.0, 2.0)

    for e in @entities_in(m, Spatial2_t )
        t = e[Spatial2_t]
        t.position = (0.0,0.0,0.0)
    end
    @test m[Spatial2_t][Entity(1)].position == (0.0, 0.0, 0.0)

end
# Issue 4: collect() and iterator length
@testset "collect" begin
    e1 = Entity(1)
    e2 = Entity(2)
    e3 = Entity(3)

    comp1 = Overseer.component_type(Test1){Test1}()
    comp2 = Overseer.component_type(Test2){Test2}()
    comp1[e1] = Test1(1)
    comp1[e2] = Test1(1); comp2[e2] = Test2(1)
    comp2[e3] = Test2(1)

    iter = @entities_in(comp1 && comp2)
    es = collect(iter)
    @test getfield.(es, :e) == [e2]
    @test eltype(es) == Overseer.EntityState{Tuple{Component{Test1},Component{Test2}}}
end

@testset "PooledComponent" begin
    @pooled_component struct Test5
        x::Int
    end

    @test Overseer.component_type(Test5) == Overseer.PooledComponent
    c5 = Overseer.component_type(Test5){Test5}()

    p1 = Entity(1)
    p2 = Entity(2)
    c5[p1] = Test5(1)
    c5[p2] = Test5(2)

    entities = [Entity(i) for i in 3:10]
    for (i, e) in enumerate(entities)
        c5[e] = (p1, p2)[mod1(i, 2)]
    end

    # check created values
    count = 0
    _sum = 0
    for e in @entities_in(c5)
        count += 1
        _sum += c5[e].x
    end
    @test count == 10
    @test _sum == 15
    @test c5.pool_size == [5, 5]

    # check for no duplication
    @test length(c5.data) == 2

    # Check some basics
    @test p1 in c5
    @test pop!(c5, p1) == Test5(1)
    @test !(p1 in c5)
    @test length(c5) == 9
    @test c5[p2] == Test5(2)
    @test !isempty(c5)
    @test c5.pool_size == [4, 5]

    count = 0
    _sum = 0
    for e in @entities_in(c5)
        count += 1
        _sum += c5[e].x
    end
    @test count == 9
    @test _sum == 14

    # adjust parent value of pool
    c5[parent(p2)] = Test5(1)
    count = 0
    _sum = 0
    for e in @entities_in(c5)
        count += 1
        _sum += c5[e].x
    end
    @test count == 9
    @test _sum == 9
    @test c5.pool_size == [4, 5]

    # adjust single value
    c5[p2] = Test5(2)
    count = 0
    _sum = 0
    for e in @entities_in(c5)
        count += 1
        _sum += c5[e].x
    end
    @test count == 9
    @test _sum == 10
    @test length(c5.data) == 3
    @test c5.pool_size == [4, 4, 1]

    Overseer.make_unique!(c5)
    @test length(c5.data) == 2
    @test c5.pool_size == [1, 8]
    @test c5.data == [Test5(2), Test5(1)]

    # remove all entites of a pool
    for i in 3:10
        pop!(c5, Entity(i))
    end
    @test length(c5.data) == 1
    @test length(c5) == 1
    @test c5.pool_size == [1]

    empty!(c5)
    @test isempty(c5)
    @test c5.pool_size == Int[]

    e1 = Entity(1)
    e2 = Entity(2)
    e3 = Entity(3)
    c5[e1] = Test5(1)
    c5[e3] = e1
    c5[e2] = Test5(1)

    @test length(collect(entity_pool(c5, 1))) == c5.pool_size[1]
    @test collect(entity_pool(c5, 1)) == [Entity(1), Entity(3)]

    c1 = Component{Test1}()
    c1[e1] = Test1(1)
    c1[e2] = Test1(2)
    c1[e3] = Test1(3)
    order = Entity[]
    for e in @entities_in(entity_pool(c5, 1) && c1)
        push!(order, e.e)
    end
    for e in @entities_in(entity_pool(c5, 2) && c1)
        push!(order, e.e)
    end
    @test order == [Entity(1), Entity(3), Entity(2)]

    @test parent(c5, 2) == Entity(2)
    @test parent(c5, 1) == Entity(1)
    @test parent(c5, Entity(3)) == Entity(1)

end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/display.jl
# Functions used to print information about solution routine to console.
# Printing is done with reference to the input problem is there is any
# ambiguity.
#############################################################################

"""
$(FUNCTIONNAME)

Prints solution information for the B&B problem. Displays first node found, solution value,
solution, and time spent solving subproblems.
"""
function print_solution!(m::GlobalOptimizer)
    if _verbosity(m) > 0
        println(" ")
        if m._end_state == GS_OPTIMAL
            println("Empty Stack: Exhaustive Search Finished")
        elseif m._end_state == GS_INFEASIBLE
            println("Empty Stack: Infeasible")
        elseif m._end_state == GS_NODE_LIMIT
            println("Node Limit Exceeded")
        elseif m._end_state == GS_ITERATION_LIMIT
            println("Maximum Iteration Exceeded")
        elseif m._end_state == GS_RELATIVE_TOL
            println("Relative Tolerance Achieved")
        elseif m._end_state == GS_ABSOLUTE_TOL
            println("Absolute Tolerance Achieved")
        elseif m._end_state == GS_TIME_LIMIT
            println("Time Limit Exceeded")
        end
        println("First Solution Found at Node $(m._first_solution_node)")
        if !_is_input_min(m)
            println("LBD = $(MOI.get(m, MOI.ObjectiveBound()))")
            println("UBD = $(MOI.get(m, MOI.ObjectiveValue()))")
        else
            println("LBD = $(MOI.get(m, MOI.ObjectiveBound()))")
            println("UBD = $(MOI.get(m, MOI.ObjectiveValue()))")
        end
        println("Solution is :")
        if m._feasible_solution_found
            k = _obj_var_slack_added(m) ? 1 : 0
            for i = 1:(m._input_problem._variable_count - k)
                println("    X[$i] = $(m._continuous_solution[i])")
            end
        end
        println(" ")
     end
     return
end

"""
$(FUNCTIONNAME)

Prints node information for the B&B problem. Node id, bound, and interval box.
"""
function print_node!(m::GlobalOptimizer)
    if _verbosity(m) >= 3
        n = m._current_node
        bound = _is_input_min(m) ? n.lower_bound : -n.lower_bound
        k = length(n) - (_obj_var_slack_added(m) ? 1 : 0)
        println(" ")
        println("Node ID: $(n.id), Lower Bound: $(bound)")
        println("Lower Variable Bounds: $(n.lower_variable_bounds[1:k])")
        println("Upper Variable Bounds: $(n.upper_variable_bounds[1:k])")
        println(" ")
    end
    return
end

"""
$(FUNCTIONNAME)

Prints the iteration information based on verbosity. The header is displayed
every `header_interval`, the iteration info is displayed every `iteration_interval`.
"""
function print_iteration!(m::GlobalOptimizer)

    if _verbosity(m) > 0

        # prints header line every B.hdr_intv times
        if mod(m._iteration_count, m._parameters.header_iterations) === 0 || m._iteration_count === 1
            println("-----------------------------------------------------------------------------------------------------------------------------")
            println("|  Iteration #  |     Nodes    | Lower Bound  |  Upper Bound  |      Gap     |     Ratio    |     Time     |    Time Left   |")
            println("-----------------------------------------------------------------------------------------------------------------------------")
        end

        # prints iteration summary every B.itr_intv times
        if mod(m._iteration_count, m._parameters.output_iterations) === 0

            print_str = "| "

            max_len = 12
            temp_str = string(m._iteration_count)
            len_str = length(temp_str)
            print_str *= (" "^(max_len - len_str))*temp_str*"  | "

            max_len = 12
            temp_str = string(m._node_count)
            len_str = length(temp_str)
            print_str *= (" "^(max_len - len_str))*temp_str*" | "

            max_len = 12
            if _is_input_min(m)
                lower = m._global_lower_bound
                upper = m._global_upper_bound
            else
                lower = m._global_lower_bound
                upper = m._global_upper_bound
            end
            #temp_str = string(round(lower, sigdigits = 5))
            #temp_str = string(lower, sigdigits = 3))
            temp_str = @sprintf "%.2E" lower
            len_str = length(temp_str)
            print_str *= (" "^(max_len - len_str))*temp_str*" | "

            #temp_str = formatted(upper, PRINTING_IOFORMAT, ndigits=4, charset=PRINTING_CHARSET)
            #temp_str = string(upper, sigdigits = 3))
            temp_str = @sprintf "%.2E" upper
            len_str = length(temp_str)
            print_str *= (" "^(max_len - len_str))*temp_str*"  |"

            max_len = 12
            #temp_str = string(round(abs(x._global_upper_bound - x._global_lower_bound), sigdigits = 3))
            temp_str = @sprintf "%.2E" abs(m._global_upper_bound - m._global_lower_bound)
            len_str = length(temp_str)
            print_str *= (" "^(max_len - len_str))*temp_str*"  | "

            max_len = 12
            #temp_str = string(round(relative_gap(x._global_lower_bound, x._global_upper_bound), sigdigits = 3))
            temp_str = @sprintf "%.2E" relative_gap(m._global_lower_bound, m._global_upper_bound)
            len_str = length(temp_str)
            print_str *= (" "^(max_len - len_str))*temp_str*" | "

            max_len = 12
            #temp_str = string(round(x._run_time, sigdigits = 3))
            temp_str = @sprintf "%.2E" m._run_time
            len_str = length(temp_str)
            print_str *= (" "^(max_len - len_str))*temp_str*" |  "

            max_len = 12
            #temp_str = string(round(x._time_left, sigdigits = 4))
            temp_str = @sprintf "%.2E" m._time_left
            len_str = length(temp_str)
            print_str *= (" "^(max_len - len_str))*temp_str*"  |"

            println(print_str)
        end
    end

    return
end

"""
$(FUNCTIONNAME)

Prints the results of a single bounding problem.
"""
function print_results!(m::GlobalOptimizer, flag::Bool)
    if _verbosity(m) > 1
        k = length(m._lower_solution) - (_obj_var_slack_added(m) ? 1 : 0)
        println(" ")
        if flag
            if _is_input_min(m)
                print("Lower Bound (First Iteration): $(m._lower_objective_value),")
            else
                print("Upper Bound (First Iteration): $(m._lower_objective_value),")
            end
            print(" Solution: $(m._lower_solution[1:k]), Feasibility: $(m._lower_feasibility)\n")
            println("Termination Status Code: $(m._lower_termination_status)")
            println("Result Code: $(m._lower_primal_status)")
        else
            if _is_input_min(m)
                print("Upper Bound: $(m._upper_objective_value), ")
            else
                print("Lower Bound: $(m._upper_objective_value), ")
            end
            print(" Solution: $(m._upper_solution[1:k]), Feasibility: $(m._upper_feasibility)\n")
            println("Termination Status Code: $(m._upper_termination_status)")
            println("Result Code: $(m._upper_result_status)")
        end
    end
    return
end

"""
$(FUNCTIONNAME)

Prints the iteration information based on verbosity. The header is displayed
every `header_interval`, the iteration info is displayed every `iteration_interval`.
"""
function print_preamble!(m::GlobalOptimizer)
    if _verbosity(m) >= 3
        if !_is_input_min(m) && isone(m._iteration_count)
            println(" ")
            println("For maximization problems a max(f) = -min(-f) transformation is applied.")
            println("Objectives values for each subproblem are the negative value of the objective")
            println("in the original problem and reconciled after branch and bound terminates.")
            println(" ")
        end
    end
    return
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
abstract type AbstractOracle{T<:AbstractTransitionSystem} end

struct Oracle{T,O,L}
    system::T
    oracle_function::O
    labelf::L
end

"""
    Oracle(system, oracle_function; label=untyped)

Create an oracle for predicting gold transitions in dependency parsing.

`system` is a transition system, defining configurations and valid transitions.

`oracle_function` is called on a paraser configuration and tree for gold predictions:

    oracle(cfg, tree, label)

`label` is a function that's called on the gold tokens for that parameters of arcs.
"""
function Oracle(system, oracle_function; label=untyped)
    Oracle(system, oracle_function, label)
end

(oracle::Oracle)(tree::DependencyTree, policy=NeverExplore()) =
    OracleSequence(oracle, tree, policy)

(oracle::Oracle)(cfg, tree::DependencyTree) =
    oracle.oracle_function(cfg, tree, oracle.labelf)

initconfig(oracle::Oracle, tree) = initconfig(oracle.system, tree)

function transitions(cfg, gold::DependencyTree, oracle::Oracle)
    A = possible_transitions(cfg, gold, oracle.labelf)
    G = oracle(cfg, gold)
    return A, G
end

"""
    OracleSequence(oracle, tree, policy=NeverExplore())

A "gold" sequence of parser configurations and transitions to build `tree`.

The sequence of transitions is performed according to

`policy` is a function that determines whether or not "incorrect" transitions
are explored. It will be called like so: `policy(
"""
struct OracleSequence{T,P}
    oracle::Oracle{T}
    tree::DependencyTree
    policy::P

    function OracleSequence(oracle::Oracle, tree::DependencyTree, policy=NeverExplore())
        T = oracle.system
        if projective_only(T) && !is_projective(tree)
            UnparsableTree(NonProjectiveGraphError(tree))
        else
            new{typeof(T),typeof(policy)}(oracle, tree, policy)
        end
    end
end

function choose_transition(ts::OracleSequence, cfg)
    A, G = transitions(cfg, ts.tree, ts.oracle)
    t = ts.policy(cfg, A, G)
    return t
end

Base.IteratorSize(ts::OracleSequence) = Base.SizeUnknown()

function Base.iterate(ts::OracleSequence, state=nothing)
    if state === nothing
        state = initconfig(ts.oracle.system, ts.tree)
    end
    if isfinal(state)
        return nothing
    else
        A, G = transitions(state, ts.tree, ts.oracle)
        isempty(A) && return nothing
        t = ts.policy(state, A, G)
        return ((state, G), t(state))
    end
end

"""
   UnparsableTree

A dependency tree that an oracle cannot parse.
"""
struct UnparsableTree
    err
end
Base.IteratorSize(pairs::UnparsableTree) = Base.HasLength()
Base.iterate(pairs::UnparsableTree, state...) = nothing
Base.length(::UnparsableTree) = 0

abstract type AbstractExplorationPolicy end

choose(rng, x::TransitionOperator) = rand(rng, [x])

function choose(rng, x)
    isempty(x) ? nothing : rand(rng, x)
end

no_model(cfg, A, G) = nothing

"""
    AlwaysExplore()

Policy for always exploring sub-optimal transitions.

If `model` predicts a legal transition, apply it. Otherwise,
sample from the possible transitions (without regard to the oracle transitions)
according to `rng`.
"""
struct AlwaysExplore{R<:AbstractRNG,M} <: AbstractExplorationPolicy
    rng::R
    model::M
end
AlwaysExplore() = AlwaysExplore(GLOBAL_RNG, no_model)
AlwaysExplore(rng::AbstractRNG) = AlwaysExplore(rng, no_model)
AlwaysExplore(model) = AlwaysExplore(GLOBAL_RNG, model)

(::AlwaysExplore)() = true
function (p::AlwaysExplore)(cfg, A::AbstractVector, G)
    t = p.model(cfg, A, G)
    is_possible(t, cfg) ?  t : choose(p.rng, A)
end

Base.show(io::IO, ::AlwaysExplore) = print(io, "AlwaysExplore")

"""
    NeverExplore()

Policy for never exploring sub-optimal transitions.

If `model` predicts a gold transition, apply it. Otherwise,
choose from the gold transitions according to `rng`.
"""
struct NeverExplore{R<:AbstractRNG,M} <: AbstractExplorationPolicy
    rng::R
    model::M
end
NeverExplore() = NeverExplore(GLOBAL_RNG, no_model)
NeverExplore(rng::AbstractRNG) = NeverExplore(rng, no_model)
NeverExplore(model, rng=GLOBAL_RNG) = NeverExplore(rng, model)

(::NeverExplore)() = false
function (p::NeverExplore)(cfg, A, G)
    t = p.model(cfg, A, G)
    is_possible(t, cfg) && t in G ? t : choose(p.rng, G)
end

Base.show(io::IO, ::NeverExplore) = print(io, "NeverExplore")

"""
    ExplorationPolicy(k, p)

Simple exploration policy from Goldberg & Nivre, 2012. Explores at rate `p`.

With rate `p`, follow `model`'s prediction if legal, or choose from the possible
transitions according to `rng` if the prediction can't be followed.
With probability 1 -`p`, choose from gold transitions according to `rng`.
"""
struct ExplorationPolicy{R<:AbstractRNG,M} <: AbstractExplorationPolicy
    p::Float64
    rng::R
    model::M

    function ExplorationPolicy(p, rng::AbstractRNG=GLOBAL_RNG, model=no_model)
        @assert 0 <= p <= 1
        new{typeof(rng),typeof(model)}(p, rng, model)
    end
    function ExplorationPolicy(p, model, rng::AbstractRNG=GLOBAL_RNG)
        @assert 0 <= p <= 1
        new{typeof(rng),typeof(model)}(p, rng, model)
    end
end

(p::ExplorationPolicy)() =
    rand(p.rng) >= 1 - p.p

(p::ExplorationPolicy)(cfg, A::AbstractVector, G::TransitionOperator) = p(cfg, A, [G])

function (p::ExplorationPolicy)(cfg, A::AbstractVector, G::AbstractVector)
    t = p.model(cfg, A, G)
    explore = rand(p.rng) >= 1 - p.p
    if explore
        return is_possible(t, cfg) ? t : rand(A)
    else
        return rand(G)
    end
end


text:
import Base.LinAlg.BLAS
using CUBLAS
using CUDArt
using Base.Test

m = 20
n = 35
k = 13

function blasabs(A)
    return abs.(real(A)) + abs.(imag(A))
end

#################
# level 1 tests #
#################

@testset "blascopy!" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        A = convert(Vector{elty}, collect(1:m))
        if elty <: Base.LinAlg.BlasComplex
            A += im*A
        end
        @test ndims(A) == 1
        n1 = length(A)
        d_A = CudaArray(A)
        d_B = CudaArray(elty, n1)
        CUBLAS.blascopy!(n,d_A,1,d_B,1)
        B = to_host(d_B)
        @test A == B
    end
end

@testset "scal!" begin
    function test_scal!{T}(alpha,A::Array{T})
        @test ndims(A) == 1
        n1 = length(A)
        d_A = CudaArray(A)
        CUBLAS.scal!(n1,alpha,d_A,1)
        A1 = to_host(d_A)
        @test alpha*A ≈ A1

        d_A = CudaArray(A)
        d_As = CUBLAS.scale(d_A, alpha)
        A1 = to_host(d_As)
        @test alpha*A ≈ A1

        CUBLAS.scale!(d_A, alpha)
        A1 = to_host(d_As)
        @test alpha*A ≈ A1
    end
    test_scal!(2.0f0,Float32[1:m;])
    test_scal!(2.0,Float64[1:m;])
    test_scal!(1.0f0+im*1.0f0,Float32[1:m;]+im*Float32[1:m;])
    test_scal!(1.0+im*1.0,Float64[1:m;]+im*Float64[1:m;])
    test_scal!(2.0f0,Float32[1:m;]+im*Float32[1:m;])
    test_scal!(2.0,Float64[1:m;]+im*Float64[1:m;])
end

@testset "dot" begin
    @testset for elty in [Float32, Float64]
        A = convert(Vector{elty}, collect(1:m))
        B = convert(Vector{elty}, collect(1:m))
        @test ndims(A) == 1
        @test ndims(B) == 1
        @test length(A) == length(B)
        n1 = length(A)
        d_A = CudaArray(A)
        d_B = CudaArray(B)
        cuda_dot1 = CUBLAS.dot(n1,d_A,1,d_B,1)
        cuda_dot2 = CUBLAS.dot(d_A,d_B)
        host_dot = dot(A,B)
        @test host_dot ≈ cuda_dot1
        @test host_dot ≈ cuda_dot2

        #d_A = CudaArray(A)
        #d_B = CudaArray(B)
        #cuda_dot3 = CUBLAS.dot(d_A, 3:5, d_B, 5:7)
        #host_dot3 = dot(A, 3:5, B, 5:7)
        #@test_approx_eq(cuda_dot3, host_dot3)
    end
end

@testset "dotu" begin
    @testset for elty in [Complex64, Complex128]
        A = rand(elty, m)
        B = rand(elty, m)
        @test ndims(A) == 1
        @test ndims(B) == 1
        @test length(A) == length(B)
        n1 = length(A)
        d_A = CudaArray(A)
        d_B = CudaArray(B)
        cuda_dot1 = CUBLAS.dotu(n1,d_A,1,d_B,1)
        cuda_dot2 = CUBLAS.dotu(d_A,d_B)
        host_dot = A.'*B
        if VERSION < v"0.6.0-dev.2074" # julia PR #19670
            @test host_dot[1] ≈ cuda_dot1
            @test host_dot[1] ≈ cuda_dot2
            @test host_dot ≈ (d_A.'*d_B)
        else
            @test host_dot ≈ cuda_dot1
            @test host_dot ≈ cuda_dot2
            @test host_dot ≈ (d_A.'*d_B)[1]
        end
    end
end

@testset "dotc" begin
    @testset for elty in [Complex64, Complex128]
        A = rand(elty, m)
        B = rand(elty, m)
        @test ndims(A) == 1
        @test ndims(B) == 1
        @test length(A) == length(B)
        n1 = length(A)
        d_A = CudaArray(A)
        d_B = CudaArray(B)
        cuda_dot1 = CUBLAS.dotc(n1,d_A,1,d_B,1)
        cuda_dot2 = CUBLAS.dotc(d_A,d_B)
        host_dot = A'*B
        if VERSION < v"0.6.0-dev.2074" # julia PR #19670
            @test host_dot[1] ≈ cuda_dot1
            @test host_dot[1] ≈ cuda_dot2
            @test host_dot ≈ (d_A'*d_B)
        else
            @test host_dot ≈ cuda_dot1
            @test host_dot ≈ cuda_dot2
            @test host_dot ≈ (d_A'*d_B)[1]
        end
    end
end

@testset "nrm2" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        A = rand(elty, m)
        @test ndims(A) == 1
        n1 = length(A)
        d_A = CudaArray(A)
        cuda_nrm2_1 = CUBLAS.nrm2(n1,d_A,1)
        cuda_nrm2_2 = CUBLAS.nrm2(d_A)
        cuda_nrm2_3 = norm(d_A)
        host_nrm2 = norm(A)
        @test host_nrm2 ≈ cuda_nrm2_1
        @test host_nrm2 ≈ cuda_nrm2_2
        @test host_nrm2 ≈ cuda_nrm2_3
    end
end

@testset "asum" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        A = rand(elty, m)
        @test ndims(A) == 1
        n1 = length(A)
        d_A = CudaArray(A)
        cuda_asum1 = CUBLAS.asum(n1,d_A,1)
        cuda_asum2 = CUBLAS.asum(d_A)
        host_asum = sum(abs.(real(A)) + abs.(imag(A)))
        @test host_asum ≈ cuda_asum1
        @test host_asum ≈ cuda_asum2
    end
end

@testset "axpy!" begin
    # test axpy!
    function test_axpy!_1(alpha,A,B)
        @test length(A) == length(B)
        n1 = length(A)
        d_A = CudaArray(A)
        d_B1 = CudaArray(B)
        CUBLAS.axpy!(n1,alpha,d_A,1,d_B1,1)
        B1 = to_host(d_B1)
        host_axpy = alpha*A + B
        @test host_axpy ≈ B1
    end
    test_axpy!_1(2.0f0,rand(Float32,m),rand(Float32,m))
    test_axpy!_1(2.0,rand(Float64,m),rand(Float64,m))
    test_axpy!_1(2.0f0+im*2.0f0,rand(Complex64,m),rand(Complex64,m))
    test_axpy!_1(2.0+im*2.0,rand(Complex128,m),rand(Complex128,m))

    function test_axpy!_2(alpha,A,B)
        @test length(A) == length(B)
        n1 = length(A)
        d_A = CudaArray(A)
        d_B1 = CudaArray(B)
        CUBLAS.axpy!(alpha,d_A,d_B1)
        B1 = to_host(d_B1)
        host_axpy = alpha*A + B
        @test host_axpy ≈ B1
    end
    test_axpy!_2(2.0f0,rand(Float32,m),rand(Float32,m))
    test_axpy!_2(2.0,rand(Float64,m),rand(Float64,m))
    test_axpy!_2(2.0f0+im*2.0f0,rand(Complex64,m),rand(Complex64,m))
    test_axpy!_2(2.0+im*2.0,rand(Complex128,m),rand(Complex128,m))

    #=function test_axpy!_3(alpha,A,B)
        @test length(A) == length(B)
        n1 = length(A)
        d_A = CudaArray(A)
        d_B1 = CudaArray(B)
        CUBLAS.axpy!(alpha,d_A,1:2:n1,d_B1,1:2:n1)
        B1 = to_host(d_B1)
        host_axpy = B
        host_axpy[1:2:n1] = alpha*A[1:2:n1] + B[1:2:n1]
        @test_approx_eq(host_axpy,B1)
    end
    test_axpy!_3(2.0f0,rand(Float32,m),rand(Float32,m))
    test_axpy!_3(2.0,rand(Float64,m),rand(Float64,m))
    test_axpy!_3(2.0f0+im*2.0f0,rand(Complex64,m),rand(Complex64,m))
    test_axpy!_3(2.0+im*2.0,rand(Complex128,m),rand(Complex128,m))

    function test_axpy!_4(alpha,A,B)
        @test length(A) == length(B)
        n1 = length(A)
        d_A = CudaArray(A)
        d_B1 = CudaArray(B)
        r = 1:div(n1,2)
        CUBLAS.axpy!(alpha,d_A,r,d_B1,r)
        B1 = to_host(d_B1)
        host_axpy = B
        host_axpy[r] = alpha*A[r] + B[r]
        @test_approx_eq(host_axpy,B1)
    end
    test_axpy!_4(2.0f0,rand(Float32,m),rand(Float32,m))
    test_axpy!_4(2.0,rand(Float64,m),rand(Float64,m))
    test_axpy!_4(2.0f0+im*2.0f0,rand(Complex64,m),rand(Complex64,m))
    test_axpy!_4(2.0+im*2.0,rand(Complex128,m),rand(Complex128,m))=#
end

@testset "iamax and iamin" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        A = rand(elty, m)
        n1 = length(A)
        d_A = CudaArray(A)
        Aabs = blasabs(A)
        imin1 = CUBLAS.iamin(n1,d_A,1)
        imax1 = CUBLAS.iamax(n1,d_A,1)
        imin2 = CUBLAS.iamin(d_A)
        imax2 = CUBLAS.iamax(d_A)
        host_imin = indmin(Aabs)
        host_imax = indmax(Aabs)
        @test imin1 == imin2 == host_imin
        @test imin1 == imin2 == host_imin
    end
end

#################
# level 2 tests #
#################

@testset "gemv!" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        alpha = convert(elty,1)
        beta = convert(elty,1)
        A = rand(elty,m,n)
        d_A = CudaArray(A)

        # test y = A*x + y
        x = rand(elty,n)
        d_x = CudaArray(x)
        y = rand(elty,m)
        d_y = CudaArray(y)
        y = A*x + y
        CUBLAS.gemv!('N',alpha,d_A,d_x,beta,d_y)
        h_y = to_host(d_y)
        @test y ≈ h_y
        A_mul_B!(d_y,d_A,d_x)
        h_y = to_host(d_y)
        @test h_y ≈ A*x

        # test x = A.'*y + x
        x = rand(elty,n)
        d_x = CudaArray(x)
        y = rand(elty,m)
        d_y = CudaArray(y)
        x = A.'*y + x
        CUBLAS.gemv!('T',alpha,d_A,d_y,beta,d_x)
        h_x = to_host(d_x)
        @test x ≈ h_x
        At_mul_B!(d_x,d_A,d_y)
        h_x = to_host(d_x)
        @test h_x ≈ A.'*y

        # test x = A'*y + x
        x = rand(elty,n)
        d_x = CudaArray(x)
        y = rand(elty,m)
        d_y = CudaArray(y)
        x = A'*y + x
        CUBLAS.gemv!('C',alpha,d_A,d_y,beta,d_x)
        h_x = to_host(d_x)
        @test x ≈ h_x
        Ac_mul_B!(d_x,d_A,d_y)
        h_x = to_host(d_x)
        @test h_x ≈ A'*y
    end
end

@testset "gemv" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        alpha = convert(elty,2)
        A = rand(elty,m,n)
        d_A = CudaArray(A)
        # test y = alpha*(A*x)
        x = rand(elty,n)
        d_x = CudaArray(x)
        y1 = alpha*(A*x)
        y2 = A*x
        d_y1 = CUBLAS.gemv('N',alpha,d_A,d_x)
        d_y2 = CUBLAS.gemv('N',d_A,d_x)
        h_y1 = to_host(d_y1)
        h_y2 = to_host(d_y2)
        @test y1 ≈ h_y1
        @test y2 ≈ h_y2
        @test y2 ≈ to_host(d_A * d_x)

        # test x = alpha*(A.'*y)
        y = rand(elty,m)
        d_y = CudaArray(y)
        x1 = alpha*(A.'*y)
        x2 = A.'*y
        d_x1 = CUBLAS.gemv('T',alpha,d_A,d_y)
        d_x2 = CUBLAS.gemv('T',d_A,d_y)
        h_x1 = to_host(d_x1)
        h_x2 = to_host(d_x2)
        @test x1 ≈ h_x1
        @test x2 ≈ h_x2
        @test x2 ≈ to_host(d_A.' * d_y)

        # test x = alpha*(A'*y)
        y = rand(elty,m)
        d_y = CudaArray(y)
        x1 = alpha*(A'*y)
        x2 = A'*y
        d_x1 = CUBLAS.gemv('C',alpha,d_A,d_y)
        d_x2 = CUBLAS.gemv('C',d_A,d_y)
        h_x1 = to_host(d_x1)
        h_x2 = to_host(d_x2)
        @test y1 ≈ h_y1
        @test y2 ≈ h_y2
        @test x2 ≈ to_host(d_A' * d_y)
    end
end

@testset "gbmv!" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # parameters
        alpha = convert(elty,2)
        beta = convert(elty,3)
        # bands
        ku = 2
        kl = 3
        # generate banded matrix
        A = rand(elty,m,n)
        A = bandex(A,kl,ku)
        # get packed format
        Ab = band(A,kl,ku)
        d_Ab = CudaArray(Ab)
        # test y = alpha*A*x + beta*y
        x = rand(elty,n)
        d_x = CudaArray(x)
        y = rand(elty,m)
        d_y = CudaArray(y)
        CUBLAS.gbmv!('N',m,kl,ku,alpha,d_Ab,d_x,beta,d_y)
        BLAS.gbmv!('N',m,kl,ku,alpha,Ab,x,beta,y)
        h_y = to_host(d_y)
        @test y ≈ h_y
        # test y = alpha*A.'*x + beta*y
        x = rand(elty,n)
        d_x = CudaArray(x)
        y = rand(elty,m)
        d_y = CudaArray(y)
        CUBLAS.gbmv!('T',m,kl,ku,alpha,d_Ab,d_y,beta,d_x)
        BLAS.gbmv!('T',m,kl,ku,alpha,Ab,y,beta,x)
        h_x = to_host(d_x)
        @test x ≈ h_x
        # test y = alpha*A'*x + beta*y
        x = rand(elty,n)
        d_x = CudaArray(x)
        y = rand(elty,m)
        d_y = CudaArray(y)
        CUBLAS.gbmv!('C',m,kl,ku,alpha,d_Ab,d_y,beta,d_x)
        BLAS.gbmv!('C',m,kl,ku,alpha,Ab,y,beta,x)
        h_x = to_host(d_x)
        @test x ≈ h_x
    end
end

@testset "gbmv" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # parameters
        alpha = convert(elty,2)
        # bands
        ku = 2
        kl = 3
        # generate banded matrix
        A = rand(elty,m,n)
        A = bandex(A,kl,ku)
        # get packed format
        Ab = band(A,kl,ku)
        d_Ab = CudaArray(Ab)
        # test y = alpha*A*x
        x = rand(elty,n)
        d_x = CudaArray(x)
        d_y = CUBLAS.gbmv('N',m,kl,ku,alpha,d_Ab,d_x)
        y = zeros(elty,m)
        y = BLAS.gbmv('N',m,kl,ku,alpha,Ab,x)
        h_y = to_host(d_y)
        @test y ≈ h_y
    end
end

@testset "symv!" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # parameters
        alpha = convert(elty,2)
        beta = convert(elty,3)
        # generate symmetric matrix
        A = rand(elty,m,m)
        A = A + A.'
        # generate vectors
        x = rand(elty,m)
        y = rand(elty,m)
        # copy to device
        d_A = CudaArray(A)
        d_x = CudaArray(x)
        d_y = CudaArray(y)
        # execute on host
        BLAS.symv!('U',alpha,A,x,beta,y)
        # execute on device
        CUBLAS.symv!('U',alpha,d_A,d_x,beta,d_y)
        # compare results
        h_y = to_host(d_y)
        @test y ≈ h_y
    end
end

@testset "symv" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate symmetric matrix
        A = rand(elty,m,m)
        A = A + A.'
        # generate vectors
        x = rand(elty,m)
        # copy to device
        d_A = CudaArray(A)
        d_x = CudaArray(x)
        # execute on host
        y = BLAS.symv('U',A,x)
        # execute on device
        d_y = CUBLAS.symv('U',d_A,d_x)
        # compare results
        h_y = to_host(d_y)
        @test y ≈ h_y
    end
end

@testset "hemv!" begin
    @testset for elty in [Complex64, Complex128]
        # parameters
        alpha = convert(elty,2)
        beta = convert(elty,3)
        # generate hermitian matrix
        A = rand(elty,m,m)
        A = A + A'
        # generate vectors
        x = rand(elty,m)
        y = rand(elty,m)
        # copy to device
        d_A = CudaArray(A)
        d_x = CudaArray(x)
        d_y = CudaArray(y)
        # execute on host
        BLAS.hemv!('U',alpha,A,x,beta,y)
        # execute on device
        CUBLAS.hemv!('U',alpha,d_A,d_x,beta,d_y)
        # compare results
        h_y = to_host(d_y)
        @test y ≈ h_y
    end
end

@testset "hemv" begin
    @testset for elty in [Complex64, Complex128]
        # generate hermitian matrix
        A = rand(elty,m,m)
        A = A + A.'
        # generate vectors
        x = rand(elty,m)
        # copy to device
        d_A = CudaArray(A)
        d_x = CudaArray(x)
        # execute on host
        y = BLAS.hemv('U',A,x)
        # execute on device
        d_y = CUBLAS.hemv('U',d_A,d_x)
        # compare results
        h_y = to_host(d_y)
        @test y ≈ h_y
    end
end

@testset "sbmv!" begin
    @testset for elty in [Float32, Float64]
        # parameters
        alpha = convert(elty,3)
        beta = convert(elty,2.5)
        # generate symmetric matrix
        A = rand(elty,m,m)
        A = A + A'
        # restrict to 3 bands
        nbands = 3
        @test m >= 1+nbands
        A = bandex(A,nbands,nbands)
        # convert to 'upper' banded storage format
        AB = band(A,0,nbands)
        # construct x and y
        x = rand(elty,m)
        y = rand(elty,m)
        # move to host
        d_AB = CudaArray(AB)
        d_x = CudaArray(x)
        d_y = CudaArray(y)
        # sbmv!
        CUBLAS.sbmv!('U',nbands,alpha,d_AB,d_x,beta,d_y)
        y = alpha*(A*x) + beta*y
        # compare
        h_y = to_host(d_y)
        @test y ≈ h_y
    end
end

@testset "sbmv" begin
    @testset for elty in [Float32, Float64]
        # parameters
        alpha = convert(elty,3)
        beta = convert(elty,2.5)
        # generate symmetric matrix
        A = rand(elty,m,m)
        A = A + A'
        # restrict to 3 bands
        nbands = 3
        @test m >= 1+nbands
        A = bandex(A,nbands,nbands)
        # convert to 'upper' banded storage format
        AB = band(A,0,nbands)
        # construct x and y
        x = rand(elty,m)
        y = rand(elty,m)
        # move to host
        d_AB = CudaArray(AB)
        d_x = CudaArray(x)
        # sbmv!
        d_y = CUBLAS.sbmv('U',nbands,d_AB,d_x)
        y = A*x
        # compare
        h_y = to_host(d_y)
        @test y ≈ h_y
    end
end

@testset "hbmv!" begin
    @testset for elty in [Complex64, Complex128]
        # parameters
        alpha = rand(elty)
        beta = rand(elty)
        # generate Hermitian matrix
        A = rand(elty,m,m)
        A = A + ctranspose(A)
        # restrict to 3 bands
        nbands = 3
        @test m >= 1+nbands
        A = bandex(A,nbands,nbands)
        # convert to 'upper' banded storage format
        AB = band(A,0,nbands)
        # construct x and y
        x = rand(elty,m)
        y = rand(elty,m)
        # move to host
        d_AB = CudaArray(AB)
        d_x = CudaArray(x)
        d_y = CudaArray(y)
        # hbmv!
        CUBLAS.hbmv!('U',nbands,alpha,d_AB,d_x,beta,d_y)
        y = alpha*(A*x) + beta*y
        # compare
        h_y = to_host(d_y)
        @test y ≈ h_y
    end
end

@testset "hbmv" begin
    @testset for elty in [Complex64, Complex128]
        # parameters
        alpha = rand(elty)
        beta = rand(elty)
        # generate Hermitian matrix
        A = rand(elty,m,m)
        A = A + ctranspose(A)
        # restrict to 3 bands
        nbands = 3
        @test m >= 1+nbands
        A = bandex(A,nbands,nbands)
        # convert to 'upper' banded storage format
        AB = band(A,0,nbands)
        # construct x and y
        x = rand(elty,m)
        y = rand(elty,m)
        # move to host
        d_AB = CudaArray(AB)
        d_x = CudaArray(x)
        # hbmv
        d_y = CUBLAS.hbmv('U',nbands,d_AB,d_x)
        y = A*x
        # compare
        h_y = to_host(d_y)
        @test y ≈ h_y
    end
end

@testset "tbmv!" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate triangular matrix
        A = rand(elty,m,m)
        # restrict to 3 bands
        nbands = 3
        @test m >= 1+nbands
        A = bandex(A,0,nbands)
        # convert to 'upper' banded storage format
        AB = band(A,0,nbands)
        # construct x and y
        x = rand(elty,m)
        # move to host
        d_AB = CudaArray(AB)
        d_x = CudaArray(x)
        # tbmv!
        CUBLAS.tbmv!('U','N','N',nbands,d_AB,d_x)
        x = A*x
        # compare
        h_x = to_host(d_x)
        @test x ≈ h_x
    end
end

@testset "tbmv" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate triangular matrix
        A = rand(elty,m,m)
        # restrict to 3 bands
        nbands = 3
        @test m >= 1+nbands
        A = bandex(A,0,nbands)
        # convert to 'upper' banded storage format
        AB = band(A,0,nbands)
        # construct x
        x = rand(elty,m)
        # move to host
        d_AB = CudaArray(AB)
        d_x = CudaArray(x)
        # tbmv!
        d_y = CUBLAS.tbmv!('U','N','N',nbands,d_AB,d_x)
        y = A*x
        # compare
        h_y = to_host(d_y)
        @test y ≈ h_y
    end
end

@testset "tbsv!" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate triangular matrix
        A = rand(elty,m,m)
        # restrict to 3 bands
        nbands = 3
        @test m >= 1+nbands
        A = bandex(A,0,nbands)
        # convert to 'upper' banded storage format
        AB = band(A,0,nbands)
        # generate vector
        x = rand(elty,m)
        # move to device
        d_AB = CudaArray(AB)
        d_x = CudaArray(x)
        #tbsv!
        CUBLAS.tbsv!('U','N','N',nbands,d_AB,d_x)
        x = A\x
        # compare
        h_x = to_host(d_x)
        @test x ≈ h_x
    end
end

@testset "tbsv" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate triangular matrix
        A = rand(elty,m,m)
        # restrict to 3 bands
        nbands = 3
        @test m >= 1+nbands
        A = bandex(A,0,nbands)
        # convert to 'upper' banded storage format
        AB = band(A,0,nbands)
        # generate vector
        x = rand(elty,m)
        # move to device
        d_AB = CudaArray(AB)
        d_x = CudaArray(x)
        #tbsv
        d_y = CUBLAS.tbsv('U','N','N',nbands,d_AB,d_x)
        y = A\x
        # compare
        h_y = to_host(d_y)
        @test y ≈ h_y
    end
end

@testset "trmv!" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate triangular matrix
        A = rand(elty,m,m)
        A = triu(A)
        # generate vector
        x = rand(elty,m)
        # move to device
        d_A = CudaArray(A)
        d_x = CudaArray(x)
        # execute trmv!
        CUBLAS.trmv!('U','N','N',d_A,d_x)
        x = A*x
        # compare
        h_x = to_host(d_x)
        @test x ≈ h_x
    end
end

@testset "trmv" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate triangular matrix
        A = rand(elty,m,m)
        A = triu(A)
        # generate vector
        x = rand(elty,m)
        # move to device
        d_A = CudaArray(A)
        d_x = CudaArray(x)
        # execute trmv!
        d_y = CUBLAS.trmv('U','N','N',d_A,d_x)
        y = A*x
        # compare
        h_y = to_host(d_y)
        @test y ≈ h_y
    end
end

@testset "trsv!" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate triangular matrix
        A = rand(elty,m,m)
        A = triu(A)
        # generate vector
        x = rand(elty,m)
        # move to device
        d_A = CudaArray(A)
        d_x = CudaArray(x)
        # execute trsv!
        CUBLAS.trsv!('U','N','N',d_A,d_x)
        x = A\x
        # compare
        h_x = to_host(d_x)
        @test x ≈ h_x
    end
end

@testset "trsv" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate triangular matrix
        A = rand(elty,m,m)
        A = triu(A)
        # generate vector
        x = rand(elty,m)
        # move to device
        d_A = CudaArray(A)
        d_x = CudaArray(x)
        # execute trsv!
        d_y = CUBLAS.trsv('U','N','N',d_A,d_x)
        y = A\x
        # compare
        h_y = to_host(d_y)
        @test y ≈ h_y
    end
end

@testset "ger!" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # construct matrix and vectors
        A = rand(elty,m,n)
        x = rand(elty,m)
        y = rand(elty,n)
        alpha = convert(elty,2)
        # move to device
        d_A = CudaArray(A)
        d_x = CudaArray(x)
        d_y = CudaArray(y)
        # perform rank one update
        CUBLAS.ger!(alpha,d_x,d_y,d_A)
        A = (alpha*x)*y' + A
        # move to host and compare
        h_A = to_host(d_A)
        @test A ≈ h_A
    end
end

@testset "syr!" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # construct matrix and vector
        A = rand(elty,m,m)
        A = A + A.'
        x = rand(elty,m)
        alpha = convert(elty,2)
        # move to device
        d_A = CudaArray(A)
        d_x = CudaArray(x)
        # perform rank one update
        CUBLAS.syr!('U',alpha,d_x,d_A)
        A = (alpha*x)*x.' + A
        # move to host and compare upper triangles
        h_A = to_host(d_A)
        A = triu(A)
        h_A = triu(h_A)
        @test A ≈ h_A
    end
end

@testset "her!" begin
    @testset for elty in [Complex64, Complex128]
        local m = 2
        # construct matrix and vector
        A = rand(elty,m,m)
        A = A + A'
        x = rand(elty,m)
        alpha = convert(elty,2)
        # move to device
        d_A = CudaArray(A)
        d_x = CudaArray(x)
        # perform rank one update
        CUBLAS.her!('U',alpha,d_x,d_A)
        A = (alpha*x)*x' + A
        # move to host and compare upper triangles
        h_A = to_host(d_A)
        A = triu(A)
        h_A = triu(h_A)
        @test A ≈ h_A
    end
end

@testset "her2!" begin
    @testset for elty in [Complex64, Complex128]
        local m = 2
        # construct matrix and vector
        A = rand(elty,m,m)
        A = A + A'
        x = rand(elty,m)
        y = rand(elty,m)
        alpha = convert(elty,2)
        # move to device
        d_A = CudaArray(A)
        d_x = CudaArray(x)
        d_y = CudaArray(y)
        # perform rank one update
        CUBLAS.her2!('U',alpha,d_x,d_y,d_A)
        A = (alpha*x)*y' + y*(alpha*x)' + A
        # move to host and compare upper triangles
        h_A = to_host(d_A)
        A = triu(A)
        h_A = triu(h_A)
        @test A ≈ h_A
    end
end

@testset "gemm!" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # parameters
        alpha = rand(elty)
        beta = rand(elty)
        # generate matrices
        A = rand(elty,m,k)
        B = rand(elty,k,n)
        C1 = rand(elty,m,n)
        C2 = copy(C1)
        # move to device
        d_A = CudaArray(A)
        d_B = CudaArray(B)
        d_C1 = CudaArray(C1)
        d_C2 = CudaArray(C2)
        # C = (alpha*A)*B + beta*C
        CUBLAS.gemm!('N','N',alpha,d_A,d_B,beta,d_C1)
        A_mul_B!(d_C2, d_A, d_B)
        h_C1 = to_host(d_C1)
        h_C2 = to_host(d_C2)
        C1 = (alpha*A)*B + beta*C1
        C2 = A*B
        # compare
        @test C1 ≈ h_C1
        @test C2 ≈ h_C2
    end
end

@testset "gemm" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate matrices
        A = rand(elty,m,k)
        B = rand(elty,k,n)
        # move to device
        d_A = CudaArray(A)
        d_B = CudaArray(B)
        # C = (alpha*A)*B + beta*C
        d_C = CUBLAS.gemm('N','N',d_A,d_B)
        C = A*B
        C2 = d_A * d_B
        # compare
        h_C = to_host(d_C)
        h_C2 = to_host(C2)
        @test C ≈ h_C
        @test C ≈ h_C2
    end
end

@testset "gemm_batched!" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # parameters
        alpha = rand(elty)
        beta = rand(elty)
        # generate matrices
        A = [rand(elty,m,k) for i in 1:10]
        B = [rand(elty,k,n) for i in 1:10]
        C = [rand(elty,m,n) for i in 1:10]
        # move to device
        d_A = CudaArray{elty, 2}[]
        d_B = CudaArray{elty, 2}[]
        d_C = CudaArray{elty, 2}[]
        for i in 1:length(A)
            push!(d_A,CudaArray(A[i]))
            push!(d_B,CudaArray(B[i]))
            push!(d_C,CudaArray(C[i]))
        end
        # C = (alpha*A)*B + beta*C
        CUBLAS.gemm_batched!('N','N',alpha,d_A,d_B,beta,d_C)
        for i in 1:length(d_C)
            C[i] = (alpha*A[i])*B[i] + beta*C[i]
            h_C = to_host(d_C[i])
            #compare
            @test C[i] ≈ h_C
        end
    end
end

@testset "gemm_batched" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate matrices
        A = [rand(elty,m,k) for i in 1:10]
        B = [rand(elty,k,n) for i in 1:10]
        # move to device
        d_A = CudaArray{elty, 2}[]
        d_B = CudaArray{elty, 2}[]
        for i in 1:length(A)
            push!(d_A, CudaArray(A[i]))
            push!(d_B, CudaArray(B[i]))
        end
        # C = A*B
        d_C = CUBLAS.gemm_batched('N','N',d_A,d_B)
        for i in 1:length(A)
            C = A[i]*B[i]
            h_C = to_host(d_C[i])
            @test C ≈ h_C
        end
    end
end

@testset "symm!" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # parameters
        alpha = rand(elty)
        beta = rand(elty)
        # generate matrices
        A = rand(elty,m,m)
        A = A + A.'
        B = rand(elty,m,n)
        C = rand(elty,m,n)
        # move to device
        d_A = CudaArray(A)
        d_B = CudaArray(B)
        d_C = CudaArray(C)
        # C = (alpha*A)*B + beta*C
        CUBLAS.symm!('L','U',alpha,d_A,d_B,beta,d_C)
        C = (alpha*A)*B + beta*C
        # compare
        h_C = to_host(d_C)
        @test C ≈ h_C
    end
end

@testset "symm" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate matrices
        A = rand(elty,m,m)
        A = A + A.'
        B = rand(elty,m,n)
        # move to device
        d_A = CudaArray(A)
        d_B = CudaArray(B)
        # C = (alpha*A)*B + beta*C
        d_C = CUBLAS.symm('L','U',d_A,d_B)
        C = A*B
        # compare
        h_C = to_host(d_C)
        @test C ≈ h_C
    end
end

@testset "syrk!" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate matrices
        A = rand(elty,m,k)
        C = rand(elty,m,m)
        C = C + C.'
        # parameters
        alpha = rand(elty)
        beta = rand(elty)
        # move to device
        d_A = CudaArray(A)
        d_C = CudaArray(C)
        # C = (alpha*A)*A.' + beta*C
        CUBLAS.syrk!('U','N',alpha,d_A,beta,d_C)
        C = (alpha*A)*A.' + beta*C
        C = triu(C)
        # move to host and compare
        h_C = to_host(d_C)
        h_C = triu(C)
        @test C ≈ h_C
    end
end

@testset "syrk" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate matrices
        A = rand(elty,m,k)
        # move to device
        d_A = CudaArray(A)
        # C = A*A.'
        d_C = CUBLAS.syrk('U','N',d_A)
        C = A*A.'
        C = triu(C)
        # move to host and compare
        h_C = to_host(d_C)
        h_C = triu(C)
        @test C ≈ h_C
    end
end

@testset "herk!" begin
    @testset for elty in [Complex64, Complex128]
        # generate matrices
        A = rand(elty,m,k)
        C = rand(elty,m,m)
        C = C + C'
        # parameters
        alpha = rand(elty)
        beta = rand(elty)
        # move to device
        d_A = CudaArray(A)
        d_C = CudaArray(C)
        CUBLAS.herk!('U','N',alpha,d_A,beta,d_C)
        C = alpha*(A*A') + beta*C
        C = triu(C)
        # move to host and compare
        h_C = to_host(d_C)
        h_C = triu(C)
        @test C ≈ h_C
    end
end

@testset "herk" begin
    @testset for elty in [Complex64, Complex128]
        # generate matrices
        A = rand(elty,m,k)
        # move to device
        d_A = CudaArray(A)
        # C = A*A'
        d_C = CUBLAS.herk('U','N',d_A)
        C = A*A'
        C = triu(C)
        # move to host and compare
        h_C = to_host(d_C)
        h_C = triu(C)
        @test C ≈ h_C
    end
end

@testset "syr2k!" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        #local m = 3
        #local k = 1
        # generate parameters
        alpha = rand(elty)
        beta = rand(elty)
        # generate matrices
        A = rand(elty,m,k)
        B = rand(elty,m,k)
        C = rand(elty,m,m)
        C = C + C.'
        # move to device
        d_A = CudaArray(A)
        d_B = CudaArray(B)
        d_C = CudaArray(C)
        # compute
        #C = alpha*(A*B.') + conj(alpha)*(B*A.') + beta*C
        C = alpha*(A*B.' + B*A.') + beta*C
        CUBLAS.syr2k!('U','N',alpha,d_A,d_B,beta,d_C)
        # move back to host and compare
        C = triu(C)
        h_C = to_host(d_C)
        h_C = triu(h_C)

        @test C ≈ h_C
    end
end

@testset "syr2k" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate parameters
        alpha = rand(elty)
        # generate matrices
        A = rand(elty,m,k)
        B = rand(elty,m,k)
        # move to device
        d_A = CudaArray(A)
        d_B = CudaArray(B)
        # compute
        #C = alpha*(A*B.') + conj(alpha)*(B*A.') + beta*C
        C = alpha*(A*B.' + B*A.')
        d_C = CUBLAS.syr2k('U','N',alpha,d_A,d_B)
        # move back to host and compare
        C = triu(C)
        h_C = to_host(d_C)
        h_C = triu(h_C)
        @test C ≈ h_C
    end
end

@testset "her2k!" begin
    @testset for (elty1, elty2) in [(Complex64, Float32), (Complex128, Float64)]
        # generate parameters
        alpha = rand(elty1)
        beta = rand(elty2)
        # generate matrices
        A = rand(elty1,m,k)
        B = rand(elty1,m,k)
        C = rand(elty1,m,m)
        C = C + C'
        # move to device
        d_A = CudaArray(A)
        d_B = CudaArray(B)
        d_C = CudaArray(C)
        # compute
        #C = alpha*(A*B') + conj(alpha)*(B*A') + beta*C
        C = alpha*(A*B') + conj(alpha)*(B*A') + beta*C
        CUBLAS.her2k!('U','N',alpha,d_A,d_B,beta,d_C)
        # move back to host and compare
        C = triu(C)
        h_C = to_host(d_C)
        h_C = triu(h_C)
        @test C ≈ h_C
    end
end

@testset "her2k" begin
    @testset for elty in [Complex64, Complex128]
        # generate matrices
        A = rand(elty,m,k)
        B = rand(elty,m,k)
        # move to device
        d_A = CudaArray(A)
        d_B = CudaArray(B)
        # compute
        C = A*B' + B*A'
        d_C = CUBLAS.her2k('U','N',d_A,d_B)
        # move back to host and compare
        C = triu(C)
        h_C = to_host(d_C)
        h_C = triu(h_C)
        @test C ≈ h_C
    end
end

@testset "trmm!" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate parameter
        alpha = rand(elty)
        # generate matrices
        A = rand(elty,m,m)
        A = triu(A)
        B = rand(elty,m,n)
        C = zeros(elty,m,n)
        # move to device
        d_A = CudaArray(A)
        d_B = CudaArray(B)
        d_C = CudaArray(C)
        # compute
        C = alpha*A*B
        CUBLAS.trmm!('L','U','N','N',alpha,d_A,d_B,d_C)
        # move to host and compare
        h_C = to_host(d_C)
        @test C ≈ h_C
    end
end

@testset "trmm" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate parameter
        alpha = rand(elty)
        # generate matrices
        A = rand(elty,m,m)
        A = triu(A)
        B = rand(elty,m,n)
        # move to device
        d_A = CudaArray(A)
        d_B = CudaArray(B)
        # compute
        C = alpha*A*B
        d_C = CUBLAS.trmm('L','U','N','N',alpha,d_A,d_B)
        # move to host and compare
        h_C = to_host(d_C)
        @test C ≈ h_C
    end
end

@testset "trsm!" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate parameter
        alpha = rand(elty)
        # generate matrices
        A = rand(elty,m,m)
        A = triu(A)
        B = rand(elty,m,n)
        # move to device
        d_A = CudaArray(A)
        d_B = CudaArray(B)
        # compute
        C = alpha*(A\B)
        CUBLAS.trsm!('L','U','N','N',alpha,d_A,d_B)
        # move to host and compare
        h_C = to_host(d_B)
        @test C ≈ h_C
    end
end

@testset "trsm" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate parameter
        alpha = rand(elty)
        # generate matrices
        A = rand(elty,m,m)
        A = triu(A)
        B = rand(elty,m,n)
        # move to device
        d_A = CudaArray(A)
        d_B = CudaArray(B)
        # compute
        C = alpha*(A\B)
        d_C = CUBLAS.trsm('L','U','N','N',alpha,d_A,d_B)
        # move to host and compare
        h_C = to_host(d_C)
        @test C ≈ h_C
    end
end

@testset "trsm_batched!" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate parameter
        alpha = rand(elty)
        # generate matrices
        A = [rand(elty,m,m) for i in 1:10]
        map!((x) -> triu(x), A, A)
        B = [rand(elty,m,n) for i in 1:10]
        # move to device
        d_A = CudaArray{elty, 2}[]
        d_B = CudaArray{elty, 2}[]
        for i in 1:length(A)
            push!(d_A,CudaArray(A[i]))
            push!(d_B,CudaArray(B[i]))
        end
        # compute
        CUBLAS.trsm_batched!('L','U','N','N',alpha,d_A,d_B)
        # move to host and compare
        for i in 1:length(d_B)
            C = alpha*(A[i]\B[i])
            h_C = to_host(d_B[i])
            #compare
            @test C ≈ h_C
        end
    end
end

@testset "trsm_batched" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate parameter
        alpha = rand(elty)
        # generate matrices
        A = [rand(elty,m,m) for i in 1:10]
        map!((x) -> triu(x), A, A)
        B = [rand(elty,m,n) for i in 1:10]
        # move to device
        d_A = CudaArray{elty, 2}[]
        d_B = CudaArray{elty, 2}[]
        for i in 1:length(A)
            push!(d_A,CudaArray(A[i]))
            push!(d_B,CudaArray(B[i]))
        end
        # compute
        d_C = CUBLAS.trsm_batched('L','U','N','N',alpha,d_A,d_B)
        # move to host and compare
        for i in 1:length(d_C)
            C = alpha*(A[i]\B[i])
            h_C = to_host(d_C[i])
            @test C ≈ h_C
        end
    end
end

@testset "hemm!" begin
    @testset for elty in [Complex64, Complex128]
        # generate parameters
        alpha = rand(elty)
        beta  = rand(elty)
        # generate matrices
        A = rand(elty,m,m)
        A = A + ctranspose(A)
        @test ishermitian(A)
        B = rand(elty,m,n)
        C = rand(elty,m,n)
        # move to device
        d_A = CudaArray(A)
        d_B = CudaArray(B)
        d_C = CudaArray(C)
        # compute
        C = alpha*(A*B) + beta*C
        CUBLAS.hemm!('L','L',alpha,d_A,d_B,beta,d_C)
        # move to host and compare
        h_C = to_host(d_C)
        @test C ≈ h_C
    end
end

@testset "hemm" begin
    @testset for elty in [Complex64, Complex128]
        # generate parameter
        alpha = rand(elty)
        # generate matrices
        A = rand(elty,m,m)
        A = A + ctranspose(A)
        @test ishermitian(A)
        B = rand(elty,m,n)
        # move to device
        d_A = CudaArray(A)
        d_B = CudaArray(B)
        # compute
        C = alpha*(A*B)
        d_C = CUBLAS.hemm('L','U',alpha,d_A,d_B)
        # move to host and compare
        h_C = to_host(d_C)
        @test C ≈ h_C
    end
end

@testset "geam!" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate parameters
        alpha = rand(elty)
        beta  = rand(elty)
        # generate matrices
        A = rand(elty,m,n)
        B = rand(elty,m,n)
        C = zeros(elty,m,n)
        # move to device
        d_A = CudaArray(A)
        d_B = CudaArray(B)
        d_C = CudaArray(C)
        # compute
        C = alpha*A + beta*B
        CUBLAS.geam!('N','N',alpha,d_A,beta,d_B,d_C)
        # move to host and compare
        h_C = to_host(d_C)
        @test C ≈ h_C

        #test in place versions too
        C = rand(elty,m,n)
        d_C = CudaArray(C)
        C = alpha*C + beta*B
        CUBLAS.geam!('N','N',alpha,d_C,beta,d_B,d_C)
        # move to host and compare
        h_C = to_host(d_C)
        @test C ≈ h_C
        C = rand(elty,m,n)
        d_C = CudaArray(C)
        C = alpha*A + beta*C
        CUBLAS.geam!('N','N',alpha,d_A,beta,d_C,d_C)
        # move to host and compare
        h_C = to_host(d_C)
        @test C ≈ h_C

        #test setting C to zero
        C = rand(elty,m,n)
        d_C = CudaArray(C)
        alpha = zero(elty)
        beta  = zero(elty)
        CUBLAS.geam!('N','N',alpha,d_A,beta,d_B,d_C)
        h_C = to_host(d_C)
        @test h_C ≈ zeros(elty,m,n)

        # bounds checking
        @test_throws DimensionMismatch CUBLAS.geam!('N','T',alpha,d_A,beta,d_B,d_C)
        @test_throws DimensionMismatch CUBLAS.geam!('T','T',alpha,d_A,beta,d_B,d_C)
        @test_throws DimensionMismatch CUBLAS.geam!('T','N',alpha,d_A,beta,d_B,d_C)
    end
end

@testset "geam" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate parameter
        alpha = rand(elty)
        beta  = rand(elty)
        # generate matrices
        A = rand(elty,m,n)
        B = rand(elty,m,n)
        # move to device
        d_A = CudaArray(A)
        d_B = CudaArray(B)
        C = zeros(elty,m,n)
        # compute
        C = alpha*A + beta*B
        d_C = CUBLAS.geam('N','N',alpha,d_A,beta,d_B)
        # move to host and compare
        h_C = to_host(d_C)
        @test C ≈ h_C
    end
end

@testset "getrf_batched!" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate matrices
        A = [rand(elty,m,m) for i in 1:10]
        # move to device
        d_A = CudaArray{elty, 2}[]
        for i in 1:length(A)
            push!(d_A,CudaArray(A[i]))
        end
        pivot, info = CUBLAS.getrf_batched!(d_A, false)
        h_info = to_host(info)
        for As in 1:length(d_A)
            C   = lufact!(copy(A[As]), Val{false}) # lufact(A[As],pivot=false)
            h_A = to_host(d_A[As])
            #reconstruct L,U
            dL = eye(elty,m)
            dU = zeros(elty,(m,m))
            k = h_info[As]
            if( k >= 0 )
                dL += tril(h_A,-k-1)
                dU += triu(h_A,k)
            end
            #compare
            @test isapprox(C[:L], dL, rtol=1e-2)
            @test isapprox(C[:U], dU, rtol=1e-2)
        end
        for i in 1:length(A)
            d_A[ i ] = CudaArray(A[i])
        end
        pivot, info = CUBLAS.getrf_batched!(d_A, true)
        h_info = to_host(info)
        h_pivot = to_host(pivot)
        for As in 1:length(d_A)
            C   = lufact(A[As])
            h_A = to_host(d_A[As])
            #reconstruct L,U
            dL = eye(elty,m)
            dU = zeros(elty,(m,m))
            k = h_info[As]
            if( k >= 0 )
                dL += tril(h_A,-k-1)
                dU += triu(h_A,k)
            end
            #compare pivots
            @test length(setdiff(h_pivot[:,As],C[:p])) == 0
            #make device pivot matrix
            P = eye(m)
            for row in 1:m
                temp = copy(P[row,:])
                P[row,:] = P[h_pivot[row,As],:]
                P[h_pivot[row,As],:] = temp
            end
            @test inv(P)*dL*dU ≈ inv(C[:P]) * C[:L] * C[:U]
        end
    end
end

@testset "getrf_batched" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate matrices
        A = [rand(elty,m,m) for i in 1:10]
        # move to device
        d_A = CudaArray{elty, 2}[]
        for i in 1:length(A)
            push!(d_A,CudaArray(A[i]))
        end
        pivot, info, d_B = CUBLAS.getrf_batched(d_A, false)
        h_info = to_host(info)
        for Bs in 1:length(d_B)
            C   = lufact!(copy(A[Bs]),Val{false}) # lufact(A[Bs],pivot=false)
            h_B = to_host(d_B[Bs])
            #reconstruct L,U
            dL = eye(elty,m)
            dU = zeros(elty,(m,m))
            k = h_info[Bs]
            if( h_info[Bs] >= 0 )
                dU += triu(h_B,k)
                dL += tril(h_B,-k-1)
            end
            #compare
            @test isapprox(C[:L], dL, rtol=1e-2)
            @test isapprox(C[:U], dU, rtol=1e-2)
        end
    end
end

@testset "getri_batched" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate matrices
        A = [rand(elty,m,m) for i in 1:10]
        # move to device
        d_A = CudaArray{elty, 2}[]
        for i in 1:length(A)
            push!(d_A,CudaArray(A[i]))
        end
        pivot, info = CUBLAS.getrf_batched!(d_A, true)
        h_info = to_host(info)
        for Cs in 1:length(h_info)
            @test h_info[Cs] == 0
        end
        pivot, info, d_C = CUBLAS.getri_batched(d_A, pivot)
        h_info = to_host(info)
        for Cs in 1:length(d_C)
            C   = inv(A[Cs])
            h_C = to_host(d_C[Cs])
            @test h_info[Cs] == 0
            @test C ≈ h_C
        end
    end
end

@testset "matinv_batched" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate matrices
        A = [rand(elty,m,m) for i in 1:10]
        # move to device
        d_A = CudaArray{elty, 2}[]
        for i in 1:length(A)
            push!(d_A,CudaArray(A[i]))
        end
        info, d_C = CUBLAS.matinv_batched(d_A)
        for Cs in 1:length(d_C)
            C   = inv(A[Cs])
            h_C = to_host(d_C[Cs])
            @test C ≈ h_C
        end
    end
end

@testset "geqrf_batched!" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate matrices
        A = [rand(elty,m,n) for i in 1:10]
        # move to device
        d_A = CudaArray{elty, 2}[]
        for i in 1:length(A)
            push!(d_A,CudaArray(A[i]))
        end
        tau, d_A = CUBLAS.geqrf_batched!(d_A)
        for As in 1:length(d_A)
            C   = qrfact(A[As])
            h_A = to_host(d_A[As])
            h_tau = to_host(tau[As])
            # build up Q
            Q = eye(elty,min(m,n))
            for i in 1:min(m,n)
                v = zeros(elty,m)
                v[i] = one(elty)
                v[i+1:m] = h_A[i+1:m,i]
                Q *= eye(elty,m) - h_tau[i] * v * v'
            end
            @test Q≈full(C[:Q])
        end
    end
end

@testset "geqrf_batched" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate matrices
        A = [rand(elty,m,n) for i in 1:10]
        # move to device
        d_A = CudaArray{elty, 2}[]
        for i in 1:length(A)
            push!(d_A,CudaArray(A[i]))
        end
        tau, d_B = CUBLAS.geqrf_batched!(d_A)
        for Bs in 1:length(d_B)
            C   = qrfact(A[Bs])
            h_B = to_host(d_B[Bs])
            h_tau = to_host(tau[Bs])
            # build up Q
            Q = eye(elty,min(m,n))
            for i in 1:min(m,n)
                v = zeros(elty,m)
                v[i] = one(elty)
                v[i+1:m] = h_B[i+1:m,i]
                Q *= eye(elty,m) - h_tau[i] * v * v'
            end
            @test Q≈full(C[:Q])
        end
    end
end

@testset "gels_batched!" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate matrices
        A = [rand(elty,n,n) for i in 1:10]
        C = [rand(elty,n,k) for i in 1:10]
        # move to device
        d_A = CudaArray{elty, 2}[]
        d_C = CudaArray{elty, 2}[]
        for i in 1:length(A)
            push!(d_A,CudaArray(A[i]))
            push!(d_C,CudaArray(C[i]))
        end
        d_A, d_C, info = CUBLAS.gels_batched!('N',d_A, d_C)
        for Cs in 1:length(d_C)
            X = A[Cs]\C[Cs]
            h_C = to_host(d_C[Cs])
            @test X≈h_C
        end
    end
end

@testset "gels_batched" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate matrices
        A = [rand(elty,n,n) for i in 1:10]
        C = [rand(elty,n,k) for i in 1:10]
        # move to device
        d_A = CudaArray{elty, 2}[]
        d_C = CudaArray{elty, 2}[]
        for i in 1:length(A)
            push!(d_A,CudaArray(A[i]))
            push!(d_C,CudaArray(C[i]))
        end
        d_B, d_D, info = CUBLAS.gels_batched('N',d_A, d_C)
        for Ds in 1:length(d_D)
            X = A[Ds]\C[Ds]
            h_D = to_host(d_D[Ds])
            @test X ≈ h_D
        end
    end
end

@testset "dgmm!" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate matrices
        A = rand(elty,m,n)
        C = rand(elty,m,n)
        X = rand(elty,m)
        # move to device
        d_A = CudaArray(A)
        d_C = CudaArray(C)
        d_X = CudaArray(X)
        # compute
        C = diagm(X) * A
        CUBLAS.dgmm!('L',d_A,d_X,d_C)
        # move to host and compare
        h_C = to_host(d_C)
        @test C ≈ h_C
        # bounds checking
        @test_throws DimensionMismatch CUBLAS.dgmm!('R',d_A,d_X,d_C)
        A = rand(elty,m,m)
        d_A = CudaArray(A)
        @test_throws DimensionMismatch CUBLAS.dgmm!('L',d_A,d_X,d_C)
    end
end

@testset "dgmm" begin
    @testset for elty in [Float32, Float64, Complex64, Complex128]
        # generate matrices
        A = rand(elty,m,n)
        X = rand(elty,m)
        # move to device
        d_A = CudaArray(A)
        d_X = CudaArray(X)
        # compute
        C = diagm(X) * A
        d_C = CUBLAS.dgmm('L',d_A,d_X)
        # move to host and compare
        h_C = to_host(d_C)
        @test C ≈ h_C
    end
end


text:
# ## Data

import CSV
import Random
import TuringModels

using DataFrames
using Statistics: mean

Random.seed!(1)

data_path = joinpath(TuringModels.project_root, "data", "Kline.csv")
df = CSV.read(data_path, DataFrame; delim=';')

df.log_pop = log.(df.population)
df.contact_high = [contact == "high" ? 1 : 0 for contact in df.contact]

# New col where we center(!) the log_pop values
mean_log_pop = mean(df.log_pop)
df.log_pop_c = map(x -> x - mean_log_pop, df.log_pop)
df

# ## Model

using Turing

@model m10_10stan_c(total_tools, log_pop_c, contact_high) = begin
    α ~ Normal(0, 100)
    βp ~ Normal(0, 1)
    βc ~ Normal(0, 1)
    βpc ~ Normal(0, 1)

    for i ∈ 1:length(total_tools)
        λ = exp(α + βp*log_pop_c[i] + βc*contact_high[i] +
            βpc*contact_high[i]*log_pop_c[i])
        total_tools[i] ~ Poisson(λ)
    end
end;

chns = sample(
    m10_10stan_c(df.total_tools, df.log_pop_c, df.contact_high),
    NUTS(),
    1000
)

# \defaultoutput{}

# ## Original output

m_10_10t_c_result = "
    mean   sd  5.5% 94.5% n_eff Rhat
 a   3.31 0.09  3.17  3.45  3671    1
 bp  0.26 0.03  0.21  0.32  5052    1
 bc  0.28 0.12  0.10  0.47  3383    1
 bcp 0.07 0.17 -0.20  0.34  4683    1
";


text:
### Initializing state and state transitions: INPUT is C, cars_per_zone, T and OUTPUT state_matrix, transition_matrix

function initializestates(C)

    state_matrix = zeros(C, T); # C X T  matrix for saving states of the system, i.e. the location of cars in each timestep. The rows stand for car ID and the columns for the timestep.
    transition_matrix = zeros(C, T, 4); # C X T X 4 matrix for saving state transitions of the system in each timestep. The entries for each car ID (C) and each time step (T) are the driving activity, destination, travel time and travel distance.
    zone = 0;

    ## assigning a zone uniformely to each car
    for i = 0:cars_per_zone:(C-cars_per_zone)
        zone = zone + 1;
        for j=1:cars_per_zone
            state_matrix[i+j, 1] = zone; # each car receives a location zone in time step 1, i.e. at 1h o'clock
        end
    end
    state_matrix = round.(Int, state_matrix);

    return state_matrix, transition_matrix

end

max_stars_repo_name: ArsamAryandoust/CarParkingMaps

max_issues_repo_name: ArsamAryandoust/CarParkingMaps

max_forks_repo_name: ArsamAryandoust/CarParkingMaps


text:
using Markdown
using InteractiveUtils
using Pkg
#Pkg.add(url = "https://github.com/hustf/Unitfu.jl")
#Pkg.add(url = "https://github.com/hustf/MechanicalUnits.jl")
#Pkg.add(url = "https://github.com/hustf/MechanicalSketch.jl")
#Pkg.develop("MechanicalSketch")
#Pkg.add("Latexify")
#Pkg.add("PlutoUI")
#Pkg.status()

Pkg.instantiate()
push!(LOAD_PATH, Base.find_package("MechanicalSketch"))
#VERSION >= v"1.5" && using Revise
import MechanicalSketch:
    text, circle, Turtle, Pencolor, Penwidth, Forward, Turn,
    HueShift, O, sethue, finish, EM, WI, m, background,
    empty_figure, preview, @svg, setline, Point, placeimage, @draw,
    @imagematrix
import Latexify: @latexify
using PlutoUI


text:
using Pkg.TOML, StaticArrays

@testset "MagneticFields.jl" begin
    t = randn()
    r = rand(SVector{3,Float64})
    G = SMatrix{3,3,Float64}([1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 -2.0])

    params = TOML.parsefile("MagneticFields.toml")

	staticHomogeneousField = MPISimulations.initialize(MPISimulations.StaticHomogeneousField,"StaticHomogeneousField",params,Float64)
    @test MPISimulations.magneticFieldStrength(staticHomogeneousField,r) ≈ SVector{3,Float64}(0.12,0.12,0.12)
    @test MPISimulations.magneticFieldStrength(staticHomogeneousField,r,t) ≈ SVector{3,Float64}(0.12,0.12,0.12)

    staticGradientField = MPISimulations.initialize(MPISimulations.StaticGradientField,"StaticGradientField",params,Float64)
    @test MPISimulations.magneticFieldStrength(staticGradientField,r) ≈ G*r
    @test MPISimulations.magneticFieldStrength(staticGradientField,r,t) ≈ G*r

    sinusoidalField1 = MPISimulations.initialize(MPISimulations.SinusoidalField,"SinusoidalField",params,Float64)
    @test MPISimulations.magneticFieldStrength(sinusoidalField1,r,t) ≈ sin(2*pi*25e3*t+π)*SVector{3,Float64}(0.12,0.12,0.12)

    spatialFieldProfile = MPISimulations.initialize(MPISimulations.StaticHomogeneousField,"SinusoidalField.spatialFieldProfile.StaticHomogeneousField",params,Float64)
    sinusoidalField2 = MPISimulations.SinusoidalField(spatialFieldProfile,2.5e6,100,1.0)
    @test sinusoidalField1.baseFrequency == sinusoidalField2.baseFrequency
    @test sinusoidalField1.divider == sinusoidalField2.divider
    @test sinusoidalField1.doubleFrequency == sinusoidalField2.doubleFrequency
    @test sinusoidalField1.relativePhase == sinusoidalField2.relativePhase
    @test sinusoidalField1.spatialFieldProfile.magneticFieldStrength == sinusoidalField2.spatialFieldProfile.magneticFieldStrength
end


text:
@testset "Linear rescaling of multi-dimensional sample" begin

    # Diagonal rescaling
    Nx = 1
    Ne = 500

    X = randn(Nx).+ randn(Nx, Ne) .* randn(Nx, Ne)
    X0 = deepcopy(X)

    # Identity transformation
    L = LinearTransform(Nx)
    @test norm(L.μ)<1e-6
    @test norm(L.L - Matrix(1.0*I, Nx, Nx))<1e-6

    TransportBasedInference.transform!(L, X)
    @test norm(X - X0)<1e-6
    TransportBasedInference.itransform!(L, X)
    @test norm(X - X0)<1e-6

    L = LinearTransform(X; diag = true)

    TransportBasedInference.transform!(L, X);

    @test norm(mean(X;dims = 2))<1e-10
    @test norm(std(X; dims = 2)[:,1] - ones(Nx))<1e-10

    TransportBasedInference.itransform!(L, X)

    @test norm(mean(X;dims = 2) - mean(X0;dims = 2))<1e-10
    @test norm(cov(X')  - cov(X0'))<1e-10

    # Diagonal rescaling
    Nx = 100
    Ne = 500
    X = zeros(Nx, Ne)

    X .= randn(Nx).+ randn(Nx, Ne) .* randn(Nx, Ne)
    X0 = deepcopy(X)

    # Identity transformation
    L = LinearTransform(Nx)
    @test norm(L.μ)<1e-6
    @test norm(L.L - Matrix(1.0*I, Nx, Nx))<1e-6

    TransportBasedInference.transform!(L, X)
    @test norm(X - X0)<1e-6
    TransportBasedInference.itransform!(L, X)
    @test norm(X - X0)<1e-6

    L = LinearTransform(X; diag = true)

    TransportBasedInference.transform!(L, X);

    @test norm(mean(X;dims = 2))<1e-10
    @test norm(std(X; dims = 2)[:,1] - ones(Nx))<1e-10

    TransportBasedInference.itransform!(L, X)

    @test norm(mean(X; dims = 2) - mean(X0; dims = 2))<1e-10
    @test norm(cov(X')  - cov(X0'))<1e-10

    # Dense rescaling
    Nx = 100
    Ne = 500
    X = zeros(Nx, Ne)

    X .= randn(Nx).+ randn(Nx, Ne) .* randn(Nx, Ne)
    X0 = deepcopy(X)
    L = LinearTransform(X; diag = false)

    TransportBasedInference.transform!(L, X);

    @test norm(mean(X;dims = 2))<1e-10
    @test norm(cov(X')  - I)<1e-10

    TransportBasedInference.itransform!(L, X)

    @test norm(mean(X;dims = 2) - mean(X0; dims = 2))<1e-10
    @test norm(cov(X')  - cov(X0'))<1e-10
end

@testset "MinMaxTransform in-place rescaling of multi-dimensional sample" begin
    # Nx = 1
    Nx = 1
    Ne = 500
    X = randn(Nx) .+ randn(Nx, Ne) .* randn(Nx, Ne)
    X0 = deepcopy(X)

    L = MinMaxTransform(X)

    transform!(L, X)

    for i=1:Nx
        @test abs(minimum(X[i,:]) - (-1.0))<1e-10
        @test abs(maximum(X[i,:]) - ( 1.0))<1e-10
    end

    itransform!(L, X)

    @test norm(X - X0)<1e-10

    # Nx = 5
    Nx = 5
    Ne = 500
    X = randn(Nx) .+ randn(Nx, Ne) .* randn(Nx, Ne)
    X0 = deepcopy(X)

    L = MinMaxTransform(X)

    transform!(L, X)

    for i=1:Nx
        @test abs(minimum(X[i,:]) - (-1.0))<1e-10
        @test abs(maximum(X[i,:]) - ( 1.0))<1e-10
    end

    itransform!(L, X)

    @test norm(X - X0)<1e-10
end

@testset "MinMaxTransform out of place rescaling of multi-dimensional sample" begin
    # Nx = 1
    Nx = 1
    Ne = 500
    Xin = randn(Nx) .+ randn(Nx, Ne) .* randn(Nx, Ne)
    Xout = zero(Xin)
    Xinverse = zero(Xin)
    X0 = deepcopy(Xin)

    L = MinMaxTransform(Xin)

    transform!(L, Xout, Xin)

    for i=1:Nx
        @test abs(minimum(Xout[i,:]) - (-1.0))<1e-10
        @test abs(maximum(Xout[i,:]) - ( 1.0))<1e-10
    end

    itransform!(L, Xinverse, Xout)

    @test norm(Xinverse - X0)<1e-10
    @test norm(Xin - X0)<1e-10

    # Nx = 5
    Nx = 5
    Ne = 500
    Xin = randn(Nx) .+ randn(Nx, Ne) .* randn(Nx, Ne)
    Xout = zero(Xin)
    Xinverse = zero(Xin)
    X0 = deepcopy(Xin)

    L = MinMaxTransform(Xin)

    transform!(L, Xout, Xin)

    for i=1:Nx
        @test abs(minimum(Xout[i,:]) - (-1.0))<1e-10
        @test abs(maximum(Xout[i,:]) - ( 1.0))<1e-10
    end

    itransform!(L, Xinverse, Xout)

    @test norm(Xinverse - X0)<1e-10
    @test norm(Xin - X0)<1e-10
end


text:
using Pkg
pkg"activate ."
pkg"add JuMP"
pkg"add Cbc"
pkg"add Juniper"
pkg"add Ipopt"

using JuMP, Cbc, Juniper
# Escolhendo o modelo
model = Model(Cbc.Optimizer)
# Definindo os parâmetros que, por sinal, coloquei aleatórios, pois não sei as entradas
m = 4
n = 4
c = rand(1:4, n,m)
a = rand(0:0, m)
b = rand(1000:1000,n)
c[pi[j],pi[i]] * x[pi[i],pi[g]]
# Implementando realmente o modelo
@variable(model,x[i = 1:n, g = 1:m], Bin)
@objective(model, Max, t)
@constraint(model, [i = 1:n], sum(x[i,g] for  g = 1:m) == 1) # 1º restrição
@constraint(model, [i = 1:n], sum(x[i,g] for i = 1:n-1, g = 1:m ) ≤ b[g]) # 2º restrição
@constraint(model, [i = 1:n], sum(x[i,g] for i = 1:n-1, g = 1:m ) ≥ a[g]) # 3º restrição
#@constraint(model, [i = 1:n - 1, g = 1:m], - sum(max(0,(c[i,j])) * x[i,g] + w[i,g] for j = i+1:n) <= 0) # 4º restrição
#@constraint(model, [i = 1:n - 1, g = 1:m], - sum(c[i,j]*x[j,g] +  sum(min(0,c[i,j]) * (1 - x[i,g]) + w[i,g]) for j = i+1:n)  <= 0) # 5º restrição

optimize!(model) # Otimizando o modelo

println("Valor da função objetivo = ", objective_value(model))
value.(x)


text:
using SpatialFields
using Base.Test
using Iterators

points = SVector{3, Float64}[[1; 0; 0], [0; 1; 0], [-1; 0; 0], [0; -1; 0], [0; 0; 1], [0; 0; -1]]
normals = SVector{3, Float64}[[1; -1; 0], [0; 1; 0], [-1; 0; 0], [0; -1; 0], [1; 0; 1], [0; 0; -1]]

field = HermiteRadialField(points, normals)

X = linspace(-2, 2)
Y = linspace(-2, 2)
Z = linspace(-2, 2)

for (i, point) in enumerate(points)
    @test isapprox(field(point), 0, atol=1e-6)
    @test isapprox(SpatialFields.gradient(field)(point), normals[i], atol=1e-6)
end


text:
function FieldODEStateAndDerivative(arg0::RealFieldElement, arg1::Vector{RealFieldElement}, arg2::Vector{RealFieldElement})
    return FieldODEStateAndDerivative((RealFieldElement, Vector{RealFieldElement}, Vector{RealFieldElement}), arg0, arg1, arg2)
end

function FieldODEStateAndDerivative(arg0::RealFieldElement, arg1::Vector{RealFieldElement}, arg2::Vector{RealFieldElement}, arg3::Vector{Vector{RealFieldElement}}, arg4::Vector{Vector{RealFieldElement}})
    return FieldODEStateAndDerivative((RealFieldElement, Vector{RealFieldElement}, Vector{RealFieldElement}, Vector{Vector{RealFieldElement}}, Vector{Vector{RealFieldElement}}), arg0, arg1, arg2, arg3, arg4)
end

function get_complete_derivative(obj::FieldODEStateAndDerivative)
    return jcall(obj, "getCompleteDerivative", Vector{RealFieldElement}, ())
end

function get_primary_derivative(obj::FieldODEStateAndDerivative)
    return jcall(obj, "getPrimaryDerivative", Vector{RealFieldElement}, ())
end

function get_secondary_derivative(obj::FieldODEStateAndDerivative, arg0::jint)
    return jcall(obj, "getSecondaryDerivative", Vector{RealFieldElement}, (jint,), arg0)
end


text:
#####
##### Navier-Stokes and tracer advection equations
#####

""" Calculate the right-hand-side of the u-momentum equation. """
function calculate_Gu!(Gu, grid, coriolis, surface_waves, closure, U, C, K, F, pHY′, parameters, time)
    @loop_xyz i j k grid begin
        @inbounds Gu[i, j, k] = ( - div_ũu(i, j, k, grid, U)
                                  - x_f_cross_U(i, j, k, grid, coriolis, U)
                                  - ∂xᶠᵃᵃ(i, j, k, grid, pHY′)
                                  + ∂ⱼ_2ν_Σ₁ⱼ(i, j, k, grid, closure, U, K)
                                  + x_curl_Uˢ_cross_U(i, j, k, grid, surface_waves, U, time)
                                  + ∂t_uˢ(i, j, k, grid, surface_waves, time)
                                  + F.u(i, j, k, grid, time, U, C, parameters)
        )
    end
    return nothing
end

""" Calculate the right-hand-side of the v-momentum equation. """
function calculate_Gv!(Gv, grid, coriolis, surface_waves, closure, U, C, K, F, pHY′, parameters, time)
    @loop_xyz i j k grid begin
        @inbounds Gv[i, j, k] = ( - div_ũv(i, j, k, grid, U)
                                  - y_f_cross_U(i, j, k, grid, coriolis, U)
                                  - ∂yᵃᶠᵃ(i, j, k, grid, pHY′)
                                  + ∂ⱼ_2ν_Σ₂ⱼ(i, j, k, grid, closure, U, K)
                                  + y_curl_Uˢ_cross_U(i, j, k, grid, surface_waves, U, time)
                                  + ∂t_vˢ(i, j, k, grid, surface_waves, time)
                                  + F.v(i, j, k, grid, time, U, C, parameters))
    end
    return nothing
end

""" Calculate the right-hand-side of the w-momentum equation. """
function calculate_Gw!(Gw, grid, coriolis, surface_waves, closure, U, C, K, F, parameters, time)
    @loop_xyz i j k grid begin
        @inbounds Gw[i, j, k] = ( - div_ũw(i, j, k, grid, U)
                                  - z_f_cross_U(i, j, k, grid, coriolis, U)
                                  + ∂ⱼ_2ν_Σ₃ⱼ(i, j, k, grid, closure, U, K)
                                  + z_curl_Uˢ_cross_U(i, j, k, grid, surface_waves, U, time)
                                  + ∂t_wˢ(i, j, k, grid, surface_waves, time)
                                  + F.w(i, j, k, grid, time, U, C, parameters))
    end
    return nothing
end

""" Calculate the right-hand-side of the tracer advection-diffusion equation. """
function calculate_Gc!(Gc, grid, c, tracer_index, closure, buoyancy, U, C, K, Fc, parameters, time)
    @loop_xyz i j k grid begin
        @inbounds Gc[i, j, k] = (- div_uc(i, j, k, grid, U, c)
                                 + ∇_κ_∇c(i, j, k, grid, closure, c, tracer_index, K, C, buoyancy)
                                 + Fc(i, j, k, grid, time, U, C, parameters))
    end
    return nothing
end

""" Store previous value of the source term and calculate current source term. """
function calculate_interior_source_terms!(G, arch, grid, coriolis, buoyancy, surface_waves, closure, U, C, pHY′, K, F, parameters, time)
    # Manually choose thread-block layout here as it's ~20% faster.
    # See: https://github.com/climate-machine/Oceananigans.jl/pull/308
    Tx, Ty = 16, 16 # CUDA threads per block
    Bx, By, Bz = floor(Int, grid.Nx/Tx), floor(Int, grid.Ny/Ty), grid.Nz  # Blocks in grid

    @launch(device(arch), threads=(Tx, Ty), blocks=(Bx, By, Bz),
            calculate_Gu!(G.u, grid, coriolis, surface_waves, closure, U, C, K, F, pHY′, parameters, time))

    @launch(device(arch), threads=(Tx, Ty), blocks=(Bx, By, Bz),
            calculate_Gv!(G.v, grid, coriolis, surface_waves, closure, U, C, K, F, pHY′, parameters, time))

    @launch(device(arch), threads=(Tx, Ty), blocks=(Bx, By, Bz),
            calculate_Gw!(G.w, grid, coriolis, surface_waves, closure, U, C, K, F, parameters, time))

    for tracer_index in 1:length(C)
        @inbounds Gc = G[tracer_index+3]
        @inbounds Fc = F[tracer_index+3]
        @inbounds  c = C[tracer_index]

        @launch(device(arch), threads=(Tx, Ty), blocks=(Bx, By, Bz),
                calculate_Gc!(Gc, grid, c, Val(tracer_index), closure, buoyancy, U, C, K, Fc, parameters, time))
    end

    return nothing
end

""" Apply boundary conditions by adding flux divergences to the right-hand-side. """
function calculate_boundary_source_terms!(Gⁿ, bcs, arch, grid, args...)

    # Velocity fields
    for i in 1:3
        ubcs = bcs[i]
        apply_z_bcs!(Gⁿ[i], arch, grid, ubcs.z.bottom, ubcs.z.top,   args...)
        apply_y_bcs!(Gⁿ[i], arch, grid, ubcs.y.left,   ubcs.y.right, args...)
    end

    # Tracer fields
    for i in 4:length(bcs)
        cbcs = bcs[i]
        apply_z_bcs!(Gⁿ[i], arch, grid, cbcs.z.bottom, cbcs.z.top,   args...)
        apply_y_bcs!(Gⁿ[i], arch, grid, cbcs.y.left,   cbcs.y.right, args...)
    end

    return nothing
end

#####
##### Pressure-related functions
#####

"Solve the Poisson equation for non-hydrostatic pressure on the CPU."
function solve_for_pressure!(pressure, ::CPU, grid, poisson_solver, ϕ)
    solve_poisson_3d!(poisson_solver, grid)
    view(pressure, 1:grid.Nx, 1:grid.Ny, 1:grid.Nz) .= real.(ϕ)
    return nothing
end

"Solve the Poisson equation for non-hydrostatic pressure on the GPU."
function solve_for_pressure!(pressure, ::GPU, grid, poisson_solver, ϕ)
    solve_poisson_3d!(poisson_solver, grid)
    @launch device(GPU()) config=launch_config(grid, :xyz) idct_permute!(pressure, grid, poisson_solver.bcs, ϕ)
    return nothing
end

"""
Update the hydrostatic pressure perturbation pHY′. This is done by integrating
the `buoyancy_perturbation` downwards:

    `pHY′ = ∫ buoyancy_perturbation dz` from `z=0` down to `z=-Lz`
"""
function update_hydrostatic_pressure!(pHY′, grid, buoyancy, C)
    @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
        @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
            @inbounds pHY′[i, j, grid.Nz] = - ℑzᵃᵃᶠ(i, j, grid.Nz, grid, buoyancy_perturbation, buoyancy, C) * ΔzF(i, j, grid.Nz, grid)
            @unroll for k in grid.Nz-1 : -1 : 1
                @inbounds pHY′[i, j, k] =
                    pHY′[i, j, k+1] - ℑzᵃᵃᶠ(i, j, k+1, grid, buoyancy_perturbation, buoyancy, C) * ΔzF(i, j, k, grid)
            end
        end
    end
    return nothing
end

"""
Calculate the right-hand-side of the Poisson equation for the non-hydrostatic
pressure

    `∇²ϕ_{NH}^{n+1} = (∇·u^n)/Δt + ∇·(Gu, Gv, Gw)`
"""
function calculate_poisson_right_hand_side!(RHS, ::CPU, grid, ::PoissonBCs, U, G, Δt)
    @loop_xyz i j k grid begin
            @inbounds RHS[i, j, k] = divᶜᶜᶜ(i, j, k, grid, U.u, U.v, U.w) / Δt +
                                     divᶜᶜᶜ(i, j, k, grid, G.u, G.v, G.w)
    end

    return nothing
end

"""
Calculate the right-hand-side of the Poisson equation for the non-hydrostatic
pressure and in the process apply the permutation

    [a, b, c, d, e, f, g, h] -> [a, c, e, g, h, f, d, b]

in the z-direction which is required by the GPU fast cosine transform algorithm for
horizontally periodic model configurations.
"""
function calculate_poisson_right_hand_side!(RHS, ::GPU, grid, ::PPN, U, G, Δt)
    Nz = grid.Nz
    @loop_xyz i j k grid begin
        if (k & 1) == 1  # isodd(k)
            k′ = convert(UInt32, CUDAnative.floor(k/2) + 1)
        else
            k′ = convert(UInt32, Nz - CUDAnative.floor((k-1)/2))
        end

        @inbounds RHS[i, j, k′] = divᶜᶜᶜ(i, j, k, grid, U.u, U.v, U.w) / Δt +
                                  divᶜᶜᶜ(i, j, k, grid, G.u, G.v, G.w)
    end
    return nothing
end

"""
Calculate the right-hand-side of the Poisson equation for the non-hydrostatic
pressure and in the process apply the permutation

    [a, b, c, d, e, f, g, h] -> [a, c, e, g, h, f, d, b]

in the y- and z-directions which is required by the GPU fast cosine transform algorithm for
reentrant channel model configurations.
"""
function calculate_poisson_right_hand_side!(RHS, ::GPU, grid, ::PNN, U, G, Δt)
    Ny, Nz = grid.Ny, grid.Nz
    @loop_xyz i j k grid begin
        if (k & 1) == 1  # isodd(k)
            k′ = convert(UInt32, CUDAnative.floor(k/2) + 1)
        else
            k′ = convert(UInt32, Nz - CUDAnative.floor((k-1)/2))
        end

        if (j & 1) == 1  # isodd(j)
            j′ = convert(UInt32, CUDAnative.floor(j/2) + 1)
        else
            j′ = convert(UInt32, Ny - CUDAnative.floor((j-1)/2))
        end

        @inbounds RHS[i, j′, k′] = divᶜᶜᶜ(i, j, k, grid, U.u, U.v, U.w) / Δt +
                                   divᶜᶜᶜ(i, j, k, grid, G.u, G.v, G.w)
    end
    return nothing
end

"""
Copy the non-hydrostatic pressure into `pNHS` and undo the permutation

    [a, b, c, d, e, f, g, h] -> [a, c, e, g, h, f, d, b]

along the z-direction.
"""
function idct_permute!(pNHS, grid, ::PPN, ϕ)
    Nz = grid.Nz
    @loop_xyz i j k grid begin
        if k <= Nz/2
            @inbounds pNHS[i, j, 2k-1] = real(ϕ[i, j, k])
        else
            @inbounds pNHS[i, j, 2(Nz-k+1)] = real(ϕ[i, j, k])
        end
    end
    return nothing
end

"""
Copy the non-hydrostatic pressure into `pNHS` and undo the permutation

    [a, b, c, d, e, f, g, h] -> [a, c, e, g, h, f, d, b]

along the y- and z-direction.
"""
function idct_permute!(pNHS, grid, ::PNN, ϕ)
    Ny, Nz = grid.Ny, grid.Nz
    @loop_xyz i j k grid begin
        if k <= Nz/2
            k′ = 2k-1
        else
            k′ = 2(Nz-k+1)
        end

        if j <= Ny/2
            j′ = 2j-1
        else
            j′ = 2(Ny-j+1)
        end

        @inbounds pNHS[i, j′, k′] = real(ϕ[i, j, k])
    end
    return nothing
end

"""
Update the horizontal velocities u and v via

    `u^{n+1} = u^n + (Gu^{n+½} - δₓp_{NH} / Δx) Δt`

Note that the vertical velocity is not explicitly time stepped.
"""
function update_velocities!(U, grid, Δt, G, pNHS)
    @loop_xyz i j k grid begin
        @inbounds U.u[i, j, k] += (G.u[i, j, k] - ∂xᶠᵃᵃ(i, j, k, grid, pNHS)) * Δt
        @inbounds U.v[i, j, k] += (G.v[i, j, k] - ∂yᵃᶠᵃ(i, j, k, grid, pNHS)) * Δt
    end
    return nothing
end

"""
Update the horizontal velocities u and v via

    `u^{n+1} = u^n + Gu^{n+½} / Δt`

Note that the vertical velocity is not explicitly time stepped.
"""
function update_velocities!(U, grid, Δt, G, ::Nothing)
    @loop_xyz i j k grid begin
        @inbounds U.u[i, j, k] += G.u[i, j, k] * Δt
        @inbounds U.v[i, j, k] += G.v[i, j, k] * Δt
    end
    return nothing
end

"""
Update tracers via

    `c^{n+1} = c^n + Gc^{n+½} Δt`
"""
function update_tracer!(c, grid, Δt, Gc)
    @loop_xyz i j k grid begin
        @inbounds c[i, j, k] += Gc[i, j, k] * Δt
    end
    return nothing
end

"Update the solution variables (velocities and tracers)."
function update_solution!(U, C, arch, grid, Δt, G, pNHS)
    @launch device(arch) config=launch_config(grid, :xyz) update_velocities!(U, grid, Δt, G, pNHS)

    for i in 1:length(C)
        @inbounds c = C[i]
        @inbounds Gc = G[i+3]
        @launch device(arch) config=launch_config(grid, :xyz) update_tracer!(c, grid, Δt, Gc)
    end

    return nothing
end

"""
Compute the vertical velocity w by integrating the continuity equation from the bottom upwards

    `w^{n+1} = -∫ [∂/∂x (u^{n+1}) + ∂/∂y (v^{n+1})] dz`
"""
function compute_w_from_continuity!(model)
    @launch(device(model.architecture), config=launch_config(model.grid, :xy),
            _compute_w_from_continuity!(datatuple(model.velocities), model.grid))
    return nothing
end

function _compute_w_from_continuity!(U, grid)
    @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
        @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
            # U.w[i, j, 1] = 0 is enforced via halo regions.
            @unroll for k in 2:grid.Nz
                @inbounds U.w[i, j, k] = U.w[i, j, k-1] - ΔzC(i, j, k, grid) * hdivᶜᶜᵃ(i, j, k-1, grid, U.u, U.v)
            end
        end
    end

    return nothing
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct FileOperationResult <: SwaggerModel
    boolean::Any # spec type: Union{ Nothing, Bool } # spec name: boolean

    function FileOperationResult(;boolean=nothing)
        o = new()
        validate_property(FileOperationResult, Symbol("boolean"), boolean)
        setfield!(o, Symbol("boolean"), boolean)
        o
    end
end # type FileOperationResult

const _property_map_FileOperationResult = Dict{Symbol,Symbol}(Symbol("boolean")=>Symbol("boolean"))
const _property_types_FileOperationResult = Dict{Symbol,String}(Symbol("boolean")=>"Bool")
Base.propertynames(::Type{ FileOperationResult }) = collect(keys(_property_map_FileOperationResult))
Swagger.property_type(::Type{ FileOperationResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_FileOperationResult[name]))}
Swagger.field_name(::Type{ FileOperationResult }, property_name::Symbol) =  _property_map_FileOperationResult[property_name]

function check_required(o::FileOperationResult)
    true
end

function validate_property(::Type{ FileOperationResult }, name::Symbol, val)
end