text:
module NullableBenchmarks

include(joinpath(dirname(@__FILE__), "..", "utils", "RandUtils.jl"))

using .RandUtils
using BenchmarkTools
using Compat

const SUITE = BenchmarkGroup()

####################
# basic operations #
####################

g = addgroup!(SUITE, "basic")

for T in (Bool, Int8, Int64, Float32, Float64, BigInt, BigFloat)
    tol = (T == BigInt || T == BigFloat) ? 0.6 : 0.3

    x = Nullable(one(T))
    g["get1", string(x)] = @benchmarkable get($x) time_tolerance=tol

    for x in (Nullable(one(T)), Nullable{T}())
        g["isnull", string(x)] = @benchmarkable isnull($x) time_tolerance=tol
        g["get2", string(x)] = @benchmarkable get($x, $(zero(T))) time_tolerance=tol

        for y in (Nullable(one(T)), Nullable(zero(T)), Nullable{T}()) time_tolerance=tol
            g["isequal", string(x), string(y)] = @benchmarkable isequal($x, $y) time_tolerance=tol
        end
    end
end

####################
# nullable array   #
####################

g = addgroup!(SUITE, "nullablearray")

immutable NullableArray{T, N} <: AbstractArray{Nullable{T}, N}
    values::Array{T, N}
    hasvalue::Array{Bool, N}
end

@inline function Base.getindex{T, N}(X::NullableArray{T, N}, I::Int...)
    if isbits(T)
        ifelse(X.hasvalue[I...], Nullable{T}(X.values[I...]), Nullable{T}())
    else
        if X.hasvalue[I...]
            Nullable{T}(X.values[I...])
        else
            Nullable{T}()
        end
    end
end

Base.size(X::NullableArray) = size(X.values)

@compat Base.IndexStyle(::Type{<:NullableArray}) = IndexLinear()

const VEC_LENGTH = 1000

function perf_sum{T<:Nullable}(X::AbstractArray{T})
    S = eltype(T)
    s = zero(S)+zero(S)
    @inbounds @simd for i in eachindex(X)
        s += get(X[i], zero(S))
    end
    s
end

function perf_countnulls{T<:Nullable}(X::AbstractArray{T})
    n = 0
    @inbounds for i in eachindex(X)
        n += isnull(X[i])
    end
    n
end

function perf_countequals{T<:Nullable}(X::AbstractArray{T}, Y::AbstractArray{T})
    n = 0
    @inbounds for i in eachindex(X, Y)
        n += isequal(X[i], Y[i])
    end
    n
end

for T in (Bool, Int8, Int64, Float32, Float64, BigInt, BigFloat, Complex{Float64})
    if T == BigInt
        S = Int128
    elseif T == BigFloat
        S = Float64
    else
        S = T
    end

    # 10% of missing values
    X = NullableArray(Vector{T}(samerand(S, VEC_LENGTH)), Vector{Bool}(samerand(VEC_LENGTH) .> .9))
    Y = NullableArray(Vector{T}(samerand(S, VEC_LENGTH)), Vector{Bool}(samerand(VEC_LENGTH) .> .9))

    for (A, X, Y) in (("NullableArray", X, Y), ("Array", collect(X), collect(Y)))
        g["perf_sum", A, string(T)] = @benchmarkable perf_sum($X)
        g["perf_countnulls", A, string(T)] = @benchmarkable perf_countnulls($X)
        g["perf_countequals", A, string(T)] = @benchmarkable perf_countequals($X, $Y)
    end
end

function perf_all(X::AbstractArray{Nullable{Bool}})
    @inbounds for i in eachindex(X)
        x = X[i]
        if isnull(x)
            return Nullable{Bool}()
        elseif !get(x)
            return Nullable(false)
        end
    end
    Nullable(true)
end

function perf_any(X::AbstractArray{Nullable{Bool}})
    allnull = true
    @inbounds for i in eachindex(X)
        x = X[i]
        if !isnull(x)
            allnull = false
            get(x) && return Nullable(true)
        end
    end
    allnull ? Nullable{Bool}() : Nullable(false)
end

# Ensure no short-circuit happens
X = NullableArray(fill(true, VEC_LENGTH), fill(true, VEC_LENGTH))
# 10% of missing values
Y = NullableArray(fill(false, VEC_LENGTH), Vector{Bool}(samerand(VEC_LENGTH) .> .1))

g["perf_all", "NullableArray"] = @benchmarkable perf_all($X)
g["perf_any", "NullableArray"] = @benchmarkable perf_any($Y)

g["perf_all", "Array"] = @benchmarkable perf_all($(collect(X)))
g["perf_any", "Array"] = @benchmarkable perf_any($(collect(Y)))

for b in values(g)
    b.params.time_tolerance = 0.50
end

end # module


text:
# ---
# title: 647. Palindromic Substrings
# id: problem647
# author: Indigo
# date: 2021-06-27
# difficulty: Medium
# categories: String, Dynamic Programming
# link: <https://leetcode.com/problems/palindromic-substrings/description/>
# hidden: true
# ---
#
# Given a string, your task is to count how many palindromic substrings in this
# string.
#
# The substrings with different start indexes or end indexes are counted as
# different substrings even they consist of same characters.
#
# **Example 1:**
#
#
#
#     Input: "abc"
#     Output: 3
#     Explanation: Three palindromic strings: "a", "b", "c".
#
#
#
#
# **Example 2:**
#
#
#
#     Input: "aaa"
#     Output: 6
#     Explanation: Six palindromic strings: "a", "a", "a", "aa", "aa", "aaa".
#
#
#
#
# **Note:**
#
#   1. The input string length won't exceed 1000.
#
#
#
#
## @lc code=start
using LeetCode

function count_substrings(s::String)
    len, res = length(s), 0
    cu = codeunits(s)
    for i in 1:len, j in (0, 1)
        l, r = i, i + j
        while l >= 1 && r <= len && cu[l] == cu[r]
            res += 1
            l -= 1
            r += 1
        end
    end
    return res
end
## @lc code=end

ext: jl


text:
let
    dict = Dict()
    GraphUtils.nested_default_dict_set!(dict,1,"my val")
    @test GraphUtils.nested_default_dict_get(dict,1) == "my val"
    GraphUtils.nested_default_dict_set!(dict,1,2,"other val")
    @test GraphUtils.nested_default_dict_get(dict,1,2) == "other val"
    # Setting a value deeper in tree should not affect previous leaf
    @test GraphUtils.nested_default_dict_get(dict,1) == "my val"
    # if partial path exists, should return the leaf of that path
    @test GraphUtils.nested_default_dict_get(dict,1,2,3) == "other val"
    # if no partial path exists, should return the default kw value
    @test GraphUtils.nested_default_dict_get(dict,2;default="default val") == "default val"
    @test GraphUtils.nested_default_dict_get(dict,2,3,4;default="default val") == "default val"
    # change an existing value
    GraphUtils.nested_default_dict_set!(dict,1,2,"yet another val")
    @test GraphUtils.nested_default_dict_get(dict,1,2) == "yet another val"
    # should return default kw from a shallow point with no value
    GraphUtils.nested_default_dict_set!(dict,3,4,5,"fifth val")
    @test GraphUtils.nested_default_dict_get(dict,3;default="no val") == "no val"

end


text:
module LocalToGlobalPosNegArraysTests

using Test
using Gridap.Arrays
using Gridap.TensorValues

p = 1
data = [2,3,-1,3,4,4,3,2,5,4,-3,4]
ptrs = [1,4,4,7,13]
gid_to_val_pos = [p, 2*p, 3*p, -p, p]
gid_to_val_neg = [3*p, -p, p]
lid_to_gid = Table(data,ptrs)
ca = LocalToGlobalPosNegArray(lid_to_gid,gid_to_val_pos,gid_to_val_neg)
a = Vector{Int64}[[2, 3, 3], [], [3, -1, -1], [3, 2, 1, -1, 1, -1]]
test_array( ca, a )

i_to_v=[[1, 2], [5, 6], [1, 5], [2, 6], [2, 3], [6, 7], [3, 7], [3, 4], [7, 8], [4, 8], [9, 10], [5, 9], [6, 10], [10, 11], [7, 11], [11, 12], [8, 12], [13, 14], [9, 13], [10, 14], [14, 15], [11, 15], [15, 16], [12, 16]]
j_to_i=Int64[]
lid_to_gid=reindex(i_to_v,j_to_i)
gid_to_val=VectorValue{2,Float64}[(0.0, 0.25), (0.25, 0.25), (0.5, 0.25), (0.75, 0.25), (0.0, 0.5), (0.25, 0.5), (0.5, 0.5), (0.75, 0.5), (0.0, 0.75), (0.25, 0.75), (0.5, 0.75), (0.75, 0.75), (0.0, 1.0), (0.25, 1.0), (0.5, 1.0), (0.75, 1.0)]
f2=LocalToGlobalPosNegArray(lid_to_gid,gid_to_val,gid_to_val)
test_array(f2,Vector{VectorValue{2,Float64}}[])

end # module


text:
using Clang, Clang.LibClang
using CEnum

include("clanghelpers.jl")

ctx = DefaultContext()
trans_unit = parse_header("../deps/include/NIDAQmx.h", args=["-fparse-all-comments"],
                          flags = CXTranslationUnit_DetailedPreprocessingRecord)

ctx.libname = "libnidaqmx"
push!(ctx.trans_units, trans_unit)
root_cursor = getcursor(trans_unit)
header = spelling(root_cursor)
ctx.children = children(root_cursor)

constants = Vector{CLMacroDefinition}()
attributes = Vector{CLMacroDefinition}()
topology = Vector{CLMacroDefinition}()
remaining = Vector{CLCursor}()

# Type aliasing map used by NIDAQmx.h
typedef_map = Dict([("int8", "Cchar"),
                    ("uInt8", "Cuchar"),
                    ("int16", "Cshort"),
                    ("uInt16", "Cushort"),
                    ("int32", "Cint"),
                    ("uInt32", "Cuint"),
                    ("float32", "Cfloat"),
                    ("float64", "Cdouble"),
                    ("int64", "Clonglong"),
                    ("uInt64", "Culonglong"),
                    ("bool32", "Cuint")
                   ])

attribute_map = Dict([("Buf", "Buffer"),
                      ("AI", "AIChannel"),
                      ("AO", "AOChannel"),
                      ("DI", "DIChannel"),
                      ("DO", "DOChannel"),
                      ("CI", "CIChannel"),
                      ("CO", "COChannel"),
                      ("Dev", "Device"),
                      ("Exported", "ExportedSignal"),
                      ("PersistedChan", "PersistedChannel"),
                      ("PersistedScale", "PersistedScale"),
                      ("PersistedTask", "PersistedTask"),
                      ("PhysicalChan", "PhysicalChannel"),
                      ("Read", "Read"),
                      ("Write", "Write"),
                      ("RealTime", "RealTime"),
                      ("Scale", "Scale"),
                      ("Task", "Task")
                     ])

blacklist = ["unix","linux", "CVIDECL", "CVITime_DECLARED", "CVICDECL",
             "DAQmxSuccess", "DAQmx_ReadWaitMode",
             "CVIAbsoluteTime_DECLARED", "TRUE", "FALSE", "NULL", "DAQmxFailed",
             "CVICALLBACK"]

neg_constants = ["DAQmx_Val_Cfg_Default", "DAQmx_Val_Default", "DAQmx_Val_Auto"]

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

for (i, child) in enumerate(ctx.children)

    # Cursor properties
    child_name = name(child)
    child_kind = kind(child)
    ctx.children_index = i

    # Skip compiler constants/garbage and type aliasing
    startswith(child_name, "_") && continue
    startswith(child_name, "DAQmx_Switch") && continue
    child_name ∈ blacklist && continue
    child_kind == CXCursor_TypedefDecl && continue
    child_name ∈ name.(constants) && continue
    child_name ∈ name.(attributes) && continue

    # Group macro constants by type (later transformed into Julia enums)
    if child_kind == CXCursor_MacroDefinition
        if startswith(child_name, "DAQmxError") || startswith(child_name, "DAQmxWarning")
            continue
        elseif startswith(child_name, "DAQmx_Val_Switch_Topology_")
            push!(topology, child)
            continue
        elseif startswith(child_name, "DAQmx_Val_")
            if child_name == "DAQmx_Val_WaitInfinitely"
                wrap!(ctx, child)
                continue
            end
            push!(constants, child)
            continue
        elseif startswith(child_name, "DAQmx_")
            push!(attributes, child)
            continue
        else
            push!(remaining, child)
            continue
        end
    end
    if child_kind == CXCursor_FunctionDecl
        wrap_dealias!(ctx, child)
        continue
    else
        #skip struct/enum and handful of type declarations (implement manually)
        push!(remaining, child)
    end
end

subattributes = Dict{String, Vector{CLMacroDefinition}}()
for prefix in values(attribute_map)
    push!(subattributes, prefix => Vector{CLMacroDefinition}())
end

sorted_attributes = sort_cursors(attributes)
sorted_constants = sort_cursors(constants)
redundant = []

for (i, cursor) in enumerate(sorted_attributes)
    cursor_name = name(cursor)
    key = split(cursor_name, "_")[2]
    if key ∈ keys(attribute_map)
        push!(subattributes[attribute_map[key]], sorted_attributes[i])
        push!(redundant, i)
    end
end

deleteat!(sorted_attributes, redundant)

for key in keys(subattributes)
    wrap_macro2enum!(ctx, subattributes[key], key)
end

wrap_macro2enum!(ctx, sorted_attributes, "DAQmxAttribute")
wrap_macro2enum!(ctx, sorted_constants, "DAQmxConstant")

api_file = joinpath(@__DIR__, "libdaqmx_api.jl")
open(api_file, "w") do f
    println(f, "# Automatically generated using Clang.jl\n")
    print_buffer(f, ctx.api_buffer)
