text:
# adapted from https://nifti.nimh.nih.gov/nifti-1/data
@testset "Orientations" begin
    @testset "Rotation matrices" begin
        R = AffineMap(RotMatrix{3,Float64}(reshape(1:9,3,3)), LinearMap((1.0, 1.0, 1.0)))
        @test NeuroCore.rotation2spatialorder(R.linear) == (:inferior, :posterior, :left)

        R = AffineMap(
            RotMatrix{3}([
                -2.0                    6.714715653593746e-19  9.081024511081715e-18
                6.714715653593746e-19  1.9737114906311035    -0.35552823543548584
                8.25548088896093e-18   0.3232076168060303     2.171081781387329
            ]),
            LinearMap((117.8551025390625, -35.72294235229492, -7.248798370361328))
        )

        @test NeuroCore.rotation2spatialorder(R.linear) == (:right, :posterior, :inferior)

        # TODO This is off and needs testing
        #Q = NeuroCore.mat2quat(R, (2.0, 2.0, 2.1999990940093994))
        #R2 = NeuroCore.quat2mat(Q, (2.0, 2.0, 2.1999990940093994))

        R = RotMatrix{3}([
            -1.0  0.0  0.0
            0.0  1.0  0.0
            0.0  0.0  1.0])

        @test NeuroCore.rotation2spatialorder(R) == (:right, :posterior, :inferior)
        @test NeuroCore.spatialorder2rotation((:right, :posterior, :inferior)) == R

        R = RotMatrix{3}([
            -1.0  0.0  0.0
            0.0  1.0  0.0
            0.0  0.0  1.0])
        @test NeuroCore.rotation2spatialorder(R) == (:right, :posterior, :inferior)
        @test NeuroCore.spatialorder2rotation((:right, :posterior, :inferior)) == R

        R = AffineMap(
                RotMatrix{3}([
                    1.000000 0.000000 0.000000
                    0.000000 -1.000000 0.000000
                    0.000000 0.000000 1.000000]),
                (1.0, 1.0, 1.0))

        @test NeuroCore.rotation2spatialorder(R) == (:left, :anterior, :inferior)
        @test NeuroCore.spatialorder2rotation((:left, :anterior, :inferior)) == R.linear

        R = AffineMap(
            RotMatrix{3}([
                1.000000 0.000000 0.000000
                0.000000 -1.000000 0.000000
                0.000000 0.000000 -1.000000]),
            LinearMap((-90.0, -126.0, -72.0))
       )

        @test NeuroCore.rotation2spatialorder(R) == (:left, :anterior, :superior)
        @test NeuroCore.spatialorder2rotation((:left, :anterior, :superior)) == R.linear

        @test NeuroCore.rotation2spatialorder(affine_map(rand(2,2,2))) ==
              NeuroCore.rotation2spatialorder(affine_map(rand(2,2))) ==
              (:left, :posterior, :inferior)
    end

    @testset "LR MNI152" begin
        Q = AffineMap(SPQuat(0.0, 0.0, 0.0),
                  ntuple(_ -> 2.0, 3))
        @test NeuroCore.rotation2spatialorder(Q) == (:left, :posterior, :inferior)

        S = AffineMap(
            RotMatrix{3}([
                -2.000000 0.000000 0.000000
                 0.000000 2.000000 0.000000
                 0.000000 0.000000 2.000000]),
            (90.0, -126.0, -72.0)
        )
        @test NeuroCore.rotation2spatialorder(S) == (:right, :posterior, :inferior)
    end

    @testset "RL MNI152" begin
        R = AffineMap(
            RotMatrix{3}([
                2.000000 0.000000 0.000000
                0.000000 2.000000 0.000000
                0.000000 0.000000 2.000000]),
            (-90.0, -126.0, -72.0)
        )
        @test NeuroCore.rotation2spatialorder(R) == (:left, :posterior, :inferior)
    end

    @testset "Z stat" begin
        Q = AffineMap(SPQuat(0.0, 0.0, 0.0), (4.0, 4.0, 6.0))
        @test NeuroCore.rotation2spatialorder(Q) == (:left, :posterior, :inferior)
    end
end

@testset "quaternion <--> matrix" begin
    R = affine_map(rand(2,2,2))
    px = pixelspacing(rand(2,2,2))
    Q = NeuroCore.mat2quat(R, px)
    R2 = NeuroCore.quat2mat(Q, px)
    @test R.linear == R2.linear

    R = AffineMap(
        RotMatrix{3}([0.14 0.0 0.0
                      0.0  0.9 0.0
                      0.0  0.0 1.1]),
        LinearMap((117.8551025390625, -35.72294235229492, -7.248798370361328))
    )
    R1 = NeuroCore.quat2mat(NeuroCore.mat2quat(R, px), px)
    @test R1 == NeuroCore.quat2mat(NeuroCore.mat2quat(R1, px), px)
end

@testset "2D orientation" begin
    Q = NeuroCore.mat2quat(rand(2,2))
    @test Q.linear.x == 0
    @test Q.linear.y == 0
    @test Q.linear.z == 0
end

@testset "Orientation errors" begin
    R = RotMatrix{3}([
        0.000000 0.000000 0.000000
        0.000000 -1.000000 0.000000
        0.000000 0.000000 -1.000000])
    @test_throws ErrorException NeuroCore.rotation2spatialorder(R)

    R = RotMatrix{3}([
        1.000000 0.000000 0.000000
        0.000000 0.000000 0.000000
        0.000000 0.000000 -1.000000])

    @test_throws ErrorException NeuroCore.rotation2spatialorder(R)
    R = RotMatrix{3}([
                -2.0                    6.714715653593746e-19  9.081024511081715e-18
                6.714715653593746e-19  1.9737114906311035    -0.35552823543548584
                8.25548088896093e-18   0.3232076168060303     2.171081781387329
            ])

    @test_throws ErrorException NeuroCore.number2dimname(9)

    @test_throws ErrorException NeuroCore.dimname2number(:foo)
end


text:
#=

  Eq10 problem in Julia ConstraintSolver.jl

  Standard benchmark problem.

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function eq20(print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))
    n = 7
    @variable(model, 0 <= x[1:n] <= 10, Int)
    x0,x1,x2,x3,x4,x5,x6 = x

    @constraint(model, -76706*x0 + 98205*x1 + 23445*x2 + 67921*x3 + 24111*x4 +
    -48614*x5 + -41906*x6 == 821228)
    @constraint(model, 87059*x0 + -29101*x1 + -5513*x2 + -21219*x3 + 22128*x4 +
    7276*x5 + 57308*x6 == 22167)
    @constraint(model, -60113*x0 + 29475*x1 + 34421*x2 + -76870*x3 + 62646*x4 +
    29278*x5 + -15212*x6 == 251591)
    @constraint(model, 49149*x0 + 52871*x1 + -7132*x2 + 56728*x3 + -33576*x4 +
    -49530*x5 + -62089*x6 == 146074)
    @constraint(model, -10343*x0 + 87758*x1 + -11782*x2 + 19346*x3 + 70072*x4 +
    -36991*x5 + 44529*x6 == 740061)
    @constraint(model, 85176*x0 + -95332*x1 + -1268*x2 + 57898*x3 + 15883*x4 +
    50547*x5 + 83287*x6 == 373854)
    @constraint(model, -85698*x0 + 29958*x1 + 57308*x2 + 48789*x3 + -78219*x4 +
    4657*x5 + 34539*x6 == 249912)
    @constraint(model, -67456*x0 + 84750*x1 + -51553*x2 + 21239*x3 + 81675*x4 +
    -99395*x5 + -4254*x6 == 277271)
    @constraint(model, 94016*x0 + -82071*x1 + 35961*x2 + 66597*x3 + -30705*x4 +
    -44404*x5 + -38304*x6 == 25334)
    @constraint(model, -60301*x0 + 31227*x1 + 93951*x2 + 73889*x3 + 81526*x4 +
    -72702*x5 + 68026*x6 == 1410723)
    @constraint(model, -16835*x0 + 47385*x1 + 97715*x2 + -12640*x3 + 69028*x4 +
    76212*x5 + -81102*x6 == 1244857)
    @constraint(model, -43277*x0 + 43525*x1 + 92298*x2 + 58630*x3 + 92590*x4 +
    -9372*x5 + -60227*x6 == 1503588)
    @constraint(model, -64919*x0 + 80460*x1 + 90840*x2 + -59624*x3 + -75542*x4 +
    25145*x5 + -47935*x6 == 18465)
    @constraint(model, -45086*x0 + 51830*x1 + -4578*x2 + 96120*x3 + 21231*x4 +
    97919*x5 + 65651*x6 == 1198280)
    @constraint(model, 85268*x0 + 54180*x1 + -18810*x2 + -48219*x3 + 6013*x4 +
    78169*x5 + -79785*x6 == 90614)
    @constraint(model, 8874*x0 + -58412*x1 + 73947*x2 + 17147*x3 + 62335*x4 +
    16005*x5 + 8632*x6 == 752447)
    @constraint(model, 71202*x0 + -11119*x1 + 73017*x2 + -38875*x3 + -14413*x4 +
    -29234*x5 + 72370*x6 == 129768)
    @constraint(model, 1671*x0 + -34121*x1 + 10763*x2 + 80609*x3 + 42532*x4 +
    93520*x5 + -33488*x6 == 915683)
    @constraint(model, 51637*x0 + 67761*x1 + 95951*x2 + 3834*x3 + -96722*x4 +
    59190*x5 + 15280*x6 == 533909)
    @constraint(model, -16105*x0 + 62397*x1 + -6704*x2 + 43340*x3 + 95100*x4 +
    -68610*x5 + 58301*x6 == 876370)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                println("x:$x_val")

            end
        end
    else
        println("status:$status")
    end

    return status
end

@time eq20()

ext: jl

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
abstract type List{T} end

struct EmptyList{T} <: List{T} end

struct ConsList{T} <: List{T}
    head::T
    tail::List{T}
end

ConsList{T}(a::T) where {T} = ConsList{T}(a, EmptyList{T}())

Base.:+(head::T, tail::List{T}) where {T} = ConsList{T}(head, tail)

Base.:+(a::List{T}, b::T) where {T} = a + ConsList{T}(b)

Base.:+(a::EmptyList{T}, b::List{T}) where {T} = b

Base.:+(a::ConsList{T}, b::List{T}) where {T} = ConsList{T}(a.head, a.tail + b)

rev(a::EmptyList{T}) where {T} = a

rev(a::ConsList{T}) where {T} = rev(a.tail) + ConsList{T}(a.head)

function Base.convert(::Type{Vector{T}}, a::EmptyList{T}) where {T}
    Vector{T}()
end

function Base.convert(::Type{Vector{T}}, a::ConsList{T}) where {T}
    unshift!(convert(Vector{T}, a.tail), a.head)
end

function Base.convert(::Type{Vector{Vector{T}}}, a::List{List{T}}) where {T}
    map((x -> convert(Vector{T}, x)), convert(Vector{List{T}}, a))
end

function zip2(a::EmptyList{T}, b::EmptyList{S}) where {T} where {S}
    EmptyList{Tuple{T, S}}()
end

function zip2(a::EmptyList{T}, b::ConsList{S}) where {T} where {S}
    EmptyList{Tuple{T, S}}()
end

function zip2(a::ConsList{T}, b::EmptyList{S}) where {T} where {S}
    EmptyList{Tuple{T, S}}()
end

function zip2(a::ConsList{T}, b::ConsList{S}) where {T} where {S}
    ConsList{Tuple{T, S}}((a.head, b.head), zip2(a.tail, b.tail))
end

function Base.map(f::Function, a::EmptyList{T}, ::Type{S}) where {T} where {S}
    EmptyList{S}()
end

function Base.map(f::Function, a::ConsList{T}, ::Type{S}) where {T} where {S}
    ConsList{S}(f(a.head), map(f, a.tail, S))
end

max_stars_repo_name: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
### Model after Gab's code - Create basic CSV files

using CSV
using DataFrames

## Canis latrans observations
# change directory
## $(homedir()) to get it in cd(~[tab])
## "../" brings you back one folder. "../../" brings you back two folders
cd("$(homedir())/Documents/UdeM/ete2019/poisot_lab/pred_prey_proj/redlist_species_data/preds")

