text:
using ..ChemistryFeaturization.Utils.BondFeatureUtils
using ..ChemistryFeaturization.Codec
using ..ChemistryFeaturization.Atoms

using MolecularGraph

include("abstractfeatures.jl")

"""
    BondFeatureDescriptor{A}

A feature object that encodes features associated with bonded pairs of atoms in a structure. For general pair features, see `PairFeatureDescriptor`.

Type parameter represents the structure representation(s) from which this feature descriptor is able to compute features.

## Fields
- `name::String`: the name of the feature
- `compute_f`: function that takes in a structure of type `A` and returns values of this feature for every bonded pair of atoms in the structure
- `encoder_decoder::AbstractCodec`: codec that encodes/decodes values of this feature
- `categorical::Bool`: flag for whether the feature is categorical or continuous-valued
- `encodable_elements::Vector{String}`: list of elements (by symbol) that can be encoded by this feature
"""
struct BondFeatureDescriptor{A,C<:AbstractCodec} <: AbstractPairFeatureDescriptor
    name::String
    compute_f::Function
    encoder_decoder::C
    categorical::Bool
    encodable_elements::Vector{String}
end

# pretty printing, long version
function Base.show(io::IO, ::MIME"text/plain", fd::BondFeatureDescriptor{A}) where {A}
    st = "$(typeof(fd)) $(fd.name):\n   categorical: $(fd.categorical)\n   works on: $(A)"
    print(io, st)
end

encodable_elements(fd::BondFeatureDescriptor) = fd.encodable_elements

"""
    BondFeatureDescriptor(name::String)

Construct a `BondFeatureDescriptor` from one of the built-in options (currently supports a variety of MolecularGraph functions, see Utils.BondFeatureUtils.bfd_names_props for details).
"""
function BondFeatureDescriptor(name::String)
    @assert name in keys(bfd_names_props) "That species feature isn't one of the built-in ones, you'll have to construct it directly. Consult ChemistryFeaturiztion.Utils.BondFeatureUtils.bfd_names_props for built-in species features."

    info = bfd_names_props[name]
    categorical = info[:categorical]
    local codec
    if categorical
        possible_vals = info[:possible_vals]
        if length(possible_vals) == 2
            @assert possible_vals == [true, false]
            codec = DirectCodec(1)
        else
            codec = OneHotOneCold(true, possible_vals)
        end
    else
        # TODO: figure out default binning situation for continuous-valued BFD's
        #codec = OneHotOneCold(false, )
    end
    BondFeatureDescriptor{info[:A],typeof(codec)}(
        name,
        info[:compute_f],
        codec,
        categorical,
        info[:encodable_elements],
    )
end
# for some reason having AbstractAtoms instead of AtomGraph doesn't get dispatched propertly, I assume this has something to do with the type parameters...
function get_value(
    bfd::BondFeatureDescriptor{GraphMol},
    a::AtomGraph{GraphMol{A,B}},
) where {A,B}
    @assert all([el in encodable_elements(bfd) for el in elements(a)]) "Feature $(bfd.name) cannot encode some element(s) in this structure!"
    vals = bfd.compute_f(a.structure)
    bonds = a.structure.edges
    mat = Matrix{Union{Missing,eltype(vals)}}(
        missing,
        length(a.structure.nodeattrs),
        length(a.structure.nodeattrs),
    )
    # force it to be symmetric...TODO: add option for asymmetric features
    for i = 1:length(bonds)
        mat[bonds[i][1], bonds[i][2]] = vals[i]
        mat[bonds[i][2], bonds[i][1]] = vals[i]
    end
    return mat
end


text:
# ---
# title: 1441. Build an Array With Stack Operations
# id: problem1441
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Stack
# link: <https://leetcode.com/problems/build-an-array-with-stack-operations/description/>
# hidden: true
# ---
#
# Given an array `target` and an integer `n`. In each iteration, you will read a
# number from  `list = {1,2,3..., n}`.
#
# Build the `target` array using the following operations:
#
#   * **Push** : Read a new element from the beginning `list`, and push it in the array.
#   * **Pop** : delete the last element of the array.
#   * If the target array is already built, stop reading more elements.
#
# You are guaranteed that the target array is strictly increasing, only
# containing numbers between 1 to `n` inclusive.
#
# Return the operations to build the target array.
#
# You are guaranteed that the answer is unique.
#
#
#
# **Example 1:**
#
#
#
#     Input: target = [1,3], n = 3
#     Output: ["Push","Push","Pop","Push"]
#     Explanation: Read number 1 and automatically push in the array -> [1]
#     Read number 2 and automatically push in the array then Pop it -> [1]
#     Read number 3 and automatically push in the array -> [1,3]
#
#
# **Example 2:**
#
#
#
#     Input: target = [1,2,3], n = 3
#     Output: ["Push","Push","Push"]
#
#
# **Example 3:**
#
#
#
#     Input: target = [1,2], n = 4
#     Output: ["Push","Push"]
#     Explanation: You only need to read the first 2 numbers and stop.
#
#
# **Example 4:**
#
#
#
#     Input: target = [2,3,4], n = 4
#     Output: ["Push","Pop","Push","Push","Push"]
#
#
#
#
# **Constraints:**
#
#   * `1 <= target.length <= 100`
#   * `1 <= target[i] <= 100`
#   * `1 <= n <= 100`
#   * `target` is strictly increasing.
#
#
## @lc code=start
using LeetCode

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

ext: jl


text:
# Julia wrapper for header: /usr/include/libavutil/mastering_display_metadata.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_mastering_display_metadata_alloc,
    av_mastering_display_metadata_create_side_data,
    av_content_light_metadata_alloc,
    av_content_light_metadata_create_side_data

function av_mastering_display_metadata_alloc()
    ccall((:av_mastering_display_metadata_alloc, libavutil), Ptr{AVMasteringDisplayMetadata}, ())
end

function av_mastering_display_metadata_create_side_data(frame)
    ccall((:av_mastering_display_metadata_create_side_data, libavutil), Ptr{AVMasteringDisplayMetadata}, (Ptr{AVFrame},), frame)
end

function av_content_light_metadata_alloc(size)
    ccall((:av_content_light_metadata_alloc, libavutil), Ptr{AVContentLightMetadata}, (Ptr{Csize_t},), size)
end

function av_content_light_metadata_create_side_data(frame)
    ccall((:av_content_light_metadata_create_side_data, libavutil), Ptr{AVContentLightMetadata}, (Ptr{AVFrame},), frame)
end


text:
# Behavior subtyping

abstract type Vehicle end

struct Car <: Vehicle end

move(v::Vehicle) = "$v has moved.";

car = Car();
move(car)
#=
julia> move(helicopter)
"Helicopter() has moved."
=#

# Inheritance over multiple levels

abstract type FlyingVehicle <: Vehicle end

liftoff(v::FlyingVehicle) = "$v has lifted off.";

struct Helicopter <: FlyingVehicle end

helicopter = Helicopter();
move(helicopter)
liftoff(helicopter)
#=
julia> move(helicopter)
"Helicopter() has moved."

julia> liftoff(helicopter)
"Helicopter() has lifted off."
=#

# Method overrides

liftoff(h::Helicopter) = "$h has lifted off vertically.";

liftoff(helicopter)
#=
julia> liftoff(h::Helicopter) = "$h has lifted off vertically.";

julia> liftoff(helicopter)
"Helicopter() has lifted off vertically."
=#

# Duck typing

abstract type Animal end

struct Horse <: Animal end

move(h::Horse) = "$h running fast.";

horse = Horse();
move(horse)


text:
using CompileBot

println("tests infer benchmark")

snoop_bench(BotConfig("TestPackage1", tmin =0.0))


text:
using DifferentialEquations, ModelingToolkit
using HoJMTK.Mechanical1D

k = 1.0
m = 1.0
b = 1.0
v = 1.0
@named spring = Spring(k=k)
@named mass = Mass(m=m)
@named damper = Damper(b=b)
@named source = ConstantVelocity(v=v)
@named fixed = Fixed()

# mech1d_eqs = [
#           connect(source.p1, damper.p1)
#           connect(damper.p2, mass.p2)
#           connect(mass.p2, source.p2, fixed.g)
#          ]
mech1d_eqs = [
          connect(source.p1, mass.p1)
          connect(mass.p2, damper.p1, spring.p1)
          connect(damper.p2, spring.p2, fixed.g)
          connect(source.p2, fixed.g)
         ]

@variables t
@named _mech1d_model = ODESystem(mech1d_eqs, t)
@named mech1d_model = compose(_mech1d_model, [spring, mass, damper, source, fixed])
sys = structural_simplify(mech1d_model)
u0 = [
    mass.p1 => 0.0
     ]
prob = ODAEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Tsit5())

using Plots
plot(sol)


text:
using Gaston
g = Gaston;

set_terminal("x11");

c = g.CurveConf();
c.plotstyle = "pm3d";

x = -15:0.33:15;
y = -15:0.33:15;
Z = g.meshgrid(x, y, (x,y) -> sin(sqrt(x .* x+y .* y))
					 / sqrt(x .* x+y .* y));
g.addcoords(x,y,Z,c);

a = g.AxesConf();
a.title = "3D: Sombrero";
g.addconf(a);
g.llplot()

max_stars_repo_path: Module 3/Chapter07/gaston-ex2.jl

max_issues_repo_path: Module 3/Chapter07/gaston-ex2.jl