end

common_file = joinpath(@__DIR__, "libdaqmx_common.jl")
open(common_file, "w") do f
    println(f, "# Automatically generated using Clang.jl\n")
    print_buffer(f, dump_to_buffer(ctx.common_buffer))
end

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


text:
function selectionsort(arr)
    N = length(arr)

    for step in 0:N-2
        m = -1
        idx = -1

        for j in 1:N-step
            if arr[j] > m
                m = arr[j]
                idx = j
            end
        end

        s = arr[N-step]
        arr[N-step] = m
        arr[idx] = s
    end
end

v = [144 89 4 9 95 12 86 25]
println("unsorted array = ",v)

selectionsort(v)

println("selection sorted array = ",v)

max_stars_repo_name: landrumb/julia

max_issues_repo_name: landrumb/julia

max_forks_repo_name: landrumb/julia


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

mutable struct LaTeX
    formula::String
end

@trigger '$' ->
function tex(stream::IO, md::MD)
    result = parse_inline_wrapper(stream, "\$", rep = true)
    return result === nothing ? nothing : LaTeX(result)
end

function blocktex(stream::IO, md::MD)
    withstream(stream) do
        ex = tex(stream, md)
        if ex ≡ nothing
            return false
        else
            push!(md, ex)
            return true
        end
    end
end

show(io::IO, tex::LaTeX) =
    print(io, '$', tex.formula, '$')

latex(io::IO, tex::LaTeX) =
    println(io, "\$\$", tex.formula, "\$\$")

latexinline(io::IO, tex::LaTeX) =
    print(io, '$', tex.formula, '$')

term(io::IO, tex::LaTeX, cols) = println_with_format(:magenta, io, tex.formula)
terminline(io::IO, tex::LaTeX) = print_with_format(:magenta, io, tex.formula)


text:
using God
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# problem definition
problem1(x) = sum( x.^2 )

# solver
solver1(problem) = problem(0:10)

# set of problem
problemSet = Array{Function, 1}([problem1])
solverSet = Array{Function, 1}([solver1])

# test cases
@test solve(problem1, solver1)    == 385
@test  385  in solve(problemSet, solver1)
@test [385] in solve(problemSet, solverSet)


text:
a = 2 + 2
println("Hello world")


text:
using FFTW, FastTransforms, LinearAlgebra, Test

FFTW.set_num_threads(ceil(Int, Sys.CPU_THREADS/2))

@testset "Nonuniform fast Fourier transforms" begin
    function nudft1(c::AbstractVector, ω::AbstractVector{T}) where {T<:AbstractFloat}
        # Nonuniform discrete Fourier transform of type I

        N = size(ω, 1)
        output = zero(c)
        for j = 1:N
        	output[j] = dot(exp.(2*T(π)*im*(j-1)/N*ω), c)
        end

        return output
    end

    function nudft2(c::AbstractVector, x::AbstractVector{T}) where {T<:AbstractFloat}
        # Nonuniform discrete Fourier transform of type II

        N = size(x, 1)
        output = zero(c)
        ω = collect(0:N-1)
        for j = 1:N
        	output[j] = dot(exp.(2*T(π)*im*x[j]*ω), c)
        end

        return output
    end

    function nudft3(c::AbstractVector, x::AbstractVector{T}, ω::AbstractVector{T}) where {T<:AbstractFloat}
        # Nonuniform discrete Fourier transform of type III

        N = size(x, 1)
        output = zero(c)
        for j = 1:N
            output[j] = dot(exp.(2*T(π)*im*x[j]*ω), c)
        end

        return output
    end

    N = round.([Int],10 .^ range(1,stop=3,length=10))

    for n in N, ϵ in (1e-4, 1e-8, 1e-12, eps(Float64))
        c = complex(rand(n))
        err_bnd = 500*ϵ*n*norm(c)

        ω = collect(0:n-1) + 0.25*rand(n)
        exact = nudft1(c, ω)
        fast = nufft1(c, ω, ϵ)
        @test norm(exact - fast, Inf) < err_bnd

        d = inufft1(fast, ω, ϵ)
        @test norm(c - d, Inf) < err_bnd

        x = (collect(0:n-1) + 0.25*rand(n))/n
        exact = nudft2(c, x)
        fast = nufft2(c, x, ϵ)
        @test norm(exact - fast, Inf) < err_bnd

        d = inufft2(fast, x, ϵ)
        @test norm(c - d, Inf) < err_bnd

        exact = nudft3(c, x, ω)
        fast = nufft3(c, x, ω, ϵ)
        @test norm(exact - fast, Inf) < err_bnd
    end

    # Check that if points/frequencies are indeed uniform, then it's equal to the fft.
    for n in (1000,), ϵ in (eps(Float64), 0.0)
        c = complex(rand(n))
        ω = collect(0.0:n-1)
        x = ω/n
        fftc = fft(c)
        if Sys.WORD_SIZE == 64
            @test_skip norm(nufft1(c, ω, ϵ) - fftc) == 0 # skip because fftw3 seems to change this
            @test norm(nufft2(c, x, ϵ) - fftc) == 0
            @test_skip norm(nufft3(c, x, ω, ϵ) - fftc) == 0 # skip because fftw3 seems to change this
        end
        err_bnd = 500*eps(Float64)*norm(c)
        @test norm(nufft1(c, ω, ϵ) - fftc) < err_bnd
        @test norm(nufft2(c, x, ϵ) - fftc) < err_bnd
        @test norm(nufft3(c, x, ω, ϵ) - fftc) < err_bnd
    end

    function nudft1(C::Matrix{Complex{T}}, ω1::AbstractVector{T}, ω2::AbstractVector{T}) where {T<:AbstractFloat}
        # Nonuniform discrete Fourier transform of type I-I

        M, N = size(C)
        output = zero(C)
        @inbounds for j1 = 1:M, j2 = 1:N
            for k1 = 1:M, k2 = 1:N
                output[j1,j2] += exp(-2*T(π)*im*((j1-1)/M*ω1[k1]+(j2-1)/N*ω2[k2]))*C[k1,k2]
            end
        end
        return output
    end

    function nudft2(C::Matrix{Complex{T}}, x::AbstractVector{T}, y::AbstractVector{T}) where {T<:AbstractFloat}
        # Nonuniform discrete Fourier transform of type II-II

        M, N = size(C)
        output = zero(C)
        @inbounds for j1 = 1:M, j2 = 1:N
            for k1 = 1:M, k2 = 1:N
                output[j1,j2] += exp(-2*T(π)*im*(x[j1]*(k1-1)+y[j2]*(k2-1)))*C[k1,k2]
            end
        end
        return output
    end

    N = round.([Int],10 .^ range(1,stop=1.7,length=5))

    for n in N, ϵ in (1e-4,1e-8,1e-12,eps(Float64))
        C = complex(rand(n,n))
        err_bnd = 500*ϵ*n*norm(C)

        x = (collect(0:n-1) + 0.25*rand(n))/n
        y = (collect(0:n-1) + 0.25*rand(n))/n
        ω1 = collect(0:n-1) + 0.25*rand(n)
        ω2 = collect(0:n-1) + 0.25*rand(n)

        exact = nudft1(C, ω1, ω2)
        fast = nufft1(C, ω1, ω2, ϵ)
        @test norm(exact - fast, Inf) < err_bnd

        exact = nudft2(C, x, y)
        fast = nufft2(C, x, y, ϵ)
        @test norm(exact - fast, Inf) < err_bnd
    end
end


text:
function get_experiment_id(uri::String)
    experiment_id = ""
    try
        experiment_id = HTTP.URIs.splitpath(uri)[2]
    catch
        throw(ArgumentError("Could not determine experiment ID in URI '$uri'"))
    end
    log_message("extracted experiment id $experiment_id")
    experiment_id
end

get_function_target(uri::String) = HTTP.URIs.splitpath(uri)[4]
get_object_id(uri::String) = HTTP.URIs.splitpath(uri)[5]

create_response(s::String) = collect(codeunits(s))


text:
# # Speeding up Quantum Mechanics - Matrix Types
# ### Carsten Bauer

# ## Free Fermions on a Chain
#
# $$\mathcal{H} = -t\sum_{\langle i,j \rangle} c_i^\dagger c_j + \mu \sum_i n_i$$
#
# Here, $t$ is the hopping amplitude, $\mu$ is the chemical potential, and $c, c^\dagger$ are creation and annihilation operators.
# For the sake of the argument of this tutorial, we'll consider **open boundary conditions**.
# Since the fermions are *not* interacting, we can work in the *single particle basis* and do not have to worry about how to construct a basis for the many-body Fock space.
# We use the canonical cartesian basis in which one uses $0$s to indicate empty sites and a $1$ for the particle's site, i.e. $|00100\rangle$ represents the basis state which has the particle exclusively on the 3rd site.
# If you aren't familiar with second quantization just think of $\mathcal{H}$ as any quantum mechanical operator that can be represented as a matrix.
#
# Let's build the Hamiltonian matrix.

using LinearAlgebra # makes Julia speak linear algebra fluently

N = 100 # number of sites
t = 1
μ = -0.5

H = diagm(0 => fill(μ, N), 1 => fill(-t, N-1), -1 => fill(-t, N-1))

# A typical thing one does in quantum mechanics is calculate the expectation value of an operator, like $\mathcal{H}$, with respect to some quantum mechanical state $\psi$.

ψ = normalize(rand(N)); # some state

#-
"""
    ev(O, ψ)

Calculates the quantum mechanical expectation value `<ψ|O|ψ>`, where `O` is an operator and `ψ` is a state.
"""
ev(O, ψ) = ψ'*O*ψ;
#-

ev(H, ψ) # calculate the energy

# Let's see how long this calculation takes.

using BenchmarkTools
@btime ev($H, $ψ);

# ### Utilizing the matrix structure
#
# Since `typeof(H) == Matrix{Float64}`, we have so far dispatched to generic matrix multiplications (and transposition) in the `ev` function.
# What we should do is exploit the *structure* of our Hamiltonian and use specialized, faster implementations of these operations. After all, most entries in $\mathcal{H}$ are equal to zero and multiplications involving those zeros are clearly unnecessary.
# To speed up our expectation value computation, we can tell Julia about the *sparsity* of our Hamiltonian. The way we do this, is by changing `H`'s type.

using SparseArrays
Hsparse = sparse(H)
#-
typeof(Hsparse)

# As we can see, `Hsparse` is now a sparse matrix. Let's see how long our computation takes after this simple type change.

@btime ev($Hsparse, $ψ);

# That's a solid **30x speedup**! Utilizing matrix structure apparently pays off a lot!
#
# Inspecting our Hamiltonian more closely, we note that it is not just sparse but actually *tridiagonal*.

# ---
# #### Interlude: Visualizing matrix structure
#
# A nice visualization of the sparsity pattern of a matrix can be obtained by calling the function `spy` of the [UnicodePlots.jl](https://github.com/Evizero/UnicodePlots.jl) package. (Note that the `canvas=DotCanvas` option is necessary due to some rendering issue in Jupyter.)

using UnicodePlots
spy(H, canvas=DotCanvas) # sparsity pattern plot

# ---

# Let's exploit this fact and benchmark the expectation value computation once again.

Htri = Tridiagonal(H)
#-
@btime ev($Htri, $ψ);

# That's almost **another 2x speedup**!

# What we learn from this is that we can tell Julia about the structure of our Hamiltonian through types. Instead of using generic (and slow) matrix operations, we then dispatch to fast special implementations that have been optimized for this particular structure.
#
# On a final note, choosing the best type (and therewith an algorithm) can be tricky and one has to play around a bit. The good thing is that it's very easy to try out different types!
# Note that apart from Julia's built-in matrix types there are also great custom types available in the ecosystem. See [JuliaMatrices](https://github.com/JuliaMatrices) for more information.
#
# # Core message of this tutorial
#
# * **Indicate the structure of a matrix**, like hermiticity or sparsity, through types. Fallback to generic types only if you run into method errors.

max_issues_repo_name: JuliaPhysics/PhysicsTutorials

max_forks_repo_name: JuliaPhysics/PhysicsTutorials


text:
#=

Julia is great because you don't have to specify types in many situations, but also, if you need or want, you can specify the types and they will be enforced.
Each time a function is used, the compiler compiles a specific version for that given signature.
To provide specialized code, the compiler needs to figure out the types of everything from the types of the input.
If that is not possible, the code is said to be "type unstable", and the performance of the code is hit.

=#
function bad_relu(x)
  if x > 0
    return x
  else
    return 0
  end
end
#=
```
@code_warntype bad_relu(1.0)
```
=#
using InteractiveUtils # hide
InteractiveUtils.@code_warntype bad_relu(1.0) # hide
#=
You should be able to see in color in your REPL the issue with this function (`::Union{Float64, Int64}`).
The compiler can't determine the type of the output given an output of type `Float64`.
If `x` is positive, then the return is of type `Float64`. Otherwise, the type will be an `Int`.

Figuring out whether there are type instabilities and fixing them are one of the main difficulties of Julia.
It is also one of the main reasons that the compiler is always being improved.
=#

max_stars_repo_name: chunggnuhc/julia-tutorial-at-nlesc

max_issues_repo_name: chunggnuhc/julia-tutorial-at-nlesc

max_forks_repo_name: chunggnuhc/julia-tutorial-at-nlesc


text:
using Test
using Yao, Yao.Blocks, QuAlgorithmZoo

@testset "pairs geometries" begin
    @test pair_ring(3) == [1=>2,2=>3,3=>1]
    ps = pair_square(2, 2)
    @test length(ps) == 8
    for item in [1=>2, 3=>4, 2=>1, 4=>3, 1=>3, 2=>4, 3=>1, 4=>2]
        @test item in ps
    end
    @test cnot_entangler(4, ps) isa ChainBlock
    @test cnot_entangler(4, ps) |> length == 8