# Select 2019 observations
canis_latrans = CSV.read("Canis_latrans.csv")
# select all that are not "NA"
c_latrans = canis_latrans[canis_latrans.year .!= "NA", :]
# convert year "string" to "integer"
c_latrans.year = parse.(Int64, c_latrans.year)
describe(c_latrans_2019)
# select all from year 2019
c_latrans_2019 = c_latrans[c_latrans.year .=== 2019, :]

#write to CSV
CSV.write("/Users/kiristern/Documents/GitHub/SDM/data/c_latrans_2019.csv", c_latrans_2019)

#test CSV file
test = CSV.read("/Users/kiristern/Documents/GitHub/SDM/data/c_latrans_2019.csv")
first(test,6)
names(test)

max_stars_repo_name: kiristern/SDM

max_issues_repo_name: kiristern/SDM

max_forks_repo_name: kiristern/SDM


text:
#=
ndgrid-oneto.jl
ndgrid for a set of "OneTo" ranges: probably the simplest possible ndgrid
=#

export ndgrid

"""
    GridOT{T,d,D} <: AbstractArray{T,D}
The `d`th component of `D`-dimensional `ndgrid(1:M, 1:N, ...)`
where `1 ≤ d ≤ D`.
"""
struct GridOT{T,d,D} <: AbstractGrid{T,d,D}
    dims::Dims{D}

    function GridOT(T::DataType, dims::Dims{D}, d::Int) where D
        1 ≤ d ≤ D || throw(ArgumentError("$d for $dims"))
        T <: Integer || throw(ArgumentError("T = $T"))
        new{T,d,D}(dims)
    end
end

Base.size(a::GridOT) = a.dims
# Base.eltype(::GridOT) = Int # default

@inline Base.@propagate_inbounds function Base.getindex(
    a::GridOT{T,d,D},
    i::Vararg{Int,D},
) where {T,d,D}
    @boundscheck checkbounds(a, i...)
    return @inbounds i[d]
end

"""
    (xg, yg, ...) = ndgrid(M, N, ...)
Shorthand for `ndgrid(1:M, 1:N, ...)`.

# Example

```jldoctest
julia> ndgrid(2,3)
([1 1 1; 2 2 2], [1 2 3; 1 2 3])
```
"""
function ndgrid(ns::Vararg{Int})
    all(>(0), ns) || throw(ArgumentError("$n ≤ 0"))
    return ndgrid(Base.OneTo.(ns)...)
end

#=
"""
    (xg, yg, ...) = ndgrid(Base.OneTo{M}, Base.OneTo{N}, ...)
Construct `ndgrid` tuple for essentially `ndgrid(1:M, 1:N, ...)`.

# Example

```jldoctest
julia> ndgrid(Base.OneTo(2), Base.OneTo(3))
([1 1 1; 2 2 2], [1 2 3; 1 2 3])
```
"""
function ndgrid(os::Vararg{Base.OneTo{Int}})
    D = length(os)
    dims = ntuple(i -> os[i].stop, D)
    return ntuple(i -> GridOT(dims, i), D)
end
=#


text:
abstract type Scalarness end
struct Scalar <: Scalarness end
struct NotScalar <: Scalarness end

isscalar(::Type{Any}) = NotScalar() # if we don't know the type we can't really know if scalar or not
isscalar(::Type{<:AbstractString}) = NotScalar() # We consider strings to be nonscalar
isscalar(::Type{<:Number}) = Scalar() # We consider Numbers to be scalar
isscalar(::Type{Char}) = Scalar() # We consider Character to be scalar
isscalar(::Type{T}) where T = hasmethod(iterate, (T,)) ? NotScalar() : Scalar()

"""
    apply_to_interior_levels(f, iter, lvls)

f is called on all levels up to `max_depth`
that have children and who's children (probably) have children.

`f` must be a function taking the current level number it's first argument,
and the level content as it's second.
That level is replaced with the return value of `f`
"""
function apply_to_interior_levels(ff, iter, lvls, max_depth)
    function inner(level_iter::T, cur_level) where T
        # This is a static Holy Trait dispatch
        this_scalar = isscalar(T)
        children_scalar = isscalar(eltype(T))
        inner(this_scalar, children_scalar, level_iter, cur_level)
    end

    # scalar, any children (Fall back if didn't manage to stop earlier)
    function inner(::Scalar,::Any, childs, cur_level)
        childs # This is just a scalare element
    end

    # not scalar, children scalar
    function inner(::NotScalar,::Scalar, childs, cur_level)
        # Only ops we care about are ones that act on nonscalar children
        childs
    end

    # not scalar, children not scalar (probably)
    function inner(::NotScalar, ::NotScalar, childs, cur_level)
        if cur_level <= max_depth
            # Only expand down if not yet at deepset level we act on
            # (This saves on allocations and time
            childs = apply(childs) do _childs
                # `apply` needs to be wrapped around all operations, including imap/map
                imap(_childs) do child
                    inner(child, cur_level+1)
                end
            end
        end

        ff(cur_level, childs)
    end

    inner(iter, 1)
end

"""
    apply_at_level(ff, iter, lvls)

Apply the given function `ff` at all levels specified in `lvls`
"""
function apply_at_level(ff, iter, lvls)
    max_depth = maximum(lvls)
    apply_to_interior_levels(iter, lvls, max_depth) do cur_level, childs
        ret = if cur_level ∈ lvls
            # If on the level where we act, then act
            apply(ff, childs)
        else
            childs
        end
        ret
    end
end

"""
    apply_at_level(iter, lvl_ops::Associative{<:Integer})

Applies a op
`ops` is a list of functions the same length as `lvls`.
`apply_pel_level` applies the given op, at the correponding level.

If only a single level and a single op is provider,
then this is identical to `apply_at_level`.
"""
function apply_at_level(iter, lvl_ops::AbstractDict{<:Integer})
    max_depth = maximum(keys(lvl_ops))
    apply_to_interior_levels(iter, lvl_ops, max_depth) do cur_level, childs
        # if level is a key then  do the thing, else no change
        if haskey(lvl_ops, cur_level)
            apply(lvl_ops[cur_level], childs)
        else
            childs
        end
    end
end


text:
module ChainedEquations
# TODO: Port the R module

import ..Concepts
import ..Concepts: MaybeMissing, type_conversion
import ..Utilities
import ..ModelFitting.ObservedOrMissing
import DataFrames

using RCall

struct MICE end

export MICE

function initialize_trackers(A::Array{MaybeMissing{Float64}}, type_assignment)
    type_tracker = Utilities.IndexTracker{Symbol}()
    if type_assignment == nothing
        Concepts.disjoint_join(type_tracker, Concepts.fit(DistributionDeduction(), A))
    else
        Concepts.disjoint_join(type_tracker, type_assignment)
    end
    Concepts.disjoint_join(type_tracker, Concepts.fit(ObservedOrMissing(), A))
    result = Concepts.groupby(type_tracker, [:Observed, :Missing])
    result[:Observed][:Total] = Base.findall(!ismissing, A)
    result[:Missing][:Total] = Base.findall(ismissing, A)
    result[:SingleView] = type_tracker
    return result
end

function Concepts.complete(model::MICE;
                           A::Array{MaybeMissing{Float64}},
                           impute_rounds = 10,
                           impute_method = "norm",
                           type_assignment)
    tracker = initialize_trackers(A, type_assignment)
    r_data = convert(DataFrames.DataFrame, A)
    r_missing = hcat(Base.map(x -> collect(x.I), tracker[:Missing][:Total])...)'
    display(r_data)
    display(r_missing)
    @rput r_data;@rput r_missing; @rput impute_rounds; @rput impute_method;
    R"""
        # print(r_missing)
        print(sessionInfo())
        library(mice)
        library(Matrix)
        data <- airquality
        data[4:10,3] <- rep(NA,7)
        data[1:5,4] <- NA
        imp.train_raw <- mice(data, m = 5, method='cart')
        # For(i in seq(1:dim(r_missing)[1]))
        # {
        #     r_data[r_missing[i,1],r_missing[i,2]] = NA
        # }
        # print(r_data)
        # imp.train_raw <- mice(r_data, m = 5, method='norm', printFlag=FALSE)
        # rImputed = as.matrix(complete(imp.train_raw))
        """
    # imputed = @rget rImputed
    # display(imputed)

    # rData = convert(Data
    # rMissing = asMatrix(findMissingEntry(size(data.sourceMatrix,1),
    #                                      size(data.sourceMatrix,2),
    #                                      data.observedEntry));
    # rImputeTimes = model.imputeTimes; rImputeMethod = model.imputeMethod;
    # # rImputeTimes = 5; rImputeMethod = "norm"
    # @rput rData;@rput rMissing;@rput rImputeTimes;@rput rImputeMethod;
    # R"""
    #     library(mice)
    #     for(i in seq(1:dim(rMissing)[1]))
    #     {
    #         rData[rMissing[i,1],rMissing[i,2]] = NA
    #     }
    #     imp.train_raw <- mice(rData, m=5, method='norm', printFlag=FALSE)
    #     rImputed = as.matrix(complete(imp.train_raw))
    #     """
    # imputed = @rget rImputed
    # ret = ImputedModel(original = data.sourceMatrix,
    #                    imputed = imputed,
    #                    missing = findMissingEntry(size(data.sourceMatrix,1),
    #                                               size(data.sourceMatrix,2),
    #                                               data.observedEntry));
    # logImputedModel(logConfig; model = ret);
    # return ret
end

end


text:
using SortingLab, SortingAlgorithms

a = rand(Float32, 1_000_000)

sortperm_sorttwo(x) = begin
	i = collect(1:length(x))
	SortingLab.sorttwo!(copy(x), i)
	i
end

@time fsortperm(a)
@time sortperm(a)
@time sortperm(a, alg=RadixSort)
@time ass = sortperm_sorttwo(a)

using BenchmarkTools
@benchmark fsortperm(a)
@benchmark sortperm(a)
@benchmark sortperm(a, alg=RadixSort)
@benchmark ass = sortperm_sorttwo(a)

# sortperm_sorttwo seems to win out!
# which is implemented in fsortperm()


text: # import TransportLayerSecurity


text:
"""
    gbrand(typeorrange, nrows, ncols, density; kwargs...)::GBMatrix
    gbrand(rng::AbstractRNG, typeorrange, nrows, ncols, density; kwargs...)::GBMatrix

Construct a random `GBMatrix`, analogous to `sprand` from SparseArrays

# Arguments
- `rng::AbstractRNG`: Random number generator for both values and indices.
- `typeorrange`: Either a type such as `Float64`, or a range such as `1:10`.
Any input which supports `eltype(typeorrange)`.
- `nrows::Integer`, `ncols::Integer`: Dimensions of the result.
- `density::AbstractFloat`: The approximate density of result.

# Keywords
- `symmetric::Bool`: The result matrix is symmetric, Aᵀ = A.
- `pattern::Bool`: The result matrix consists solely of `one(eltype(typeorrange))`.
- `skewsymmetric::Bool`: The result matrix is skew-symmetric, Aᵀ = -A.
- `hermitian::Bool`: The result matrix is hermitian, aᵢⱼ = āⱼᵢ.
- `nodiagonal::Bool`: The result matrix has no values on the diagonal.

# Returns
- `GBMatrix`
"""
function gbrand(
    rng::AbstractRNG, typeorrange, nrows::Integer, ncols::Integer, density::AbstractFloat;
    symmetric=false, pattern=false, skewsymmetric=false, hermitian=false, nodiagonal=false
)
    type = eltype(typeorrange)
    A = GBMatrix{type}(nrows, ncols)
    (type <: Unsigned || type == Bool) && (skewsymmetric = false)
    if nrows != ncols
        symmetric = false
        skewsymmetric = false
        hermitian = false
    end
    if pattern || symmetric
        skewsymmetric = false
        hermitian = false
    end
    if skewsymmetric
        hermitian = false
        nodiagonal = false
    end
    if !(type <: Complex)
        hermitian = false
    end

    for _ ∈ 1:round(Int64, nrows * ncols * density)
        i = rand(rng, 1:nrows)
        j = rand(rng, 1:ncols)
        nodiagonal && (i == j) && continue
        if pattern
            x = one(type)
        else
            x = rand(rng, typeorrange)
        end
        A[i, j] = x
        symmetric && (A[j, i] = x)
        skewsymmetric && (A[j, i] = -x)
        hermitian && (A[j, i] = conj(x))
    end

    return A