max_forks_repo_path: Module 3/Chapter07/gaston-ex2.jl


text:
module EBM

include("Computation.jl")
using .Computation

include("Format.jl")
using .Format

include("FC.jl")
using .FC

export ebm

# Determines upper boudn alpha with Easy Buddy Match, optionally outputs proof
# Author: Antara Hebbar
function ebm(m::Int64, s::Int64; output::Int64 = 2)
    output > 0 && printHeader(center("EASY BUDDY MATCH"))

    if !(m > s > 0)
        if output > 0
            printf("Easy Buddy Match does not apply", line = true)
            printEnd()
        end
        return 1 // 1
    elseif m % s == 0
        if output > 0
            printfT("Easy Buddy Match", "Since $m % $s = 0, muffins($m,$s) = 1")
            printEnd()
        end
        return 1 // 1
    end

    (V, W, sV, sW) = sv(m, s)

    alpha = 1
    if V == 3
        d = m - s

        i = 0
        while 3 * d * i < s
            i += 1
        end

        k = i - 1
        a = s - (3 * d * k)

        if 2 * d + 1 <= a <= 3 * d
            alpha = (d * k + a // 3) // (3 * d * k + a)
        elseif a == 2 * d
            alpha = fc(m, s, output = 0)
        elseif 1 <= a <= d
            alpha = (d * k + min(a // 2)) // (3 * d * k + a)
        elseif d <= a <= 2 * d - 1
            alpha = (d * k + min((a + d) // 4)) // (3 * d * k + a)
        end

        if output > 0 && alpha != 1
            if output == 2
                ebmproof(m, s, alpha)
            else
                printfT(
                    "Easy Buddy Match",
                    "Upper bound of muffins($m,$s) is $(formatFrac(alpha))",
                )
                printEnd()
            end
        end
        return alpha
    end

    if output > 0
        printf("Easy Buddy Match inconclusive", line = true)
        printEnd()
    end
    alpha
end

# Output proof of Easy Buddy Match
function ebmproof(m, s, alpha)

    #finding V and W
    V, W, Vnum, Wnum = sv(m, s)
    Wshares = W * Wnum
    Vshares = V * Vnum
    abuddy = 1 - alpha

    #formatting fractions into strings
    den = lcm(s, denominator(alpha))
    total = m // s
    alphaS = formatFrac(alpha, den)
    alphaden = formatFrac(alpha, denominator(alpha))
    totalS = formatFrac(total, den)
    aB = formatFrac(abuddy, den)

    d = m - s

    i = 0
    while 3 * d * i < s
        i += 1
    end

    k = i - 1
    a = s - (3 * d * k)
    comden = 3 * d * k + a

    if V == 3 #ebm only works if there are 2-shares
        firstparameter = formatFrac(a // 3)
        case1 = (3 * d * k + a + d) // (comden)
        case1str = formatFrac(case1)
        V_case1 = formatFrac(1 // (V + 1))
        v_conj1 = formatFrac(case1 * (1 // (V + 1)))
        total = formatFrac(m // s)

        #claim
        printHeader("CLAIM:")
        printf("There is a ($m, $s) procedure where the smallest piece is ≥ $alphaS.")
        printLine()
        printfT(
            "Note",
            "Let the common denominator $den. Therefore, alpha will be referred to as $alphaden.",
        )

        printHeader("DEFINING TERMS FOR EBM ALGORITHM: ")
        printfT(
            "Variables",
            "We will use the following variables for the EBM algorithm (where m=muffins and s=students):",
            "",
            "d = m-s = $d",
            "",
            "k is the largest number >= 0 such that 3dk > s. k = $k.",
            "",
            "a = s - 3dk = $a",
        )
        printfT(
            "Common Terms",
            "A buddy of share size x is equal to 1-x. A match of share size y is equal to $m/$s-y.",
        )

        if a == 2 * d
            a = fc(m, s)
            printf("Alpha equals FC($m, $s) because a=2d.")
        else

            #assumptions
            printHeader("ASSUMPTIONS FOR EBM ALGORITHM:")
            printfT(
                "Important note",
                "The proof is conducted using parameter X. Assume that there is a ($m, $s) procedure where X ≥ $(firstparameter) so alpha ≥ $firstparameter.",
                "",
                "Hence, there is a procedure where every muffin is cut into 2 pieces. There are total $(6*d*k + 2*a +2*d) pieces.",
            )

            #CASEWORK, checking if V+1 or V-2 cases are impossible
            printHeader("CASEWORK:")

            if case1 * (1 // (V + 1)) < alpha
                printfT(
                    "Case 1",
                    "If Alice gets ≥ $(V+1) shares then some share is ≤ ",
                    "",
                    "$case1str * $V_case1 ≤ $v_conj1, which is < $alphaden",
                    "",
                    "Note: the inequality follows X > a/3 > a/4 and k ≥ 1",
                )
            end

            if (V - 2) == 1 && (3 * d * k + a + d) // (comden) > 1 // 2
                printfT(
                    "Case 2",
                    "If Alice gets $(V-2) share, since each share is ≤ 1/2, she has:",
                    "",
                    "≤ 1/2 ≤ $(3*d*k+a+d)/$(comden), which is impossible",
                )
            end

            printfT(
                "Case 3",
                "Everyone is a $V or $(V-1) student. Let s_2 & s_3 be the number of 2-students and 3-students. Since there are $(6*d*k + 2*a +2*d) shares: ",
                "",
                "",
                "$(V)s_$V + $(V-1)s_$(V-1) = $(6*d*k + 2*a +2*d)",
                "",
                "s_$V + s_$(V-1) = $(3*d*k +a)",
                "",
                "Hence, s_$(V-1) = $(3*d*k+a-2*d), s_$(V) = $(2*d). So there are $(6*d*k+2*a-4*d) $(V-1) shares and $(6*d) $V-shares.",
            )

            printHeader("CONTINUING CASEWORK")
            printfT(
                "Case 3a",
                "There is a $V-share ≥ ($(d*k+a+d)-2X)/$comden. The remaining $(V-1) $V-shares add up to ",
                "",
                "≤ $m/$s - ($(d*k+a+d)-2X)/$comden = ($(2*d*k)+2X)/$comden.",
                "",
                "Hence there is a $V-share ($(2*d*k)+2X)/$comden * 1/2 = ($(d*k)+X)/$comden",
            )
            printfT(
                "Case 3b",
                "There is a $(V-1) share ≥ ($(2*d*k+a)-X)/$comden. Its buddy is ",
                "",
                "≤ 1 - ($(2*d*k+a)-X)/$comden  = ($(d*k)+X)/$comden",
            )
            printfT(
                "Case 3c",
                "There is a $(V-1) share ≤ ($(d*k+d)+X)/$comden. Its match is ",
                "",
                "≥ $total - ($(d*k+d)+X)/$comden  = ($(2*d*k+a)-X)/$comden",
                "",
                "Its buddy is ≤ 1 - ($(2*d*k+a)-X)/$comden  = ($(d*k)+X)/$comden",
            )
            printfT(
                "Case 3d",
                "There is a $V-share ≤ ($(d*k)+X)/$comden. This case is obvious.",
            )
            printfT(
                "Case 3e",
                "All $V-shares are in (($(d*k)+X)/$comden, ($(d*k+a+d)-2X)/$comden))",
                "",
                "All $(V-1) shares are in (($(d*k+d)+X)/$comden, ($(2*d*k+a)-X)/$comden)",
            )

            printHeader("INTERVAL DIAGRAM:")
            printf("The following diagram depicts what we know so far: ")
            println(
                "\n",
                interval(
                    ["(", "($(d*k)+X)/$comden"],
                    [")[", "($(d*k+a+d)-2X)/$comden)"],
                    ["](", "($(d*k+d)+X)/$comden"],
                    [")", "($(2*d*k+a)-X)/$comden"],
                    labels = ["$(V)s_$V $V-shs", "0", "$(V-1)s_$(V-1) $(V-1)-shs"],
                ),
            )
            bound = formatFrac(a // 3)
            printf("The assumption that X ≥ $bound ensures that the interval of $V-shares and $(V-1) shares do not intersect.")
            printLine()

            j = k - 1
            #Define bijections
            printHeader("BUDDY/MATCH ANALYSIS:")

            printfT(
                "Buddy/Match Bijections",
                "We define the following (where B signifies buddying and M signifies matching):",
                "",
                "M₀ = (($(d*k)+X)/$comden, ($(d*k+a+d)-2X)/$comden))",
                "B₀ = B(M₀) = (($(2*d*k+a-d)-X)/$comden), ($(2*d*k+a)-X)/$comden)",
                "",
                "∀ ≤ i ≤ $(k-1), [Mᵢ = M(Bᵢ₋₁) = (($(d*k)+$(d)i+X)/$comden, $(d*k)+$(d)(i+1)+X)/$comden)]",
                "",
                "∀ ≤ i ≤ $(k-1), [Bᵢ = B(Mᵢ) = (($(2*d*k+a)-$(d)(i-1)-X)/$comden, ($(2*d*k+a)-$(d)(i)-X)/$comden)]",
            ) #printing bijections

            printfT(
                "Important note",
                "The last B (M) is B𝒌 (M𝒌₋₁)).",
                "",
                "Since we might have an interval that overlaps the 3-share region, and can't apply M to it, we stop here.",
                "",
                "If k ≥ 2, B𝒌₋₂ is in the same $(V-1)-share region and M() can be applied to it. This is good since M𝒌₋₁ needs B𝒌₋₂ to be in the $(V-1)-shr region.",
            ) #defining restrictions

            printfT(
                "Translating into endpoints",
                "Therefore, we want the left endpoint of B𝒌₋₂ ≥ the right endpoint of the $V-shares.",
                "",
                "($(d*k+a+d)-2X)/$comden) ≤ ($(2*d*k+a-(k-1)*d)-X)/$comden",
                "This is equal to 1 ≤ 2, a true statement.",
            )
            printfT(
                "Final note",
                "|B₀| = $(6*d). Since B and M are bijections: ",
                "",
                "1. every Bᵢ and Mᵢ have $(6*d) shares",
                "",
                "2. if x is an endpoint of Bᵢ or Mᵢ, there are no shares size x",
            ) #final point

            contra = true

            printHeader("CONTRADICTION: ")
            if 2 * d + 1 <= a <= 3 * d #where X=a/3
                printfT(
                    "Case 3e.1",
                    "Assume by contradiction that X > $a/3.",
                    "",
                    "We want the left endpoint of B𝒌₋₁ to be ≤ the right endpoint of the two-shares. Hence:",
                    "",
                    "($(2*d*k+a - k*d)-X)/$comden ≤ ($(d*k+a+d)-2X)/$comden",
                    "",
                    "This is X<d, which is true since X=$(a)/3 and a ≤ 3d which equals $a ≤ $(3*d)",
                    "",
                    "B₀∪...∪B𝒌₋₁ ⊇ (($(d*k+d)+X)/$comden, ($(2*d*k+a)-X)/$comden)",
                    "",
                    "(Right hand set is $(V-1) shares)",
                    "",
                    "|B₀∪...∪B𝒌₋₁| ≥(($(d*k+d)+X)/$comden, ($(2*d*k+a)-X)/$comden)",
                    "",
                    "6*d*k ≥ 6*d*k + 2*a - 4*d",
                    "",
                    "Simplifying: $(2*d)≥$a. This is our contradiction.",
                )

            elseif 1 <= a <= d #where X=a/2
                printfT(
                    "Case 3e.2",
                    "Assume by way of contradiction that X>$a/2.",
                    "",
                    "We want the left endpoint of B𝒌₋₁ to be ≤ the right endpoint of the three-shares. Hence:",
                    "",
                    "($(2*d*k+a - k*d)-X)/$comden ≤ ($(d*k)+X)/$comden, which is X= $a/2.",
                    "",
                    "Since X=$a/2 is true by assumption",
                    "",
                    "B₀∪...∪B𝒌₋₁ ⊇ (($(d*k)+X)/$comden, ($(2*d*k+a)-X)/$comden)",
                    "",
                    "(Right hand set is $(V-1) shares)",
                    "",
                    "|B₀∪...∪B𝒌₋₁| ≥(($(d*k)+X)/$comden, ($(2*d*k+a)-X)/$comden)",
                    "",
                    "6*d*k ≥ 6*d*k + 2*a + 2*d",
                    "",
                    "Simplifying: a+d=0 or $a+$d=0. This is our contradiction, since d≥1.",
                )

            elseif d <= a <= 2 * d - 1 #where X=(a+d)/4 this is unfinished - have not tested yet.
                printfT(
                    "Case 3e.3",
                    "Assume by way of contradiction that X > $(a+d)/4 and X ≤ $a/2.",
                    "",
                    "We will use a buddy-match method to find intervals that cover the entire inverval, causing a contradiction.",
                    "",
                    "M₀ = (($(d*k)+X)/$comden, ($(d*k+d)-2X)/$comden))",
                    "B₀ = B(M₀) = (($(2*d*k+a-d)-X)/$comden), ($(2*d*k+a)-X)/$comden)",
                    "",
                    "∀ ≤ i ≤ $(k), [Mᵢ = M(Bᵢ₋₁) = (($(d*k)+$(d)i+X)/$comden, $(d*k)+$(d)(i+1)+X)/$comden)]",
                    "",
                    "∀ ≤ i ≤ $(k), [Bᵢ = B(Mᵢ) = (($(2*d*k+a)-$(d)(i-1)-X)/$comden, ($(2*d*k+a)-$(d)(i)-X)/$comden)]",
                    "",
                    "Note: B𝒌₋₁ = (($(d*k+a)-X)/$comden, ($(d*k+a+d)-X)/$comden",
                    "",
                    "We need ($(d*k+a)-X)/$comden ≤ ($(d*k+a+d)-X. This is true.",
                    "",
                    "We also want M₀∪...∪M𝒌 to cover the entire interval. We need ($(2*d*k+a)-X)/$comden ≤ ($(d*k)+$(d*(k+1))+X)/$comden, which is X ≥ $(a-d)/2.",
                    "",
                    "We have this since X > $(a+d)/4. Since |B₀|=$(6*d) and all Bᵢ and Mᵢ have the same shares, no shares have the same size as the endpoints of any Bᵢ and Mᵢ. Therefore:",
                    "",
                    "|M₀∪...∪M𝒌| = $(6*d)k +2a + 2d = $(6(k-1)),  which implies a=2d, or $a=$(2*d). This is our contradiction.",
                )
            else
                printf("No contradiction found, ebm failed")
                contra = false
            end

            if contra
                printHeader("CONCLUSION:")
                printfT(
                    "Conclude alpha",
                    "Each possible case verifies a lower bound for α that contradicts the assumption",
                    "",
                    "All possible cases contradict the assumption iff. α = $alpha",
                )
            end

        end
    end
end

end


text:
N = 8000
#N = 100

srand(1)

A = complex(sprand(N, N, 0.1))

v1 = rand(Complex128, N)
result = rand(Complex128, N)

function mv1(alpha::Complex128, M::SparseMatrixCSC{Complex128}, v::Vector{Complex128}, beta::Complex128, result::Vector{Complex128})
    @inbounds for i=1:length(result)
        result[i] *= beta
    end
    @inbounds for colindex = 1:M.n::Int
        vj = alpha*v[colindex]
        for i=M.colptr[colindex]:M.colptr[colindex+1]-1
            result[M.rowval[i]] += M.nzval[i]*vj
        end
    end
    nothing
end

function mv1(alpha::Complex128, v::Vector{Complex128}, M::SparseMatrixCSC{Complex128}, beta::Complex128, result::Vector{Complex128})
    @inbounds for i=1:length(result)
        result[i] *= beta
    end
    @inbounds for colindex = 1:M.n::Int
        for i=M.colptr[colindex]:M.colptr[colindex+1]-1
            result[colindex] += M.nzval[i]*alpha*v[M.rowval[i]]
        end
    end
    nothing
end

@time mv1(Complex(1.), A, v1, Complex(0.), result)
@time mv1(Complex(1.), A, v1, Complex(0.), result)

println(norm(result-A*v1))

@time mv1(Complex(1.), v1, A, Complex(0.), result)
@time mv1(Complex(1.), v1, A, Complex(0.), result)

println(norm(transpose(result)-transpose(v1)*A))


text:
#============================================

Benchmarks for generating and performing initial insertions into
a skip list

============================================#

using BenchmarkTools
using Random
using SkipLists
using UUIDs

function generate_data(N)
    collect(uuid4() for ii = 1:N)
end

function benchmark_generate_vector(; N = 100_000, kws...)
    function gen_list(X::Vector{T}) where T
        list = Vector{T}(undef, 0)
        for ii in X
            insert_idx = searchsorted(list, ii)
            insert!(list, first(insert_idx), ii)
        end
    end
    @benchmark $gen_list(X) setup=(X = generate_data($N))
end

function benchmark_generate_skiplist(::Type{SkipList}; N = 100_000, kws...)
    function gen_list(X::Vector{T}) where T
        list = SkipList{T}(; kws...)
        for ii in X
            insert!(list, ii)
        end
    end
    @benchmark $gen_list(X) setup=(X = generate_data($N))
end


text:
include("taxes.jl")
include("utility_revenue.jl")

function monthly_bill(energy_dict::Dict, consumer::Residential; print_output = false, tax_cutoff = 280, SD=SD)

    savings = Dict{String, Any}()
    bill = Dict{String,Any}()

    grid_energy_cost = 0.0
    counterfactual_cost = 0.0

    for block in consumer.tariff.e_cost
        energy_dict["utility_supplied_energy"] > block[1][end] ? grid_energy_cost += block[1][end]*block[2] : grid_energy_cost += (energy_dict["utility_supplied_energy"] - block[1][1])*block[2]*(energy_dict["utility_supplied_energy"] > block[1][1])
    end

    for block in consumer.tariff.e_cost
        energy_dict["consumer_energy"] > block[1][end] ?  counterfactual_cost += block[1][end]*block[2] :  counterfactual_cost += (energy_dict["consumer_energy"] - block[1][1])*block[2]*(energy_dict["consumer_energy"] > block[1][1])
    end

    bill["Balance"] = energy_dict;

    bill["grid_cost"] = grid_energy_cost < 0.01 ? consumer.tariff.min_cost : grid_energy_cost
    bill["withdrawn_energy_cost"] = energy_dict["withdrawn_energy"]*consumer.tariff.access

    bill["total_energy_cost"] = grid_energy_cost + bill["withdrawn_energy_cost"]
    savings["grid_energy_savings"] = counterfactual_cost - grid_energy_cost

    bill["demand_charge"] = consumer.tariff.p_cost[1][2]*energy_dict["peak_demand"]
    bill["counterfactual_demand_charge"] = consumer.tariff.p_cost[1][2]*energy_dict["peak_power"]

    savings["demand_charges"] = bill["counterfactual_demand_charge"] - bill["demand_charge"]

    bill["street_light"] = street_light(energy_dict["grid_energy"], consumer.tariff.street_light)
    savings["street_light"] = street_light(energy_dict["consumer_energy"], consumer.tariff.street_light) - bill["street_light"]

    bill["VAT"] = VAT(energy_dict["grid_energy"], grid_energy_cost + bill["demand_charge"] + bill["withdrawn_energy_cost"], cutoff = tax_cutoff)
    savings["VAT"] = VAT(energy_dict["consumer_energy"], counterfactual_cost + bill["counterfactual_demand_charge"], cutoff = tax_cutoff) -  bill["VAT"]

    bill["firefighters"] = firefighters(energy_dict["grid_energy"], bill["total_energy_cost"] + bill["demand_charge"])
    savings["firefighters"] = firefighters(energy_dict["consumer_energy"], counterfactual_cost + bill["counterfactual_demand_charge"]) - bill["firefighters"]

    bill["total_cost"] = bill["total_energy_cost"] + bill["street_light"] + bill["VAT"] + bill["firefighters"] + bill["demand_charge"]

    bill["counterfactual_cost"] = counterfactual_cost + bill["counterfactual_demand_charge"] +
                                  street_light(energy_dict["consumer_energy"], consumer.tariff.street_light) +
                                  VAT(energy_dict["consumer_energy"], counterfactual_cost + bill["counterfactual_demand_charge"], cutoff = tax_cutoff) +
                                  firefighters(energy_dict["consumer_energy"], counterfactual_cost + bill["counterfactual_demand_charge"])

    bill["total_savings"] = bill["counterfactual_cost"] - bill["total_cost"]

    bill["savings"] = savings

	utility_revenue_change(bill, SD)

    if print_output
            println(round(bill["grid_cost"], digits=2), " ",
                    round(bill["withdrawn_energy_cost"], digits =2), " ",
                    round(bill["demand_charge"], digits = 2), " ",
                    round(bill["total_cost"], digits = 2), " ",
                    " | ",
                    round(bill["counterfactual_cost"], digits=2)," ",
                    round( bill["total_savings"], digits=2),)
    end

    return bill

end

function monthly_bill(energy_dict::Dict, consumer::CommIndus; print_output = false, SD=SD)

    savings = Dict{String, Any}()
    bill = Dict{String,Any}()

    grid_energy_cost = 0.0
    demand_charges = 0.0
    counterfactual_energy_cost = 0.0
    counterfactual_demand_cost = 0.0

    #energy costs

    if energy_dict["utility_supplied_energy"] <= consumer.tariff.e_cost[1][1][end]
        grid_energy_cost = energy_dict["utility_supplied_energy"]*consumer.tariff.e_cost[1][2]
    elseif energy_dict["utility_supplied_energy"] > consumer.tariff.e_cost[1][1][end]
        grid_energy_cost = energy_dict["utility_supplied_energy"]*consumer.tariff.e_cost[2][2]
    end

    if energy_dict["consumer_energy"] <= consumer.tariff.e_cost[1][1][end]
        counterfactual_energy_cost = energy_dict["consumer_energy"]*consumer.tariff.e_cost[1][2]
    elseif energy_dict["consumer_energy"] > consumer.tariff.e_cost[1][1][end]
        counterfactual_energy_cost = energy_dict["consumer_energy"]*consumer.tariff.e_cost[2][2]
    end

    #demand charges

    if !(energy_dict["grid_energy"] <= consumer.tariff.e_cost[1][1][end])

        if energy_dict["peak_demand"] <= consumer.tariff.p_cost[1][1][end]
            demand_charges = consumer.tariff.p_cost[1][2]
        elseif energy_dict["peak_demand"] > consumer.tariff.p_cost[1][1][end]
            demand_charges = consumer.tariff.p_cost[2][2]*energy_dict["peak_demand"]
        end
    end

    if !(energy_dict["consumer_energy"] <= consumer.tariff.e_cost[1][1][end])

        if energy_dict["peak_power"] <= consumer.tariff.p_cost[1][1][end]
            counterfactual_demand_cost = consumer.tariff.p_cost[1][2]
        elseif energy_dict["peak_power"] > consumer.tariff.p_cost[1][1][end]
            counterfactual_demand_cost = consumer.tariff.p_cost[2][2]*energy_dict["peak_power"]
        end

    end

    bill["Balance"] = energy_dict;
    bill["counterfactual_energy_cost"] = counterfactual_energy_cost
    bill["grid_cost"] = grid_energy_cost < 0.01 ? consumer.tariff.min_cost : grid_energy_cost

    savings["grid_energy_savings"] = counterfactual_energy_cost - grid_energy_cost
    bill["withdrawn_energy_cost"] = energy_dict["withdrawn_energy"]*consumer.tariff.access
    bill["total_energy_cost"] = grid_energy_cost + bill["withdrawn_energy_cost"]

    bill["demand_charges"] = demand_charges
    bill["counterfactual_demand_charge"] = counterfactual_demand_cost
    savings["demand_charges"] = bill["counterfactual_demand_charge"] - bill["demand_charges"]

    bill["street_light"] = street_light(energy_dict["grid_energy"], consumer.tariff.street_light)
    savings["street_light"] = street_light(energy_dict["consumer_energy"], consumer.tariff.street_light) - bill["street_light"]

    bill["VAT"] = VAT(energy_dict["grid_energy"], grid_energy_cost+bill["withdrawn_energy_cost"]+bill["demand_charges"])
    savings["VAT"] = VAT(energy_dict["consumer_energy"], counterfactual_energy_cost+counterfactual_demand_cost) - bill["VAT"]

    bill["firefighters"] = firefighters(energy_dict["grid_energy"], bill["total_energy_cost"]+bill["demand_charges"])
    savings["firefighters"] = firefighters(energy_dict["consumer_energy"], counterfactual_energy_cost+counterfactual_demand_cost) - bill["firefighters"]

    bill["total_cost"] = bill["total_energy_cost"] +  bill["demand_charges"] + bill["street_light"] + bill["VAT"] + bill["firefighters"]

    bill["counterfactual_cost"] = counterfactual_energy_cost + counterfactual_demand_cost +
                                  street_light(energy_dict["consumer_energy"], consumer.tariff.street_light) +
                                  VAT(energy_dict["consumer_energy"], counterfactual_energy_cost+counterfactual_demand_cost) +
                                  firefighters(energy_dict["consumer_energy"], counterfactual_energy_cost+counterfactual_demand_cost)

    bill["total_savings"] = bill["counterfactual_cost"] - bill["total_cost"]

    bill["savings"] = savings

	utility_revenue_change(bill, SD)

    if print_output
            println(round(bill["grid_cost"], digits=2), " ",
                    round(bill["withdrawn_energy_cost"], digits =2), " ",
                    round(bill["demand_charges"], digits = 2), " ",
                    round(bill["total_cost"], digits = 2), " ",
                    " | ",
                    round(bill["counterfactual_energy_cost"], digits=2)," ",
                    round(bill["counterfactual_demand_charge"], digits=2)," ",
                    round(bill["counterfactual_cost"], digits=2)," ",
                     " | ",
                    round( bill["total_savings"], digits=2),)
    end

    return bill

end

function monthly_bill(energy_dict::Dict, consumer::TMT; print_output = false, SD=SD)

	savings = Dict{String, Any}()
    bill = Dict{String,Any}()

	grid_energy = energy_dict["grid_energy"]

    grid_energy_cost_peak = energy_dict["utility_supplied_energy_peak"]*consumer.tariff.e_cost[1][2]
    grid_energy_cost_valley = energy_dict["utility_supplied_energy_valley"]*consumer.tariff.e_cost[2][2]
    grid_energy_cost_night = energy_dict["utility_supplied_energy_night"]*consumer.tariff.e_cost[3][2]

    counterfactual_energy_cost_peak = energy_dict["consumer_energy_peak"]*consumer.tariff.e_cost[1][2]
    counterfactual_energy_cost_valley = energy_dict["consumer_energy_valley"]*consumer.tariff.e_cost[2][2]
    counterfactual_energy_cost_night = energy_dict["consumer_energy_night"]*consumer.tariff.e_cost[3][2]

    grid_demand_cost_peak = energy_dict["peak_demand_peak"]*consumer.tariff.p_cost[1][2]
    grid_demand_cost_valley = energy_dict["peak_demand_valley"]*consumer.tariff.p_cost[2][2]
    grid_demand_cost_night = energy_dict["peak_demand_night"]*consumer.tariff.p_cost[3][2]

    counterfactual_demand_cost_peak = energy_dict["peak_power_peak"]*consumer.tariff.p_cost[1][2]
    counterfactual_demand_cost_valley = energy_dict["peak_power_valley"]*consumer.tariff.p_cost[2][2]
    counterfactual_demand_cost_night = energy_dict["peak_power_night"]*consumer.tariff.p_cost[3][2]

    bill["Balance"] = energy_dict;
    #Energy Costs
    bill["grid_energy_cost_peak"] = grid_energy_cost_peak
    bill["grid_energy_cost_valley"] = grid_energy_cost_valley
    bill["grid_energy_cost_night"] = grid_energy_cost_night
    bill["withdrawn_energy_cost"] = energy_dict["withdrawn_energy"]*consumer.tariff.access

	#Demand Costs
    bill["grid_demand_charge_peak"] = grid_demand_cost_peak
    bill["grid_demand_charge_valley"] = grid_demand_cost_valley
    bill["grid_demand_charge_night"] = grid_demand_cost_night

	grid_taxable_energy_cost = (grid_energy_cost_peak + grid_energy_cost_valley + grid_energy_cost_night) < 0.01 ? consumer.tariff.min_cost : (grid_energy_cost_peak + grid_energy_cost_valley + grid_energy_cost_night)

    grid_taxable_cost = grid_taxable_energy_cost + grid_demand_cost_peak + grid_demand_cost_valley + grid_demand_cost_night + bill["withdrawn_energy_cost"]

	 bill["counterfactual_energy_cost_peak"] = counterfactual_energy_cost_peak
    bill["counterfactual_energy_cost_valley"] = counterfactual_energy_cost_valley
    bill["counterfactual_energy_cost_night"] = counterfactual_energy_cost_night

    bill["counterfactual_demand_charge_peak"] = counterfactual_demand_cost_peak
    bill["counterfactual_demand_charge_valley"] = counterfactual_demand_cost_valley
    bill["counterfactual_demand_charge_night"] = counterfactual_demand_cost_night

    counterfactual_grid_taxable_cost = counterfactual_energy_cost_peak + counterfactual_energy_cost_valley + counterfactual_energy_cost_night + counterfactual_demand_cost_peak + counterfactual_demand_cost_valley + counterfactual_demand_cost_night

    counterfactual_total_grid_energy=energy_dict["consumer_energy_peak"]+energy_dict["consumer_energy_valley"]+energy_dict["consumer_energy_night"]

    bill["street_light"] = street_light(grid_energy,  consumer.tariff.street_light)
	savings["street_light"] =  street_light(counterfactual_total_grid_energy,  consumer.tariff.street_light) - bill["street_light"]

    bill["VAT"] = VAT(grid_energy, grid_taxable_cost)
	savings["VAT"] = VAT(counterfactual_total_grid_energy,  counterfactual_grid_taxable_cost) - bill["VAT"]

    bill["firefighters"] = firefighters(grid_energy, grid_taxable_cost)
	savings["firefighters"] = firefighters(counterfactual_total_grid_energy,counterfactual_grid_taxable_cost) - bill["firefighters"]

    bill["total_cost"] = grid_taxable_cost +
                                  bill["street_light"] +
                                  bill["VAT"] +
                                  bill["firefighters"]

    bill["counterfactual_cost"] = counterfactual_grid_taxable_cost +
                                  street_light(counterfactual_total_grid_energy,  consumer.tariff.street_light) +
                                  VAT(counterfactual_total_grid_energy,  counterfactual_grid_taxable_cost) +
                                  firefighters(counterfactual_total_grid_energy,  counterfactual_grid_taxable_cost)

    bill["savings_energy_cost_peak"]     = counterfactual_energy_cost_peak - grid_energy_cost_peak
    bill["savings_energy_cost_valley"]   = counterfactual_energy_cost_valley - grid_energy_cost_valley
    bill["savings_energy_cost_night"]    = counterfactual_energy_cost_night - grid_energy_cost_night
    bill["savings_demand_charge_peak"]   =  counterfactual_demand_cost_peak - grid_demand_cost_peak
    bill["savings_demand_charge_valley"] = counterfactual_demand_cost_valley - grid_demand_cost_valley
    bill["savings_demand_charge_night"]  = counterfactual_demand_cost_night - grid_demand_cost_night

	bill["total_savings"] = bill["counterfactual_cost"] -   bill["total_cost"]

    savings["savings_energy_cost_peak"]     = bill["savings_energy_cost_peak"]
	savings["savings_energy_cost_valley"] = bill["savings_energy_cost_valley"]
	savings["savings_energy_cost_night"]  = bill["savings_energy_cost_night"]

	savings["savings_demand_charge_peak"] = bill["savings_demand_charge_peak"]
	savings["savings_demand_charge_valley"]= bill["savings_demand_charge_valley"]
	savings["savings_demand_charge_night"]  = bill["savings_demand_charge_night"]

	bill["savings"] = savings

	utility_revenue_change(bill, SD)

    if print_output
            println(round(bill["grid_energy_cost_peak"], digits=2), " ",
                    round(bill["grid_energy_cost_valley"], digits=2), " ",
                    round(bill["grid_energy_cost_night"], digits=2), " ",
                    " | ",
                    round(bill["withdrawn_energy_cost"], digits=2), " ",
                    " | ",
                    round(bill["counterfactual_energy_cost_peak"], digits=2)," ",
                    round(bill["counterfactual_energy_cost_valley"], digits=2)," ",
                    round(bill["counterfactual_energy_cost_night"],  digits=2)," ",
                     " | ",
                    round(bill["grid_demand_charge_peak"], digits=2), " ",
                    round(bill["grid_demand_charge_valley"], digits=2), " ",
                    round(bill["grid_demand_charge_night"], digits=2), " ",
                    " | ",
                    round(bill["counterfactual_demand_charge_peak"], digits=2)," ",
                    round(bill["counterfactual_demand_charge_valley"], digits=2)," ",
                    round(bill["counterfactual_demand_charge_night"],  digits=2)," ",
                     " | ",
                    round(bill["total_savings"], digits=2)," ",)

    end

    return bill

end

max_stars_repo_name: arnavgautam/PVCR

max_issues_repo_name: arnavgautam/PVCR

max_forks_repo_name: arnavgautam/PVCR


text:
# Script for processing BATSRUS implicit shocktube test weak scaling runs.
#
# The runlogs are stored in Google Drive, ScalingTests.
#
# Hongyang Zhou, hyzhou@umich.edu 05/09/2019, last modified 09/03/2020

include("Utility.jl")

using PyPlot

Machine = "Frontera"

if Machine == "Bluewaters"
   # New results after modifying the message passing
   result = postprocess("results/gnu_Impl")

   result_0  = sort(result[result.Thread.==0, :],(:Proc))
   result_1  = sort(result[result.Thread.==1, :],(:Proc))
   result_2  = sort(result[result.Thread.==2, :],(:Proc))
   result_4  = sort(result[result.Thread.==4, :],(:Proc))
   result_8  = sort(result[result.Thread.==8, :],(:Proc))
   result_16 = sort(result[result.Thread.==16,:],(:Proc))

   time_0 = copy(result_0.Time)
   time_1 = copy(result_1.Time)
   time_2 = copy(result_2.Time)
   time_4 = copy(result_4.Time)
   time_8 = copy(result_8.Time)
   time_16 = copy(result_16.Time)

   append!(time_0,fill(NaN,size(time_16)[1]-size(time_0)[1]))
   append!(time_1,fill(NaN,size(time_16)[1]-size(time_1)[1]))
   append!(time_2,fill(NaN,size(time_16)[1]-size(time_2)[1]))
   append!(time_4,fill(NaN,size(time_16)[1]-size(time_4)[1]))
   append!(time_8,fill(NaN,size(time_16)[1]-size(time_8)[1]))

   # Calculate the number of cell updates per second for 2nd order scheme
   cell_size = sort(8^3 .* 256 .* result_16.Thread .* result_16.Proc .* 2)

   speed_0 = cell_size ./ time_0
   speed_1 = cell_size ./ time_1
   speed_2 = cell_size ./ time_2
   speed_4 = cell_size ./ time_4
   speed_8 = cell_size ./ time_8
   speed_16 = cell_size ./ time_16

   proc = result_16.Proc .* result_16.Thread
   speed_ref = speed_0[1] ./ proc[1] .* proc

   plt.figure(figsize=(18,5))

   speed = hcat(speed_0, speed_1, speed_2, speed_4, speed_8, speed_16)
   ax1 = plt.subplot(121)
   plt.plot(proc[1:11], speed[1:11,:]/1e9, "-o", markersize=2, alpha=0.8)
   plt.plot(proc[1:11], speed_ref[1:11]/1e9, ":")
   plt.legend(loc="upper left",
         ("pure MPI", "1 thread", "2 threads", "4 threads", "8 threads",
          "16 threads", "linear ref."))
   plt.grid(true)
   ax1.xaxis.set_major_locator(plt.MaxNLocator(5))
   plt.title("BATSRUS Weak Scaling, Implicit Scheme", y=1.03, fontsize=16)
   plt.xlabel("number of cores", fontsize=16)
   plt.ylabel(L"giga cell ($10^9$) updates per second", fontsize=16)

   speed = hcat(speed_0, speed_1, speed_2, speed_4, speed_8, speed_16)
   ax2 = plt.subplot(122)
   plt.plot(proc, speed/1e9, "-o", markersize=3, alpha=0.6)
   plt.plot(proc, speed_ref/1e9, ":")
   plt.legend(loc="upper left",
         ("pure MPI", "1 thread", "2 threads", "4 threads", "8 threads",
          "16 threads", "linear ref."))
   plt.grid(true)
   ax2.xaxis.set_major_locator(plt.MaxNLocator(5))
   plt.title("BATSRUS Weak Scaling, Implicit Scheme", y=1.03, fontsize=16)
   plt.xlabel("number of cores", fontsize=16)
   plt.ylabel(L"giga cell ($10^9$) updates per second", fontsize=16)
   #plt.text(0.95,-0.08, "number of cores", transform=ax2.transAxes)

   #=
   mylabel = ["1 thread", "2 thread", "4 thread", "8 thread", "16 threads",
      "", "", "", "", ""]
   mylabel = reshape(mylabel,1,:)

   p1 = plot(proc,speednew, seriestype=[:line, :scatter],
        markersize=3, lw=1.5, markercolor=[1 2 3 4 5],
        markerstrokewidth = 0.1,
        linecolor=[1 2 3 4 5],
        xscale=:identity, yscale=:identity, legend=:best,
        label=mylabel,
        xlabel="number of cores",
        ylabel="cell updates per second",
        title="BATSRUS Hybrid Weak Scaling for Implicit Scheme",
        xtickfontsize=16, ytickfontsize=16,
        legendfontsize=16, axis=:equal)

   speed_ref = speednew[1] ./ proc[1] .* proc
   plot!(proc,speed_ref, linestyle=:dashdot, linecolor=:black,
         label="linear ref.")
   =#

elseif Machine == "Frontera"
   res = postprocess("RESULT_gnu_Impl")

   result_1 = sort( res[res.Thread.==1,:],(:Proc))
   result_2 = sort( res[res.Thread.==2,:],(:Proc))
   result_4 = sort( res[res.Thread.==4,:],(:Proc))
   result_7 = sort( res[res.Thread.==7,:],(:Proc))
   result_14 = sort( res[res.Thread.==14,:],(:Proc))
   result_28 = sort( res[res.Thread.==28,:],(:Proc))
   result_56 = sort(res[res.Thread.==56,:],(:Proc))

   time_1 = copy(result_1.Time)
   time_2 = copy(result_2.Time)
   time_4 = copy(result_4.Time)
   time_7 = copy(result_7.Time)
   time_14 = copy(result_14.Time)
   time_28 = copy(result_28.Time)
   time_56 = copy(result_56.Time)

   append!(time_1,fill(NaN,size(time_56)[1]-size(time_1)[1]))
   append!(time_2,fill(NaN,size(time_56)[1]-size(time_2)[1]))
   append!(time_4,fill(NaN,size(time_56)[1]-size(time_4)[1]))
   append!(time_7,fill(NaN,size(time_56)[1]-size(time_7)[1]))
   append!(time_14,fill(NaN,size(time_56)[1]-size(time_14)[1]))
   append!(time_28,fill(NaN,size(time_56)[1]-size(time_28)[1]))
   append!(time_56,fill(NaN,size(time_56)[1]-size(time_56)[1]))

   # Calculate the number of cell updates per second for 2nd order scheme
   cell_size = sort(8^3 .* 224 .* result_56.Thread .* result_56.Proc .* 2)

   speed_1 = cell_size ./ time_1
   speed_2 = cell_size ./ time_2
   speed_4 = cell_size ./ time_4
   speed_7 = cell_size ./ time_7
   speed_14 = cell_size ./ time_14
   speed_28 = cell_size ./ time_28
   speed_56 = cell_size ./ time_56

   proc = result_56.Proc .* result_56.Thread

   speed = hcat(speed_1, speed_2, speed_4, speed_7, speed_14, speed_28, speed_56)
   speed_ref = speed_1[1] ./ proc[1] .* proc

   plt.figure(figsize=(6,5))

   ax1 = plt.subplot(111)
   plt.plot(proc, speed/1e9, "-o", markersize=3, alpha=0.6)
   plt.plot(proc, speed_ref/1e9, ":")
   plt.legend(loc="upper left",
         ("1 thread", "2 threads", "4 threads", "7 threads",
          "14 threads", "28 threads", "56 threads", "linear ref."))
   plt.grid(true)
   ax1.xaxis.set_major_locator(plt.MaxNLocator(5))
   plt.title("BATSRUS Weak Scaling, Implicit Scheme", y=1.03, fontsize=16)
   plt.xlabel("number of cores", fontsize=16)
   plt.ylabel(L"giga cell ($10^9$) updates per second", fontsize=16)
end


text:
function valid_perturbation(x; ϵ=0.0)
    n = length(x)
    p = zeros(Float64, n)
    j = rand(1:n)
    for i ∈ 1:n
        if i ≠ j
            p[i] = randn()
            if x[i] < ϵ+1e-5
                p[i] = max(0, p[i])
            elseif x[i] > 1-1e-5
                p[i] = min(0, p[i])
            end
        end
    end
    sp = sum(p)
    if ≈(sp, 0; atol=1e-4) || (x[j] < ϵ+1e-5 && sp > 0) || (x[j] > 1-1e-5 && sp < 0)
        p .= 0
        return p
    end
    p[j] = -sum(p)
    p = p ./ norm(p)
    return p
end

@testset "equilibrium conditions" begin
    sol_tol = 1e-9
    ϵ = 1e-3
    for N ∈ [2,3,4]
        for min_prob ∈ [0.0, 0.05]
            d = 3*ones(Int,N)
            cost_tensors = [randn(d...) for n ∈ 1:N]
            sol = compute_equilibrium(cost_tensors; ϵ=min_prob, convergence_tolerance=sol_tol)
            xvec = vcat(sol.x...)
            @test all(xvec .≥ min_prob)
            x = sol.x

            for n ∈ 1:N
                @test sum(sol.x[n]) ≈ 1.0
                @test all(sol.x[n] .≥ 0.0)
                true_cost = expected_cost(cost_tensors[n], xvec, sol._deriv_info.tensor_indices, sol._deriv_info.primal_indices)
                for j ∈ 1:100
                    p = valid_perturbation(sol.x[n]; ϵ=min_prob)
                    x2 = deepcopy(sol.x)
                    x2[n] += 1e-4*p
                    xvec2 = vcat(x2...)
                    other_cost = expected_cost(cost_tensors[n], xvec2, sol._deriv_info.tensor_indices, sol._deriv_info.primal_indices)
                    @test other_cost ≥ true_cost - ϵ
                end
            end
        end
    end
end


text:
module Sonic

using Reexport

@reexport using Statistical_Analysis
@reexport using Optimisation
@reexport using Performance_measurement
@reexport using Classification
@reexport using Market_risk

include("The_Hedge_Fund_Industry.jl")
using LinearAlgebra
using DataStructures
using Distributions
using ForwardDiff
using Plots

# Statistical Analysis
export VAI # Value Added Index
export μ  # MEAN
export σ  # standard deviation
export ̃x  # median
export m3  # Skewness
export m4  # Kurtosis
# COVARIANCE AND CORRELATION
export R2  # coefficient of determination = the square of the correlation
# A residual plot
export SE  # standard error
export JBTest  #Jarque-Bera Test

#  Mean-Variance Optimisation
export minVar # MININIMISE VARIANCE

# Bound Constraints
# Linear Equality Constraints
# Sharpe Ratio Maximisation

#Classification
export my_f, derivate_of_my_f

export my_f1
export fVAI, getXLData
#Performance Measurement
export VAMI #diversified managed futures
export Risk-Adjusted_Returns

# Hedge Fund Classification
export dendrogram

# Market Risk Management
export VaR,  # Value at Risk
export MVaR, # Modified Value at Risk
export ES,   # Expected Shortfall
export EVT,  # Extreme Value Theory
export BM,   # Block Maxima
export POT,  # Peaks Over Threshold

end # module


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 #
#################

# test blascopy!
function test_blascopy!{T}(A::Array{T})
    @test ndims(A) == 1
    n1 = length(A)
    d_A = CudaArray(A)
    d_B = CudaArray(T,n1)
    CUBLAS.blascopy!(n,d_A,1,d_B,1)
    B = to_host(d_B)
    @test A == B
end
test_blascopy!(Float32[1:m;])
test_blascopy!(Float64[1:m;])
test_blascopy!(Float32[1:m;]+im*Float32[1:m;])
test_blascopy!(Float64[1:m;]+im*Float64[1:m;])

# test scal!
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_approx_eq(alpha*A,A1)

    d_A = CudaArray(A)
    d_As = scale(d_A, alpha)
    A1 = to_host(d_As)
    @test_approx_eq(alpha*A,A1)

    scale!(d_A, alpha)
    A1 = to_host(d_As)
    @test_approx_eq(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;])

# test dot
function test_dot(A,B)
    @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_approx_eq(cuda_dot1,host_dot)
    @test_approx_eq(cuda_dot2,host_dot)

    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
test_dot(Float32[1:m;],Float32[1:m;])
test_dot(Float64[1:m;],Float64[1:m;])

# test dotu
function test_dotu(A,B)
    @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
    @test_approx_eq(cuda_dot1,host_dot)
    @test_approx_eq(cuda_dot2,host_dot)
    @test_approx_eq(d_A.'*d_B, host_dot)
end
test_dotu(rand(Complex64,m),rand(Complex64,m))
test_dotu(rand(Complex128,m),rand(Complex128,m))

# test dotc
function test_dotc(A,B)
    @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
    @test_approx_eq(cuda_dot1,host_dot)
    @test_approx_eq(cuda_dot2,host_dot)
    @test_approx_eq(d_A'*d_B, host_dot)
end
test_dotc(rand(Complex64,m),rand(Complex64,m))
test_dotc(rand(Complex128,m),rand(Complex128,m))

# test nrm2
function test_nrm2(A)
    @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_approx_eq(cuda_nrm2_1,host_nrm2)
    @test_approx_eq(cuda_nrm2_2,host_nrm2)
    @test_approx_eq(cuda_nrm2_3,host_nrm2)
end
test_nrm2(rand(Float32,m))
test_nrm2(rand(Float64,m))
test_nrm2(rand(Complex64,m))
test_nrm2(rand(Complex128,m))

# test asum
function test_asum(A)
    @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_approx_eq(cuda_asum1,host_asum)
    @test_approx_eq(cuda_asum2,host_asum)
end
test_asum(Float32[1:m;])
test_asum(Float64[1:m;])
test_asum(rand(Complex64,m))
test_asum(rand(Complex128,m))

# 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_approx_eq(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_approx_eq(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))

# test iamax & iamin
function test_iamaxmin(A)
    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
test_iamaxmin(rand(Float32,m))
test_iamaxmin(rand(Float64,m))
test_iamaxmin(rand(Complex64,m))
test_iamaxmin(rand(Complex128,m))

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

# test gemv!
function test_gemv!(elty)
    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_approx_eq(y,h_y)
    A_mul_B!(d_y,d_A,d_x)
    h_y = to_host(d_y)
    @test_approx_eq(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_approx_eq(x,h_x)
    At_mul_B!(d_x,d_A,d_y)
    h_x = to_host(d_x)
    @test_approx_eq(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_approx_eq(x,h_x)
    Ac_mul_B!(d_x,d_A,d_y)
    h_x = to_host(d_x)
    @test_approx_eq(h_x,A'*y)
end
test_gemv!(Float32)
test_gemv!(Float64)
test_gemv!(Complex64)
test_gemv!(Complex128)

# test gemv
function test_gemv(elty)
    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_approx_eq(y1,h_y1)
    @test_approx_eq(y2,h_y2)
    @test_approx_eq(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_approx_eq(x1,h_x1)
    @test_approx_eq(x2,h_x2)
    @test_approx_eq(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_approx_eq(x1,h_x1)
    @test_approx_eq(x2,h_x2)
    @test_approx_eq(x2, to_host(d_A' * d_y))
end
test_gemv(Float32)
test_gemv(Float64)
test_gemv(Complex64)
test_gemv(Complex128)

##############
# test gbmv! #
##############

function test_gbmv!(elty)
    # 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_approx_eq(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_approx_eq(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_approx_eq(x,h_x)
end
test_gbmv!(Float32)
test_gbmv!(Float64)
test_gbmv!(Complex64)
test_gbmv!(Complex128)

# TODO: complete this test when BLAS.gbmv is fixed
function test_gbmv(elty)
    # 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_approx_eq(y,h_y)
end
test_gbmv(Float32)
test_gbmv(Float64)
test_gbmv(Complex64)
test_gbmv(Complex128)

#############
# test symv #
#############

function test_symv!(elty)
    # 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_approx_eq(y,h_y)
end
test_symv!(Float32)
test_symv!(Float64)
test_symv!(Complex64)
test_symv!(Complex128)

function test_symv(elty)
    # 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_approx_eq(y,h_y)
end
test_symv(Float32)
test_symv(Float64)
test_symv(Complex64)
test_symv(Complex128)

##############
# test hemv! #
##############

function test_hemv!(elty)
    # 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_approx_eq(y,h_y)
end
# BLAS.hemv! is not in julia v0.3.0
test_hemv!(Complex64)
test_hemv!(Complex128)

function test_hemv(elty)
    # 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_approx_eq(y,h_y)
end
test_hemv(Complex64)
test_hemv(Complex128)

##############
# test sbmv! #
##############

function test_sbmv!(elty)
    # 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_approx_eq(y,h_y)
end
test_sbmv!(Float32)
test_sbmv!(Float64)

function test_sbmv(elty)
    # 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_approx_eq(y,h_y)
end
test_sbmv(Float32)
test_sbmv(Float64)

##############
# test hbmv! #
##############

function test_hbmv!(elty)
    # 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_approx_eq(y,h_y)
end
test_hbmv!(Complex64)
test_hbmv!(Complex128)

function test_hbmv(elty)
    # 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_approx_eq(y,h_y)
end
test_hbmv(Complex128)
test_hbmv(Complex64)

##############
# test tbmv! #
##############

function test_tbmv!(elty)
    # 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_approx_eq(x,h_x)
end
test_tbmv!(Float64)
test_tbmv!(Float32)
test_tbmv!(Complex64)
test_tbmv!(Complex128)

function test_tbmv(elty)
    # 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_approx_eq(y,h_y)
end
test_tbmv(Float64)
test_tbmv(Float32)
test_tbmv(Complex64)
test_tbmv(Complex128)

##############
# test tbsv! #
##############

function test_tbsv!(elty)
    # 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_approx_eq(x,h_x)
end
test_tbsv!(Float32)
test_tbsv!(Float64)
test_tbsv!(Complex64)
test_tbsv!(Complex128)

function test_tbsv(elty)
    # 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_approx_eq(y,h_y)
end
test_tbsv(Float32)
test_tbsv(Float64)
test_tbsv(Complex64)
test_tbsv(Complex128)

##############
# test trmv! #
##############

function test_trmv!(elty)
    # 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_approx_eq(x,h_x)
end
test_trmv!(Float32)
test_trmv!(Float64)
test_trmv!(Complex64)
test_trmv!(Complex128)

function test_trmv(elty)
    # 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_approx_eq(y,h_y)
end
test_trmv(Float32)
test_trmv(Float64)
test_trmv(Complex64)
test_trmv(Complex128)

##############
# test trsv! #
##############

function test_trsv!(elty)
    # 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_approx_eq(x,h_x)
end
test_trsv!(Float32)
test_trsv!(Float64)
test_trsv!(Complex64)
test_trsv!(Complex128)

function test_trsv(elty)
    # 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_approx_eq(y,h_y)
end
test_trsv(Float32)
test_trsv(Float64)
test_trsv(Complex64)
test_trsv(Complex128)

#############
# test ger! #
#############

function test_ger!(elty)
    # 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_approx_eq(A,h_A)
end
test_ger!(Float32)
test_ger!(Float64)
test_ger!(Complex64)
test_ger!(Complex128)

#############
# test syr! #
#############

function test_syr!(elty)
    # 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_approx_eq(A,h_A)
end
test_syr!(Float32)
test_syr!(Float64)
test_syr!(Complex64)
test_syr!(Complex128)

#############
# test her! #
#############

function test_her!(elty)
    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_approx_eq(A,h_A)
end
test_her!(Complex64)
test_her!(Complex128)

##############
# test her2! #
##############

function test_her2!(elty)
    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_approx_eq(A,h_A)
end
test_her2!(Complex64)
test_her2!(Complex128)

##############
# test gemm! #
##############

function test_gemm!(elty)
    # 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_approx_eq(C1,h_C1)
    @test_approx_eq(C2,h_C2)

end
test_gemm!(Float32)
test_gemm!(Float64)
test_gemm!(Complex64)
test_gemm!(Complex128)

function test_gemm(elty)
    # 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(d_C)
    @test_approx_eq(C,h_C)
    @test_approx_eq(C,h_C2)
end
test_gemm(Float32)
test_gemm(Float64)
test_gemm(Complex64)
test_gemm(Complex128)

######################
# test gemm_batched! #
######################

function test_gemm_batched!(elty)
    # 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_approx_eq(C[i],h_C)
    end
end
test_gemm_batched!(Float32)
test_gemm_batched!(Float64)
test_gemm_batched!(Complex64)
test_gemm_batched!(Complex128)

function test_gemm_batched(elty)
    # 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_approx_eq(C,h_C)
    end
end
test_gemm_batched(Float32)
test_gemm_batched(Float64)
test_gemm_batched(Complex64)
test_gemm_batched(Complex128)

##############
# test symm! #
##############

function test_symm!(elty)
    # 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_approx_eq(C,h_C)
end
test_symm!(Float32)
test_symm!(Float64)
test_symm!(Complex64)
test_symm!(Complex128)

function test_symm(elty)
    # 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_approx_eq(C,h_C)
end
test_symm(Float32)
test_symm(Float64)
test_symm(Complex64)
test_symm(Complex128)

##############
# test syrk! #
##############

function test_syrk!(elty)
    # 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_approx_eq(C,h_C)
end
test_syrk!(Float32)
test_syrk!(Float64)
test_syrk!(Complex64)
test_syrk!(Complex128)

function test_syrk(elty)
    # 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_approx_eq(C,h_C)
end
test_syrk(Float32)
test_syrk(Float64)
test_syrk(Complex64)
test_syrk(Complex128)

##############
# test herk! #
##############

function test_herk!(elty)
    # 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.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_approx_eq(C,h_C)
end
test_herk!(Complex64)
test_herk!(Complex128)

function test_herk(elty)
    # 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_approx_eq(C,h_C)
end
test_herk(Complex64)
test_herk(Complex128)

###############
# test syr2k! #
###############

function test_syr2k!(elty)
    #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)

    #println("elty: $elty")
    #println("C:")
    #show(C)
    #println()
    #println("h_C:")
    #show(h_C)
    #println()

    @test_approx_eq(C,h_C)
end
test_syr2k!(Float32)
test_syr2k!(Float64)
test_syr2k!(Complex64)
test_syr2k!(Complex128)

function test_syr2k(elty)
    # 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_approx_eq(C,h_C)
end
test_syr2k(Float32)
test_syr2k(Float64)
test_syr2k(Complex64)
test_syr2k(Complex128)

###############
# test her2k! #
###############

function test_her2k!(elty1,elty2)
    # 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_approx_eq(C,h_C)
end
test_her2k!(Complex64,Float32)
test_her2k!(Complex128,Float64)

function test_her2k(elty1)
    # generate matrices
    A = rand(elty1,m,k)
    B = rand(elty1,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_approx_eq(C,h_C)
end
test_her2k(Complex64)
test_her2k(Complex128)

##############
# test trmm! #
##############

function test_trmm!(elty)
    # 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_approx_eq(C,h_C)
end
test_trmm!(Float32)
test_trmm!(Float64)
test_trmm!(Complex64)
test_trmm!(Complex128)

function test_trmm(elty)
    # 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_approx_eq(C,h_C)
end
test_trmm(Float32)
test_trmm(Float64)
test_trmm(Complex64)
test_trmm(Complex128)

##############
# test trsm! #
##############

function test_trsm!(elty)
    # 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
    B = alpha*(A\B)
    CUBLAS.trsm!('L','U','N','N',alpha,d_A,d_B)
    # move to host and compare
    h_B = to_host(d_B)
    @test_approx_eq(B,h_B)
end
test_trsm!(Float32)
test_trsm!(Float64)
test_trsm!(Complex64)
test_trsm!(Complex128)

function test_trsm(elty)
    # 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_approx_eq(C,h_C)
end
test_trsm(Float32)
test_trsm(Float64)
test_trsm(Complex64)
test_trsm(Complex128)

######################
# test trsm_batched! #
######################

function test_trsm_batched!(elty)
    # generate parameter
    alpha = rand(elty)
    # generate matrices
    A = [rand(elty,m,m) for i in 1:10]
    map!((x) -> triu(x), 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)
        B = alpha*(A\B)
        h_B = to_host(d_B[i])
        #compare
        @test_approx_eq(B,h_B)
    end
end
test_trsm!(Float32)
test_trsm!(Float64)
test_trsm!(Complex64)
test_trsm!(Complex128)

function test_trsm(elty)
    # generate parameter
    alpha = rand(elty)
    # generate matrices
    A = [rand(elty,m,m) for i in 1:10]
    map!((x) -> triu(x), 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_approx_eq(C,h_C)
    end
end
test_trsm(Float32)
test_trsm(Float64)
test_trsm(Complex64)
test_trsm(Complex128)

##############
# test hemm! #
##############

function test_hemm!(elty)
    # 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_approx_eq(C,h_C)
end

test_hemm!(Complex64)
test_hemm!(Complex128)

function test_hemm(elty)
    # 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_approx_eq(C,h_C)
end
test_hemm(Complex64)
test_hemm(Complex128)

##############
# test geam! #
##############

function test_geam!(elty)
    # 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_approx_eq(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_approx_eq(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_approx_eq(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_approx_eq(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

test_geam!(Float32)
test_geam!(Float64)
test_geam!(Complex64)
test_geam!(Complex128)

function test_geam(elty)
    # 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_approx_eq(C,h_C)
end

test_geam(Float32)
test_geam(Float64)
test_geam(Complex64)
test_geam(Complex128)

######################
# test getrf_batched #
######################

function test_getrf_batched!(elty)
    # 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_approx_eq(C[:L],dL)
        @test_approx_eq(C[:U],dU)
    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_approx_eq(inv(P)*dL*dU, inv(C[:P]) * C[:L] * C[:U])
    end
end
test_getrf_batched!(Float32)
test_getrf_batched!(Float64)
test_getrf_batched!(Complex64)
test_getrf_batched!(Complex128)

function test_getrf_batched(elty)
    # 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_approx_eq(C[:L],dL)
        @test_approx_eq(C[:U],dU)
    end
end
test_getrf_batched(Float32)
test_getrf_batched(Float64)
test_getrf_batched(Complex64)
test_getrf_batched(Complex128)

######################
# test getri_batched #
######################

function test_getri_batched(elty)
    # 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_approx_eq(C,h_C)
    end
end
test_getri_batched(Float32)
test_getri_batched(Float64)
test_getri_batched(Complex64)
test_getri_batched(Complex128)

#######################
# test matinv_batched #
#######################

function test_matinv_batched(elty)
    # 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_approx_eq(C,h_C)
    end
end
test_matinv_batched(Float32)
test_matinv_batched(Float64)
test_matinv_batched(Complex64)
test_matinv_batched(Complex128)

######################
# test geqrf_batched #
######################

function test_geqrf_batched!(elty)
    # 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_approx_eq(Q,full(C[:Q]))
    end
end
test_geqrf_batched!(Float32)
test_geqrf_batched!(Float64)
test_geqrf_batched!(Complex64)
test_geqrf_batched!(Complex128)

function test_geqrf_batched(elty)
    # 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_approx_eq(Q,full(C[:Q]))
    end
end
test_geqrf_batched(Float32)
test_geqrf_batched(Float64)
test_geqrf_batched(Complex64)
test_geqrf_batched(Complex128)

#####################
# test gels_batched #
#####################

function test_gels_batched!(elty)
    # 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_approx_eq(X,h_C)
    end
end
test_gels_batched!(Float32)
test_gels_batched!(Float64)
test_gels_batched!(Complex64)
test_gels_batched!(Complex128)

function test_gels_batched(elty)
    # 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_approx_eq(X,h_D)
    end
end
test_gels_batched(Float32)
test_gels_batched(Float64)
test_gels_batched(Complex64)
test_gels_batched(Complex128)

##############
# test dgmm! #
##############

function test_dgmm!(elty)
    # 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_approx_eq(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

test_dgmm!(Float32)
test_dgmm!(Float64)
test_dgmm!(Complex64)
test_dgmm!(Complex128)

function test_dgmm(elty)
    # 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_approx_eq(C,h_C)
end

test_dgmm(Float32)
test_dgmm(Float64)
test_dgmm(Complex64)
test_dgmm(Complex128)


text:
# Type that keeps mode coupling analysis data
# 2014 murozturk@ku.edu.tr

include("helper.jl")

# fields are kept in the order they are used in analysis
immutable mcAnalysis
	md::protein
	R2r::Matrix{Float64}
	r::Matrix{Float64}
	coef::Matrix{Float64}
	rf1::Matrix{Float64}
	r2R::Matrix{Float64}
	Rf1::Matrix{Float64}
	kl::Vector{Float64}
	js::Vector{Float64}
end

# Constructor performs mode coupling analysis for a protein
# takes md coordinates R as input
# ith cloumn is all samples for mode i
# nModes = 3* nResidues

function mcAnalysis(R::Matrix{Float64}, index=1:size(R,2))

	nSamples, nModes = size(R)
	println("Beginning mode coupling analysis")
	println("nSamples, nModes :", nSamples, " ", nModes)

	# convert R to onur's convention (for debugging purposes)
	#R= [ R[:,1:3:end] R[:,2:3:end] R[:,3:3:end] ]

	md = protein(R, index)

	R2r, r2R = modeSpaceTransform (md)

	println ("Transforming into modespace...")

	r=(R2r*md.R')'

	r=removeRows(r, findOutliers(r))
	nSamples = size(r,1)

	# calculating hermite coefficients
	coef = coefHermit(r)

	# sampling f1 distribution
	rf1 = sampleHermit(nSamples, coef)

	println ("Transforming back into cartesian space...")

	Rf1 = (r2R*rf1')'

	println ("Calculating Kullback–Leibler distances...")
	kl = perResidue( klDist(md.R, Rf1))

	println ("Calculating Jensen-Shannon distances...")
	js = perResidue( jsDist(md.R, Rf1))

	mcAnalysis(md, R2r, r, coef, rf1, r2R, Rf1, kl, js )
end

# high level recall analysis that also handles indexing and normalizing
function recall(mcA::mcAnalysis, targetSet::Vector{Int})
	normalizedJs= mcA.js ./ mcA.md.stdPerRes
	predictPerm = sortperm(normalizedJs, rev=true)
	recall(mcA.md.index[predictPerm], targetSet  )
end

println ("Making dummy run for JIT compiler...")
tmp = mcAnalysis(randn(9,9));
recall(tmp,[1])
println ("Ready to roll.")


text:
include("spectral.jl")

function compute_accuracy(nearest, k)
    count = zeros(k, k)
    for (i, c_algo) in nearest
        c_truth = i <= 100 ? 1 : 2
        count[c_truth, c_algo] += 1
    end
    count
end

function trial(x, cluster_num; neighbor_num=10, sigma=1)
    nearest = spectral_clustering(x, cluster_num; neighbor_num=neighbor_num, sigma=sigma)
    count = compute_accuracy(nearest, cluster_num)

    println("k=$neighbor_num\tsigma=$sigma\t$(collect(count))")
end

function main()
    x = include("spiral_data.jl");
    cluster_num  = 2;

    # fix sigma = 1, neighbor_num varies
    sigma = 1
    neighbor_nums = [5, 10, 20, 30, 50, 200]
    for neighbor_num in neighbor_nums
        trial(x, cluster_num, neighbor_num=neighbor_num, sigma=sigma)
    end

    # fix neighbor_num = 10/20/30, sigma varies
    neighbor_num = 30
    sigmas = [0.01, 0.1, 0.5, 1, 2, 5, 10, 15]
    for sigma in sigmas
        trial(x, cluster_num, neighbor_num=neighbor_num, sigma=sigma)
    end

end

main()

max_stars_repo_name: zxteloiv/julia-ex

max_issues_repo_name: zxteloiv/julia-ex

max_forks_repo_name: zxteloiv/julia-ex