end

@testset "random circuit" begin
    c = rand_circuit(1)
    @test c isa ChainBlock
    @test length(c) == 5
    c = rand_circuit(9)
    @test c isa ChainBlock
    @test length(c) == 45
end

@testset "rotter, collect_rotblocks, num_gradient, opgrad" begin
    @test  merged_rotor(true, true) == Rx(0)
    @test  merged_rotor(false, false) == merged_rotor() == chain(Rz(0), Rx(0), Rz(0))
    @test  merged_rotor(false, true) == chain(Rz(0), Rx(0))
    @test  merged_rotor(true, false) == chain(Rx(0), Rz(0))
    @test collect(rotorset(:Merged, 5, true, false), RotationGate) |> length == 10

    @test rotor(5, 2, true, true) isa ChainBlock
    @test rotor(5, 2, true, true) |> length == 1
    @test rotor(5, 2, true, true) |> nqubits == 5
    @test collect(rotorset(:Split, 5, true, false), PutBlock{<:Any, <:Any, <:RotationGate}) |> length == 10
end

@testset "random diff circuit" begin
    c = random_diff_circuit(4, 3, [1=>3, 2=>4, 2=>3, 4=>1])
    rots = collect(c, RotationGate)
    @test length(rots) == nparameters(c) == 40
    @test dispatch!(+, c, ones(40)*0.1) |> parameters == ones(40)*0.1
    @test dispatch!(+, c, :random) |> parameters != ones(40)*0.1

    nbit = 4
    c = random_diff_circuit(nbit, 1, pair_ring(nbit), mode=:Split) |> autodiff(:BP)
    reg = rand_state(4)
    dispatch!(c, randn(nparameters(c)))

    dbs = collect(c, BPDiff)
    op = kron(4, 1=>Z, 2=>X)
    loss1z() = expect(op, copy(reg) |> c)  # return loss please

    # back propagation
    ψ = copy(reg) |> c
    δ = copy(ψ) |> op
    backward!(δ, c)
    bd = gradient(c)

    # get num gradient
    nd = numdiff.(loss1z, dbs)
    ed = opdiff.(()->copy(reg)|>c, dbs, Ref(op))

    @test isapprox.(nd, ed, atol=1e-4) |> all
    @test isapprox.(nd, bd, atol=1e-4) |> all
end


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

using PDMesh