end

function gbrand(
    typeorrange, nrows::Integer, ncols::Integer, density::AbstractFloat;
    symmetric=false, pattern=false, skewsymmetric=false, hermitian=false, nodiagonal=false
)
    return gbrand(
        default_rng(), typeorrange, nrows, ncols, density;
        symmetric, pattern, skewsymmetric, hermitian, nodiagonal
    )
end


text:
# LibTOAST.jl: interface to the TOAST++ library
# Copyright (C) 2019 Samuel Powell

# Imports
import Base: string, print, show

# Export types
export Mesh

# Export methods
export delete, string, print, show, load, save,
  nodecount, elemcount, dimensions, boundingbox, fullsize, data, maxnodes, surface

# Type definitions
"""
    Mesh

A reference to a TOAST mesh, with cached sparsity structure for re-use in the
generation of, e.g., system matrices.
"""
mutable struct Mesh

  ptr::Cxx.CxxCore.CppPtr
  rowptr::Array{Cint,1}
  colidx::Array{Cint,1}

  function Mesh(ptr::Cxx.CxxCore.CppPtr)

    # Remove unused nodes
    if @cxx ptr->Shrink()
      warn("Removing unused nodes.");
    end

    # Perform initialisation
    _info("Initiliasing mesh.")
    @cxx ptr->Setup()

    # Compute sparsity pattern
    _info("Calculating sparsity pattern.")
    rowptr, colidx = _sparsity!(ptr)

    mesh = new(ptr, rowptr, colidx)
    finalizer(_mesh_delete, mesh)

    return mesh

  end

end

"""
    Mesh(vtxs, elem, eltp)

Construct a new mesh given a list of vertices, elements, and element types.

# Arguments
* `vtxs::Matrix{Float64}`: matrix of vertices (nvtx x dim)
* `elem::Matrix{Integer}`: element connectivity matrix (nels x mnnd)
* `eltp::Vector{Integer}`: a vector of element types (nels x 1)

...where `nvtx` is the number of vertices in the mesh, `dim` is the spatial
dimension, `nels` is the number of elements, `mnnd` is the maximium number of
nodes per element.
"""
function Mesh(vtxs::Matrix{Float64}, elem::Matrix{I}, eltp::Vector{I}) where {I <: Integer}
  meshptr = _mesh_new()
  _make!(meshptr,vtxs,elem,eltp)
  return Mesh(meshptr)
end

"""
    Mesh(filename)

Construct a new mesh from a specified file.
"""
function Mesh(fn::String)
  meshptr = _mesh_new()
  _load!(meshptr,fn)
  return Mesh(meshptr)
end

# Create new mesh pointer
_mesh_new() = @cxxnew TOAST_Mesh()

# Delete and deallocate mesh pointer
# TODO: Does Cxx.jl do this for me?
_mesh_delete(mesh) = finalize(mesh.ptr) # icxx"""delete $(mesh.ptr);"""

# Compute sparsity structure in zero-based CSR format and 1-based CSC.
function _sparsity!(ptr::Cxx.CxxCore.CppPtr)

  rowptr = Ptr{Int32}[0]
  colidx = Ptr{Int32}[0]
  nnzero = Int32[0]

  icxx"""
    $(ptr)->SparseRowStructure($(pointer(rowptr))[0], $(pointer(colidx))[0], $(pointer(nnzero))[0]);
  """

  nnd = @cxx ptr->nlen()

  _info("Number of vertices: $nnd, number of nonzeros: $(nnzero[1])")

  # Zero based rowptr and colidx, for resuse when building a system matrix
  rowptr = unsafe_wrap(Array, rowptr[1], nnd+1, own = true)
  colidx = unsafe_wrap(Array, colidx[1], nnzero[1], own = true)

  return rowptr, colidx

end

# Build a mesh from a matrix of vertices, elements, and element types
function _make!(meshptr::Cxx.CxxCore.CppPtr, node, elem, eltp)

  jnvtx = size(node,1)
  jnel = size(elem,1)
  jdim = size(node,2)
  jnnd0 = size(elem,2)

  # The following code fragment is modified from mtMesh.cc
  icxx"""
      int i, j, k;
      int nvtx = $jnvtx;
      int nel  = $jnel;
      int dim  = $jdim;
      int nnd0 = $jnnd0;
      double *vtx = $(pointer(node));
      int *idx = $(pointer(elem));
      int *etp = $(pointer(eltp));

      Mesh *mesh = $(meshptr);

      // Create node list
      mesh->nlist.New (nvtx);
      for (i = 0; i < nvtx; i++) {
        mesh->nlist[i].New(dim);
        mesh->nlist[i].SetBndTp (BND_NONE); // Unknown (yet)
      }
      for (j = k = 0; j < dim; j++) {
        for (i = 0; i < nvtx; i++) {
          mesh->nlist[i][j] = vtx[k++];
        }
      }

      // Create element list
      Element *el, **list = new Element*[nel];
      for (i = 0; i < nel; i++) {
        int eltp = (int)(etp[i]+0.5);
        switch (eltp) {
          case ELID_TRI3OLD:
            list[i] = new Triangle3old;
            break;
          case ELID_TET4:
            list[i] = new Tetrahedron4;
            break;
          case ELID_WDG6:
            list[i] = new Wedge6;
            break;
          case ELID_VOX8:
            list[i] = new Voxel8;
            break;
          case ELID_TRI6:
            list[i] = new Triangle6;
            break;
          case ELID_TET10:
            list[i] = new Tetrahedron10;
            break;
          case ELID_TRI6_IP:
            list[i] = new Triangle6_ip;
            break;
          case ELID_TRI10:
            list[i] = new Triangle10;
            break;
          case ELID_TRI10_IP:
            list[i] = new Triangle10_ip;
            break;
          case ELID_TET10_IP:
            list[i] = new Tetrahedron10_ip;
            break;
          case ELID_PIX4:
            list[i] = new Pixel4;
            break;
          case ELID_TRI3:
            list[i] = new Triangle3;
            break;
          case ELID_TRI3D3:
            list[i] = new Triangle3D3;
            break;
          case ELID_TRI3D6:
            list[i] = new Triangle3D6;
            break;
          default:
            printf("Element type not supported!\n");
            list[i] = 0;
            break;
        }
      }
      mesh->elist.SetList (nel, list);
      delete []list;

      for (j = k = 0; j < nnd0; j++) {
        for (i = 0; i < nel; i++) {
          if ((el = mesh->elist[i])) {
            if (j < el->nNode())
              el->Node[j] = (int)(idx[k]-0.5);
            }
          k++;
        }
      }
    """

  return nothing

end

# Method definitions
# data
# surf
# maxnodes

# String conversion
string(mesh::Mesh) = ("")

# Print method
print(io::IO, mesh::Mesh) = print(io, string(mesh))

# Show method
show(io::IO, mesh::Mesh) = print(io, mesh)

# Load a mesh from the specified filename
function _load!(ptr::Cxx.CxxCore.CppPtr, fn::String)

  # TODO: We would like to do something like this:
  # ifs = @cxx std::ifstream(pointer(filename))
  # icxx"""$(ifs) >> *$(obj.meshptr);"""

  # @show ptr
  # Check if file exists
  if !isfile(fn)
    error("File not found ", fn)
  end

  # Load the mesh
  icxx"""
    ifstream ifs($(pointer(fn)));
    ifs >> *($(ptr));
    ifs.close();
  """

  _info("Loaded $(stat(fn).size) bytes from $(fn).")

  return

end

"""
    save(mesh, filename)

Save mesh to the specified filename.
"""
function save(mesh::Mesh, fn::String)

  # Output the mesh
  icxx"""
    ofstream ofs($(pointer(fn)));
    ofs << *($(mesh.ptr));
    ofs.close();
  """

  _info("Mesh written to $(fn), $(stat(fn).size) bytes.")

  return nothing

end

"""
    nodecount(mesh)

Return the number of nodes in the mesh.
"""
nodecount(mesh::Mesh) = @cxx mesh.ptr->nlen()

"""
    elemcount(mesh)

Return the number of elements in the mesh.
"""
elemcount(mesh::Mesh) = @cxx mesh.ptr->elen()

"""
    dimensions(mesh)

Return the number of spatial dimensions of the mesh.
"""
dimensions(mesh::Mesh) = @cxx mesh.ptr->Dimension()

"""
    maxnodes(mesh)

Return the maximum number of nodes in any element of the mesh.
"""
function maxnodes(mesh::Mesh)

  meshptr = mesh.ptr

  mnnd = icxx"""
    int elen = $meshptr->elen();
    int mnnd = $meshptr->elist[0]->nNode();
    for (int i = 0; i < elen; i++)
      mnnd = ::max (mnnd, $meshptr->elist[i]->nNode());
    return mnnd;
  """

  return mnnd

end

"""
  data(mesh)

Return matrices defining the vertices, element connectivity, and element types
of the specified mesh.
"""
function data(mesh::Mesh)

  meshptr = mesh.ptr

  ndim = dimensions(mesh)
  nvtx = nodecount(mesh)
  vtx = Array{Float64,2}(undef, nvtx, ndim)
  vtxptr = pointer(vtx)

  # for i in 0:(ndim-1), j in 0:(nvtx-1)
  #   icxx"""
  #     *$(vtxptr)++ = $(meshptr)->nlist[$j][$i];
  #     """
  # end

  icxx"""
    // vertex coordinate list
    for (int i = 0; i < $ndim; i++)
      for (int j = 0; j < $nvtx; j++)
            *$(vtxptr)++ = $(meshptr)->nlist[j][i];
  """

  mnnd = maxnodes(mesh)
  nele = elemcount(mesh)
  ele = Array{Cint}(undef, nele, mnnd)
  eleptr = pointer(ele)

  icxx"""
    // element index list
    // (1-based; value 0 indicates unused matrix entry)
    for (int i = 0; i < $mnnd; i++)
        for (int j = 0; j < $nele; j++)
            if (i < $(meshptr)->elist[j]->nNode())
                *$(eleptr)++ = $(meshptr)->elist[j]->Node[i]+1;
            else
                *$(eleptr)++ = 0;
  """

  elt = Vector{Cint}(undef, nele)
  eltptr = pointer(elt)

  icxx"""
    // element type list
    for (int i = 0; i < $nele; i++)
        *$(eltptr)++ = $(meshptr)->elist[i]->Type();
  """

  return vtx, ele, elt

end

"""
    surface(mesh)

Extract vertices and element connectivity of the surface of the mesh.
"""
function surface(mesh::Mesh)

  meshptr = mesh.ptr

  nlen = nodecount(mesh)
  ndim = dimensions(mesh)
  nbnd = icxx"""$(meshptr)->nlist.NumberOf(BND_ANY);"""

  vtx = Array{Float64}(undef, nbnd, ndim)

  icxx"""
    for (int i = 0; i < $(ndim); i++)
      for (int j = 0; j < $(nlen); j++)
        if ($(meshptr)->nlist[j].isBnd())
          *$(pointer(vtx))++ = $(meshptr)->nlist[j][i];
  """

  bndidx = Array{Cint}(undef, nlen)
  nnd = [Cint(0)]
  nface = [Cint(0)]

  icxx"""
    int  *bndellist, *bndsdlist;

    for (int j = 0, k = 0; j < $(nlen); j++)
      $(pointer(bndidx))[j] = ($(meshptr)->nlist[j].isBnd() ? k++ : -1);

    // NB: This assumes all elements contain the same number of vertices
    *$(pointer(nnd)) = 0;
    *$(pointer(nface)) = $(meshptr)->BoundaryList (&bndellist, &bndsdlist);
    for (int j = 0; j < *$(pointer(nface)); j++)
      *$(pointer(nnd)) = ::max (*$(pointer(nnd)), $(meshptr)->elist[bndellist[j]]->nSideNode(bndsdlist[j]));

    delete[] bndellist;
    delete[] bndsdlist;
  """

  idx = Array{Cint}(undef, nnd[1], nface[1])

  icxx"""
    int sd, nn, nd, bn, *bndellist, *bndsdlist;

    $(meshptr)->BoundaryList (&bndellist, &bndsdlist);

    for (int i = 0; i < *$(pointer(nnd)); i++)
      for (int j = 0; j < *$(pointer(nface)); j++) {
        Element *pel = $(meshptr)->elist[bndellist[j]];
        sd = bndsdlist[j];
        nn = pel->nSideNode (sd);
        if (i < nn) {
          nd = pel->Node[pel->SideNode (sd, i)];
          bn = $(pointer(bndidx))[nd]+1;
        } else bn = 0;
        *$(pointer(idx))++ = bn;
      }
  """

  return vtx, idx

end

"""
    boundingbox(mesh)

Return the bounding box of the mesh in a ``2 \times dim`` matrix.
"""
function boundingbox(mesh::Mesh)

  dim = dimensions(mesh)
  bb = Array{Cdouble}(undef, 2,dim)

  icxx"""
    double *jptr = $(pointer(bb));

    Point pmin($dim), pmax($dim);
    $(mesh.ptr)->BoundingBox(pmin, pmax);

    for(int i=0; i<$dim; i++) {
      *jptr++ = pmin[i];
      *jptr++ = pmax[i];
    }
  """

  return bb

end

"""
    size(mesh)

Compute the total area (2D) for volume (3D) of the mesh.
"""
fullsize(mesh::Mesh) = @cxx mesh.ptr->FullSize()


text:
using .TurbulenceClosures: ∂z_aaf, ▶z_aaf

const g_Earth = 9.80665

#=
Supported buoyancy types:

- Nothing
- SeawaterBuoyancy
=#

#####
##### Functions for buoyancy = nothing
#####

@inline buoyancy_perturbation(i, j, k, grid::AbstractGrid{T}, ::Nothing, C) where T = zero(T)
@inline buoyancy_frequency_squared(i, j, k, grid::AbstractGrid{T}, ::Nothing, C) where T = zero(T)

#####
##### Seawater buoyancy for buoyancy determined by temperature and salinity
#####

"""
    BuoyancyTracer <: AbstractBuoyancy{Nothing}

Type indicating that the tracer `T` represents buoyancy.
"""
struct BuoyancyTracer <: AbstractBuoyancy{Nothing} end

@inline buoyancy_perturbation(i, j, k, grid, ::BuoyancyTracer, C) = @inbounds C.T[i, j, k]
@inline buoyancy_frequency_squared(i, j, k, grid, ::BuoyancyTracer, C) = ∂z_aaf(i, j, k, grid, C.T)

"""
    SeawaterBuoyancy{G, EOS} <: AbstractBuoyancy{EOS}

Buoyancy model for temperature- and salt-stratified seawater.
"""
struct SeawaterBuoyancy{G, EOS} <: AbstractBuoyancy{EOS}
    gravitational_acceleration :: G
    equation_of_state :: EOS
end

"""
    SeawaterBuoyancy([T=Float64;] gravitational_acceleration = g_Earth,
                                  equation_of_state = LinearEquationOfState(T))

Returns parameters for a temperature- and salt-stratified seawater buoyancy model
with a `gravitational_acceleration` constant (typically called 'g'), and an
`equation_of_state` that related temperature and salinity (or conservative temperature
and absolute salinity) to density anomalies and buoyancy.
"""
function SeawaterBuoyancy(T=Float64;
                          gravitational_acceleration = g_Earth,
                          equation_of_state = LinearEquationOfState(T))
    return SeawaterBuoyancy{T, typeof(equation_of_state)}(gravitational_acceleration, equation_of_state)
end

"""
    buoyancy_frequency_squared(i, j, k, grid, b::SeawaterBuoyancy, C)

Returns the buoyancy frequency squared for temperature and salt-stratified water,

```math
N^2 = g \\left ( \\alpha \\partial_z T - \\beta \\partial_z S \\right ) \\, ,
```

where ``\$ g \$`` is gravitational acceleration, ``\$ \\alpha \$`` is the thermal expansion
coefficient, ``\$ \\beta \$`` is the haline contraction coefficient, ``\$ T \$`` is
temperature or conservative temperature, where applicable, and ``\$ S \$`` is the
salinity or absolute salinity, where applicable.
"""
@inline buoyancy_frequency_squared(i, j, k, grid, b::SeawaterBuoyancy, C) =
    b.gravitational_acceleration * (
           thermal_expansion(i, j, k, grid, b.equation_of_state, C) * ∂z_aaf(i, j, k, grid, C.T)
        - haline_contraction(i, j, k, grid, b.equation_of_state, C) * ∂z_aaf(i, j, k, grid, C.S) )

#####
##### Linear equation of state
#####

"""
    LinearEquationOfState{T} <: AbstractEquationOfState

Linear equation of state for seawater.
"""
struct LinearEquationOfState{T} <: AbstractEquationOfState
    α :: T
    β :: T
end

"""
    LinearEquationOfState([T=Float64;] α=1.67e-4, β=7.80e-4)

Returns parameters for a linear equation of state for seawater with
thermal expansion coefficient `α` [K⁻¹] and haline contraction coefficient
`β` [ppt⁻¹]. The buoyancy perturbation associated with a linear equation of state is

    `b = g * α * T - g * β * S`

Default constants are taken from Table 1.2 (page 33) of Vallis, "Atmospheric and Oceanic Fluid
Dynamics: Fundamentals and Large-Scale Circulation" (2ed, 2017).
"""
LinearEquationOfState(T=Float64; α=1.67e-4, β=7.80e-4) =
    LinearEquationOfState{T}(α, β)

const LinearSeawaterBuoyancy = SeawaterBuoyancy{FT, <:LinearEquationOfState} where FT

@inline buoyancy_perturbation(i, j, k, grid, b::LinearSeawaterBuoyancy, C) =
    return @inbounds b.gravitational_acceleration * (   b.equation_of_state.α * C.T[i, j, k]
                                                      - b.equation_of_state.β * C.S[i, j, k] )

@inline  thermal_expansion(i, j, k, grid, eos::LinearEquationOfState, C) = eos.α
@inline haline_contraction(i, j, k, grid, eos::LinearEquationOfState, C) = eos.β

#####
##### Nonlinear equations of state
#####

@inline buoyancy_perturbation(i, j, k, grid, b::AbstractBuoyancy{<:AbstractNonlinearEquationOfState}, C) =
    - b.gravitational_acceleration * ρ′(i, j, k, grid, b.equation_of_state, C) / b.ρ₀

#####
##### Roquet et al 2015 idealized nonlinear equations of state
#####

# Reference: Table 3 in Roquet et el, "Defining a Simplified yet 'Realistic' Equation of State for Seawater", (JPO, 2015)
optimized_roquet_coeffs = Dict(
                  :linear => (R₀₁₀ = -1.775e-1, R₁₀₀ = 7.718e-1, R₀₂₀ = 0,         R₀₁₁ = 0,          R₂₀₀ = 0,         R₁₀₁ = 0,         R₁₁₀ = 0),
               :cabbeling => (R₀₁₀ = -0.844e-1, R₁₀₀ = 7.718e-1, R₀₂₀ = -4.561e-3, R₀₁₁ = 0,          R₂₀₀ = 0,         R₁₀₁ = 0,         R₁₁₀ = 0),
:cabbeling_thermobaricity => (R₀₁₀ = -0.651e-1, R₁₀₀ = 7.718e-1, R₀₂₀ = -5.027e-3, R₀₁₁ = -2.5681e-5, R₂₀₀ = 0,         R₁₀₁ = 0,         R₁₁₀ = 0),
                :freezing => (R₀₁₀ = -0.491e-1, R₁₀₀ = 7.718e-1, R₀₂₀ = -5.027e-3, R₀₁₁ = -2.5681e-5, R₂₀₀ = 0,         R₁₀₁ = 0,         R₁₁₀ = 0),
            :second_order => (R₀₁₀ =  0.182e-1, R₁₀₀ = 8.078e-1, R₀₂₀ = -4.937e-3, R₀₁₁ = -2.4677e-5, R₂₀₀ = -1.115e-4, R₁₀₁ = -8.241e-6, R₁₁₀ = -2.446e-3)
)

"""
    RoquetIdealizedNonlinearEquationOfState{F, C, T} <: AbstractNonlinearEquationOfState

Parameters associated with the idealized nonlinear equation of state proposed by
Roquet et al., "Defining a Simplified yet 'Realistic' Equation of State for Seawater",
Journal of Physical Oceanography (2015).
"""
struct RoquetIdealizedNonlinearEquationOfState{F, C, T} <: AbstractNonlinearEquationOfState
                   ρ₀ :: T
    polynomial_coeffs :: C
end

type_convert_roquet_coeffs(T, coeffs) = NamedTuple{propertynames(coeffs)}(Tuple(T(R) for R in coeffs))

"""
    RoquetIdealizedNonlinearEquationOfState([T=Float64,] flavor, ρ₀=1024.6,
                                            polynomial_coeffs=optimized_roquet_coeffs[flavor])

Returns parameters for the idealized polynomial nonlinear equation of state with
reference density `ρ₀` and `polynomial_coeffs` proposed by Roquet et al., "Defining a
Simplified yet 'Realistic' Equation of State for Seawater", Journal of Physical
Oceanography (2015). The default reference density is `ρ₀ = 1024.6 kg m⁻³`, the average
surface density of seawater in the world ocean.

The `flavor` of the nonlinear equation of state is a symbol that selects a set of optimized
polynomial coefficients defined in Table 2 of Roquet et al., "Defining a Simplified yet
'Realistic' Equation of State for Seawater", Journal of Physical Oceanography (2015), and
further discussed in the text surrounding equations (12)--(15). The optimization minimizes
errors in estimated horizontal density gradient estiamted from climatological temperature
and salinity distributions between the 5 simplified forms chosen by Roquet et. al
and the full-fledged [TEOS-10](http://www.teos-10.org) equation of state.

The equations of state define the density anomaly `ρ′`, and have
the polynomial form

    `ρ′(T, S, D) = Σᵢⱼₐ Rᵢⱼₐ Tⁱ Sʲ Dᵃ`,

where `T` is conservative temperature, `S` is absolute salinity, and `D` is the
geopotential depth, currently just `D = -z`. The `Rᵢⱼₐ` are constant coefficients.

Flavors of idealized nonlinear equations of state
=================================================

    - `:linear`: a linear equation of state, `ρ′ = R₁₀₀ * T + R₀₁₀ * S`

    - `:cabbeling`: includes quadratic temperature term,
                    `ρ′ = R₁₀₀ * T + R₀₁₀ * S + R₀₂₀ * T^2`

    - `:cabbeling_thermobaricity`: includes 'thermobaricity' term,
                                   `ρ′ = R₁₀₀ * T + R₀₁₀ * S + R₀₂₀ * T^2 + R₀₁₁ * T * D`

    - `:freezing`: same as `:cabbeling_thermobaricity` with modified constants to increase
                   accuracy near freezing

    - `:second_order`: includes quadratic salinity, halibaricity, and thermohaline term,
                       `ρ′ = R₁₀₀ * T + R₀₁₀ * S + R₀₂₀ * T^2 + R₀₁₁ * T * D`
                             + R₂₀₀ * S^2 + R₁₀₁ * S * D + R₁₁₀ * S * T`

Example
=======

julia> using Oceananigans

julia> eos = Oceananigans.RoquetIdealizedNonlinearEquationOfState(:cabbeling);

julia> eos.polynomial_coeffs
(R₀₁₀ = -0.0844, R₁₀₀ = 0.7718, R₀₂₀ = -0.004561, R₀₁₁ = 0.0, R₂₀₀ = 0.0, R₁₀₁ = 0.0, R₁₁₀ = 0.0)

References
==========

    - Roquet et al., "Defining a Simplified yet 'Realistic' Equation of State for
      Seawater", Journal of Physical Oceanography (2015).

    - "Thermodynamic Equation of State for Seawater" (TEOS-10), http://www.teos-10.org
"""
function RoquetIdealizedNonlinearEquationOfState(T, flavor=:cabbeling_thermobaricity;
                                                 polynomial_coeffs=optimized_roquet_coeffs[flavor], ρ₀=1024.6)
    typed_coeffs = type_convert_roquet_coeffs(T, polynomial_coeffs)
    return RoquetIdealizedNonlinearEquationOfState{flavor, typeof(typed_coeffs), T}(ρ₀, typed_coeffs)