x, v, y, vol, type_ = unpack(create(Cone(10.0, 20.0), resolution=0.5, rand_=0.0)

write_data("./examples/ovito.data", x, vol, ones(Int64, size(x)[2]))

max_stars_repo_name: sajidmannan/PDMesh

max_issues_repo_name: sajidmannan/PDMesh

max_forks_repo_name: sajidmannan/PDMesh


text:
using NbodyGradient, DelimitedFiles
include("../src/integrator/convert.jl")

fname = "elements.txt"
el = readdlm(fname,',')

ic = ElementsIC(el,8,0.0)
s = State(ic)

elements = get_orbital_elements(s,ic)

println(isapprox(elements,ic.elements))


text:
include(string(PROJECT_HOME,"/simulation_code/updatefuns.jl"))
include(string(PROJECT_HOME,"/simulation_code/create_hydraulics_layered.jl"))

function make_layers(_totaldepth)
	linear_layers = collect(0:0.1:0.8);
	toadd = [0,0,1,2,3,4,5,6,7];
	quad_layers = cumsum(toadd)/sum(toadd);

	if _totaldepth < 0.8
		lower_layer_thick = (_totaldepth-0.1)/Int(floor((_totaldepth-0.1) / 0.1));
		layers = vcat(0, 0.1:lower_layer_thick:_totaldepth);
	end

	if _totaldepth >= 0.8
		layers = linear_layers + quad_layers*(_totaldepth-0.8);
	end
	return convert(Array{FT}, -1*layers)
end

function create_spac(
            sm::AbstractStomatalModel = OSMWang{FT}(),
            vcmax::FT = FT(30),
            kmax::FT = FT(0.3),
		z_soil::FT = FT(1000),
     chl::FT = FT(40),
	 N_slice::Int = 1,
	alpha::FT = FT(1.368),
	n::FT = FT(2.6257)
) where {FT<:AbstractFloat}
(

	#_soil_hs = create_soil_VC(VanGenuchten{FT}(), "Silt Loam");

	_soil_hs = VanGenuchten{FT}(stype = "Ozark",
									α = alpha,
									n = n,
								   Θs = 0.45,
								   Θr = 0.067);

	#_soil_hs = BrooksCorey{FT}(stype= "Ozark", Θs = FT(porosity), Θr = FT(0.05), ϕs = psi_sat, b = b_soil);

	_totaldepth = FT(z_soil/1000);
	_rootdepth = FT(z_soil/1000);
	#_soil_bounds = collect(FT,[0,-0.1,-0.25,-0.45,-0.7]);
	#_soil_bounds = collect(FT,[0,-0.1,-0.2,-0.3,-0.4]);
	#_soil_bounds = collect(FT,[0,-0.1,-0.2,-0.3,-0.4,-0.5,-0.6,-0.7]);

	#_soil_bounds = collect(FT,[0,-0.1,-0.15,-0.2,-0.25,-0.3,-0.35,-0.4]);
	#_soil_bounds[3:end] *= _totaldepth/FT(0.4);

	#_soil_bounds = collect(FT,[0,-0.1,-0.25,-0.5,-0.75,-1.0,-1.5,-2]);
	#_soil_bounds *= _totaldepth/FT(2);

	#_soil_bounds = collect(FT,[0,-0.1,-0.2, -0.3, -0.4, -0.5, -0.6, -0.7]);
	#_soil_bounds *= _totaldepth/FT(0.7);

	_soil_bounds = make_layers(_totaldepth);

	_totaldepth = -1*_soil_bounds[end];
	_rootdepth = -1*_soil_bounds[end];

	_air_bounds = convert(Array{FT},[0,5,9,12,16,18.5,20] );

	_tree_hs = create_tree2(FT(-1*_rootdepth), FT(9), FT(18.5), _soil_bounds,
						   _air_bounds, N_slice);
	for _root in _tree_hs.roots
		_root.sh = deepcopy(_soil_hs);
	end;

	_node = SPACMono{FT}(soil_bounds = _soil_bounds,
						  air_bounds = _air_bounds,
							  z_root = FT(-1*_rootdepth),
							z_canopy = 18.5,
							plant_hs = _tree_hs,
								  ga = 413.223,
								  la = 1735.537,
							latitude = 38.74,
						   longitude = -92.2,
						   elevation = 219.4,
					   stomata_model = sm);
	_node.swc = collect(FT,0.4*ones(length(_soil_bounds)-1));

	n_canopy = length(_node.plant_hs.canopy_index_in_air);

	#=
	incl_grid =  collect(FT,[22.5,45+22.5]);
	azi_grid = collect(FT,[90,270]);
	incl_grid_bnd = [collect(FT,[0,45]) collect(FT,[45,90])];
	=#

	incl_grid =  collect(FT,[45]);
	azi_grid = collect(FT,[180]);
	incl_grid_bnd = [collect(FT,[0]) collect(FT,[90])];

	nleaf = length(incl_grid)*length(azi_grid) + 1;

	_node.plant_ps = [CanopyLayer{FT}(n_leaf=nleaf) for
										  i in 1:n_canopy];

	_node.canopy_rt = Canopy4RT{FT}(nLayer=n_canopy, LAI=_node.la/_node.ga, litab  = incl_grid, litab_bnd = incl_grid_bnd, lazitab = azi_grid);

	#"RT dimensions"
	_node.rt_dim = create_rt_dims(_node.canopy_rt, _node.wl_set);
	#"CanopyRads container"
	_node.can_rad = create_canopy_rads(FT, _node.rt_dim);
	#"CanopyOpticals container"
	_node.can_opt = create_canopy_opticals(FT, _node.rt_dim);
	#"Array of LeafBios container"
	_node.leaves_rt = [create_leaf_bios(FT, _node.rt_dim) for
										 i in 1:n_canopy];

	#"RT container"
	_node.rt_con = create_rt_cache(FT, _node.rt_dim);
	#"Container for sunlit leaf area fraction in each layer"
	_node.f_SL = repeat(_node.canopy_rt.lidf, outer=[ _node.canopy_rt.nAzi ]) /
						_node.canopy_rt.nAzi;

	#=
	for _iPS in _node.plant_ps
		_iPS.g_min   = 0.001;
		_iPS.g_min25 = 0.001;
		_iPS.g_max   = 0.08; #0.50;
		_iPS.g_max25 = 0.08; #0.50;
	end;
	=#

	for _iPS in _node.plant_ps
		_iPS.g_min   = 0.001;
		_iPS.g_min25 = 0.001;
		_iPS.g_max   = 0.12; #0.50;
		_iPS.g_max25 = 0.12; #0.50;
	end;

	_node.canopy_rt.Ω       = 0.69;
	_node.canopy_rt.clump_a = 0.69;
	update_LAI!(_node, FT(4.2));
	update_Weibull!(_node, FT(5.703), FT(0.953));

	#update_Weibull!(_node, FT(4.0), FT(3.0));
	#for _leaf in _node.plant_hs.leaves
        #_leaf.vc.b = 1.0;
        #_leaf.vc.c = 3.0;
    #end;

	#for _leaf in _node.plant_hs.leaves
    #    _leaf.vc.b = 2.0;
        #_leaf.vc.c = c;
    #end;

    # update fitted Vcmax, Kmax, and Chlrophyll content
    update_VJRWW!(_node, vcmax);
    update_Kmax!(_node, kmax);
    update_Cab!(_node, chl);
    initialize_spac_canopy!(_node);

    return _node
)
end


text:
# using Pkg
# Pkg.add("Match")
# using Match

function nl(al)
    templ = [];
    for a in al
        if typeof(a) <: Array
            append!(templ, nl(a))
        else
            push!(templ, a)
        end
    end
    templ
end

function nl2(al)
    function innerfunc(nl, last) :: Array
        if length(nl) == 0
            return last
        end

        if typeof(nl[1]) <: Array
            innerfunc(nl[2:end], innerfunc(nl[1],last))
        else
            innerfunc(nl[2:end], push!(last, nl[1]))
        end
    end

    innerfunc(al,[])
end

function nlFromLisp(nl)
    if length(nl) == 0
        []
    elseif typeof(nl[1]) <: Array
        append!(nlFromLisp(nl[1]), nlFromLisp(nl[2:end]))
    else
        append!([nl[1]], nlFromLisp(nl[2:end]))
    end
end

using Test

@test nl([1, [2, 3, [4, 5], [6, 7, [8]], [9]], 10]) == [1,2,3,4,5,6,7,8,9,10]
@test nl2([1, [2, 3, [4, 5], [6, 7, [8]], [9]], 10]) == [1,2,3,4,5,6,7,8,9,10]
@test nlFromLisp([1, [2, 3, [4, 5], [6, 7, [8]], [9]], 10]) == [1,2,3,4,5,6,7,8,9,10]

max_stars_repo_name: ccqpein/Arithmetic-Exercises

max_issues_repo_name: ccqpein/Arithmetic-Exercises

max_forks_repo_name: ccqpein/Arithmetic-Exercises


text:
# Note on genome representation - p_int is a vector of integers of length n_gene,encoding gene type.
# p_float is an n_gene x 8 array encoding the real values for each gene (each row is a gene)
# Hence p_int[1] is the associated type of gene p_float[1,:]

function initialise_grn(p_float::Array{Float64,2},p_int::Vector{Int64},grn_config::FGRNConfig)

    n_rec = 0
    n_beh = 0
    n_reg = 0
    n_env = 0

    n = grn_config.protein_size

    for i in 1:size(p_int,1)

        if p_int[i] == 1
            n_rec += 1
        elseif p_int[i] == 2
            n_beh += 1
        elseif p_int[i] == 3
            n_reg += 1
        else
            n_env +=1
        end
    end

    grn = ( rec_coding = ones(Int,n,n,n_rec),
            beh_coding = ones(Int,n,n,n_beh),
            reg_coding = ones(Int,n,n,n_reg),
            env_coding = ones(Int,n,n,n_env),

            rec_prom = ones(Int,n,n,n_rec),
            beh_prom = ones(Int,n,n,n_beh),
            reg_prom = ones(Int,n,n,n_reg),
            env_prom = ones(Int,n,n,n_env),

            rec_cod_thresh = zeros(Float64,n_rec),
            rec_aff_thresh = zeros(Float64,n_rec),
            rec_id = ones(Int,n_rec),

            beh_cod_thresh  = zeros(Float64,n_beh),
            beh_aff_thresh = zeros(Float64,n_beh),
            beh_id = ones(Int,n_beh),

            reg_cod_thresh  = zeros(Float64,n_reg),
            reg_aff_thresh = zeros(Float64,n_reg),
            reg_id = ones(Int,n_reg),

            env_cod_thresh = zeros(Float64,n_env),
            env_aff_thresh = zeros(Float64,n_env),
            env_id = ones(Int,n_env),

            Ct = grn_config.Ct, # assume average affinity = 0.3 ish, average affinity th - 0.5 as uniform
            Cs = grn_config.Cs, # assume average affinity = 0.3 ish, average affinity th - 0.5 as uniform
            Cw = grn_config.Cw,
            Cp = grn_config.Cp,
            Ci = grn_config.Ci,
           )

    n_rec = 0
    n_beh = 0
    n_reg = 0
    n_env = 0

    for i in 1:size(p_int,1)

        coding = mandelbrot_protein(p_float[i,1:3],grn_config.mandel_thresh,grn_config.max_iter,n)
        promoter = mandelbrot_protein(p_float[i,4:6],grn_config.mandel_thresh,grn_config.max_iter,n)
        threshold = p_float[i,7]
        affinity = p_float[i,8]

        if p_int[i] == 1
            n_rec += 1
            grn.rec_coding[:,:,n_rec] = coding
            grn.rec_prom[:,:,n_rec] = promoter
            grn.rec_cod_thresh[n_rec] = threshold
            grn.rec_aff_thresh[n_rec] = affinity
            grn.rec_id[n_rec] = i

        elseif p_int[i] == 2
            n_beh += 1
            grn.beh_coding[:,:,n_beh] = coding
            grn.beh_prom[:,:,n_beh] = promoter
            grn.beh_cod_thresh[n_beh] = threshold
            grn.beh_aff_thresh[n_beh] = affinity
            grn.beh_id[n_beh] = i

        elseif p_int[i] == 3
            n_reg += 1
            grn.reg_coding[:,:,n_reg] = coding
            grn.reg_prom[:,:,n_reg] = promoter
            grn.reg_cod_thresh[n_reg] = threshold
            grn.reg_aff_thresh[n_reg] = affinity
            grn.reg_id[n_reg] = i

        else
            n_env += 1
            grn.env_coding[:,:,n_env] = coding
            grn.env_prom[:,:,n_env] = promoter
            grn.env_cod_thresh[n_env] = threshold
            grn.env_aff_thresh[n_env] = affinity
            grn.env_id[n_env] = i
        end
    end

    return grn

end

function calculate_affinity(protein_a::Array{Int,2}, protein_b::Array{Int,2})

    return convert(Float64,sum(@. abs(protein_a - protein_b))) / (size(protein_a,1)*size(protein_a,1)*255)

end

function bind_probability(affinity::Float64,affinity_th::Float64,config::FGRNConfig)

    return (1 + tanh((affinity_th - affinity - config.Ct)/config.Cs))/2

end

function update_rule(product_conc::Float64,original_conc::Float64,config::FGRNConfig)

    return product_conc*tanh((product_conc - 0.5)/config.Cw)/config.Ci - original_conc/config.Cp

end

function updateConcentrations!(dc::Array{Float64,1},grn::NamedTuple, concentrations::Array{Float64,1},grn_config::FGRNConfig)

    product = zeros(Int,grn_config.protein_size,grn_config.protein_size)
    product_conc = zeros(Float64,grn_config.protein_size,grn_config.protein_size)

    fill!(dc, zero(Float64))

    nProteins::Int = 0
    affinity = 0

    @inbounds for j in 1:grn_config.protein_size
        @inbounds for i in 1:grn_config.protein_size

            for k in 1:size(grn.rec_coding,3)

                if (concentrations[grn.rec_id[k]] > 0.5*grn.rec_cod_thresh[k]) && (grn.rec_coding[i,j,k] > product[i,j])
                    product[i,j] = grn.rec_coding[i,j,k]
                    product_conc[i,j] = concentrations[grn.rec_id[k]]

                    nProteins+=1

                end
            end

            for k in 1:size(grn.reg_coding,3)

                if (concentrations[grn.reg_id[k]] > 0.5*grn.reg_cod_thresh[k]) && (grn.reg_coding[i,j,k] > product[i,j])
                    product[i,j] = grn.reg_coding[i,j,k]
                    product_conc[i,j] = concentrations[grn.reg_id[k]]

                    nProteins+=1

                end
            end

            for k in 1:size(grn.env_coding,3)

                if (concentrations[grn.env_id[k]] > 0.5*grn.env_cod_thresh[k]) && (grn.env_coding[i,j,k] > product[i,j])
                    product[i,j] = grn.env_coding[i,j,k]
                    product_conc[i,j] = concentrations[grn.env_id[k]]

                    nProteins+=1

                end

            end
        end
    end

    product_conc = sum(product_conc)/(grn_config.protein_size*grn_config.protein_size)

    if nProteins != 0

        @inbounds for k in 1:size(grn.rec_prom,3)

            dc[grn.rec_id[k]] = 0

        end

        @inbounds for k in 1:size(grn.beh_prom,3)

            affinity = calculate_affinity(grn.beh_prom[:,:,k],product)

            if affinity > grn.beh_aff_thresh[k]

                dc[grn.beh_id[k]] = update_rule((1 - affinity)*product_conc,concentrations[grn.beh_id[k]],grn_config)

            else

                dc[grn.beh_id[k]] = 0

            end
        end

        @inbounds for k in 1:size(grn.reg_prom,3)

            affinity = calculate_affinity(grn.reg_prom[:,:,k],product)

            if affinity > grn.reg_aff_thresh[k]

                dc[grn.reg_id[k]] = update_rule((1 - affinity)*product_conc,concentrations[grn.reg_id[k]],grn_config)

            else

                dc[grn.reg_id[k]] = 0
            end
        end

        @inbounds for k in 1:size(grn.env_prom,3)

            affinity = calculate_affinity(grn.env_prom[:,:,k],product)

            if affinity > grn.env_aff_thresh[k]

                dc[grn.env_id[k]] = update_rule((1 - affinity)*product_conc,concentrations[grn.env_id[k]],grn_config)

            else

                dc[grn.env_id[k]] = 0
            end
        end

    else
        fill!(dc, zero(Float64))

    end

    @. concentrations = concentrations + dc;

end

max_stars_repo_name: harrybooth/FractalDevelopment

max_issues_repo_name: harrybooth/FractalDevelopment

max_forks_repo_name: harrybooth/FractalDevelopment


text:
using Test
using CLIMA.Mesh.Topologies
using Combinatorics, MPI

@testset "cubedshellwarp tests" begin
  import CLIMA.Mesh.Topologies: cubedshellwarp

  @testset "check radius" begin
    @test hypot(cubedshellwarp( 3.0,-2.2, 1.3)...) ≈ 3.0 rtol=eps()
    @test hypot(cubedshellwarp(-3.0,-2.2, 1.3)...) ≈ 3.0 rtol=eps()
    @test hypot(cubedshellwarp( 1.1,-2.2, 3.0)...) ≈ 3.0 rtol=eps()
    @test hypot(cubedshellwarp( 1.1,-2.2,-3.0)...) ≈ 3.0 rtol=eps()
    @test hypot(cubedshellwarp( 1.1, 3.0, 0.0)...) ≈ 3.0 rtol=eps()
    @test hypot(cubedshellwarp( 1.1,-3.0, 0.0)...) ≈ 3.0 rtol=eps()
  end

  @testset "check sign" begin
    @test sign.(cubedshellwarp( 3.0,-2.2, 1.3)) == sign.(( 3.0,-2.2, 1.3))
    @test sign.(cubedshellwarp(-3.0,-2.2, 1.3)) == sign.((-3.0,-2.2, 1.3))
    @test sign.(cubedshellwarp( 1.1,-2.2, 3.0)) == sign.(( 1.1,-2.2, 3.0))
    @test sign.(cubedshellwarp( 1.1,-2.2,-3.0)) == sign.(( 1.1,-2.2,-3.0))
    @test sign.(cubedshellwarp( 1.1, 3.0, 0.0)) == sign.(( 1.1, 3.0, 0.0))
    @test sign.(cubedshellwarp( 1.1,-3.0, 0.0)) == sign.(( 1.1,-3.0, 0.0))
  end

  @testset "check continuity" begin
    for (u,v) in zip(permutations([3.0, 2.999999999, 1.3]), permutations([2.999999999, 3.0, 1.3]))
      @test all(cubedshellwarp(u...) .≈ cubedshellwarp(v...))
    end
    for (u,v) in zip(permutations([3.0, -2.999999999, 1.3]), permutations([2.999999999, -3.0, 1.3]))
      @test all(cubedshellwarp(u...) .≈ cubedshellwarp(v...))
    end
    for (u,v) in zip(permutations([-3.0, 2.999999999, 1.3]), permutations([-2.999999999, 3.0, 1.3]))
      @test all(cubedshellwarp(u...) .≈ cubedshellwarp(v...))
    end
    for (u,v) in zip(permutations([-3.0, -2.999999999, 1.3]), permutations([-2.999999999, -3.0, 1.3]))
      @test all(cubedshellwarp(u...) .≈ cubedshellwarp(v...))
    end
  end
end

@testset "cubedshellunwarp" begin
  import CLIMA.Mesh.Topologies: cubedshellwarp, cubedshellunwarp

  for u in permutations([3.0, 2.999999999, 1.3])
    @test all(cubedshellunwarp(cubedshellwarp(u...)...) .≈ u)
  end
  for u in permutations([3.0, -2.999999999, 1.3])
    @test all(cubedshellunwarp(cubedshellwarp(u...)...) .≈ u)
  end
  for u in permutations([-3.0, 2.999999999, 1.3])
    @test all(cubedshellunwarp(cubedshellwarp(u...)...) .≈ u)
  end
  for u in permutations([-3.0, -2.999999999, 1.3])
    @test all(cubedshellunwarp(cubedshellwarp(u...)...) .≈ u)
  end
end

@testset "grid1d" begin
  g = grid1d(0,10,nelem=10)
  @test eltype(g) == Float64
  @test length(g) == 11
  @test g[1] == 0
  @test g[end] == 10

  g = grid1d(10f0,20f0,elemsize=0.1)
  @test eltype(g) == Float32
  @test length(g) == 101
  @test g[1] == 10
  @test g[end] == 20

  g = grid1d(10f0,20f0,InteriorStretching(0),elemsize=0.1)
  @test eltype(g) == Float32
  @test length(g) == 101
  @test g[1] == 10
  @test g[end] == 20

  g = grid1d(10f0,20f0,SingleExponentialStretching(2.5f0),elemsize=0.1)
  @test eltype(g) == Float32
  @test length(g) == 101
  @test g[1] == 10
  @test g[end] == 20
end

@testset "BrickTopology tests" begin

  let
    comm = MPI.COMM_SELF

    elemrange = (0:10,)
    periodicity=(true,)

    topology = BrickTopology(comm, elemrange, periodicity=periodicity)

    nelem = length(elemrange[1])-1

    for e = 1:nelem
      @test topology.elemtocoord[:,:, e] == [e-1  e]
    end

    @test topology.elemtoelem == [nelem collect(1:(nelem-1))'; collect(2:nelem)' 1]
    @test topology.elemtoface == repeat(2:-1:1, outer=(1,nelem))

    @test topology.elemtoordr == ones(Int, size(topology.elemtoordr))
    @test topology.elemtobndy == zeros(Int, size(topology.elemtoordr))

    @test topology.elems == 1:nelem
    @test topology.realelems == 1:nelem
    @test topology.ghostelems == nelem.+(1:0)

    @test length(topology.sendelems) == 0

    @test topology.nabrtorank == Int[]
    @test topology.nabrtorecv == UnitRange{Int}[]
    @test topology.nabrtosend == UnitRange{Int}[]
  end

  let
    comm = MPI.COMM_SELF
    topology = BrickTopology(comm, (0:4,5:9), periodicity=(false,true))

    nelem = 16

    @test topology.elemtocoord[:,:, 1] == [0 1 0 1; 5 5 6 6]
    @test topology.elemtocoord[:,:, 2] == [1 2 1 2; 5 5 6 6]
    @test topology.elemtocoord[:,:, 3] == [1 2 1 2; 6 6 7 7]
    @test topology.elemtocoord[:,:, 4] == [0 1 0 1; 6 6 7 7]
    @test topology.elemtocoord[:,:, 5] == [0 1 0 1; 7 7 8 8]
    @test topology.elemtocoord[:,:, 6] == [0 1 0 1; 8 8 9 9]
    @test topology.elemtocoord[:,:, 7] == [1 2 1 2; 8 8 9 9]
    @test topology.elemtocoord[:,:, 8] == [1 2 1 2; 7 7 8 8]
    @test topology.elemtocoord[:,:, 9] == [2 3 2 3; 7 7 8 8]
    @test topology.elemtocoord[:,:,10] == [2 3 2 3; 8 8 9 9]
    @test topology.elemtocoord[:,:,11] == [3 4 3 4; 8 8 9 9]
    @test topology.elemtocoord[:,:,12] == [3 4 3 4; 7 7 8 8]
    @test topology.elemtocoord[:,:,13] == [3 4 3 4; 6 6 7 7]
    @test topology.elemtocoord[:,:,14] == [2 3 2 3; 6 6 7 7]
    @test topology.elemtocoord[:,:,15] == [2 3 2 3; 5 5 6 6]
    @test topology.elemtocoord[:,:,16] == [3 4 3 4; 5 5 6 6]

    @test topology.elemtoelem ==
      [1   1   4  4  5  6   6  5   8   7  10   9  14   3   2  15
       2  15  14  3  8  7  10  9  12  11  11  12  13  13  16  16
       6   7   2  1  4  5   8  3  14   9  12  13  16  15  10  11
       4   3   8  5  6  1   2  7  10  15  16  11  12   9  14  13]

    @test topology.elemtoface ==
      [1  2  2  1  1  1  2  2  2  2  2  2  2  2  2  2
       1  1  1  1  1  1  1  1  1  1  2  2  2  1  1  2
       4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4
       3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3]

    @test topology.elemtoordr == ones(Int, size(topology.elemtoordr))

    @test topology.elems == 1:nelem
    @test topology.realelems == 1:nelem
    @test topology.ghostelems == nelem.+(1:0)

    @test length(topology.sendelems) == 0

    @test topology.nabrtorank == Int[]
    @test topology.nabrtorecv == UnitRange{Int}[]
    @test topology.nabrtosend == UnitRange{Int}[]
  end

  let
    comm = MPI.COMM_SELF
    for px in (true, false)
      topology = BrickTopology(comm, (0:10,), periodicity=(px,))
      @test Topologies.hasboundary(topology) == !px
    end
    for py = (true, false), px in (true, false)
      topology = BrickTopology(comm, (0:10,0:3), periodicity=(px,py))
      @test Topologies.hasboundary(topology) == !(px && py)
    end
    for pz = (true, false), py = (true, false), px in (true, false)
      topology = BrickTopology(comm, (0:10,0:3,-1:3), periodicity=(px,py,pz))
      @test Topologies.hasboundary(topology) == !(px && py && pz)
    end
  end
end

@testset "StackedBrickTopology tests" begin
  let
    comm = MPI.COMM_SELF
    topology = StackedBrickTopology(comm, (2:5,4:6), periodicity=(false,true),
                                    boundary=((1,2),(3,4)))

    nelem = 6

    @test topology.elemtocoord[:,:, 1] == [2 3 2 3; 4 4 5 5]
    @test topology.elemtocoord[:,:, 2] == [2 3 2 3; 5 5 6 6]
    @test topology.elemtocoord[:,:, 3] == [3 4 3 4; 4 4 5 5]
    @test topology.elemtocoord[:,:, 4] == [3 4 3 4; 5 5 6 6]
    @test topology.elemtocoord[:,:, 5] == [4 5 4 5; 4 4 5 5]
    @test topology.elemtocoord[:,:, 6] == [4 5 4 5; 5 5 6 6]

    @test topology.elemtoelem ==
      [1 2 1 2 3 4
       3 4 5 6 5 6
       2 1 4 3 6 5
       2 1 4 3 6 5]

    @test topology.elemtoface ==
      [1 1 2 2 2 2
       1 1 1 1 2 2
       4 4 4 4 4 4
       3 3 3 3 3 3]

    @test topology.elemtoordr == ones(Int, size(topology.elemtoordr))

    @test topology.elemtobndy ==
      [1 1 0 0 0 0
       0 0 0 0 2 2
       0 0 0 0 0 0
       0 0 0 0 0 0]

    @test topology.elems == 1:nelem
    @test topology.realelems == 1:nelem
    @test topology.ghostelems == nelem.+(1:0)

    @test length(topology.sendelems) == 0

    @test topology.nabrtorank == Int[]
    @test topology.nabrtorecv == UnitRange{Int}[]
    @test topology.nabrtosend == UnitRange{Int}[]
  end
  let
    comm = MPI.COMM_SELF
    for py = (true, false), px in (true, false)
      topology = StackedBrickTopology(comm, (0:10,0:3), periodicity=(px,py))
      @test Topologies.hasboundary(topology) == !(px && py)
    end
    for pz = (true, false), py = (true, false), px in (true, false)
      topology = StackedBrickTopology(comm, (0:10,0:3,-1:3),
                                      periodicity=(px,py,pz))
      @test Topologies.hasboundary(topology) == !(px && py && pz)
    end
  end
end

@testset "StackedCubedSphereTopology tests" begin
  topology = StackedCubedSphereTopology(MPI.COMM_SELF, 3, 1.0:3.0)
  @test Topologies.hasboundary(topology)
end

@testset "CubedShellTopology tests" begin
  topology = CubedShellTopology(MPI.COMM_SELF, 3, Float64)
  @test !Topologies.hasboundary(topology)
end

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

using Base.Test

using Base.LinAlg: BlasComplex, BlasFloat, BlasReal

n = 10

# Split n into 2 parts for tests needing two matrices
n1 = div(n, 2)
n2 = 2*n1

srand(1234321)

areal = randn(n,n)/2
aimg  = randn(n,n)/2
a2real = randn(n,n)/2
a2img  = randn(n,n)/2
breal = randn(n,2)/2
bimg  = randn(n,2)/2

@testset for eltya in (Float32, Float64, Complex64, Complex128)
    a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
    a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real)
    asym = a'+a                  # symmetric indefinite
    apd  = a'*a                 # symmetric positive-definite
    ε = εa = eps(abs(float(one(eltya))))

    @testset for eltyb in (Float32, Float64, Complex64, Complex128, Int)
        b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex.(breal, bimg) : breal)
        εb = eps(abs(float(one(eltyb))))
        ε = max(εa,εb)

        α = rand(eltya)
        aα = fill(α,1,1)
        @test lqfact(α)[:L]*lqfact(α)[:Q] ≈ lqfact(aα)[:L]*lqfact(aα)[:Q]
        @test lq(α)[1]*lq(α)[2] ≈ lqfact(aα)[:L]*lqfact(aα)[:Q]
        @test abs(lqfact(α)[:Q][1,1]) ≈ one(eltya)
        tab = promote_type(eltya,eltyb)

        for i = 1:2
            let a = i == 1 ? a : view(a, 1:n - 1, 1:n - 1), b = i == 1 ? b : view(b, 1:n - 1), n = i == 1 ? n : n - 1
                lqa   = lqfact(a)
                l,q   = lqa[:L], lqa[:Q]
                qra   = qrfact(a)
                @testset "Basic ops" begin
                    @test size(lqa,1) == size(a,1)
                    @test size(lqa,3) == 1
                    @test size(lqa[:Q],3) == 1
                    @test Base.LinAlg.getq(lqa) == lqa[:Q]
                    @test_throws KeyError lqa[:Z]
                    @test full(lqa') ≈ a'
                    @test lqa * lqa' ≈ a * a'
                    @test lqa' * lqa ≈ a' * a
                    @test q*full(q, thin = false)' ≈ eye(eltya,n)
                    @test l*q ≈ a
                    @test full(lqa) ≈ a
                    @test full(copy(lqa)) ≈ a
                    lstring = sprint(show,l)
                    qstring = sprint(show,q)
                    @test sprint(show,lqa) == "$(typeof(lqa)) with factors L and Q:\n$lstring\n$qstring"
                end
                @testset "Binary ops" begin
                    @test a*(lqa\b) ≈ b atol=3000ε
                    @test lqa*b ≈ qra[:Q]*qra[:R]*b atol=3000ε
                    @test A_mul_Bc(eye(eltyb,size(q.factors,2)),q)*full(q,thin=false) ≈ eye(n) atol=5000ε
                    if eltya != Int
                        @test eye(eltyb,n)*q ≈ convert(AbstractMatrix{tab},q)
                    end
                    @test q*b ≈ full(q,thin=false)*b atol=100ε
                    @test q.'*b ≈ full(q,thin=false).'*b atol=100ε
                    @test q'*b ≈ full(q,thin=false)'*b atol=100ε
                    @test a*q ≈ a*full(q,thin=false) atol=100ε
                    @test a*q.' ≈ a*full(q,thin=false).' atol=100ε
                    @test a*q' ≈ a*full(q,thin=false)' atol=100ε
                    @test a'*q ≈ a'*full(q,thin=false) atol=100ε
                    @test a'*q' ≈ a'*full(q,thin=false)' atol=100ε
                    @test_throws DimensionMismatch q*b[1:n1 + 1]
                    @test_throws DimensionMismatch Ac_mul_B(q,ones(eltya,n+2,n+2))
                    @test_throws DimensionMismatch ones(eltyb,n+2,n+2)*q
                    if isa(a, DenseArray) && isa(b, DenseArray)
                        #use this to test 2nd branch in mult code
                        pad_a = vcat(eye(a), a)
                        pad_b = hcat(eye(b), b)
                        @test pad_a*q ≈ pad_a*full(q,thin=false) atol=100ε
                        @test q.'*pad_b ≈ full(q,thin=false).'*pad_b atol=100ε
                        @test q'*pad_b ≈ full(q,thin=false)'*pad_b atol=100ε
                    end
                end
            end
        end

        @testset "Matmul with LQ factorizations" begin
            lqa = lqfact(a[:,1:n1])
            l,q = lqa[:L], lqa[:Q]
            @test full(q)*full(q)' ≈ eye(eltya,n1)
            @test (full(q,thin=false)'*full(q,thin=false))[1:n1,:] ≈ eye(eltya,n1,n)
            @test_throws DimensionMismatch A_mul_B!(eye(eltya,n+1),q)
            @test Ac_mul_B!(q,full(q)) ≈ eye(eltya,n1)
            @test_throws DimensionMismatch A_mul_Bc!(eye(eltya,n+1),q)
            @test_throws BoundsError size(q,-1)
        end
    end
end


text:
module AtomShell

using Sockets

abstract type Shell end

_shell = nothing

include("install.jl")
include("process.jl")
include("window.jl")

end


text:
function get_jacobian(obj::DenseWeightedEvaluation)
    return jcall(obj, "getJacobian", RealMatrix, ())
end

function get_point(obj::DenseWeightedEvaluation)
    return jcall(obj, "getPoint", RealVector, ())
end

function get_residuals(obj::DenseWeightedEvaluation)
    return jcall(obj, "getResiduals", RealVector, ())
end


text:
rtol = 1e-15

# Physical grid tests
testdx(g) = abs(sum(g.x[2:end].-g.x[1:end-1]) - (g.nx-1)*g.dx) < 1e-15*(g.nx-1)
testdy(g) = abs(sum(g.y[2:end].-g.y[1:end-1]) - (g.ny-1)*g.dy) < 1e-15*(g.nx-1)

testx(g) = isapprox(g.x[end]-g.x[1], g.Lx-g.dx, rtol=rtol)
testy(g) = isapprox(g.y[end]-g.y[1], g.Ly-g.dy, rtol=rtol)

testX(g) = sum(g.X[end, :].-g.X[1, :]) - (g.Lx-g.dx)*g.ny < 1e-15*g.ny
testY(g) = sum(g.Y[:, end].-g.Y[:, 1]) - (g.Ly-g.dy)*g.nx < 1e-15*g.nx

# Test proper arrangement of fft wavenumbers
testk(g) = sum(g.k[2:g.nkr-1] .+ reverse(g.k[g.nkr+1:end], dims=1)) == 0.0
testkr(g) = sum(g.k[1:g.nkr] .- g.kr) == 0.0
testl(g) = sum(g.l[:, 2:Int(g.ny/2)] .+ reverse(g.l[:, Int(g.ny/2+2):end], dims=2)) == 0.0

testK(g) = sum(g.K[2:g.nkr-1, :] .+ reverse(g.K[g.nkr+1:end, :], dims=1)) == 0.0
testL(g) = sum(g.L[:, 2:Int(g.ny/2)] .+ reverse(g.L[:, Int(g.ny/2+2):end], dims=2)) == 0.0
testKr(g) = sum(g.K[1:g.nkr, :] .- g.Kr) == 0.0
testLr(g) = sum(g.L[1:g.nkr, :] .- g.Lr) == 0.0

# Test 1D grid
g1 = OneDGrid(32, 2π)

@test testdx(g1)
@test testx(g1)
@test testk(g1)
@test testkr(g1)

# Test 2D rectangular grid
g2 = TwoDGrid(32, 2π, 24, 4π)

@test testdx(g2)
@test testdy(g2)
@test testx(g2)
@test testy(g2)
@test testX(g2)
@test testY(g2)
@test testk(g2)
@test testkr(g2)
@test testl(g2)
@test testK(g2)
@test testL(g2)
@test testKr(g2)
@test testLr(g2)


text:
@testgroup "run_testcase" begin

function get_outcome(func)
    tcs = Jute.collect_testobjs() do
        @testcase "tc" begin
            func()
        end
    end
    tc = tcs[1]
    Jute.run_testcase(tc, [])
end

function get_results(func)
    return get_outcome(func).results
end

# Check that if a testcase calls no assertions, a single Pass gets added to the results.
@testcase "no assertions" begin
    results = get_results() do
    end
    @test length(results) == 1
    @test isa(results[1], Jute.Test.Pass)
end

@testcase "report value" begin
    results = get_results() do
        @test 1 == 1
        @test_result 1
        @test_result "a"
    end

    @test length(results) == 3
    @test isa(results[1], Jute.Test.Pass)
    @test isa(results[2], Jute.ReturnValue) && results[2].value == 1
    @test isa(results[3], Jute.ReturnValue) && results[3].value == "a"
end

@testcase "multiple fails" begin
    results = get_results() do
        @test 1 == 1
        @test 1 == 2 # first fail, execution should not stop

        # exception thrown, but since it happened inside an assertion,
        # execution should not stop
        @test error("caught")

        error("uncaught") # uncaught exception, execution should stop
        @test 1 == 1
    end

    @test length(results) == 4
    @test isa(results[1], Jute.Test.Pass)
    @test isa(results[2], Jute.Test.Fail)
    @test isa(results[3], Jute.Test.Error)
    @test isa(results[4], Jute.Test.Error)
end

@testcase "@test_throws" begin
    results = get_results() do
        @test 1 == 1
        @test_throws ErrorException error("Caught exception")
        @test 1 == 1
    end

    @test length(results) == 3
    @test isa(results[1], Jute.Test.Pass)
    @test isa(results[2], Jute.Test.Pass)
    @test isa(results[3], Jute.Test.Pass)
end

@testcase "@test_skip" begin
    results = get_results() do
        @test 1 == 1
        @test_skip 1 == 2
        @test 1 == 1
    end

    @test length(results) == 3
    @test isa(results[1], Jute.Test.Pass)
    @test isa(results[2], Jute.Test.Broken)
    @test isa(results[3], Jute.Test.Pass)
end

@testcase "@test_broken" begin
    results = get_results() do
        @test 1 == 1
        @test_broken 1 == 2
        @test 1 == 1
    end

    @test length(results) == 3
    @test isa(results[1], Jute.Test.Pass)
    @test isa(results[2], Jute.Test.Broken)
    @test isa(results[3], Jute.Test.Pass)
end

@testcase "unexpected pass" begin
    results = get_results() do
        @test 1 == 1
        @test_broken 1 == 1 # marked as broken, but succeeds
        @test 1 == 1
    end

    @test length(results) == 3
    @test isa(results[1], Jute.Test.Pass)
    @test isa(results[2], Jute.Test.Error)
    @test isa(results[3], Jute.Test.Pass)
end

@testcase "@inferred" begin

    f(a,b,c) = b > 1 ? 1 : 1.0

    results = get_results() do
        @inferred f(1,2,3)
    end

    @test length(results) == 1
    @test isa(results[1], Jute.Test.Error)

    results = get_results() do
        @inferred max(1, 2)
    end

    @test length(results) == 1
    @test isa(results[1], Jute.Test.Pass)
end

@testcase "@test_warn" begin
    results = get_results() do
        @test_warn "warn" println(stderr, "warn")
        @test_warn "warn" println(stderr, "foo")
    end

    @test length(results) == 2
    @test isa(results[1], Jute.Test.Pass)
    @test isa(results[2], Jute.Test.Fail)
end

@testcase "@test_nowarn" begin
    results = get_results() do
        @test_nowarn println(stderr, "warn")
        @test_nowarn sin(1.0)
    end

    @test length(results) == 2
    @test isa(results[1], Jute.Test.Fail)
    @test isa(results[2], Jute.Test.Pass)
end

function critical_fail()
    i = 1 # to check that an external variable can be used in the assertion
    @critical @test i == 2 # critical failure, execution should stop
end

@testcase "@critical" begin
    results = get_results() do
        @test 1 == 1
        critical_fail() # test that @critical works even inside a function
        @test 1 == 1
    end

    @test length(results) == 2
    @test isa(results[1], Jute.Test.Pass)
    @test isa(results[2], Jute.Test.Fail)
end

@testcase "@critical on an unsupported expression" begin
    @test_throws LoadError eval(:( @critical 1 == 1 ))
end

@testcase "is_failed()" begin
    outcome = get_outcome() do
        @test 1 == 1
        @test 2 == 2
    end
    @test !Jute.is_failed(outcome)

    outcome = get_outcome() do
        @test 1 == 1
        @test 1 == 2
    end
    @test Jute.is_failed(outcome)

    outcome = get_outcome() do
        @test 1 == 1
        error("uncaught")
    end
    @test Jute.is_failed(outcome)

    outcome = get_outcome() do
        @test 1 == 1
        @test_broken 1 == 2
    end
    @test !Jute.is_failed(outcome)

    outcome = get_outcome() do
        @test 1 == 1
        @test_broken 1 == 1
    end
    @test Jute.is_failed(outcome)

    outcome = get_outcome() do
        @test 1 == 1
        @test_result 1
    end
    @test !Jute.is_failed(outcome)
end

end


text:
__precompile__(true)
module SatelliteDynamics

# Usings
using Dates
using Printf
using LinearAlgebra
using SOFA

using StaticArrays: SVector, SMatrix

# Define constants for array indexing
const idx1t3 = SVector(1, 2, 3)
const idx4t6 = SVector(4, 5, 6)
const idx1t6 = SVector(1, 2, 3, 4, 5, 6)

# Includes
include("constants.jl")
include("universe.jl")
include("time.jl")
include("reference_systems.jl")
include("attitude.jl")
include("coordinates.jl")
include("astrodynamics.jl")
include("orbit_dynamics.jl")
include("sgp_models.jl")

# Export EarthEnvironment submodule
include(joinpath(".", "earth_environment", "earth_environment.jl"))
include(joinpath(".", "simulation", "simulation.jl"))

end # module


text:
using ..Ocean: coriolis_parameter
using ..HydrostaticBoussinesq
using ...DGMethods

import ...BalanceLaws
import ..HydrostaticBoussinesq:
    viscosity_tensor,
    coriolis_force!,
    velocity_gradient_argument!,
    velocity_gradient_flux!,
    hydrostatic_pressure!,
    compute_flow_deviation!

@inline function velocity_gradient_argument!(m::HBModel, ::Coupled, G, Q, A, t)
    G.∇u = Q.u
    G.∇uᵈ = A.uᵈ

    return nothing
end

@inline function velocity_gradient_flux!(m::HBModel, ::Coupled, D, G, Q, A, t)
    ν = viscosity_tensor(m)
    ∇u = @SMatrix [
        G.∇uᵈ[1, 1] G.∇uᵈ[1, 2]
        G.∇uᵈ[2, 1] G.∇uᵈ[2, 2]
        G.∇u[3, 1] G.∇u[3, 2]
    ]
    D.ν∇u = -ν * ∇u

    return nothing
end

@inline hydrostatic_pressure!(::HBModel, ::Coupled, _...) = nothing

@inline function coriolis_force!(m::HBModel, ::Coupled, S, Q, A, t)
    # f × u
    f = coriolis_parameter(m, m.problem, A.y)
    uᵈ, vᵈ = A.uᵈ # Horizontal components of velocity
    S.u -= @SVector [-f * vᵈ, f * uᵈ]

    return nothing
end

# Compute Horizontal Flow deviation from vertical mean
@inline function compute_flow_deviation!(dg, m::HBModel, ::Coupled, Q, t)
    FT = eltype(Q)
    Nq, Nqk, _, _, nelemv, nelemh, nrealelemh, _ = basic_grid_info(dg)

    #### integrate the tendency
    model_int = dg.modeldata.integral_model
    integral = model_int.balance_law
    update_auxiliary_state!(model_int, integral, Q, 0)

    ### properly shape MPIStateArrays
    num_int = number_states(integral, Auxiliary())
    data_int = model_int.state_auxiliary.data
    data_int = reshape(data_int, Nq^2, Nqk, num_int, nelemv, nelemh)

    num_aux = number_states(m, Auxiliary())
    data_aux = dg.state_auxiliary.data
    data_aux = reshape(data_aux, Nq^2, Nqk, num_aux, nelemv, nelemh)

    num_state = number_states(m, Prognostic())
    data_state = reshape(Q.data, Nq^2, Nqk, num_state, nelemv, nelemh)

    ### get vars indices
    index_∫u = varsindex(vars_state(integral, Auxiliary(), FT), :(∫x))
    index_uᵈ = varsindex(vars_state(m, Auxiliary(), FT), :uᵈ)
    index_u = varsindex(vars_state(m, Prognostic(), FT), :u)

    ### get top value (=integral over full depth)
    ∫u = @view data_int[:, end:end, index_∫u, end:end, 1:nrealelemh]
    uᵈ = @view data_aux[:, :, index_uᵈ, :, 1:nrealelemh]
    u = @view data_state[:, :, index_u, :, 1:nrealelemh]

    ## make a copy of horizontal velocity
    ## and remove vertical mean velocity
    uᵈ .= u
    uᵈ .-= ∫u / m.problem.H

    return nothing
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
module SystemStructures

using DataStructures
using SymbolicUtils: istree, operation, arguments, Symbolic
using ..ModelingToolkit
import ..ModelingToolkit: isdiffeq, var_from_nested_derivative, vars!, flatten,
    value, InvalidSystemException, isdifferential, _iszero, isparameter
using ..BipartiteGraphs
using UnPack
using Setfield
using SparseArrays

#=
When we don't do subsitution, variable information is split into two different
places, i.e. `states` and the right-hand-side of `observed`.

eqs = [0 ~ z + x; 0 ~ y + z^2]
states = [y, z]
observed = [x ~ sin(y) + z]
struct Reduced
    var
    expr
    idxs
end
fullvars = [Reduced(x, sin(y) + z, [2, 3]), y, z]
active_𝑣vertices = [false, true, true]
      x   y   z
eq1:  1       1
eq2:      1   1

      x   y   z
eq1:      1   1
eq2:      1   1

for v in 𝑣vertices(graph); active_𝑣vertices[v] || continue

end
=#

export SystemStructure, SystemPartition
export initialize_system_structure, find_linear_equations
export isdiffvar, isdervar, isalgvar, isdiffeq, isalgeq
export dervars_range, diffvars_range, algvars_range

@enum VariableType::Int8 DIFFERENTIAL_VARIABLE ALGEBRAIC_VARIABLE DERIVATIVE_VARIABLE

Base.@kwdef struct SystemPartition
    e_solved::Vector{Int}
    v_solved::Vector{Int}
    e_residual::Vector{Int}
    v_residual::Vector{Int}
end

function Base.:(==)(s1::SystemPartition, s2::SystemPartition)
    tup1 = (s1.e_solved, s1.v_solved, s1.e_residual, s1.v_residual)
    tup2 = (s2.e_solved, s2.v_solved, s2.e_residual, s2.v_residual)
    tup1 == tup2
end

Base.@kwdef struct SystemStructure
    fullvars::Vector
    vartype::Vector{VariableType}
    varassoc::Vector{Int}
    inv_varassoc::Vector{Int}
    varmask::BitVector # `true` if the variable has the highest order derivative
    algeqs::BitVector
    graph::BipartiteGraph{Int,Vector{Vector{Int}},Int,Nothing}
    solvable_graph::BipartiteGraph{Int,Vector{Vector{Int}},Int,Nothing}
    assign::Vector{Int}
    inv_assign::Vector{Int}
    scc::Vector{Vector{Int}}
    partitions::Vector{SystemPartition}
end

isdervar(s::SystemStructure, var::Integer) = s.vartype[var] === DERIVATIVE_VARIABLE
isdiffvar(s::SystemStructure, var::Integer) = s.vartype[var] === DIFFERENTIAL_VARIABLE
isalgvar(s::SystemStructure, var::Integer) = s.vartype[var] === ALGEBRAIC_VARIABLE

dervars_range(s::SystemStructure) = Iterators.filter(Base.Fix1(isdervar, s), eachindex(s.vartype))
diffvars_range(s::SystemStructure) = Iterators.filter(Base.Fix1(isdiffvar, s), eachindex(s.vartype))
algvars_range(s::SystemStructure) = Iterators.filter(Base.Fix1(isalgvar, s), eachindex(s.vartype))

isalgeq(s::SystemStructure, eq::Integer) = s.algeqs[eq]
isdiffeq(s::SystemStructure, eq::Integer) = !isalgeq(s, eq)

function initialize_system_structure(sys)
    sys = flatten(sys)

    iv = independent_variable(sys)
    eqs = copy(equations(sys))
    neqs = length(eqs)
    algeqs = trues(neqs)
    dervaridxs = OrderedSet{Int}()
    var2idx = Dict{Any,Int}()
    symbolic_incidence = []
    fullvars = []
    var_counter = Ref(0)
    addvar! = let fullvars=fullvars, var_counter=var_counter
        var -> begin
            get!(var2idx, var) do
                push!(fullvars, var)
                var_counter[] += 1
            end
        end
    end

    vars = OrderedSet()
    for (i, eq′) in enumerate(eqs)
        if _iszero(eq′.lhs)
            eq = eq′
        else
            eq = 0 ~ eq′.rhs - eq′.lhs
        end
        vars!(vars, eq.rhs)
        isalgeq = true
        statevars = []
        for var in vars
            isequal(var, iv) && continue
            if isparameter(var) || (istree(var) && isparameter(operation(var)))
                continue
            end
            varidx = addvar!(var)
            push!(statevars, var)

            dvar = var
            idx = varidx
            while isdifferential(dvar)
                if !(idx in dervaridxs)
                    push!(dervaridxs, idx)
                end
                isalgeq = false
                dvar = arguments(dvar)[1]
                idx = addvar!(dvar)
            end
        end
        push!(symbolic_incidence, copy(statevars))
        empty!(statevars)
        empty!(vars)
        algeqs[i] = isalgeq
        if isalgeq
            eqs[i] = eq
        end
    end

    # sort `fullvars` such that the mass matrix is as diagonal as possible.
    dervaridxs = collect(dervaridxs)
    sorted_fullvars = OrderedSet(fullvars[dervaridxs])
    for dervaridx in dervaridxs
        dervar = fullvars[dervaridx]
        diffvar = arguments(dervar)[1]
        if !(diffvar in sorted_fullvars)
            push!(sorted_fullvars, diffvar)
        end
    end
    for v in fullvars
        if !(v in sorted_fullvars)
            push!(sorted_fullvars, v)
        end
    end
    fullvars = collect(sorted_fullvars)
    var2idx = Dict(fullvars .=> eachindex(fullvars))
    dervaridxs = 1:length(dervaridxs)

    nvars = length(fullvars)
    diffvars = []
    vartype = fill(DIFFERENTIAL_VARIABLE, nvars)
    varassoc = zeros(Int, nvars)
    inv_varassoc = zeros(Int, nvars)
    for dervaridx in dervaridxs
        vartype[dervaridx] = DERIVATIVE_VARIABLE
        dervar = fullvars[dervaridx]
        diffvar = arguments(dervar)[1]
        diffvaridx = var2idx[diffvar]
        push!(diffvars, diffvar)
        varassoc[diffvaridx] = dervaridx
        inv_varassoc[dervaridx] = diffvaridx
    end

    algvars = setdiff(states(sys), diffvars)
    for algvar in algvars
        # it could be that a variable appeared in the states, but never appeared
        # in the equations.
        algvaridx = var2idx[algvar]
        vartype[algvaridx] = ALGEBRAIC_VARIABLE
    end

    graph = BipartiteGraph(neqs, nvars, Val(false))
    for (ie, vars) in enumerate(symbolic_incidence), v in vars
        jv = var2idx[v]
        add_edge!(graph, ie, jv)
    end

    @set! sys.eqs = eqs
    @set! sys.structure = SystemStructure(
        fullvars = fullvars,
        vartype = vartype,
        varassoc = varassoc,
        inv_varassoc = inv_varassoc,
        varmask = iszero.(varassoc),
        algeqs = algeqs,
        graph = graph,
        solvable_graph = BipartiteGraph(nsrcs(graph), ndsts(graph), Val(false)),
        assign = Int[],
        inv_assign = Int[],
        scc = Vector{Int}[],
        partitions = SystemPartition[],
    )
    return sys
end

function find_linear_equations(sys)
    s = structure(sys)
    @unpack fullvars, graph = s
    is_linear_equations = falses(nsrcs(graph))
    eqs = equations(sys)
    eadj = Vector{Int}[]
    cadj = Vector{Int}[]
    coeffs = Int[]
    for (i, eq) in enumerate(eqs); isdiffeq(eq) && continue
        empty!(coeffs)
        linear_term = 0
        all_int_vars = true

        term = value(eq.rhs - eq.lhs)
        for j in 𝑠neighbors(graph, i)
            var = fullvars[j]
            c = expand_derivatives(Differential(var)(term), false)
            # test if `var` is linear in `eq`.
            if !(c isa Symbolic) && c isa Number
                if c == 1 || c == -1
                    c = convert(Integer, c)
                    linear_term += c * var
                    push!(coeffs, c)
                else
                    all_int_vars = false
                end
            end
        end

        # Check if all states in the equation is both linear and homogeneous,
        # i.e. it is in the form of
        #
        #       ``∑ c_i * v_i = 0``,
        #
        # where ``c_i`` ∈ ℤ and ``v_i`` denotes states.
        if all_int_vars && isequal(linear_term, term)
            is_linear_equations[i] = true
            push!(eadj, copy(𝑠neighbors(graph, i)))
            push!(cadj, copy(coeffs))
        else
            is_linear_equations[i] = false
        end
    end

    return is_linear_equations, eadj, cadj
end

function Base.show(io::IO, s::SystemStructure)
    @unpack graph = s
    S = incidence_matrix(graph, Num(Sym{Real}(:×)))
    print(io, "Incidence matrix:")
    show(io, S)
end

end # module


text:
#=
    Following the description in
        _Exploring corner transfer matrices and corner tensors for the classical
        simulation of quantum lattice systems_
    by Orus

    this implementation follows the description in III.B, specifically
    the simplified one-directional 1D-method as in e.g. Fig 12, 13
    DOI: 10.1103/PhysRevB.85.205117
=#
struct transconjCTMIterable{T,
        TA<:AbstractTensor{T,4},
        TC<:AbstractTensor{T,2},
        TT<:AbstractTensor{T,3}} <: AbstractCTMIterable
    A::TA
    χ::Int
    Cinit::Union{TC,Nothing}
    Tsinit::Union{Tuple{TT,TT},Nothing}
end

function transconjctmiterable(A::DTensor{T,4}, χ::Int,
        Cinit = nothing,
        Tsinit = nothing) where T
    transconjCTMIterable{T,DTensor{T,4},DTensor{T,2},DTensor{T,3}}(A, χ, Cinit, Tsinit)
end

function transconjctmiterable(A::DASTensor{T,4,SYM,CHS,SS,CH}, χ::Int,
        Cinit = nothing,
        Tsinit = nothing) where {T,N,SYM,CHS,SS,CH}
    TC = DASTensor{T,2,SYM,CHS,SS,CH}
    TT = DASTensor{T,3,SYM,CHS,SS,CH}
    TA = DASTensor{T,4,SYM,CHS,SS,CH}
    transconjCTMIterable{T,TA,TC,TT}(A, χ, Cinit, Tsinit)
end

struct transconjCTMState{S,
        TA <: AbstractTensor,
        TC <: AbstractTensor,
        TT <: AbstractTensor} <: AbstractCTMState
    C::TC
    Ts::Tuple{TT,TT}
    oldsvdvals::Vector{S}
    diffs::Vector{S}
    n_it::Ref{Int}
end

function iterate(iter::transconjCTMIterable{S,TA,TC,TT}) where {S,TA,TC,TT}
    @unpack A, χ, Cinit, Tsinit = iter
    C  = isnothing(Cinit)  ? initializeC(A,χ) : Cinit
    Ts = isnothing(Tsinit) ? (initializeT(A,χ),initializeT(A,χ)) : Tsinit

    l = ifelse(C isa DASTensor, 2χ, χ)
    oldsvdvals = zeros(real(S),l)
    state = transconjCTMState{real(S),TA,TC,TT}(C, Ts, oldsvdvals, [], Ref(0))
    return state, state
end

function iterate(iter::transconjCTMIterable, state::transconjCTMState)
    @unpack A, χ = iter
    @unpack C, Ts, oldsvdvals, diffs, n_it = state
    T1, T2 = Ts

    xmove!((C,T2),(T1,A), χ)
    ymove!((C,T1),(T2,A), χ)

    vals = diag(tensorsvd(C)[2])
    maxval = maximum(vals)
    apply!(C, x -> x .= x ./ maxval)
    apply!(T1, x -> x .= x ./ norm(T1))
    apply!(T2, x -> x .= x ./ norm(T2))
    normalize!(vals,1)

    #compare
    push!(diffs, sum(abs, oldsvdvals - vals))
    oldsvdvals[:] = vals
    n_it[] += 1

    return state, state
end

function xmove!((C,T2),(T1,A), χ)
    # grow
    @tensor begin
        Cp[1,2,3]      := C[1,-1]    * T1[-1,2,3]
        T2p[1,2,3,4,5] := T2[1,-1,5] * A[2,3,4,-1]
    end

    #renormalize
    U, = tensorsvd(Cp, ((1,2),(3,)), svdtrunc = svdtrunc_maxχ(χ))
    @tensor begin
        C[1,2]    = U'[-1,-2,1] * Cp[-1,-2,2]
        T2[1,2,3] = U'[-1,-2,1] * T2p[-1,-2,2,-3,-4] * U[-4,-3,3]
    end
end

function ymove!((C,T1),(T2,A), χ)
    # grow
    @tensor begin
        Cp[1,2,3]      := C[-1,3]    * T2[1,2,-1]
        T1p[1,2,3,4,5] := T1[1,-1,5] * A[3,4,-1,2]
    end

    #renormalize
    u,s,Vd = tensorsvd(Cp, ((1,),(2,3)), svdtrunc = svdtrunc_maxχ(χ))
    @tensor begin
        C[1,2]    = Cp[1,-1,-2] * Vd'[2,-1,-2]
        T1[1,2,3] = Vd[1,-1,-2] * T1p[-2,-1,2,-3,-4] * Vd'[3,-3,-4]
    end
end


text:
"""Biharmonic diffusion operator with constant viscosity coefficient
    Viscosity = ν∇⁴ ⃗u. Although constant, the coefficient is actually inside
    Viscosity = ∇⋅ν∇∇² ⃗u."""
function diffusion_constant!(dudx,dudy,dvdx,dvdy,DS,DS_q,DT,νSmag,νSmag_q,Lu,Lv,dLudx,dLudy,dLvdx,dLvdy,
        S11,S12,S21,S22,LLu1,LLu2,LLv1,LLv2,u,v)
    stress_tensor!(dLudx,dLudy,dLvdx,dLvdy,Lu,Lv,u,v)
    viscous_tensor_Constant!(S11,S12,S21,S22,dLudx,dLudy,dLvdx,dLvdy)

    ∂x!(LLu1,S11)
    ∂y!(LLu2,S12)
    ∂x!(LLv1,S21)
    ∂y!(LLv2,S22)
end

""" Smagorinsky-like biharmonic viscosity
    Viscosity = ∇ ⋅ (cSmag Δ⁴ |D| ∇∇² ⃗u)
The Δ⁴-scaling is omitted as gradient operators are dimensionless."""
function diffusion_Smagorinsky!(dudx,dudy,dvdx,dvdy,DS,DS_q,DT,νSmag,νSmag_q,Lu,Lv,dLudx,dLudy,dLvdx,dLvdy,
        S11,S12,S21,S22,LLu1,LLu2,LLv1,LLv2,u,v)
    ∂x!(dudx,u)
    ∂y!(dvdy,v)
    ∂x!(dvdx,v)
    ∂y!(dudy,u)

    # biharmonic diffusion
    stress_tensor!(dLudx,dLudy,dLvdx,dLvdy,Lu,Lv,u,v)
    Smagorinsky_coeff!(νSmag,νSmag_q,DS,DS_q,DT,dudx,dvdy,dudy,dvdx)
    viscous_tensor_Smagorinsky!(S11,S12,S21,S22,νSmag,νSmag_q,dLudx,dLudy,dLvdx,dLvdy)

    ∂x!(LLu1,S11)
    ∂y!(LLu2,S12)
    ∂x!(LLv1,S21)
    ∂y!(LLv2,S22)
end

"""νSmag = cSmag * |D|, where deformation rate |D| = √((∂u/∂x - ∂v/∂y)^2 + (∂u/∂y + ∂v/∂x)^2).
The grid spacing Δ is omitted here as the operators are dimensionless."""
function Smagorinsky_coeff!(νSmag,νSmag_q,DS,DS_q,DT,dudx,dvdy,dudy,dvdx)
    # horizontal shearing strain squared
    m,n = size(DS_q)
    @boundscheck (m+ep,n) == size(dudy) || throw(BoundsError())
    @boundscheck (m,n) == size(dvdx) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            DS_q[i,j] = (dudy[i+ep,j] + dvdx[i,j])^2
        end
    end

    Ixy!(DS,DS_q)

    # horizontal tension squared
    m,n = size(DT)
    @boundscheck (m+ep,n+2) == size(dudx) || throw(BoundsError())
    @boundscheck (m+2,n) == size(dvdy) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            DT[i,j] = (dudx[i+ep,j+1] + dvdy[i+1,j])^2
        end
    end

    # viscosity = Smagorinsky coefficient times deformation rate
    m,n = size(νSmag)
    @boundscheck (m,n) == size(DS) || throw(BoundsError())
    @boundscheck (m,n) == size(DT) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            νSmag[i,j] = cSmag*sqrt(DS[i,j] + DT[i,j])
        end
    end

    Ixy!(νSmag_q,νSmag)
end

"""Biharmonic stress tensor ∇∇²(u,v) = (∂/∂x(∇²u), ∂/∂y(∇²u); ∂/∂x(∇²v), ∂/∂y(∇²v))"""
function stress_tensor!(dLudx,dLudy,dLvdx,dLvdy,Lu,Lv,u,v)
    ∇²!(Lu,u)
    ∇²!(Lv,v)
    ∂x!(dLudx,Lu)
    ∂y!(dLudy,Lu)
    ∂x!(dLvdx,Lv)
    ∂y!(dLvdy,Lv)
end

"""Biharmonic stress tensor times Smagorinsky coefficient
νSmag * ∇∇² ⃗u = (S11, S12; S21, S22)."""
function viscous_tensor_Smagorinsky!(S11,S12,S21,S22,νSmag,νSmag_q,dLudx,dLudy,dLvdx,dLvdy)
    m,n = size(S11)
    @boundscheck (m+2-ep,n) == size(νSmag) || throw(BoundsError())
    @boundscheck (m,n) == size(dLudx) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            S11[i,j] = νSmag[i+1-ep,j] * dLudx[i,j]
        end
    end

    m,n = size(S12)
    @boundscheck (m,n) == size(νSmag_q) || throw(BoundsError())
    @boundscheck (m+ep,n) == size(dLudy) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            S12[i,j] = νSmag_q[i,j] * dLudy[i+ep,j]
        end
    end

    m,n = size(S21)
    @boundscheck (m,n) == size(νSmag_q) || throw(BoundsError())
    @boundscheck (m,n) == size(dLvdx) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            S21[i,j] = νSmag_q[i,j] * dLvdx[i,j]
        end
    end

    m,n = size(S22)
    @boundscheck (m,n+2) == size(νSmag) || throw(BoundsError())
    @boundscheck (m,n) == size(dLvdy) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            S22[i,j] = νSmag[i,j+1] * dLvdy[i,j]
        end
    end
end

"""Biharmonic stress tensor times Smagorinsky coefficient
νSmag * ∇∇² ⃗u = (S11, S12; S21, S22)"""
function viscous_tensor_Constant!(S11,S12,S21,S22,dLudx,dLudy,dLvdx,dLvdy)
    m,n = size(S11)
    @boundscheck (m,n) == size(dLudx) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            S11[i,j] = νB * dLudx[i,j]
        end
    end

    m,n = size(S12)
    @boundscheck (m+ep,n) == size(dLudy) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            S12[i,j] = νB * dLudy[i+ep,j]
        end
    end

    m,n = size(S21)
    @boundscheck (m,n) == size(dLvdx) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            S21[i,j] = νB * dLvdx[i,j]
        end
    end

    m,n = size(S22)
    @boundscheck (m,n) == size(dLvdy) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            S22[i,j] = νB * dLvdy[i,j]
        end
    end
end

"""Update u with bottom friction tendency (Bu,Bv) and biharmonic viscosity."""
function add_drag_diff_tendencies!(u,v,Bu,Bv,LLu1,LLu2,LLv1,LLv2)
    m,n = size(u) .- (2*halo,2*halo)
    @boundscheck (m+2-ep,n+2) == size(Bu) || throw(BoundsError())
    @boundscheck (m,n+2) == size(LLu1) || throw(BoundsError())
    @boundscheck (m+2-ep,n) == size(LLu2) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            u[i+2,j+2] += nΔt_diff*(Bu[i+1-ep,j+1] + LLu1[i,j+1] + LLu2[i+1-ep,j])
        end
    end

    m,n = size(v) .- (2*halo,2*halo) # cut off the halo
    @boundscheck (m+2,n+2) == size(Bv) || throw(BoundsError())
    @boundscheck (m,n+2) == size(LLv1) || throw(BoundsError())
    @boundscheck (m+2,n) == size(LLv2) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
             v[i+2,j+2] += nΔt_diff*(Bv[i+1,j+1] + LLv1[i,j+1] + LLv2[i+1,j])
        end
    end
end

if diffusion == "Constant"
    diffusive! = diffusion_constant!
elseif diffusion == "Smagorinsky"
    diffusive! = diffusion_Smagorinsky!
else
    throw(error("Diffusion not correctly specified. Only Constant or Smagorinsky allowed."))
end

max_stars_repo_name: eschnett/juls

max_issues_repo_name: eschnett/juls

max_forks_repo_name: eschnett/juls


text:
using DFWannier
T= Float32
x = WannierModel{T}("/home/ponet/Documents/PhD/GeTe/NSOC/paperxsf/","/home/ponet/Documents/PhD/GeTe/SOC/GeTe_bands.out",[[Atom(T[0.0,0.0,-0.0239129,-2*0.155854]...) for i=1:4]...,[Atom(T[0.0,0.0,5.5540692,2*0.318205]...) for i=1:4]...],true);
x2 = WannierModel{T}("/home/ponet/Documents/PhD/GeTe/NSOC/paperxsf/","/home/ponet/Documents/PhD/GeTe/SOC/GeTe_bands.out",[[Atom(T[0.0,0.0,-0.0239129,-0.155854]...) for i=1:4]...,[Atom(T[0.0,0.0,5.5540692,0.318205]...) for i=1:4]...]);
using Plots
@time test = calculate_eig_angmom_soc_bloch_gpu(x,90:0.1:110.);
@time test1 = calculate_eig_angmom_soc_bloch_gpu(x,90:0.1:110.);
@time test1 = calculate_eig_angmom_soc_bloch(x2,90:0.1:110.);
plot(plot(test1[8],:angmom2_x),plot(test[8],:angmom2_x))

@time benchmark = construct_bloch_sum(x.wfcs[1],x.kpoints[1]);

begin
test1 = construct_bloch_sum(x.wfcs[1],x.kpoints[1])
assert(Array(test1.values)==Array(benchmark.values))
end

test2 = calculate_eig_angmom_soc_bloch(x2,90:110.);
using Plots
test = construct_bloch_sum(x.wfcs[1],x.kpoints[1])
Array(test.values)
test = construct_bloch_sum(x2.wfcs[1],x2.kpoints[1])
T=Float64
test_wfc1 =  DFWannier.host2gpu(read_xsf_file("/home/ponet/Documents/PhD/GeTe/NSOC/paperxsf/wan_00003.xsf",Atom{T}(0.0,0.0,0.1,0.1),T))
test_wfc2=  DFWannier.host2gpu(read_xsf_file("/home/ponet/Documents/PhD/GeTe/NSOC/paperxsf/wan_00004.xsf",Atom{T}(0.0,0.0,0.1,0.1),T))
dims = size(test_wfc1.values)
cu_dims = CuArray(UInt32[dims...])
Lx = CuArray(zeros(Complex{T},dims))
Ly = CuArray(zeros(Complex{T},dims))
Lz = CuArray(zeros(Complex{T},dims))
n1 = CuArray(zeros(Complex{T},dims))
n2 = CuArray(zeros(Complex{T},dims))
grid = Array(test_wfc1.grid)
origin = grid[1,1,1]
a = grid[2,1,1] .- origin
b = grid[1,2,1] .- origin
c = grid[1,1,2] .- origin
V = CuArray(inv([[a...] [b...] [c...]]))
begin
  # indices = Array{Tuple{Tuple{Int32,Int32,Int32},Tuple{Int32,Int32,Int32}},1}()
  # coeff = Array{Complex{T},1}()
  indices = Array{Tuple{CuArray{Int32,1},CuArray{Int32,1}},1}()
  coeff = Array{CuArray{Complex{T},1},1}()
  k=T[0.2,0.2,0.3]
  for R1=-1:0,R2=-1:0,R3=-1:0
    if R1+R2+R3 == 0
      continue
    end
    R= R1*test_wfc1.cell[1]+R2*test_wfc1.cell[2]+R3*test_wfc1.cell[3]
    ind1,ind2 = DFWannier.find_start(test_wfc1,R,27)
    # push!(indices,((Int32(ind1[1]),Int32(ind1[2]),Int32(ind1[3])),(Int32(ind1[1]-ind2[1]),Int32(ind1[2]-ind2[2]),Int32(ind1[3]-ind2[3]))))
    push!(indices,(CuArray{Int32}([ind1...]),CuArray{Int32}([ind2...])))
    # push!(coeff,Complex{T}(exp(dot(-2*pi*k,[R1,R2,R3])*1im)))
    push!(coeff,CuArray(Complex{T}(exp(dot(-2*pi*k,[R1,R2,R3])*1im))))
  end
  dfprintln(length(indices),length(coeff))
  k_wfcs = Array{Wfc3D_gpu{T},1}()
  #optimize so we dont redo the zero index ...
  for wfc in [test_wfc1,test_wfc2]
    push!(k_wfcs,Wfc3D_gpu(wfc.grid,CuArray(zeros(Complex{T},size(wfc.values))),wfc.cell,wfc.atom))
  end
end
@time  for i =1:2000
  construct_bloch_sums([test_wfc1,test_wfc2],k_wfcs,k,CuArray(indices),CuArray(coeff))
end
calculate_angmom(test_wfc1,test_wfc2,V,CuArray([test_wfc1.atom.center.x,test_wfc1.atom.center.y,test_wfc1.atom.center.z]),dims,Lx,Ly,Lz,n2,n2)

using CUDAnative, CUDAdrv

function haversine_cpu(lat1::Float32, lon1::Float32, lat2::Float32, lon2::Float32, radius::Float32)
  c1 = cospi(lat1 / 180.0f0)
  c2 = cospi(lat2 / 180.0f0)
  dlat = lat2 - lat1
  dlon = lon2 - lon1
  d1 = sinpi(dlat / 360.0f0)
  d2 = sinpi(dlon / 360.0f0)
  t = d2 * d2 * c1 * c2
  a = d1 * d1 + t
  c = 2.0f0 * asin(min(1.0f0, sqrt(a)))
  return radius * c
end

function pairwise_dist_cpu(lat::Vector{Float32}, lon::Vector{Float32})
  # allocate
  n = length(lat)
  rowresult = Array{Float32}(n, n)

  # brute force fill in each cell
  for i in 1:n, j in 1:n
    @inbounds rowresult[i, j] = haversine_cpu(lat[i], lon[i], lat[j], lon[j] , 6372.8f0)
  end

  return rowresult
end

# from https://devblogs.nvidia.com/parallelforall/fast-great-circle-distance-calculation-cuda-c/
function haversine_gpu(lat1::Float32, lon1::Float32, lat2::Float32, lon2::Float32, radius::Float32)
  # XXX: need to prefix math intrinsics with CUDAnative
  c1 = CUDAnative.cospi(lat1 / 180.0f0)
  c2 = CUDAnative.cospi(lat2 / 180.0f0)
  dlat = lat2 - lat1
  dlon = lon2 - lon1
  d1 = CUDAnative.sinpi(dlat / 360.0f0)
  d2 = CUDAnative.sinpi(dlon / 360.0f0)
  t = d2 * d2 * c1 * c2
  a = d1 * d1 + t
  c = 2.0f0 * CUDAnative.asin(CUDAnative.min(1.0f0, CUDAnative.sqrt(a)))
  return radius * c
end

# pairwise distance calculation kernel
function pairwise_dist_kernel(lat::CuDeviceVector{Float32}, lon::CuDeviceVector{Float32},
  rowresult::CuDeviceMatrix{Float32}, n)
  i = (blockIdx().x-1) * blockDim().x + threadIdx().x
  j = (blockIdx().y-1) * blockDim().y + threadIdx().y

  if i <= n && j <= n
    # store to shared memory
    shmem = @cuDynamicSharedMem(Float32, 2*blockDim().x + 2*blockDim().y)
    if threadIdx().y == 1
      shmem[threadIdx().x] = lat[i]
      shmem[blockDim().x + threadIdx().x] = lon[i]
    end
    if threadIdx().x == 1
      shmem[2*blockDim().x + threadIdx().y] = lat[j]
      shmem[2*blockDim().x + blockDim().y + threadIdx().y] = lon[j]
    end
    sync_threads()

    # load from shared memory
    lat_i = shmem[threadIdx().x]
    lon_i = shmem[blockDim().x + threadIdx().x]
    lat_j = shmem[2*blockDim().x + threadIdx().y]
    lon_j = shmem[2*blockDim().x + blockDim().y + threadIdx().y]

    @inbounds rowresult[i, j] = haversine_gpu(lat_i, lon_i, lat_j, lon_j, 6372.8f0)
  end
end

function pairwise_dist_gpu(lat::Vector{Float32}, lon::Vector{Float32})
  # upload
  lat_gpu = CuArray(lat)
  lon_gpu = CuArray(lon)

  # allocate
  n = length(lat)
  rowresult_gpu = CuArray{Float32}(n, n)

  # calculate launch configuration
  # NOTE: we want our launch configuration to be as square as possible,
  #       because that minimizes shared memory usage
  ctx = CuCurrentContext()
  dev = device(ctx)
  total_threads = min(n, attribute(dev, CUDAdrv.MAX_THREADS_PER_BLOCK))
  threads_x = floor(Int, sqrt(total_threads))
  threads_y = total_threads ÷ threads_x
  threads = (threads_x, threads_y)
  blocks = ceil.(Int, n ./ threads)

  # calculate size of dynamic shared memory
  shmem = 2 * sum(threads) * sizeof(Float32)
  dfprintln(shmem)
  @cuda (blocks, threads, shmem) pairwise_dist_kernel(lat_gpu, lon_gpu, rowresult_gpu, n)

  return Array(rowresult_gpu)
end

# generate reasonable data
const n = 10000
const lat = rand(Float32, n) .* 45
const lon = rand(Float32, n) .* -120

@test pairwise_dist_cpu(lat, lon) ≈ pairwise_dist_gpu(lat, lon)


text:
using TiledViews

function reading_writing()
    flat = collect(reshape(1:49,(7,7)))
    tiled = TiledView(flat, (4, 4),(0, 0));
    tileW, matching_window = TiledWindowView(flat, (5, 5);verbose=true);
    @show size(tileW)
    windowed = collect(tileW .* matching_window);
    tileW[:,:,:,:] .= 0  # cleares the original array
    @show tileW.parent
    tileW .+= windowed  # writes the windowed data back into the array
    @show flat
end

data = ones(10,10).+0.0;
myview, matching_window = TiledWindowView(data, (5, 5);verbose=true);
size(myview)
windowed = collect(myview .* matching_window);
myview[:,:,:,:].=0  # cleares the original array


text:
using CxxWrap
using Test

if !@isdefined libjlcxx_containers

  depfile = joinpath(@__DIR__, "deps.jl")
  if !isfile(depfile)
    include(joinpath(@__DIR__, "build.jl"))
  end

  include(depfile)
end


text:
@testset "argmax" begin
    arg, val = GeometryTypes.argmax(sin, range(0, stop=π, length=101))
    @test arg ≈ π/2
    @test val ≈ 1

    arg, val = GeometryTypes.argmax(cos, range(0, stop=π, length=101))
    @test arg ≈ 0  # 0, not π, because argmax requires a subsequent value to be *strictly better* than the current best
    @test val ≈ 1.0

    @test @allocated(GeometryTypes.argmax(sin, range(0, stop=π, length=101))) == 0

    @test_throws ArgumentError GeometryTypes.argmax(identity, [])
end


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

mutable struct IoK8sApiFlowcontrolV1alpha1ServiceAccountSubject <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    namespace::Any # spec type: Union{ Nothing, String } # spec name: namespace

    function IoK8sApiFlowcontrolV1alpha1ServiceAccountSubject(;name=nothing, namespace=nothing)
        o = new()
        validate_property(IoK8sApiFlowcontrolV1alpha1ServiceAccountSubject, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiFlowcontrolV1alpha1ServiceAccountSubject, Symbol("namespace"), namespace)
        setfield!(o, Symbol("namespace"), namespace)
        o
    end
end # type IoK8sApiFlowcontrolV1alpha1ServiceAccountSubject

const _property_map_IoK8sApiFlowcontrolV1alpha1ServiceAccountSubject = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("namespace")=>Symbol("namespace"))
const _property_types_IoK8sApiFlowcontrolV1alpha1ServiceAccountSubject = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("namespace")=>"String")
Base.propertynames(::Type{ IoK8sApiFlowcontrolV1alpha1ServiceAccountSubject }) = collect(keys(_property_map_IoK8sApiFlowcontrolV1alpha1ServiceAccountSubject))
Swagger.property_type(::Type{ IoK8sApiFlowcontrolV1alpha1ServiceAccountSubject }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiFlowcontrolV1alpha1ServiceAccountSubject[name]))}
Swagger.field_name(::Type{ IoK8sApiFlowcontrolV1alpha1ServiceAccountSubject }, property_name::Symbol) =  _property_map_IoK8sApiFlowcontrolV1alpha1ServiceAccountSubject[property_name]

function check_required(o::IoK8sApiFlowcontrolV1alpha1ServiceAccountSubject)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("namespace")) === nothing) && (return false)
    true
end

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


text:
using   Pkg

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath( @__DIR__, "..","deps","deps.jl")

if !isfile(depsjl_path)
error(
      "Package not installed properly, run Pkg.build(\"Package\"), restart Julia and try again",
     )
end
include(depsjl_path)

# Module initialization function
function __init__()
# Always check your dependencies from `deps.jl`
check_deps(   )
end

include(joinpath("..", "deps", "build.jl"))

const MaybeFloat64 = Union{Float64, Nothing }

    function foo(
        x ::Union{Nothing,Vector{UInt}} = nothing,
        y:: Union{Nothing,UInt64} = nothing;
        z :: Bool = true,
        a::AbstractString="John Doe",
        b::Integer= 12345,
        c::AbstractString ="XYZ\n\nABC",
    )

        for i in 1:b

            println( x[i] )
        end

    end