end

RoquetIdealizedNonlinearEquationOfState(flavor::Symbol=:cabbeling_thermobaricity; kwargs...) =
    RoquetIdealizedNonlinearEquationOfState(Float64, flavor; kwargs...)

""" Return the geopotential depth at `i, j, k` at cell centers. """
@inline D_aac(i, j, k, grid) = @inbounds -grid.zC[k]
const D = D_aac

""" Return the geopotential depth at `i, j, k` at cell z-interfaces. """
@inline D_aaf(i, j, k, grid) = @inbounds -grid.zF[k]

@inline ρ′(i, j, k, grid, eos::RoquetIdealizedNonlinearEquationOfState, C) =
    @inbounds (   eos.polynomial_coeffs.R₁₀₀ * C.S[i, j, k]
                + eos.polynomial_coeffs.R₀₁₀ * C.T[i, j, k]
                + eos.polynomial_coeffs.R₀₂₀ * C.T[i, j, k]^2
                + eos.polynomial_coeffs.R₀₁₁ * C.T[i, j, k] * D(i, j, k, grid)
                + eos.polynomial_coeffs.R₂₀₀ * C.S[i, j, k]^2
                + eos.polynomial_coeffs.R₁₀₁ * C.S[i, j, k] * D(i, j, k, grid)
                + eos.polynomial_coeffs.R₁₁₀ * C.S[i, j, k] * C.T[i, j, k] )

@inline thermal_expansion(i, j, k, grid, eos::RoquetIdealizedNonlinearEquationOfState, C) =
    @inbounds (   eos.polynomial_coeffs.R₀₁₀
                + 2 * eos.polynomial_coeffs.R₀₂₀ * ▶z_aaf(i, j, k, grid, C.T)
                + eos.polynomial_coeffs.R₀₁₁ * D_aaf(i, j, k, grid)
                + eos.polynomial_coeffs.R₁₁₀ * ▶z_aaf(i, j, k, grid, C.T) )

@inline haline_contraction(i, j, k, grid, eos::RoquetIdealizedNonlinearEquationOfState, C) =
    @inbounds (   eos.polynomial_coeffs.R₁₀₀
                + 2 * eos.polynomial_coeffs.R₂₀₀ * ▶z_aaf(i, j, k, grid, C.S)
                + eos.polynomial_coeffs.R₁₀₁ * D_aaf(i, j, k, grid)
                + eos.polynomial_coeffs.R₁₁₀ * ▶z_aaf(i, j, k, grid, C.S) )


text:
# order 2
h = [33743097329453577701//Int128(10)^20,
    97759682018833491296//Int128(10)^20,
    93278808104030343530//Int128(10)^20,
    1, 1]
e = [1//1, 0, 0, 0, 0]
H = Diagonal(h)

q1 = [0,
    55932483188770411252//Int128(10)^20,
    -5932483188770411252//Int128(10)^20,
    0, 0]'
q2 = [-q1[2] 0 55932483188770411252//Int128(10)^20 0 0]
q3 = [-q1[3] -q2[3] 0 1//2 0]
q4 = [0 0 -1//2 0 1//2]
q5 = [0 0 0 -1//2 0]
Q = vcat(q1, q2, q3, q4, q5)

D1 = H \ (Q - 1//2*e*e'); display(D1)

# order 4
h = [21427296612044126417//Int128(10)^20,
    1123759588488739348//Int128(10)^18,
    1434458792494126//Int128(10)^15,
    10917323021736130836//Int128(10)^19,
    09883816115129601975//Int128(10)^19,
    1, 1, 1]
e = [1//1, 0, 0, 0, 0, 0, 0, 0]
H = Diagonal(h)

q1 = [0,
     66884898686930380508//Int128(10)^20,
    -25171531878753856238//Int128(10)^20,
    10997619816825822803//Int128(10)^20,
    -27109866250023470592//Int128(10)^21,
    0, 0, 0]'
q2 = [-q1[2], 0,
     92214436948640491071//Int128(10)^20,
    -32412368653542520402//Int128(10)^20,
    70828303918324098284//Int128(10)^21,
    0, 0, 0]'
q3 = [-q1[3], -q2[3], 0,
    8180378089216779335//Int128(10)^19,
    -14760875822281158529//Int128(10)^20,
    0, 0, 0]'
q4 = [-q1[4], -q2[4], -q3[4], 0,
    68722365388784429092//Int128(10)^20,
    -1//12, 0, 0]'
q5 = [-q1[5] -q2[5] -q3[5] -q4[5] 0 2//3 -1//12 0]
q6 = [0 0 0 1//12 -2//3 0 2//3 -1//12]
q7 = [0 0 0 0 1//12 -2//3 0 2//3]
q8 = [0 0 0 0 0 1//12 -2//3 0]
Q = vcat(q1, q2, q3, q4, q5, q6, q7, q8)

D1 = H \ (Q - 1//2*e*e'); display(D1)

# order 6
h = [15109714532036117328//Int128(10)^20,
    80967585357107013003//Int128(10)^20,
    10911427148079254850//Int128(10)^19,
    10435269041571577756//Int128(10)^19,
    98680905919946100728//Int128(10)^20,
    10037581831426163456//Int128(10)^19,
    99943556356761752125//Int128(10)^20,
    1, 1, 1, 1]
e = [1//1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
H = Diagonal(h)

q1 = [0,
    66670790901888837033//Int128(10)^20,
    -23418791580399147484//Int128(10)^20,
    84251264588860596867//Int128(10)^21,
    -15923290838179674350//Int128(10)^21,
    -15653772860347171721//Int128(10)^22,
    071741032045689717567//Int128(10)^23,
    0, 0, 0, 0]'
q2 = [-q1[2], 0,
    89405599296515541581//Int128(10)^20,
    -28597427787314667440//Int128(10)^20,
    57056178538117177397//Int128(10)^21,
    41320613074890940489//Int128(10)^22,
    -25620459187266476645//Int128(10)^22,
    0, 0, 0, 0]'
q3 = [-q1[3], -q2[3], 0,
    82961715259707113283//Int128(10)^20,
    -18233747042994439227//Int128(10)^20,
    83784382166533084621//Int128(10)^22,
    42099567773838744673//Int128(10)^22,
    0, 0, 0, 0]'
q4 = [-q1[4], -q2[4], -q3[4], 0,
    75419218459746682761//Int128(10)^20,
    -14034899831339963049//Int128(10)^20,
    14050953028717865444//Int128(10)^21,
    0, 0, 0, 0]'
q5 = [-q1[5], -q2[5], -q3[5], -q4[5], 0,
    74751473989919011204//Int128(10)^20,
    -15119380469839682133//Int128(10)^20,
    1//60, 0, 0, 0]'
q6 = [-q1[6], -q2[6], -q3[6], -q4[6], -q5[6], 0,
    75144419715723149872//Int128(10)^20,
    -3//20, 1//60, 0, 0]'
q7 = [-q1[7] -q2[7] -q3[7] -q4[7] -q5[7] -q6[7] 0 3//4 -3//20 1//60 0]
q8 = [0 0 0 0 -1//60 3//20 -3//4 0 3//4 -3//20 1//60]
q9 = [0 0 0 0 0 -1//60 3//20 -3//4 0 3//4 -3//20]
q10 = [0 0 0 0 0 0 -1//60 3//20 -3//4 0 3//4]
q11 = [0 0 0 0 0 0 0 -1//60 3//20 -3//4 0]
Q = vcat(q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11)

D1 = H \ (Q - 1//2*e*e'); display(D1)

# order 8
h = [12163222110707502878//Int128(10)^20,
    65235832636546639982//Int128(10)^20,
    87730414198101010954//Int128(10)^20,
    97388951771079542799//Int128(10)^20,
    10072514376844677230//Int128(10)^19,
    99768726657776478834//Int128(10)^20,
    10005302998791085514//Int128(10)^19,
    99994066100338390832//Int128(10)^20,
    1, 1, 1, 1, 1]
e = [1//1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
H = Diagonal(h)

q1 = [0,
    66670790901888837033//Int128(10)^20,
    -21994030190635039046//Int128(10)^20,
    61752567584332553851//Int128(10)^21,
    -32312350944133128873//Int128(10)^22,
    -33934980320003350186//Int128(10)^22,
    15157027970563223705//Int128(10)^24,
    32350133072942893419//Int128(10)^23,
    0, 0, 0, 0, 0]'
q2 = [-q1[2], 0,
    86688767821045233147//Int128(10)^20,
    -24298087640343350527//Int128(10)^20,
    39549469619698650847//Int128(10)^21,
    20763528371484737510//Int128(10)^22,
    -65045489396961912976//Int128(10)^23,
    -40836028016484152445//Int128(10)^23,
    0, 0, 0, 0, 0]'
q3 = [-q1[3], -q2[3], 0,
    82065092584472835146//Int128(10)^20,
    -21014872891771196683//Int128(10)^20,
    38572177503610408523//Int128(10)^21,
    -16637807199883547459//Int128(10)^23,
    -46321740653650899692//Int128(10)^24,
    0, 0, 0, 0, 0]'
q4 = [-q1[4], -q2[4], -q3[4], 0,
    81102837324727866266//Int128(10)^20,
    -20423896795865859484//Int128(10)^20,
    34947121761434371005//Int128(10)^21,
    23139100244270367378//Int128(10)^22,
    0, 0, 0, 0, 0]'
q5 = [-q1[5], -q2[5], -q3[5], -q4[5], 0,
    80054065093594950025//Int128(10)^20,
    -19699167992690472530//Int128(10)^20,
    37220336417235830241//Int128(10)^21,
    -1//280, 0, 0, 0, 0]'
q6 = [-q1[6], -q2[6], -q3[6], -q4[6], -q5[6], 0,
    79903079313046586260//Int128(10)^20,
    -19999788736822592697//Int128(10)^20,
    4//105, -1//280, 0, 0, 0]'
q7 = [-q1[7], -q2[7], -q3[7], -q4[7], -q5[7], -q6[7], 0,
    80016334685519857774//Int128(10)^20,
    -1//5, 4//105, -1//280, 0, 0]'
q8 = [-q1[8], -q2[8], -q3[8], -q4[8], -q5[8], -q6[8], -q7[8], 0, 4//5, -1//5, 4//105, -1//280, 0]'
q9 = [0, 0, 0, 0, 1//280, -4//105, 1//5, -4//5, 0, 4//5, -1//5, 4//105, -1//280]'
q10 = [0, 0, 0, 0, 0, 1//280, -4//105, 1//5, -4//5, 0, 4//5, -1//5, 4//105]'
q11 = [0, 0, 0, 0, 0, 0, 1//280, -4//105, 1//5, -4//5, 0, 4//5, -1//5]'
q12 = [0, 0, 0, 0, 0, 0, 0, 1//280, -4//105, 1//5, -4//5, 0, 4//5]'
q13 = [0, 0, 0, 0, 0, 0, 0, 0, 1//280, -4//105, 1//5, -4//5, 0]'
Q = vcat(q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13)

D1 = H \ (Q - 1//2*e*e'); display(D1)


text:
function save_api_key(api_key::String)
    io = open(joinpath(homedir(), ".fred_api_key"), "w")
    println(io, api_key)
    close(io)
    return nothing
end


text:
module FinancialSeries

# documentation
using Docile

if VERSION < v"0.4-"
  using Dates
else
  using Base.Dates
end
using Reexport
@reexport using TimeSeries

import Base: show, getindex

@document

export AbstractInstrument, AbstractCurrency, AbstractFinancialID,
       Stock, Currency, CurrencyPair,
       Ticker, CUSIP, BloombergID, ReutersID,
       FinancialTimeSeries,
       timestamp, values, colnames, instrument,
       TickData,
       USD, EUR, GBP, AUD, JPY, EURUSD, EURJPY, USDGBP, USDAUD, USDNZD, USDJPY, F, G, H, J, K, M, N, Q, U, V, X, Z,
       merge, parsedatetime, parsedatetime1, parsedatetime_from_TOS, makedatetime, datetolastsecond, discretesignal
       # yahoo, fred

include("instruments.jl")
include("financialtimeseries.jl")

#include("tickdata.jl")
#include("getindex.jl")
#include("utilities.jl")
#include("readwrite.jl")

end


text:
@testset "NNlib" begin

@testset "batched_mul" begin
    using NNlib: batched_mul, batched_adjoint, batched_transpose

    A = randn(Float32, 3,3,2);
    B = randn(Float32, 3,3,2);

    C = batched_mul(A, B)
    @test cu(C) ≈ batched_mul(cu(A), cu(B))

    Ct = batched_mul(batched_transpose(A), B)
    @test cu(Ct) ≈ batched_mul(batched_transpose(cu(A)), cu(B))

    Ca = batched_mul(A, batched_adjoint(B))
    @test cu(Ca) ≈ batched_mul(cu(A), batched_adjoint(cu(B)))
end

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "GLMNet"
version = v"5.0.0"

# Collection of sources required to build GLMNet
sources = [
    "https://github.com/JuliaStats/GLMNet.jl.git" =>
    "3030d4f6e034dca06f25b1343a33b9c799cda7f7",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/GLMNet.jl/deps/

flags="-fdefault-real-8 -ffixed-form -shared -O3"
if [[ ${target} != *mingw* ]]; then
     flags="${flags} -fPIC"
fi
if [[ ${target} != aarch64* ]] && [[ ${target} != arm* ]]; then
     flags="${flags} -m${nbits}"
fi

libdir=$prefix/lib
if [[ ${target} == *mingw* ]]; then
     libdir=$prefix/bin
fi

${FC} ${LDFLAGS} ${flags} glmnet5.f90 -o libglmnet.${dlext}

mkdir -p $libdir
mv libglmnet.${dlext} $libdir/

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf),
    MacOS(:x86_64),
    FreeBSD(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]
platforms = expand_gcc_versions(platforms)

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

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

]

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

max_stars_repo_name: JackDunnNZ/GLMNetBuilder

max_issues_repo_name: JackDunnNZ/GLMNetBuilder

max_forks_repo_name: JackDunnNZ/GLMNetBuilder


text:
using ApproxFun

function smoothness(f::Fun)
	s = f.space
	∫ = DefiniteIntegral(s)
	∂² = Derivative(s,2)
	return Number(-∫ * (∂² * f)^2)
end

s = Chebyshev()
ℓ(x) = @inbounds smoothness(Fun(s,x)) - x[1]^2 - x[2]^2

using ForwardDiff

ForwardDiff.gradient(ℓ, randn(10))

using TransformVariables
t = as((x=as(Vector, 10),))

using SampleChainsDynamicHMC
chain = initialize!(DynamicHMCChain, ℓ, t)

smoothness(f)

julia> using ApproxFun

julia> function smoothness(f)
           s = f.space
           ∫ = DefiniteIntegral(s)
           -∫ * (Derivative(s,2) * f)^2
       end
smoothness (generic function with 1 method)

julia> x = Fun(identity, 0..1)
Fun(Chebyshev(0..1),[0.5, 0.5])

julia> f = exp(x);

julia> smoothness(f)
-3.1945280494653217 on ApproxFunBase.AnyDomain()

function b(i)
    Derivative(2) * Fun(s,[zeros(i-1);1.0])
end

b(4)

function Sᵢⱼ(i,j)
    (DefiniteIntegral(s) * Multiplication(b(i)) * b(j)).coefficients[1]
end

Sᵢⱼ(7,5)

function makeS(n, space=Chebyshev())
    z =  Fun(Ultraspherical(2), zeros(n))
    b = [ z; z;[ Derivative(2) * Fun(space,[zeros(j-1);1.0; zeros(n-j+2)]) for j in 3:n]]
    S = Matrix{Float64}(undef, n, n)
    for j in 1:n
        for i in 1:n
            S[i,j] = prod((DefiniteIntegral() * Multiplication(b[i],space) * b[j]).coefficients)
        end
    end
    return S
end

n = 10
makeS(20)


text:
"""
    surface(cmd0::String="", arg1=nothing; kwargs...)

Reads randomly-spaced (x,y,z) triples and produces a binary grid file of gridded values z(x,y) by solving:

		(1 - T) * L (L (z)) + T * L (z) = 0

Full option list at [`surface`]($(GMTdoc)surface.html)

Parameters
----------

- $(GMT.opt_R)
- **I** | **inc** :: [Type => Str | Number]

    *x_inc* [and optionally *y_inc*] is the grid spacing.
    ($(GMTdoc)surface.html#i)
- **A** | **aspect_ratio** :: [Type => Number]

    Aspect ratio. If desired, grid anisotropy can be added to the equations.
    ($(GMTdoc)surface.html#a)
- **C** | **convergence** :: [Type => Number]

    Convergence limit. Iteration is assumed to have converged when the maximum absolute change in any
    grid value is less than convergence_limit.
    ($(GMTdoc)surface.html#c)
- **D** | **breakline** :: [Type => String | NamedTuple]     `Arg = breakline[+z[level]] | (data=Array, [zlevel=x])`

    Use xyz data in the breakline file as a ‘soft breakline’, that is a line whose vertices will be used to
    constrain the nearest grid nodes without any further interpolation.
    ($(GMTdoc)surface.html#d)
- **G** | **outgrid** :: [Type => Str]

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = surface(....) form.
    ($(GMTdoc)surface.html#g)
- **Ll** | **lower** :: [Type => Str | Number]

    Impose limits on the output solution. lower sets the lower bound. lower can be the name of a grid
    file with lower bound values, a fixed value, d to set to minimum input value,
    ($(GMTdoc)surface.html#l)
- **Lu** | **upper** :: [Type => Str | Number]

    ($(GMTdoc)surface.html#l)
- **N** | **max_iter** :: [Type => Number]

    Number of iterations. Iteration will cease when convergence_limit is reached or when number of
    iterations reaches max_iterations.
    ($(GMTdoc)surface.html#n)
- **Q** | **suggest** :: [Type => Bool]

    Suggest grid dimensions which have a highly composite greatest common factor.
    ($(GMTdoc)surface.html#q)
- **S** | **search_radius** :: [Type => Number | Str]

    Sets the resolution of the projected grid that will be created.
    ($(GMTdoc)surface.html#s)
- **T** | **tension** :: [Type => Number | Str]

    Tension factor[s]. These must be between 0 and 1.
    ($(GMTdoc)surface.html#t)
- $(GMT.opt_V)
- **Z** | **over_relaxation** :: [Type => Str | GMTgrid]

    Over-relaxation factor. This parameter is used to accelerate the convergence; it is a number between 1 and 2.
    ($(GMTdoc)surface.html#z)
- $(GMT.opt_a)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_r)
- $(GMT.opt_swap_xy)
"""
function surface(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && return monolitic("surface", cmd0, arg1)

	d = KW(kwargs);     arg2 = nothing
	help_show_options(d)		# Check if user wants ONLY the HELP mode
	cmd, = parse_common_opts(d, "", [:R :I :V_params :a :bi :di :e :f :h :i :r :yx])
	cmd  = parse_these_opts(cmd, d, [[:A :aspect_ratio], [:C :convergence], [:G :grid :outgrid],
	                                 [:Ll :lower], [:Lu :upper], [:N :max_iter], [:Q :suggest], [:S :search_radius], [:T :tension], [:Z :over_relaxation]])
	cmd, args, n, = add_opt(d, cmd, 'D', [:D :breakline], :data, Array{Any,1}([arg1, arg2]), (zlevel="+z",))
	if (n > 0)  arg1, arg2 = args[:]  end

	common_grd(d, cmd0, cmd, "surface ", arg1, arg2)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
surface(arg1, cmd0::String=""; kw...) = surface(cmd0, arg1; kw...)


text:
# Code to asses the effect of moment-compensation on the diffusion signal

## Read
#Waveform
using MRIsim, PlotlyJS
b = [0, 50, 200, 400, 600, 790]
E = Array{Any, 2}(undef,6,4)
FID = false
path_file = "/home/ccp/Documents/MRIsim.jl/"
#Timings
f = readlines(path_file*"qte_timings.txt")
timings = tryparse.(Float64, split(f[2]," "))
δ1, rf180, δ2 = timings[timings .≠ nothing]*1e-3
#Normal DTI
G, n1, n2 = MRIsim.read_diff_fwf(path_file*"QTE_Waveforms/MX_MC2_b790.txt")
Gmax = 21e-3
g1, g2 = G[1:20:n1,1], G[1:20:n2,4]
seq_dti = Sequence(Gmax*[Grad(g1,δ1) delay(rf180) -Grad(g2,δ2)])
#ACQ
if FID
    seq_ACQ = Sequence([Grad(0,0)], [RF(0,0)], DAC(1,0))
else
    FOV = 4e-2 #10cm
    Nxy = 20
    ΔT = 20e-6 #100um
    Gmax_acq = 60e-3
    seq_ACQ,_,_,_ = PulseDesigner.EPI_base(FOV, Nxy, ΔT, Gmax_acq)
end
#Phantom
Ns = 10_000; Nsx = floor(Int64, sqrt(Ns) ); Ns = Nsx*Nsx
ADC = 2e-9 #Normal diffusion coefficient
v = 6e-2 # 2mm/s Capilar flow
Δx = .5e-2 #3mm
x = range(-Δx/2,Δx/2,length=Nsx) #spin coordinates
y = range(-Δx/2,Δx/2,length=Nsx)  #spin coordinates
xx, yy = x .+ y'*0, x*0 .+ y' #grid points
T = dur(seq_dti)+dur(seq_ACQ)/2
phantom  = Phantom(x=xx[:], y=yy[:], Dλ1=ADC*ones(Ns), Dλ2=ADC*ones(Ns))
phantom2 = Phantom(x=xx[:], y=yy[:], ux=(x,y,z,t)->v*t, Dλ1=ADC*ones(Ns), Dλ2=ADC*ones(Ns))
#Iter
for idx = 1:5
#Seq
G, n1, n2 = MRIsim.read_diff_fwf(path_file*"QTE_Waveforms/MX_MC2_b$(b[idx+1]).txt")
g1, g2 = G[1:20:n1,1], G[1:20:n2,4]
seq = Sequence(Gmax*[Grad(g1,δ1) delay(rf180) -Grad(g2,δ2)])
seq += seq_ACQ
#Simulate
f = sqrt(b[idx+1]/790)
simParams = Dict(:Δt => 10e-6, :Nblocks=> 1)
if FID
    recParams = Dict(:skip_rec => true)
else
    recParams = Dict(:epi => true, :recon => :fft, :Nx => Nxy)
end
S_MM = MRIsim.simulate(phantom, seq, simParams, recParams)
S_DTI = MRIsim.simulate(phantom, f*seq_dti+seq_ACQ, simParams, recParams)
S_MM2 = MRIsim.simulate(phantom2, seq, simParams, recParams)
S_DTI2 = MRIsim.simulate(phantom2, f*seq_dti+seq_ACQ, simParams, recParams)
S0 = MRIsim.simulate(phantom, 0*seq, simParams, recParams)
S02 = MRIsim.simulate(phantom2, 0*seq, simParams, recParams)
#Saving results
E[idx+1,1] = S_MM
E[idx+1,2] = S_DTI
E[idx+1,3] = S_MM2
E[idx+1,4] = S_DTI2
E[1,1] = S0
E[1,2] = S0
E[1,3] = S02
E[1,4] = S02
end

##Plot
i = 4
zmax = 12
plot(heatmap(z=abs.(E[i,1]), transpose = false, zmin=0, zmax=zmax), Layout(;title="No mov. MM-DTI b$(b[i])"))
plot(heatmap(z=abs.(E[i,2]), transpose = false, zmin=0, zmax=zmax), Layout(;title="No mov. DTI b$(b[i])"))
plot(heatmap(z=abs.(E[i,3]), transpose = false, zmin=0, zmax=zmax), Layout(;title="Mov. MM-DTI b$(b[i])"))
plot(heatmap(z=abs.(E[i,4]), transpose = false, zmin=0, zmax=zmax), Layout(;title="Mov. DTI b$(b[i])"))

p = MRIsim.plot_grads_moments(seq_dti+seq_ACQ)


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dcc_upload

"""
    dcc_upload(;kwargs...)
    dcc_upload(children::Any;kwargs...)
    dcc_upload(children_maker::Function;kwargs...)

An Upload component.
Upload components allow your app to accept user-uploaded files via drag'n'drop
Keyword arguments:
- `children` (a list of or a singular dash component, string or number | String; optional): Contents of the upload component
- `id` (String; optional): The ID of this component, used to identify dash components
in callbacks. The ID needs to be unique across all of the
components in an app.
- `contents` (String | Array of Strings; optional): The contents of the uploaded file as a binary string
- `filename` (String | Array of Strings; optional): The name of the file(s) that was(were) uploaded.
Note that this does not include the path of the file
(for security reasons).
- `last_modified` (Real | Array of Reals; optional): The last modified date of the file that was uploaded in unix time
(seconds since 1970).
- `accept` (String; optional): Allow specific types of files.
See https://github.com/okonet/attr-accept for more information.
Keep in mind that mime type determination is not reliable across
platforms. CSV files, for example, are reported as text/plain
under macOS but as application/vnd.ms-excel under Windows.
In some cases there might not be a mime type set at all.
See: https://github.com/react-dropzone/react-dropzone/issues/276
- `disabled` (Bool; optional): Enable/disable the upload component entirely
- `disable_click` (Bool; optional): Disallow clicking on the component to open the file dialog
- `max_size` (Real; optional): Maximum file size. If `-1`, then infinite
- `min_size` (Real; optional): Minimum file size
- `multiple` (Bool; optional): Allow dropping multiple files
- `className` (String; optional): HTML class name of the component
- `className_active` (String; optional): HTML class name of the component while active
- `className_reject` (String; optional): HTML class name of the component if rejected
- `className_disabled` (String; optional): HTML class name of the component if disabled
- `style` (Dict; optional): CSS styles to apply
- `style_active` (Dict; optional): CSS styles to apply while active
- `style_reject` (Dict; optional): CSS styles if rejected
- `style_disabled` (Dict; optional): CSS styles if disabled
- `loading_state` (optional): Object that holds the loading state object coming from dash-renderer. loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
"""
function dcc_upload(; kwargs...)
        available_props = Symbol[:children, :id, :contents, :filename, :last_modified, :accept, :disabled, :disable_click, :max_size, :min_size, :multiple, :className, :className_active, :className_reject, :className_disabled, :style, :style_active, :style_reject, :style_disabled, :loading_state]
        wild_props = Symbol[]
        return Component("dcc_upload", "Upload", "dash_core_components", available_props, wild_props; kwargs...)
end

dcc_upload(children::Any; kwargs...) = dcc_upload(;kwargs..., children = children)
dcc_upload(children_maker::Function; kwargs...) = dcc_upload(children_maker(); kwargs...)

max_stars_repo_name: asc-csa/dash-core-components

max_issues_repo_name: asc-csa/dash-core-components

max_forks_repo_name: asc-csa/dash-core-components


text:
"""
    function RefElemData(elementType::Line, approxType::SBP, N)
    function RefElemData(elementType::Quad, approxType::SBP, N)
    function RefElemData(elementType::Hex,  approxType::SBP, N)
    function RefElemData(elementType::Tri,  approxType::SBP, N)

SBP reference element data for `Quad()`, `Hex()`, and `Tri()` elements.

For `Line()`, `Quad()`, and `Hex()`, `approxType` is `SBP{TensorProductLobatto}`.

For `Tri()`, `approxType` can be `SBP{Kubatko{LobattoFaceNodes}}`, `SBP{Kubatko{LegendreFaceNodes}}`, or `SBP{Hicken}`.
"""
function RefElemData(elementType::Line, approxType::SBP{TensorProductLobatto}, N; tol = 100*eps())

    rd = RefElemData(elementType, N, quad_rule_vol = gauss_lobatto_quad(0,0,N))

    rd = @set rd.Vf = droptol!(sparse(rd.Vf), tol)
    rd = @set rd.LIFT = Diagonal(rd.wq) \ (rd.Vf' * Diagonal(rd.wf)) # TODO: make this more efficient with LinearMaps?

    return _convert_RefElemData_fields_to_SBP(rd, approxType)
end

function RefElemData(elementType::Quad, approxType::SBP{TensorProductLobatto}, N; tol = 100*eps())

    # make 2D SBP nodes/weights
    r1D, w1D = gauss_lobatto_quad(0, 0, N)
    sq, rq = vec.(NodesAndModes.meshgrid(r1D)) # this is to match ordering of nrstJ
    wr, ws = vec.(NodesAndModes.meshgrid(w1D))
    wq = wr .* ws
    quad_rule_vol = (rq, sq, wq)
    quad_rule_face = (r1D, w1D)

    rd = RefElemData(elementType, N, quad_rule_vol = quad_rule_vol, quad_rule_face = quad_rule_face)

    rd = @set rd.Vf = droptol!(sparse(rd.Vf), tol)
    rd = @set rd.LIFT = Diagonal(rd.wq) \ (rd.Vf' * Diagonal(rd.wf)) # TODO: make this more efficient with LinearMaps?

    return _convert_RefElemData_fields_to_SBP(rd, approxType)
end

function RefElemData(elementType::Hex, approxType::SBP{TensorProductLobatto}, N; tol = 100*eps())

    # make 2D SBP nodes/weights
    r1D, w1D = gauss_lobatto_quad(0, 0, N)
    rf, sf = vec.(NodesAndModes.meshgrid(r1D, r1D))
    wr, ws = vec.(NodesAndModes.meshgrid(w1D, w1D))
    wf = wr .* ws
    sq, rq, tq = vec.(NodesAndModes.meshgrid(r1D, r1D, r1D)) # this is to match ordering of nrstJ
    wr, ws, wt = vec.(NodesAndModes.meshgrid(w1D, w1D, w1D))
    wq = wr .* ws .* wt
    quad_rule_vol = (rq, sq, tq, wq)
    quad_rule_face = (rf, sf, wf)

    rd = RefElemData(elementType, N, quad_rule_vol = quad_rule_vol, quad_rule_face = quad_rule_face)

    rd = @set rd.Vf = droptol!(sparse(rd.Vf), tol)
    rd = @set rd.LIFT = Diagonal(rd.wq) \ (rd.Vf' * Diagonal(rd.wf)) # TODO: make this more efficient with LinearMaps?

    return _convert_RefElemData_fields_to_SBP(rd, approxType)
end

function RefElemData(elementType::Tri, approxType::SBP, N; tol = 100*eps())

    quad_rule_vol, quad_rule_face = diagE_sbp_nodes(elementType, approxType, N)

    # build polynomial reference element using quad rules; will be modified to create SBP RefElemData
    rd = RefElemData(elementType, Polynomial(), N; quad_rule_vol=quad_rule_vol, quad_rule_face=quad_rule_face)

    # determine Fmask = indices of face nodes among volume nodes
    Ef,Fmask = build_Ef_Fmask(rd)

    # Build traditional SBP operators from hybridized operators. See Section 3.2 of
    # [High-order entropy stable dG methods for the SWE](https://arxiv.org/pdf/2005.02516.pdf)
    # by Wu and Chan 2021. [DOI](https://doi.org/10.1016/j.camwa.2020.11.006)
    (Qrh, Qsh), _ = hybridized_SBP_operators(rd)
    Nq = length(rd.wq)
    Vh_sbp = [I(Nq); Ef]
    Qr = Vh_sbp'*Qrh*Vh_sbp
    Qs = Vh_sbp'*Qsh*Vh_sbp
    Dr,Ds = (x->diagm(1 ./ rd.wq) * x).((Qr,Qs))

    rd = @set rd.rst = quad_rule_vol[1:2]   # set nodes = SBP nodes
    rd = @set rd.rstq = quad_rule_vol[1:2]  # set quad nodes = SBP nodes
    rd = @set rd.Drst = (Dr,Ds)
    rd = @set rd.Fmask = vec(Fmask)

    # TODO: make these more efficient with custom operator?
    rd = @set rd.Vf = droptol!(sparse(Ef),tol)
    rd = @set rd.LIFT = Diagonal(rd.wq)\(rd.Vf'*Diagonal(rd.wf))

    # make V1 the interpolation matrix from triangle vertices to SBP nodal points
    rd = @set rd.V1 = vandermonde(elementType,N,rd.rst...)/rd.VDM * rd.V1

    # Vp operator = projects SBP nodal vector onto degree N polynomial, then interpolates to plotting points
    rd = @set rd.Vp = vandermonde(elementType,N,rd.rstp...)/rd.VDM * rd.Pq

    return _convert_RefElemData_fields_to_SBP(rd, approxType)
end

#####
##### Utilities for SBP
#####

function diagE_sbp_nodes(elem::Tri, approxType::SBP{Kubatko{LobattoFaceNodes}}, N)

    if N==6
        @warn "N=6 SBP operators with quadrature strength 2N-1 and Lobatto face nodes may require very small timesteps."
    end

    # from Ethan Kubatko, private communication
    vars = matread((@__DIR__)*"/data/sbp_nodes/KubatkoQuadratureRules.mat");
    rs = vars["Q_GaussLobatto"][N]["Points"]
    r,s = (rs[:,i] for i = 1:size(rs,2))
    w = vec(vars["Q_GaussLobatto"][N]["Weights"])
    quad_rule_face = gauss_lobatto_quad(0,0,N+1)

    return (r,s,w), quad_rule_face
end

function diagE_sbp_nodes(elem::Tri, approxType::SBP{Kubatko{LegendreFaceNodes}}, N)

    vars = matread((@__DIR__)*"/data/sbp_nodes/KubatkoQuadratureRules.mat");
    rs = vars["Q_GaussLegendre"][N]["Points"]
    r,s = (rs[:,i] for i = 1:size(rs,2))
    w = vec(vars["Q_GaussLegendre"][N]["Weights"])
    quad_rule_face = gauss_quad(0,0,N)

    return (r,s,w), quad_rule_face
end

parsevec(type, str) = str |>
  (x -> split(x, ", ")) |>
  (x -> map(y -> parse(type, y), x))

function diagE_sbp_nodes(elem::Tri, approxType::SBP{Hicken}, N)

    # from Jason Hicken https://github.com/OptimalDesignLab/SummationByParts.jl/tree/work
    lines = readlines((@__DIR__)*"/data/sbp_nodes/tri_diage_p$N.dat")
    r = parsevec(Float64,lines[11])
    s = parsevec(Float64,lines[12])
    w = parsevec(Float64,lines[13])

    # convert Hicken format to biunit right triangle
    r = @. 2*r-1
    s = @. 2*s-1
    w = 2.0 * w/sum(w)

    quad_rule_face = gauss_lobatto_quad(0,0,N+1)

    return (r,s,w), quad_rule_face
end

function build_Ef_Fmask(rd_sbp::RefElemData; tol = 100*eps())

    @unpack rq, sq, rf, sf, Nfaces = rd_sbp
    rf, sf = (x->reshape(x, length(rf) ÷ Nfaces, Nfaces)).((rf, sf))
    Fmask = zeros(Int, length(rf) ÷ Nfaces, Nfaces) #
    Ef = zeros(length(rf), length(rq)) # extraction matrix
    for i in eachindex(rq)
        for f = 1:rd_sbp.Nfaces
            id = findall(@. abs(rq[i]-rf[:,f]) + abs(sq[i]-sf[:,f]) .< tol)
            Fmask[id,f] .= i
            Ef[id .+ (f-1) * size(rf,1), i] .= 1
        end
    end
    return Ef, Fmask
end

get_face_nodes(x::AbstractVector, Fmask) = view(x, Fmask)
get_face_nodes(x::AbstractMatrix, Fmask) = view(x, Fmask, :)

function _convert_RefElemData_fields_to_SBP(rd, approxType::SBP)
    rd = @set rd.M = Diagonal(rd.wq)
    rd = @set rd.Pq = I
    rd = @set rd.Vq = I
    rd = @set rd.approximationType = approxType
    return rd
end

"""
    function hybridized_SBP_operators(rd::RefElemData{DIMS})

Constructs hybridized SBP operators given a `RefElemData`. Returns operators `Qrsth..., VhP, Ph`.
"""
function hybridized_SBP_operators(rd)
    @unpack M, Vq, Pq, Vf, wf, Drst, nrstJ = rd
    Qrst = (D->Pq' * M * D * Pq).(Drst)
    Ef = Vf * Pq
    Brst = (nJ->diagm(wf .* nJ)).(nrstJ)
    Qrsth = ((Q, B)->.5*[Q-Q' Ef'*B; -B*Ef B]).(Qrst, Brst)
    Vh = [Vq; Vf]
    Ph = M \ transpose(Vh)
    VhP = Vh * Pq
    return Qrsth, VhP, Ph, Vh
end


text:
module PINN

using Gmsh
using Gridap
#using GridapGmsh
using SparseArrays
#using Gridap.Geometry
using ChainRules
using Zygote
using LinearAlgebra
import ChainRules: rrule
import Gmsh: gmsh

export MeshGenerator
export MatrixA
export MatrixB
export pf_p
export u_pf
export g_u
export g_p
export Threshold
export p_vec
export extract_design
export Filter
export a_f

include("FilterAndThreshold.jl")
include("Helper.jl")
include("MeshGenerator.jl")
include("PML.jl")
include("Objective.jl")
include("Model.jl")

end

max_stars_repo_name: Tiamat-Tech/hpinn

max_issues_repo_name: louhz/hpinn

max_forks_repo_name: louhz/hpinn


text:
include("../TTM/ttm21f.jl")
include("/home/heindelj/dev/julia_development/wally/src/molecule_tools/read_xyz.jl")
using BenchmarkTools, DelimitedFiles, Test

header, labels, geoms = read_xyz("/home/heindelj/dev/julia_development/MolecularPotentials/Reference_Implementation/ttm-water-models-04-Jul-2013/hexamers/W20_fused_prisms.xyz")
ttm21f = TTM21F(20)
grads = zeros(3, 60)
energy = evaluate!(ttm21f, geoms[1], grads)
@btime evaluate!(ttm21f, geoms[1], grads)

correct_energy = -215.22312443016003
correct_grads  = readdlm("/home/heindelj/OneDrive/Documents/Coding_Projects/julia_development/MolecularPotentials/MolecularPotentials/src/tests/w20_grads_correct.txt", Float64)
Test.@test energy ≈ correct_energy
Test.@test grads  ≈ correct_grads
println(energy)
#display(grads)

max_stars_repo_name: heindelj/MolecularPotentials

max_issues_repo_name: heindelj/MolecularPotentials

max_forks_repo_name: heindelj/MolecularPotentials


text:
module Tracker

using MacroTools
using MacroTools: @q, @forward

import Base: ==

export TrackedArray, TrackedVector, TrackedMatrix, Params, param, back!

tracker(x) = nothing

istracked(x) = tracker(x) ≠ nothing
isleaf(x) = !istracked(x) || isleaf(tracker(x))
grad(x) = grad(tracker(x))
grad(::Nothing) = nothing
data(x) = x

struct Call{F,As<:Tuple}
  func::F
  args::As
end

Call(f::F, args::T) where {F,T} = Call{F,T}(f, args)
Call() = Call(nothing, ())

# When deserialising, the object_id changes
a::Call == b::Call = a.func == b.func && a.args == b.args

@inline (c::Call)() = c.func(data.(c.args)...)

mutable struct Tracked{T}
  ref::UInt32
  f::Call
  isleaf::Bool
  grad::T
  Tracked{T}(f::Call) where T = new(0, f, false)
  Tracked{T}(f::Call, grad::T) where T = new(0, f, false, grad)
  Tracked{T}(f::Call{Nothing}, grad::T) where T = new(0, f, true, grad)
end

istracked(x::Tracked) = true
isleaf(x::Tracked) = x.f == Call()
grad(x::Tracked) = x.grad

track(f::Call, x) = Tracked{typeof(x)}(f)

function _forward end

function track(f::F, xs...; kw...) where F
  y, back = _forward(f, xs...; kw...)
  track(Call(back, tracker.(xs)), y)
end

macro grad(ex)
  @capture(shortdef(ex), (name_(args__) = body_) |
                         (name_(args__) where {T__} = body_)) || error("Need a function definition")
  T == nothing && (T = [])
  isexpr(name, :(::)) || (name = :(::typeof($name)))
  insert!(args, 1+isexpr(args[1], :parameters) , name)
  @q(Tracker._forward($(args...)) where $(T...) = $body) |> esc
end

function update!(x, Δ)
  x.data .+= data(Δ)
  tracker(x).grad .= 0
  return x
end

include("idset.jl")
include("back.jl")
include("scalar.jl")
include("array.jl")
include("numeric.jl")

"""
    hook(f, x) -> x′

Hook into gradient backpropagation. `x` is unmodified, but when backpropagating
`f` will be applied to the incoming gradient. For example, `hook(-, x)` will reverse
the sign of the gradient applied to `x`."""
hook(f, x) = istracked(x) ? track(hook, f, x) : x
@grad hook(f, x) = data(x), Δ -> (nothing, f(Δ))

"""
    checkpoint(f, args...)

Behaves like `f(args...)`, but avoids storing the intermediate values needed for
calculating gradients. Instead, `f(args...)` will be called again during the
backward pass. This can be used to save memory in larger models.
"""
checkpoint(f, args...) = track(checkpoint, f, args...)

@grad function checkpoint(f, args...)
  data(f(args...)), function (Δ)
    y, back = forward(f, args...)
    (nothing, back(Δ)...)
  end
end

nobacksies(f, x) = track(nobacksies, f, x)
nobacksies(f, xs::Tuple) = map(x -> nobacksies(f, x), xs)
@grad nobacksies(f, x) = data(x), Δ -> error("Nested AD not defined for $f")

param(x::Number) = TrackedReal(float(x))
param(xs::AbstractArray) = TrackedArray(float.(xs))

@grad identity(x) = data(x), Δ -> (Δ,)
param(x::TrackedReal) = track(identity, x)
param(x::TrackedArray) = track(identity, x)

import Adapt.adapt

adapt(T, xs::TrackedArray) = param(adapt(T, data(xs)))

end


text:
function areaJl!{T <: FloatingPoint}(direction :: Array{T, 2}, dout :: Array{Number4{T}, 2},  LONG :: Real)
    D1, D2 = size(direction)
    SHORT = one(T)
    PI_4 = pi/4
    PI = pi

    for i in 1:length(direction)
            dir = direction[i];

            area1 = LONG*SHORT/2 * (   PI_4/2-dir - atan((SHORT-LONG)*sin(2*(  PI_4/2-dir)) / (SHORT+LONG + (SHORT-LONG)*cos(2*(  PI_4/2-dir)))));
            area2 = LONG*SHORT/2 * ( 3*PI_4/2-dir - atan((SHORT-LONG)*sin(2*(3*PI_4/2-dir)) / (SHORT+LONG + (SHORT-LONG)*cos(2*(3*PI_4/2-dir)))));
            area3 = LONG*SHORT/2 * ( 5*PI_4/2-dir - atan((SHORT-LONG)*sin(2*(5*PI_4/2-dir)) / (SHORT+LONG + (SHORT-LONG)*cos(2*(5*PI_4/2-dir)))));
            area4 = LONG*SHORT/2 * ( 7*PI_4/2-dir - atan((SHORT-LONG)*sin(2*(7*PI_4/2-dir)) / (SHORT+LONG + (SHORT-LONG)*cos(2*(7*PI_4/2-dir)))));
            area5 = LONG*SHORT/2 * ( 9*PI_4/2-dir - atan((SHORT-LONG)*sin(2*(9*PI_4/2-dir)) / (SHORT+LONG + (SHORT-LONG)*cos(2*(9*PI_4/2-dir)))));

            lps = SHORT * LONG * PI;

            s0 = (area2-area1)/lps;
            s1 = (area3-area2)/lps;
            s2 = (area4-area3)/lps;
            s3 = (area5-area4)/lps;

            dout[i] = Number4(s0, s1, s2, s3);
    end
    return dout
end


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 5
n = 5
k = 10
blockdim = 5

@testset "geam" begin
    @testset for elty in [Float32,Float64,Complex64,Complex128]
        A = sparse(rand(elty,m,n))
        B = sparse(rand(elty,m,n))
        alpha = rand(elty)
        beta = rand(elty)
        @testset "csr" begin
            d_A = CudaSparseMatrixCSR(A)
            d_B = CudaSparseMatrixCSR(B)
            C = alpha * A + beta * B
            d_C = CUSPARSE.geam(alpha,d_A,beta,d_B,'O','O','O')
            h_C = to_host(d_C)
            @test C ≈ h_C
            d_C = CUSPARSE.geam(d_A,beta,d_B,'O','O','O')
            h_C = to_host(d_C)
            C = A + beta * B
            @test C ≈ h_C
            d_C = CUSPARSE.geam(alpha,d_A,d_B,'O','O','O')
            h_C = to_host(d_C)
            C = alpha * A + B
            @test C ≈ h_C
            d_C = CUSPARSE.geam(d_A,d_B,'O','O','O')
            h_C = to_host(d_C)
            C = A + B
            @test C ≈ h_C
            d_C = d_A + d_B
            h_C = to_host(d_C)
            C = A + B
            @test C ≈ h_C
            d_C = d_A - d_B
            h_C = to_host(d_C)
            C = A - B
            @test C ≈ h_C
            B_ = sparse(rand(elty,k,n))
            d_B = CudaSparseMatrixCSR(B_)
            @test_throws DimensionMismatch CUSPARSE.geam(d_B,d_A,'O','O','O')
        end
        @testset "csc" begin
            d_A = CudaSparseMatrixCSC(A)
            d_B = CudaSparseMatrixCSC(B)
            C = alpha * A + beta * B
            d_C = CUSPARSE.geam(alpha,d_A,beta,d_B,'O','O','O')
            h_C = to_host(d_C)
            @test C ≈ h_C
            d_C = CUSPARSE.geam(d_A,beta,d_B,'O','O','O')
            h_C = to_host(d_C)
            C = A + beta * B
            @test C ≈ h_C
            d_C = CUSPARSE.geam(alpha,d_A,d_B,'O','O','O')
            h_C = to_host(d_C)
            C = alpha * A + B
            @test C ≈ h_C
            d_C = CUSPARSE.geam(d_A,d_B,'O','O','O')
            h_C = to_host(d_C)
            C = A + B
            @test C ≈ h_C
            d_C = d_A + d_B
            h_C = to_host(d_C)
            C = A + B
            @test C ≈ h_C
            d_C = d_A - d_B
            h_C = to_host(d_C)
            C = A - B
            @test C ≈ h_C
            B_ = sparse(rand(elty,k,n))
            d_B = CudaSparseMatrixCSC(B_)
            @test_throws DimensionMismatch CUSPARSE.geam(d_B,d_A,'O','O','O')
        end
        @testset "mixed" begin
            A = spdiagm(rand(elty,m))
            B = spdiagm(rand(elty,m)) + sparse(diagm(rand(elty,m-1),1))
            d_A = CudaSparseMatrixCSR(A)
            d_B = CudaSparseMatrixCSC(B)
            d_C = d_B + d_A
            h_C = to_host(d_C)
            @test h_C ≈ A + B
            d_C = d_A + d_B
            h_C = to_host(d_C)
            @test h_C ≈ A + B
            d_C = d_A - d_B
            h_C = to_host(d_C)
            @test h_C ≈ A - B
            d_C = d_B - d_A
            h_C = to_host(d_C)
            @test h_C ≈ B - A
        end
    end
end