text:
module SimpleTrees

using ..Factories: Inference,Argument,Parameter,FunctionFactory,TypeFactory
using ..Factories: addfields!,addparams!,addargs!,addwhereparams!,extendbody!
using ..Factories: MixEscaped,Escaped,UnEscaped,rawexpr
using ..TypeTraits: efficientoperations

export simpletreedepth,simpletreewidth
export AbstractSimpleTree
export root,parent,children
export addnode!,deletenode!
export isleaf,level
export ancestor,descendants,siblings,leaves
export subtree,move!
export SimpleTreeCore,@simpletree,SimpleTree

abstract type SimpleTreeIteration end
struct SimpleTreeDepth <: SimpleTreeIteration end
struct SimpleTreeWidth <: SimpleTreeIteration end

"""
    simpletreedepth

Indicate that the iteration over a tree is depth-first.
"""
const simpletreedepth=SimpleTreeDepth()

"""
    simpletreewidth

Indicate that the iteration over a tree is width-first.
"""
const simpletreewidth=SimpleTreeWidth()

"""
    AbstractSimpleTree{Node,Data}

Abstract type for all concrete trees.
"""
abstract type AbstractSimpleTree{N,D} end

"""
    eltype(tree::AbstractSimpleTree)
    eltype(::Type{<:AbstractSimpleTree{N,D}}) where {N,D}

Get the eltype of a tree.
"""
Base.eltype(tree::AbstractSimpleTree)=tree|>typeof|>eltype
Base.eltype(::Type{<:AbstractSimpleTree{N,D}}) where {N,D}=Pair{N,D}

"""
    root(tree::AbstractSimpleTree) -> Union{keytype(tree),Nothing}

Get a tree's root node.
"""
root(tree::AbstractSimpleTree)=tree.TREECORE.root

"""
    haskey(tree::AbstractSimpleTree{N},node::N) where N -> Bool

Check whether a node is in a tree.
"""
Base.haskey(tree::AbstractSimpleTree{N},node::N) where N=haskey(tree.TREECORE.contents,node)

"""
    length(tree::AbstractSimpleTree) -> Int

Get the number of a tree's nodes.
"""
Base.length(tree::AbstractSimpleTree)=length(tree.TREECORE.contents)

"""
    parent(tree::AbstractSimpleTree{N},node::N,superparent::Union{N,Nothing}=nothing) where N -> Union{N,Nothing}

Get the parent of a tree's node. When `node` is the tree's root, return `superparent`.
"""
parent(tree::AbstractSimpleTree{N},node::N,superparent::Union{N,Nothing}=nothing) where N=(node==tree|>root ? superparent : tree.TREECORE.parent[node])

"""
    children(tree::AbstractSimpleTree) -> Vector{keytype(tree)}
    children(tree::AbstractSimpleTree,::Nothing) -> Vector{keytype(tree)}
    children(tree::AbstractSimpleTree{N},node::N) where N -> Vector{N}

Get the children of a tree's node.
"""
children(tree::AbstractSimpleTree)=children(tree,nothing)
children(tree::AbstractSimpleTree,::Nothing)=(tree|>root===nothing ? error("children error: empty tree!") : [tree|>root])
children(tree::AbstractSimpleTree{N},node::N) where N=tree.TREECORE.children[node]

"""
    addnode!(tree::AbstractSimpleTree{N},node::N) where N} -> typeof(tree)
    addnode!(tree::AbstractSimpleTree{N},::Nothing,node::N) where N -> typeof(tree)
    addnode!(tree::AbstractSimpleTree{N},parent::N,node::N) where N -> typeof(tree)

Update the structure of a tree by adding a node. When the parent is `nothing`, the input tree must be empty and the input node becomes the tree's root.
"""
addnode!(tree::AbstractSimpleTree{N},node::N) where N=addnode!(tree,nothing,node)
function addnode!(tree::AbstractSimpleTree{N},::Nothing,node::N) where N
    @assert tree|>root===nothing "addnode! error: not empty tree."
    tree.TREECORE.root=node
    tree.TREECORE.children[node]=N[]
    return tree
end
function addnode!(tree::AbstractSimpleTree{N},parent::N,node::N) where N
    @assert haskey(tree,parent) "addnode! error: parent($parent) not in tree."
    @assert !haskey(tree,node) "addnode! error: node($node) already in tree."
    push!(tree.TREECORE.children[parent],node)
    tree.TREECORE.parent[node]=parent
    tree.TREECORE.children[node]=N[]
    return tree
end

"""
    deletenode!(tree::AbstractSimpleTree{N},node::N) where N -> typeof(tree)

Update the structure of a tree by deleting a node.
"""
function deletenode!(tree::AbstractSimpleTree{N},node::N) where N
    @assert haskey(tree,node) "deletenode! error: node($node) not in tree."
    if node==tree|>root
        tree.TREECORE.root=nothing
    else
        pnode=parent(tree,node)
        haskey(tree.TREECORE.children,pnode) && filter!(x->x!=node,tree.TREECORE.children[pnode])
    end
    pop!(tree.TREECORE.contents,node)
    pop!(tree.TREECORE.parent,node,nothing)
    pop!(tree.TREECORE.children,node)
    return tree
end

"""
    getindex(tree::AbstractSimpleTree{N},node::N) where N -> N

Get the data of a tree's node.
"""
Base.getindex(tree::AbstractSimpleTree{N},node::N) where N=tree.TREECORE.contents[node]

"""
    setindex!(tree::AbstractSimpleTree{N,D},data::D,node::N) where {N,D}

Set the data of a tree's node.
"""
Base.setindex!(tree::AbstractSimpleTree{N,D},data::D,node::N) where {N,D}=(tree.TREECORE.contents[node]=data)

"""
    empty(tree::AbstractSimpleTree)

Construct an empty tree of the same type with the input one.
"""
@generated function Base.empty(tree::AbstractSimpleTree{N,D}) where {N,D}
    contents=Expr[:(getfield(tree,$i)) for i=1:length(tree|>fieldnames)-1]
    return :(($tree)($(contents...),SimpleTreeCore{$N,$D}()))
end

"""
    keytype(tree::AbstractSimpleTree)
    keytype(::Type{<:AbstractSimpleTree{N}}) where N

Get a tree's node type.
"""
Base.keytype(tree::AbstractSimpleTree)=tree|>typeof|>keytype
Base.keytype(::Type{<:AbstractSimpleTree{N}}) where N=N

"""
    valtype(tree::AbstractSimpleTree)
    valtype(::Type{<:AbstractSimpleTree{N,D} where N}) where D=D

Get a tree's data type.
"""
Base.valtype(tree::AbstractSimpleTree)=tree|>typeof|>valtype
Base.valtype(::Type{<:AbstractSimpleTree{N,D} where N}) where D=D

"""
    ==(t1::T,t2::T) where T<:AbstractSimpleTree -> Bool
    isequal(t1::T,t2::T) where T<:AbstractSimpleTree -> Bool

Overloaded equivalent operator.
"""
Base.:(==)(t1::T,t2::T) where T<:AbstractSimpleTree = ==(efficientoperations,t1,t2)
Base.isequal(t1::T,t2::T) where T<:AbstractSimpleTree=isequal(efficientoperations,t1,t2)

"""
    keys(tree::AbstractSimpleTree{N},::SimpleTreeDepth,node::Union{N,Nothing}=tree|>root) where N
    keys(tree::AbstractSimpleTree{N},::SimpleTreeWidth,node::Union{N,Nothing}=tree|>root) where N

Iterate over a tree's nodes starting from a certain `node` by depth first search or width first search.
"""
Base.keys(tree::AbstractSimpleTree{N},ti::SimpleTreeIteration,node::Union{N,Nothing}=tree|>root) where N=TreeKeys{typeof(ti),N,valtype(tree),typeof(tree)}(tree,node)
struct TreeKeys{P<:SimpleTreeIteration,N,D,T<:AbstractSimpleTree{N,D}}
    tree::T
    node::Union{N,Nothing}
end
Base.eltype(::Type{TreeKeys{<:SimpleTreeIteration,N,D,<:AbstractSimpleTree{N,D}} where D}) where N=N
Base.IteratorSize(::Type{<:TreeKeys})=Base.SizeUnknown()
Base.iterate(tk::TreeKeys)=tk.tree|>root==nothing ? nothing : (tk.node,copy(children(tk.tree,tk.node)))
function Base.iterate(tk::TreeKeys{SimpleTreeDepth,N},state::Vector{N}) where N
    length(state)==0 ? nothing : (node=popfirst!(state);prepend!(state,children(tk.tree,node));(node,state))
end
function Base.iterate(tk::TreeKeys{SimpleTreeWidth,N},state::Vector{N}) where N
    length(state)==0 ? nothing : (node=popfirst!(state);append!(state,children(tk.tree,node));(node,state))
end

"""
    values(tree::AbstractSimpleTree{N},::SimpleTreeDepth,node::Union{N,Nothing}=tree|>root) where N
    values(tree::AbstractSimpleTree{N},::SimpleTreeWidth,node::Union{N,Nothing}=tree|>root) where N

Iterate over a tree's data starting from a certain `node` by depth first search or width first search.
"""
Base.values(tree::AbstractSimpleTree{N},ti::SimpleTreeIteration,node::Union{N,Nothing}=tree|>root) where N=(tree[key] for key in keys(tree,ti,node))

"""
    pairs(tree::AbstractSimpleTree{N},::SimpleTreeDepth,node::Union{N,Nothing}=tree|>root) where N
    pairs(tree::AbstractSimpleTree{N},::SimpleTreeWidth,node::Union{N,Nothing}=tree|>root) where N

Iterate over a tree's (node,data) pairs starting from a certain `node` by depth first search or width first search.
"""
Base.pairs(tree::AbstractSimpleTree{N},ti::SimpleTreeIteration,node::Union{N,Nothing}=tree|>root) where N=((key,tree[key]) for key in keys(tree,ti,node))

"""
    isleaf(tree::AbstractSimpleTree{N},node::N) where N -> Bool

Judge whether a tree's node is a leaf (a node without children) or not.
"""
isleaf(tree::AbstractSimpleTree{N},node::N) where N=length(children(tree,node))==0

"""
    level(tree::AbstractSimpleTree{N},node::N) where N -> Int

Get the level of tree's node.
"""
function level(tree::AbstractSimpleTree{N},node::N) where N
    result=1
    while node!=tree|>root
        result+=1
        node=parent(tree,node)
    end
    return result
end

"""
    ancestor(tree::AbstractSimpleTree{N},node::N,generation::Int=1) where N -> N

Get the ancestor of a tree's node of the n-th generation.
"""
function ancestor(tree::AbstractSimpleTree{N},node::N,generation::Int=1) where N
    @assert generation>=0 "ancestor error: generation($generation) must be non-negative."
    result=node
    for i=1:generation
        result=parent(tree,result)
    end
    result
end

"""
    descendants(tree::AbstractSimpleTree{N},node::N,generation::Int=1) where N -> Vector{N}

Get the descendants of a tree's node of the n-th generation.
"""
function descendants(tree::AbstractSimpleTree{N},node::N,generation::Int=1) where N
    @assert generation>=0 "descendants error: generation($generation) must be non-negative."
    result=N[node]
    for i=1:generation
        result=vcat((children(tree,node) for node in result)...)
    end
    result
end

"""
    siblings(tree::AbstractSimpleTree{N},node::N) where N -> Vector{N}

Get the siblings (other nodes sharing the same parent) of a tree's node.
"""
siblings(tree::AbstractSimpleTree{N},node::N) where N=filter(x->x!=node,children(tree,parent(tree,node)))

"""
    leaves(tree::AbstractSimpleTree) -> Vector{keytype(tree)}

Get a tree's leaves.
"""
leaves(tree::AbstractSimpleTree)=keytype(tree)[node for node in keys(tree,simpletreedepth) if isleaf(tree,node)]

"""
    push!(tree::AbstractSimpleTree{N,D},node::N,data::D) where {N,D} -> typeof(tree)
    push!(tree::AbstractSimpleTree{N,D},parent::Union{N,Nothing},node::N,data::D) where {N,D} -> typeof(tree)

Push a new node to a tree. When `parent` is `nothing`, this function set the root node of an empty tree.
"""
Base.push!(tree::AbstractSimpleTree{N,D},node::N,data::D) where {N,D}=push!(tree,nothing,node,data)
function Base.push!(tree::AbstractSimpleTree{N,D},parent::Union{N,Nothing},node::N,data::D) where {N,D}
    addnode!(tree,parent,node)
    tree[node]=data
    return tree
end

"""
    append!(tree::AbstractSimpleTree{N,D},subtree::AbstractSimpleTree{N,D}) where {N,D} -> typeof(tree)
    append!(tree::AbstractSimpleTree{N,D},node::Union{N,Nothing},subtree::AbstractSimpleTree{N,D}) where {N,D} -> typeof(tree)

Append a subtree to a tree.
"""
Base.append!(tree::AbstractSimpleTree{N,D},subtree::AbstractSimpleTree{N,D}) where {N,D}=append!(tree,nothing,subtree)
function Base.append!(tree::AbstractSimpleTree{N,D},node::Union{N,Nothing},subtree::AbstractSimpleTree{N,D}) where {N,D}
    for (key,value) in pairs(subtree,simpletreewidth)
        push!(tree,parent(subtree,key,node),key,value)
    end
    return tree
end

"""
    delete!(tree::AbstractSimpleTree{N},node::N) where N -> typeof(tree)

Delete a node and all its descendants from a tree.
"""
function Base.delete!(tree::AbstractSimpleTree{N},node::N) where N
    for key in collect(N,keys(tree,simpletreedepth,node))
        deletenode!(tree,key)
    end
    return tree
end

"""
    empty!(tree::AbstractSimpleTree) -> typeof(tree)

Empty a tree.
"""
Base.empty!(tree::AbstractSimpleTree)=delete!(tree,tree|>root)

"""
    subtree(tree::AbstractSimpleTree{N},node::N) where N -> typeof(tree)

Get a subtree whose root is `node`.
"""
function subtree(tree::AbstractSimpleTree{N},node::N) where N
    result=empty(tree)
    for (i,(key,value)) in enumerate(pairs(tree,simpletreedepth,node))
        push!(result,i==1 ? nothing : parent(tree,key),key,value)
    end
    return result
end

"""
    move!(tree::AbstractSimpleTree{N},node::N,parent::N) where N -> typeof(tree)

Move a subtree to a new position.
"""
function move!(tree::AbstractSimpleTree{N},node::N,parent::N) where N
    sub=subtree(tree,node)
    delete!(tree,node)
    append!(tree,parent,sub)
    return tree
end

"""
    SimpleTreeCore()

The core of a tree.
"""
mutable struct SimpleTreeCore{N,D}
    root::Union{N,Nothing}
    contents::Dict{N,D}
    parent::Dict{N,N}
    children::Dict{N,Vector{N}}
    SimpleTreeCore{N,D}() where {N,D}=new{N,D}(nothing,Dict{N,D}(),Dict{N,N}(),Dict{N,Vector{N}}())
end

"""
    ==(tc1::TC,tc2::TC) where TC<:SimpleTreeCore -> Bool
    isequal(tc1::TC,tc2::TC) where TC<:SimpleTreeCore -> Bool

Overloaded equivalent operator.
"""
Base.:(==)(tc1::TC,tc2::TC) where TC<:SimpleTreeCore = ==(efficientoperations,tc1,tc2)
Base.isequal(tc1::TC,tc2::TC) where TC<:SimpleTreeCore=isequal(efficientoperations,tc1,tc2)

"""
    @simpletree structdef treeparams::Union{Expr,Nothing}=nothing

Decorate a "raw" struct to be a subtype of `AbstractSimpleTree`.
!!! note
    1. A "raw" struct means:
       - It has no explicit supertype;
       - It has no inner constructor;
       - It has no attribute `:TREECORE`.
    2. The keytype and valtype can be assigned by the argument `treeparams` in the form `{keytype,valtype}`.
       - When the formal argument names of keytype and valtype are not assigned, they can be automatically generated by the functioin `gensym`.
         For example, the struct `Tree` after the decoration by the following code
         ```julia
         @simpletree(struct Tree{N} info::Vector{N} end)
         ```
         will have three type parameters.
       - When the keytype or valtype is a definite type, use symbol `:(::)` before that type.
         For example, the struct `Tree` decorated by the following code
         ```julia
         @simpletree(struct Tree{N} info::Vector{N} end,{::String,D<:Number})
         ```
         only has two type parameters `N` and `D<:Number`.
       - When the formal argument names of keytype and valtype overlap with those of the raw struct type parameters, the duplicates will be considered as the same.
         For example, the decorated struct `Tree` by the following code
         ```julia
         @simpletree(struct Tree{D} info::Vector{D} end,{::String,D<:Number})
         ```
         has only one type parameter `D<:Number`, where the `D` in the `info::Vector{D}` is the same `D` with that in the decorated attribute `TREECORE::SimpleTreeCore{String,D}`.
       - When the formal argument names of keytype and valtype have no intersection with those of the raw struct type parameters,
         the type parameters of the decorated struct will be just extended by keytype and valtype.
         For example, the decorated struct `Tree` by the following code
         ```julia
         @simpletree(struct Tree{T} info::Vector{T} end,{N<:Tuple,D<:Number})
         ```
         has 3 type parameters, `T`, `N<:Tuple` and `D<:Number`.
"""
macro simpletree(structdef,treeparams::Union{Expr,Nothing}=nothing)
    tf=TypeFactory(structdef)
    fieldnames=[field.name for field in tf.fields]
    paramnames=[param.name for param in tf.params]
    @assert tf.supertype==Inference(:Any) "@simpletree error: no explicit supertype except `Any` is allowed."
    @assert length(tf.constructors)==0 "@simpletree error: no inner constructor is allowed."
    @assert :TREECORE ∉ fieldnames "@simpletree error: :TREECORE is a reserved attribute name."
    if treeparams==nothing
        treeparamnames=(gensym(),gensym())
        treeparams=[Parameter(treeparamnames[1]),Parameter(treeparamnames[2])]
        truetreeparamnames=treeparamnames
        escapednames=()
    else
        @assert (treeparams.head==:braces && length(treeparams.args)==2) "@simpletree error: not supported treeparams."
        treeparams=[Parameter(arg) for arg in treeparams.args]
        treeparams[1].name===nothing && treeparams[1].head≠:(::) && (treeparams[1].name=gensym())
        treeparams[2].name===nothing && treeparams[2].head≠:(::) && (treeparams[2].name=gensym())
        treeparamnames=tuple((tp.name===nothing ? tp(rawexpr) : tp.name for tp in treeparams)...)
        truetreeparamnames=tuple((tp.name for tp in treeparams if tp.name≠nothing)...)
        escapednames=vcat((names(tp(rawexpr)) for tp in treeparams if tp.name===nothing)...)
        filter!(tp->tp.name∉paramnames && tp.name≠nothing,treeparams)
    end
    tf.supertype=Inference(:(AbstractSimpleTree{$(treeparamnames...)}))
    append!(paramnames,(tp.name for tp in treeparams))
    addfields!(tf,:(TREECORE::SimpleTreeCore{$(treeparamnames...)}))
    addparams!(tf,treeparams...)
    outer=FunctionFactory(name=tf.name)
    addparams!(outer,truetreeparamnames...)
    addargs!(outer,(Argument(name=field.name,type=field.type) for field in tf.fields[1:end-1])...)
    addwhereparams!(outer,tf.params...)
    extendbody!(outer,:($(tf.name)($(fieldnames...),SimpleTreeCore{$(treeparamnames...)}())))
    sbtreedef=tf(MixEscaped(Escaped(tf.name,escapednames...),UnEscaped(paramnames...,:AbstractSimpleTree,:SimpleTreeCore)))
    outerdef=outer(MixEscaped(Escaped(tf.name,escapednames...),UnEscaped(paramnames...)))
    return Expr(:block,:(Base.@__doc__($sbtreedef)),outerdef)
end

"""
    SimpleTree{N,D}() where {N,D}

The minimum tree structure that implements all the default tree methods.
"""
struct SimpleTree{N,D} <: AbstractSimpleTree{N,D}
    TREECORE::SimpleTreeCore{N,D}
end
SimpleTree{N,D}() where {N,D}=SimpleTree(SimpleTreeCore{N,D}())

end #module


text:
"""
  This function tests GLS for a given state loaded in the eqn object.
  This function is not a test function, but is called by test functions.
"""
function test_GLS(mesh::AbstractMesh{Tmsh}, sbp, eqn::AbstractSolutionData{Tsol, Tres}, opts) where {Tsol, Tres, Tmsh}

  eqn.params.tau_type = 2
  Dxi = diagm(1./sbp.w)*sbp.Q[:, :, 1]
  Deta = diagm(1./sbp.w)*sbp.Q[:, :, 2]

  # create indices
  idx_range = Array{UnitRange{Int64}}(mesh.numNodesPerElement)
  for i=1:mesh.numNodesPerElement
    start_idx = (i-1)*mesh.numDofPerNode + 1
    end_idx = i*mesh.numDofPerNode
    idx_range[i] = copy(start_idx:end_idx)
  end

#    println("idx_range = ", idx_range)

    size_block = mesh.numNodesPerElement*mesh.numDofPerNode
#    println("size_block = ", size_block)

    # get GLS from the code
    fill!(eqn.res, 0.0)
    EulerEquationMod.applyGLS3(mesh, sbp, eqn, opts)

  # testing: only do one element
  for el =1:mesh.numEl
#    println("testing element ", el)

    # constant mapping elements only
    dxidx = zeros(2,2)
    dxidx_hat_el = sview(mesh.dxidx, :, :, 1, el)
    jac_el = sview(mesh.jac, :, el)
    q_el = reshape(copy(eqn.q[:, :, el]), size_block)

    for i=1:2
      for j=1:2
        dxidx[i,j] = dxidx_hat_el[i, j, 1]*jac_el[1]
      end
    end

    # calculate Dx, Dy
    Dx = dxidx[1,1]*Dxi + dxidx[2, 1]*Deta
    Dy = dxidx[1,2]*Dxi + dxidx[2, 2]*Deta

#    println("Dx = \n", Dx)
#    println("Dy = \n", Dy)

    # create block Dx, Dy
    Dx_tilde = zeros(Tmsh, size_block, size_block)
    Dy_tilde = zeros(Dx_tilde)

    for i=1:mesh.numNodesPerElement
      idx_i = idx_range[i]
#      println("idx_i = ", idx_i)
      for j=1:mesh.numNodesPerElement
        idx_j = idx_range[j]
#        println("  idx_j = ", idx_j)
        Dx_tilde[idx_i, idx_j] = Dx[i, j]*eye(mesh.numDofPerNode)
        Dy_tilde[idx_i, idx_j] = Dy[i, j]*eye(mesh.numDofPerNode)
      end
    end

#    println("Dx_tilde = \n", Dx_tilde)
#    println("Dy_tilde = \n", Dy_tilde)

    # create A1 tilde and A2 tilde
    A1_tilde = zeros(Tsol, size_block, size_block)
    A2_tilde = zeros(Tsol, size_block, size_block)

    for i=1:mesh.numNodesPerElement
      idx_i = idx_range[i]
      q_i = q_el[idx_i]

      A1 = sview(A1_tilde, idx_i, idx_i)
      EulerEquationMod.calcA1(eqn.params, q_i, A1)

      A2 = sview(A2_tilde, idx_i, idx_i)
      EulerEquationMod.calcA2(eqn.params, q_i, A2)
    end

#    println("A1 = \n", A1_tilde)
#    println("A2 = \n", A2_tilde)

    # create middle terms, including tau
    middle_tilde = zeros(Tres, size_block, size_block)

    for i=1:mesh.numNodesPerElement

      idx_i = idx_range[i]
      tau = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode)
      EulerEquationMod.getTau(eqn.params, jac_el[i], tau)  # tau number 2

      middle_tilde[idx_i, idx_i] = (sbp.w[i]/jac_el[i])*tau
    end

    # create the operator
    fancy_L  = A1_tilde*Dx_tilde + A2_tilde*Dy_tilde
    gls_operator = fancy_L.'*middle_tilde*fancy_L

    @test ( isSymmetric(gls_operator, 1e-12) )== true
    #println("max asymmetry = ", maximum(abs.(gls_operator - gls_operator.')))

    gls_test = -gls_operator*q_el

    gls_code = reshape(copy(eqn.res[:, :, el]), size_block)

    @test isapprox( gls_code, gls_test) atol=1e-12
#    println("gls_test = \n", gls_test)
#    println("gls_code = \n", gls_code)

    trial_term = fancy_L*q_el
    middle_term = middle_tilde*trial_term
    gls_term = -fancy_L.'*middle_term
#    println("element $el trial term = \n", trial_term)
#    println("elemetn $el middle_term = \n", middle_term)
#    println("element $el gls_term = \n", gls_term)

    @test isapprox( gls_term, gls_test) atol=1e-12
    # test matrix transpose
    tmp1 = A1_tilde.'*middle_term
    tmp2 = A2_tilde.'*middle_term
    gls_test2 = -(Dx_tilde.'*tmp1 + Dy_tilde.'*tmp2)
    @test isapprox( gls_test2, gls_test) atol=1e-12

  end  # end loop over elements

  return nothing

end  # end function

# run the tests
"""
  Test GLS on a channel flow
"""
function test_gls_channel(mesh, sbp, eqn, opts)
  #=
  include("input_vals_channel.jl")
  arg_dict["solve"] = false
  arg_dict["variable_type"] = :entropy
  f = open("input_vals_channel_gls.jl", "w")
  println(f, "arg_dict = ")
  println(f, arg_dict)
  close(f)
  =#

  @testset "----- Testing GLS3 channel -----" begin
 #   fname = "input_vals_channel_gls.jl"
 #   mesh, sbp, eqn, opts = solvePDE(fname)
    array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)
    test_GLS(mesh, sbp, eqn, opts)
  end  # end facts block

  return nothing
end

test_gls_channel_inputname = "input_vals_channel.jl"
test_gls_channel_moddict = Dict{String, Any}("solve" => false, "variable_type" => :entropy, "new_fname" => "input_vals_channel_gls")
#test_gls_channel(mesh, sbp, eqn, opts)
add_func3!(EulerTests, test_gls_channel, test_gls_channel_inputname, test_gls_channel_moddict, [TAG_ENTROPYVARS, TAG_SHORTTEST])

"""
  Test GLS on the isentropic vortex for all degree operators
"""
function test_gls_vortex()
  for p = 1:4
    if true
      @testset "----- Testing GLS3 p$p  on isentropic vortex -----" begin
        # test on isentropic vortex
        arg_dict = read_input_file("input_vals_vortex3.jl")
        arg_dict["order"] = p
        arg_dict["solve"] = false
        arg_dict["variable_type"] = :entropy
        #TODO: use make_input function instead
        make_input(arg_dict, "input_vals_vortex3_gls")

        fname = "input_vals_vortex3_gls.jl"
        mesh, sbp, eqn, opts = solvePDE(fname)
        array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)
        test_GLS(mesh, sbp, eqn ,opts)
      end  # end facts block
    end
  end

  return nothing
end

#test_gls_vortex(mesh, sbp, eqn, opts)
#add_func1!(EulerTests, test_gls_vortex, [TAG_ENTROPYVARS, TAG_SHORTTEST])

"""
  Test finite differencing and complex stepping of the GLS term
"""
function test_gls_fd()
  for p = 1:4
    if true
      @testset "----- Performing GLS3 p$p finite difference checks -----" begin
        fname = "input_vals_vortex3_gls.jl"
        mesh, sbp, eqn, opts = solvePDE(fname)

        arg_dict["order"] = p
        f = open("input_vals_vortex3_gls.jl", "w")
        print(f, "arg_dict = ")
        println(f, arg_dict)
        close(f)
        fname = "input_vals_vortex3_gls.jl"
        mesh, sbp, eqn, opts = solvePDE(fname)

        # rescale the problem
        for i = 1:length(eqn.q)
          eqn.q[i] = 1000*eqn.q[i]
        end

        # trick code into only doing the first element
        mesh.numEl = 1

        len = mesh.numDofPerNode*mesh.numNodesPerElement
        array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)
        test_GLS(mesh, sbp, eqn, opts)
        jac_fd = zeros(len, len)
        eps_fd = 1e-8
        res0 = copy(reshape(eqn.res[:, :, 1], len))  # use res from previous run
        println("doing finite differences")
        for j=1:mesh.numNodesPerElement
          for i=1:mesh.numDofPerNode
            pos = (j-1)*mesh.numDofPerNode + i
            eqn.q[i, j, 1] += eps_fd
            test_GLS(mesh, sbp, eqn, opts)
            res_ij = copy(reshape(eqn.res[:, :, 1], len))
            jac_fd[:, pos] = (res_ij - res0)/eps_fd
            eqn.q[i, j, 1] -= eps_fd  # undo perturbation
          end
        end

        # do complex step
        println("doing complex step")
        arg_dict["run_type"] = 5
        arg_dict["jac_method"] = 2
        f = open("input_vals_vortex3c_gls.jl", "w")
        print(f, "arg_dict = ")
        println(f, arg_dict)
        close(f)
        fname = "input_vals_vortex3c_gls.jl"
        mesh, sbp, eqn, opts = solvePDE(fname)

        # rescale the problem
        for i = 1:length(eqn.q)
          eqn.q[i] = 1000*eqn.q[i]
        end

        # trick code into only doing 1 element
        mesh.numEl = 1

        array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)
        eqn.params.tau_type = p

        eps_c = 1e-20
        jac_c = zeros(len, len)
        for j=1:mesh.numNodesPerElement
          for i=1:mesh.numDofPerNode
            pos = (j-1)*mesh.numDofPerNode + i
            eqn.q[i, j, 1] += complex(0, eps_c)
            test_GLS(mesh, sbp, eqn, opts)
            res_ij = copy(reshape(eqn.res[:, :, 1], len))
            jac_c[:, pos] = imag(res_ij)/eps_c
            eqn.q[i, j, 1] -= complex(0, eps_c)  # undo perturbatino
          end
        end

        for j=1:len
          tol = 0.003
          for k = 1:len
            if abs(jac_fd[k,j]) > 1e-4
              @test  abs((jac_c[k, j] - jac_fd[k, j])/jac_c[k,j])  < tol
            end
          end
        end

    #=
        for j=1:12
          println("column $j")
          println("jac_c = \n", jac_c[:, j])
          println("jac_fd = \n", jac_fd[:, j])
          println("jac diff = \n", jac_c[:,j] - jac_fd[:, j])
        end
    =#
      end  # end facts block
    end  # end if statement
  end  # end p=1:4

  return nothing

end  # end function

#test_gls_fd()
#add_func1!(EulerTests, test_gls_fd, [TAG_COMPLEX, TAG_SHORTTEST])

ext: jl


text:
# This file was generated, do not modify it. # hide
Q = MatrixDatabase(randn(2, 100))
k = 30
I, D = searchbatch(G, Q, k)
@show typeof(I) => size(I)
@show typeof(D) => size(D)

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
module StabilizingDissipations

using StaticArrays
using LinearAlgebra

using ClimateMachine.Mesh.Grids: min_node_distance, HorizontalDirection

import ClimateMachine.Ocean.HydrostaticBoussinesq: diffusivity_tensor, viscosity_tensor

using ClimateMachine.Ocean.HydrostaticBoussinesq: HydrostaticBoussinesqModel

struct StabilizingDissipation{T}
    κʰ_min :: T
    κʰ_max :: T
    νʰ_min :: T
    νʰ_max :: T
    smoothness_exponent :: T
    minimum_node_spacing :: T
    Δu :: T
    Δθ :: T
end

function StabilizingDissipation(;
                                time_step,
                                Δu,
                                Δθ,
                                minimum_node_spacing,
                                diffusive_cfl = 0.1,
                                κʰ_min = 0,
                                νʰ_min = 0,
                                smoothness_exponent = 2)

    FT = Float64 #eltype(grid)

    κʰ_max = diffusive_cfl * minimum_node_spacing^2 / time_step
    νʰ_max = diffusive_cfl * minimum_node_spacing^2 / time_step

    return StabilizingDissipation{FT}(FT(κʰ_min),
                                      FT(κʰ_max),
                                      FT(νʰ_min),
                                      FT(νʰ_max),
                                      FT(smoothness_exponent),
                                      FT(minimum_node_spacing),
                                      FT(Δu),
                                      FT(Δθ))
end

@inline function viscosity_tensor(m::HydrostaticBoussinesqModel{<:StabilizingDissipation}, ∇u)
    ϰ = m.stabilizing_dissipation
    Δh = ϰ.minimum_node_spacing

    # ∇u is a 2 × 3 tensor (eg ∇uʰ, where ∇ is 3D)
    @inbounds begin
        ux = ∇u[1, 1]
        uy = ∇u[1, 2]
        vx = ∇u[2, 1]
        vy = ∇u[2, 2]
    end

    ∇u² = ux^2 + uy^2 + vx^2 + vy^2

    arg = (sqrt(∇u²) * Δh / ϰ.Δu)^ϰ.smoothness_exponent

    νʰ = m.νʰ + ϰ.νʰ_min + (ϰ.νʰ_max - ϰ.νʰ_min) * tanh(arg)

    return Diagonal(@SVector [νʰ, νʰ, m.νᶻ])
end

@inline function diffusivity_tensor(m::HydrostaticBoussinesqModel{<:StabilizingDissipation}, ∇θ)

    @inbounds begin
        θx = ∇θ[1]
        θy = ∇θ[2]
        θz = ∇θ[3]
    end

    # Compute convective adjustement diffusivity
    θz < 0 ? κᶻ = m.κᶜ : κᶻ = m.κᶻ

    # Compute stabilizing dissipation
    ϰ = m.stabilizing_dissipation
    Δh = ϰ.minimum_node_spacing

    ∇θ² = θx^2 + θy^2 + θz^2

    arg = (sqrt(∇θ²) * Δh / ϰ.Δθ)^ϰ.smoothness_exponent

    κʰ = m.κʰ + ϰ.κʰ_min + (ϰ.κʰ_max - ϰ.κʰ_min) * tanh(arg)

    return Diagonal(@SVector [κʰ, κʰ, κᶻ])
end

end # module


text:
using LoopFieldCalc

# Define current loop with radius of 50 cm, current of 1 A, centered at the origin
loop = LoopFieldCalc.CurrentLoop(
  radius = 0.5,         # radius in meters
  current = 1.0,        # current in Amperes
  center = LoopFieldCalc.CartesianPoint(0.0, 0.0, 0.0)    # cartesian coordinates of the loop center
)

# Define grid points along x, y, and z axes. We have a 2D grid of points in the y-z plane.
xs = 0.0:0.0
ys = -1.0:0.01:1.0
zs = -2.0:0.01:2.0

# Define current loops. Here, we only have the single loop, defined above
loops = [loop]

# Compute the magnetic field components
Bx, By, Bz = LoopFieldCalc.field_on_grid(loops, xs, ys, zs)

LoopFieldCalc.write_field("examples/example_1.dat", ys, zs, By[1, :, :], Bz[1, :, :])


text:
module SfPhysics

using Reexport

include("units.jl")

@reexport using .SfUnits

include("kinematics.jl")
include("fluid-dynamics.jl")

@reexport using .SfFluidDynamics

include("geometry.jl")

@reexport using .SfGeometry

include("gravity.jl")
include("blackhole.jl")

@reexport using .SfGravity
@reexport using .SfBlackHole

include("rocketry.jl")

@reexport using .SfRocketry

include("relativity.jl")

@reexport using .SfRelativity

include("matter.jl")

@reexport using .SfMatter

include("materials.jl")

@reexport using .SfMaterials

include("objects.jl")

@reexport using .SfObjects

include("astronomy.jl")

@reexport using .SfAstronomy

include("thermo.jl")

@reexport using .SfThermo

include("planetary.jl")
include("solar.jl")

@reexport using .SfPlanetary

export SfSolarSystem

include("atmosphere.jl")

@reexport using .SfAtmosphere

include("laser.jl")

@reexport using .SfLaser

include("coriolis.jl")

@reexport using .SfCoriolis

include("electro.jl")

@reexport using .SfElectro

include("weaponry.jl")

@reexport using .SfWeaponry

end

max_stars_repo_name: icantbelieveitsnotphysics/SfPhysics

max_issues_repo_name: icantbelieveitsnotphysics/SfPhysics

max_forks_repo_name: icantbelieveitsnotphysics/SfPhysics


text:
using Test
using DFTK
using ForwardDiff
import FiniteDiff
using MPI
include("testcases.jl")

# Hellmann-Feynman stress
# via ForwardDiff & custom FFTW overloads on ForwardDiff.Dual

@testset "ForwardDiff stresses on silicon" begin
    function make_basis(lattice, symmetries)
        model = model_PBE(lattice, silicon.atoms, silicon.positions; symmetries)
        PlaneWaveBasis(model; Ecut=7, kgrid=(3, 3, 3))
    end

    function recompute_energy(lattice, symmetries)
        basis = make_basis(lattice, symmetries)
        scfres = self_consistent_field(basis; is_converged=DFTK.ScfConvergenceDensity(1e-13))
        energies, H = energy_hamiltonian(basis, scfres.ψ, scfres.occupation; ρ=scfres.ρ)
        energies.total
    end

    function hellmann_feynman_energy(scfres, lattice, symmetries)
        basis = make_basis(lattice, symmetries)
        ρ = DFTK.compute_density(basis, scfres.ψ, scfres.occupation)
        energies, H = energy_hamiltonian(basis, scfres.ψ, scfres.occupation; ρ=ρ)
        energies.total
    end

    a = 10.0  # slightly compressed and twisted
    lattice = a / 2 * [[0 1 1.1];
                       [1 0 1.];
                       [1 1 0.]]
    is_converged = DFTK.ScfConvergenceDensity(1e-11)
    scfres = self_consistent_field(make_basis(lattice, true); is_converged)
    scfres_nosym = self_consistent_field(make_basis(lattice, false); is_converged)
    stresses = compute_stresses_cart(scfres)
    @test isapprox(stresses, compute_stresses_cart(scfres_nosym), atol=1e-10)

    dir = MPI.bcast(randn(3, 3), 0, MPI.COMM_WORLD)

    dE_stresses = dot(dir, stresses) * scfres.basis.model.unit_cell_volume
    ref_recompute = FiniteDiff.finite_difference_derivative(0.0) do ε
        recompute_energy(lattice + ε*dir*lattice, false)
    end
    ref_HF = FiniteDiff.finite_difference_derivative(0.0) do ε
        hellmann_feynman_energy(scfres_nosym, lattice+ε*dir*lattice, false)
    end
    FD_HF = ForwardDiff.derivative(0.0) do ε
        hellmann_feynman_energy(scfres_nosym, lattice+ε*(dir*lattice), false)
    end

    @test isapprox(ref_HF, ref_recompute, atol=1e-5)
    @test isapprox(ref_HF, FD_HF, atol=1e-5)
    @test isapprox(dE_stresses, ref_recompute, atol=1e-5)
end


text:
using SelectTP
using Test

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


text:
function getConfigDescriptor()

    return Dict(

        "MODEL_MISC" => [

            ConfigEntry(
                "init_file",
                :optional,
                [nothing, String],
                nothing
            ),

            ConfigEntry(
                "rpointer_file",
                :optional,
                [String,],
                "rpointer.iom";
                desc="If read_restart of the init function is set true, then this entry has to contain a valid rpointer filename.",
            ),

            ConfigEntry(
                "enable_archive",
                :required,
                [Bool,],
            ),

            ConfigEntry(
                "daily_record",
                :optional,
                [AbstractArray,],
                [],
            ),

            ConfigEntry(
                "monthly_record",
                :optional,
                [AbstractArray,],
                [],
            ),

        ],

        "MODEL_CORE" => [

            ConfigEntry(
                "timetype",
                :optional,
                [String,],
                "DateTimeNoLeap",
            ),

            ConfigEntry(
                "domain_file",
                :required,
                [String,],
            ),

            ConfigEntry(
                "topo_file",
                :optional,
                [String, Nothing],
                nothing,
            ),

            ConfigEntry(
                "cdata_var_file_map",
                :optional,
                [Dict, Nothing],
                nothing,
            ),

            ConfigEntry(
                "cdata_beg_time",
                :optional,
                [Any,],
            ),

            ConfigEntry(
                "cdata_end_time",
                :optional,
                [Any,],
            ),

            ConfigEntry(
                "cdata_align_time",
                :optional,
                [Any,],
            ),

            ConfigEntry(
                "substeps",
                :optional,
                [Integer,],
                8;
            ),

            ConfigEntry(
                "advection_scheme",
                :required,
                ["static", "ekman_KSC2018", "ekman_CO2012", "ekman_AGA2020", "ekman_AGA2020_allowU"],
            ),

            ConfigEntry(
                "MLD_scheme",
                :required,
                ["prognostic", "datastream", "static"],
            ),

            ConfigEntry(
                "Qflx",
                :optional,
                ["on", "off"],
                "off",
            ),

            ConfigEntry(
                "Qflx_finding",
                :optional,
                ["on", "off"],
                "off",
            ),

            ConfigEntry(
                "weak_restoring",
                :optional,
                ["on", "off"],
                "off",
            ),

            ConfigEntry(
                "convective_adjustment",
                :optional,
                ["on", "off"],
                "off",
            ),

            ConfigEntry(
                "z_w",
                :optional,
                [AbstractArray{Float64, 1}],
                nothing;
                desc = "Will be overwritten if `init_file` is used.",
            ),

            ConfigEntry(
                "Ks_H",
                :optional,
                [Float64,],
                1e3;
                desc = "Horizontal tracer diffusivity. Will be overwritten if `init_file` is used.",
            ),

            ConfigEntry(
                "Ks_V",
                :optional,
                [Float64,],
                1e-4;
                desc = "Vertical tracer diffusivity. Will be overwritten if `init_file` is used.",
            ),

            ConfigEntry(
                "Ks_V_cva",
                :optional,
                [Float64,],
                1.0;
                desc = "Convective adjustment diffusivity. Will be overwritten if `init_file` is used.",
            ),

            ConfigEntry(
                "τwk_TEMP",
                :optional,
                [Float64,],
                Inf;
                desc = "Timescale of weak-restoring of temperature. Will be overwritten if `init_file` is used",
            ),

            ConfigEntry(
                "τwk_SALT",
                :optional,
                [Float64,],
                Inf;
                desc = "Timescale of weak-restoring of salinity. Will be overwritten if `init_file` is used",
            ),

            ConfigEntry(
                "τfrz",
                :optional,
                [Float64,],
                3600.0;
                desc = "Timescale of surface water below freezing point.",
            ),

            ConfigEntry(
                "rad_R",
                :optional,
                [Float64,],
                0.58;
                desc = "Fast absorption portion of sunlight as described in Paulson & Simpson (1977). Will be overwritten if `init_file` is used",
            ),

            ConfigEntry(
                "rad_ζ1",
                :optional,
                [Float64,],
                0.15;
                desc = "Light penetration length scale as described in Paulson & Simpson (1977). Will be overwritten if `init_file` is used",
            ),

            ConfigEntry(
                "rad_ζ2",
                :optional,
                [Float64,],
                23.0;
                desc = "Light penetration length scale as described in Paulson & Simpson (1977). Will be overwritten if `init_file` is used",
            ),

            ConfigEntry(
                "ϵ",
                :optional,
                [Float64,],
                1.4e-5;
                desc = "Rayleigh friction of momentum. Will be overwritten if `init_file` is used",
            ),

            ConfigEntry(
                "γ",
                :optional,
                [Float64,],
                0.06;
                desc = "Rayleigh friction of momentum. Will be overwritten if `init_file` is used",
            ),

            ConfigEntry(
                "Ekman_layers",
                :optional,
                [Integer,],
                5;
                desc = "Number of Ekman layers. Will be overwritten if `init_file` is used",
            ),

            ConfigEntry(
                "Returnflow_layers",
                :optional,
                [Integer,],
                5;
                desc = "Number of Ekman return flows. Will be overwritten if `init_file` is used",
            ),

            ConfigEntry(
                "transform_vector_field",
                :optional,
                [Bool,],
                true;
                desc = "If this is set true, then TAUX_east and TAUY_north are considered pointing to true east and north. Thus, a vector transformation onto grid-x and -y direction will be performed. If this is set false, TAUX_east and TAUY_north are considered aligned with grid-x and -y so no transformation will be performed.",
            ),

        ],
    )
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
# Test coherent state
o2 = osc(15)
v = coherentvec(o2, 1)
@test v ⋅ v ≈ 1
@test expectvec(v, o2('n')) ≈ 1 atol=1e-11
@test expectvec(v, o2('n')) ≈ weakvaluevec(v, v, o2('n'))
@test expectvec(v, o2('d')) ≈ weakvaluevec(v, groundvec(o2), o2('d')) atol=1e-11
@test weakvaluevec(v, groundvec(o2), o2('n')) ≈ 0

vp = projector(v)
@test vp ≈ coherent(o2, 1)

@test tr(vp * vp) ≈ 1
@test expect(vp, o2('n')) ≈ 1 atol=1e-11
@test expect(vp, o2('n')) ≈ weakvalue(vp, vp, o2('n'))
@test expect(vp, o2('d')) ≈ weakvalue(vp, ground(o2), o2('d')) atol=1e-11
@test weakvalue(vp, ground(o2), o2('n')) ≈ 0


text:
module _SpecEn
export SpecEn
using FFTW: fft
using DSP: conv
    """
        Spec, BandEn = SpecEn(Sig)

    Returns the spectral entropy estimate of the full spectrum (`Spec`)
    and the within-band entropy (`BandEn`) estimated from the data sequence (`Sig`)
    using the default  parameters:
    N-point FFT = 2*len(`Sig`) + 1, normalised band edge frequencies = [0 1],
    logarithm = base 2, normalisation = w.r.t # of spectrum/band frequency values.

        Spec, BandEn = SpecEn(Sig::AbstractArray{T,1} where T<:Real; N::Int=1 + (2*size(Sig,1)), Freqs::Tuple{Real,Real}=(0,1), Logx::Real=exp(1), Norm::Bool=true)

    Returns the spectral entropy (`Spec`) and the within-band entropy (`BandEn`)
    estimate for the data sequence (`Sig`) using the specified 'keyword' arguments:

    # Arguments:
    `N`     - Resolution of spectrum (N-point FFT), an integer > 1  \n
    `Freqs` - Normalised band edge frequencies, a 2 element tuple with values \n
              in range [0 1] where 1 corresponds to the Nyquist frequency (Fs/2).
              Note: When no band frequencies are entered, BandEn == SpecEn
    `Logx`  - Logarithm base, a positive scalar (enter 0 for natural log) \n
    `Norm`  - Normalisation of `Spec` value:\n
              [false]  no normalisation.
              [true]   normalises w.r.t # of spectrum/band frequency values - default.

    For more info, see the EntropyHub guide.

    # See also `XSpecEn`, `fft`, `MSEn`,  `XMSEn`

    # References:
        [1] G.E. Powell and I.C. Percival,
            "A spectral entropy method for distinguishing regular and
            irregular motion of Hamiltonian systems."
            Journal of Physics A: Mathematical and General
            12.11 (1979): 2053.

        [2] Tsuyoshi Inouye, et al.,
            "Quantification of EEG irregularity by use of the entropy of
            the power spectrum."
            Electroencephalography and clinical neurophysiology
            79.3 (1991): 204-210.

    """
    function SpecEn(Sig::AbstractArray{T,1} where T<:Real; N::Int=1 + (2*size(Sig,1)),
        Freqs::Tuple{Real,Real}=(0,1), Logx::Real=exp(1), Norm::Bool=true)

    (size(Sig)[1] > 4) ? nothing : error("Sig:   must be a numeric vector")
    (N > 1) ? nothing :  error("N:     must be an integer > 1")
    (0<=Freqs[1]<1 && 0<Freqs[2]<=1 && Freqs[1]<Freqs[2]) ? nothing :
        error("Freq:    must be a two element tuple with values in range [0 1].
                The values must be in increasing order.")
    (Logx>0) ? nothing :   error("Logx:     must be a positive number > 0")

    Freqs = collect(Freqs)
    Fx = Int(ceil(N/2))
    Freqs = Int.(round.(Freqs.*Fx))
    Freqs[Freqs.==0] .= 1

    if Freqs[1] > Freqs[2]
        error("Lower band frequency must come first.")
    elseif Freqs[2]-Freqs[1]<1
        error("Spectrum resoution too low to determine bandwidth.")
    elseif minimum(Freqs)<0 || maximum(Freqs)>Fx
        error("Freqs must be normalized w.r.t sampling frequency [0 1].")
    end

    Temp = conv(Sig,Sig)
    N <= size(Temp,1) ? Temp = Temp[1:N]  : Temp = vcat(Temp,zeros(N-size(Temp)[1]))

    Pt = abs.(fft(Temp))
    Pxx = Pt[1:Fx]/sum(Pt[1:Fx])
    Spec = -(transpose(Pxx)*log.(Logx, Pxx))
    Pband = (Pxx[Freqs[1]:Freqs[2]])/sum(Pxx[Freqs[1]:Freqs[2]])
    BandEn = -(transpose(Pband)*log.(Logx, Pband))

    if Norm
        Spec = Spec/(log(Logx, Fx));
        BandEn = BandEn/(log(Logx, Freqs[2]-Freqs[1]+1));
    end

    return Spec, BandEn
    end

end

"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

     http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub

"""


text:
function runcompare()
    println("Running deep dag comparison")
    println("Node Meta: ", ENV["NODE_META_IMPL"])
    println("Cluster Meta: ", ENV["CLUSTER_META_IMPL"])
    J = joinpath(JULIA_HOME, "julia")
    D = dirname(@__FILE__)
    F = map(x->joinpath(D,x), ("dagger.jl", "dagscheduler.jl", "dagscheduler_2node.jl"))
    #F = map(x->joinpath(D,x), ("dagger.jl", "dagscheduler_2node.jl"))

    for f in F
        command = `$J $f`
        println("- $(basename(f))")
        run(command)
    end
end

ENV["NODE_META_IMPL"] = "DagScheduler.ShmemMeta.ShmemExecutorMeta"
ENV["CLUSTER_META_IMPL"] = "DagScheduler.SimpleMeta.SimpleExecutorMeta"
runcompare()

ENV["NODE_META_IMPL"] = "DagScheduler.ShmemMeta.ShmemExecutorMeta"
ENV["CLUSTER_META_IMPL"] = "DagScheduler.FdbMeta.FdbExecutorMeta"
runcompare()


text:
function simple_call(fname, args...; kwargs...)
    try
        if endswith(fname, ".")
            fname = chop(fname);
            # f = eval(Main, parse(fname));
            f = funcfind(fname);
            r = broadcast(f, args...);
        else
            # f = eval(Main, parse(fname));
            f = funcfind(fname);
            r = f(args...; kwargs...);
        end
        sexp(r)
    catch e
        Rerror(e, stacktrace(catch_backtrace())).p;
    end;
end

# function simple_call(fname, arg1, arg2)
#     try
#         if endswith(fname, ".")
#             fname = chop(fname);
#             # f = eval(Main, parse(fname));
#             f = funcfind(fname);
#             r = broadcast(f, arg1, arg2);
#         else
#             # f = eval(Main, parse(fname));
#             f = funcfind(fname);
#             r = f(arg1, arg2);
#         end
#         RObject(r).p
#     catch e
#         Rerror(e, stacktrace(catch_backtrace())).p;
#     end;
# end

max_stars_repo_name: KevCaz/JuliaCall

max_issues_repo_name: tinagui/JuliaCall

max_forks_repo_name: tinagui/JuliaCall


text:
mutable struct Core

    wksp      :: Workspace

#    mask_sT   :: AbstractArray{Float64, 3}
#    mask_T    :: AbstractArray{Float64, 3}

    amo_slab  :: Union{AdvancedMatrixOperators, Nothing}
    amo       :: Union{AdvancedMatrixOperators, Nothing}

    mtx       :: Dict

    vd        :: VerticalDiffusion

    cdatam     :: Union{CyclicDataManager, Nothing}

    function Core(
        ev :: Env,
        tmpfi :: TempField,
    )

        cfg = ev.config

        wksp = Workspace(Nx=ev.Nx, Ny=ev.Ny, Nz=ev.Nz)

        gf = PolelikeCoordinate.CurvilinearSphericalGridFile(
            cfg["domain_file"];
            R  = Re,
            Ω  = Ω,
        )

        gd      = ev.gd
        gd_slab = ev.gd_slab

        amo_slab = AdvancedMatrixOperators(;
            gd = gd_slab,
            mask_T     = ev.topo.sfcmask_sT,
            deepmask_T = ev.topo.sfcmask_sT,
        )

        amo = AdvancedMatrixOperators(;
            gd = gd,
            mask_T = ev.topo.mask_T,
            deepmask_T = ev.topo.deepmask_T,
        )

        # Build Advection Matrix
        function build!(id_mtx, idx)
            local result
            # using transpose speeds up by 100 times
            tp = transpose(id_mtx) |> sparse
            result = transpose(tp[:, view(idx, :)]) |> sparse
            dropzeros!(result)

            idx .= 0 # clean so that debug is easir when some girds are not assigned
            return result
        end

        # Build the matrix to broadcast sW to W grid
        # Notice that we cannot use W_pts because gd_slab has two W layers
        num_sT = reshape( collect(1:amo_slab.bmo.T_pts), amo_slab.bmo.T_dim...)
        mapping_T = repeat(num_sT, outer=(ev.Nz+1, 1, 1))
        W_broadcast_sT = build!(amo_slab.bmo.T_I_T, mapping_T)

        # Build Radiation Matrix
        swflx_factor_W =  cfg["rad_R"]  * exp.(gd.z_W / cfg["rad_ζ1"]) + (1.0 - cfg["rad_R"]) * exp.(gd.z_W / cfg["rad_ζ2"])

        Nz_bot = view(ev.topo.Nz_bot_sT, 1, :, :)
        #println("Nz_bot: ", Nz_bot[20, 114])
        #println("Before: ", swflx_factor_W[:, 20, 114])
        for j=1:ev.Ny, i=1:ev.Nx
            swflx_factor_W[(Nz_bot[i, j]+1):end, i, j] .= 0.0 # Bottom absorbs everything
        end
        #println("After: ", swflx_factor_W[:, 20, 114])

        #swflx_factor_W .= 0
        #swflx_factor_W[1, :, :] .= 1.0

        # Surface flux is immediately absorbed at the top layer
        sfcflx_factor_W = 0.0 * gd.z_W
        sfcflx_factor_W[1, :, :] .= 1.0

        #Ks_H matrix
#        Ks_H_U = zeros(Float64, amo.bmo.V_dim...)

        # This design is to prevent the equator develop wind-Ekman-SST feedback
        # that creates spurious rainband in Central Pacific
        Ks_H_func = (ϕ, z) -> 500.0 + (20000.0 - 500.0) * exp( - 0.5 * (ϕ/deg2rad(10.0))^2 ) * exp(z/100.0)

#        println(size(gd.ϕ_U))
#        println(size(gd.ϕ_V))
#        println(size(gd.ϕ_UV))
#        println(size(gd.z_V))

        Ks_H_U = reshape(Ks_H_func.(gd.ϕ_U, gd.z_U), amo.bmo.U_dim...)
        Ks_H_V = reshape(Ks_H_func.(gd.ϕ_V, gd.z_V), amo.bmo.V_dim...)
        Ks_H_T = reshape(Ks_H_func.(gd.ϕ_T, gd.z_T), amo.bmo.T_dim...)

        # f and ϵ matrices
        f_sT = 2 * gd.Ω * sin.(gd_slab.ϕ_T)
        β_sT = (2 * gd.Ω / gd.R) * cos.(gd_slab.ϕ_T)
        ϵ_sT = f_sT * 0 .+ cfg["ϵ"]
        D_sT = f_sT.^2 + ϵ_sT.^2
        invD_sT = D_sT.^(-1.0)

        mtx = Dict(
            :ones_T          => ones(Float64, amo.bmo.T_pts),
            :T_swflxConv_sT  => - amo.T_mask_T * amo.T_DIVz_W * spdiagm(0 => view(swflx_factor_W, :)) * W_broadcast_sT,
            :T_sfcflxConv_sT => - amo.T_mask_T * amo.T_DIVz_W * spdiagm(0 => view(sfcflx_factor_W, :)) * W_broadcast_sT,
            :invD_sT         => invD_sT,
            :f_sT            => f_sT,
            :β_sT            => β_sT,
            :ϵ_sT            => ϵ_sT,
            :D_sT            => D_sT,
            :Ks_H_T          => Ks_H_T,
            :Ks_H_U          => Ks_H_U,
            :Ks_H_V          => Ks_H_V,
        )

        tmpfi.check_usage[:Ks_H_U] .= Ks_H_U
        tmpfi.check_usage[:Ks_H_V] .= Ks_H_V

        vd = VerticalDiffusion(
            amo;
            K_iso=cfg["Ks_V"],
            K_cva=(cfg["convective_adjustment"] == "on") ? cfg["Ks_V_cva"] : cfg["Ks_V"],
        )

        # surface mask but is in 3D T grid. This one is different from the topo.sfcmask_sT.
        sfcmask_T = zeros(Float64, amo.bmo.T_dim...)
        sfcmask_T[1, :, :] .= 1
        mtx[:T_sfcmask_T] = spdiagm(0 => amo.T_mask_T * reshape(sfcmask_T, :))

        if length(ev.cdata_varnames) == 0
            cdatam = nothing
        else

            if cfg["cdata_var_file_map"] == nothing
                throw(ErrorException("Some config require cyclic data forcing file: $(ev.cdata_varnames)"))
            else
                function parseDateTime(timetype, str)
                    m = match(r"(?<year>[0-9]+)-(?<month>[0-9]{2})-(?<day>[0-9]{2})\s+(?<hour>[0-9]{2}):(?<min>[0-9]{2}):(?<sec>[0-9]{2})", str)
                    if m == nothing
                        throw(ErrorException("Unknown time format: " * (str)))
                    end

                    return timetype(
                        parse(Int64, m[:year]),
                        parse(Int64, m[:month]),
                        parse(Int64, m[:day]),
                        parse(Int64, m[:hour]),
                        parse(Int64, m[:min]),
                        parse(Int64, m[:sec]),
                    )
                end

                timetype   = getproperty(CFTime, Symbol(cfg["timetype"]))

                var_file_map = Dict()
                for varname in ev.cdata_varnames
                    var_file_map[varname] = cfg["cdata_var_file_map"][varname]
                end

                cdatam = CyclicDataManager(;
                    timetype     = timetype,
                    var_file_map = var_file_map,
                    beg_time     = parseDateTime(timetype, cfg["cdata_beg_time"]),
                    end_time     = parseDateTime(timetype, cfg["cdata_end_time"]),
                    align_time   = parseDateTime(timetype, cfg["cdata_align_time"]),
                    sub_yrng     = ev.sub_yrng,
                )

                tmpfi.datastream = makeDataContainer(cdatam)
            end
        end

        if cfg["weak_restoring"] == "on"
            mtx[:T_invτwk_TEMP_T] = - amo.T_mask_T * spdiagm(0 => ones(Float64, amo.bmo.T_pts)) / cfg["τwk_TEMP"]
            mtx[:T_invτwk_SALT_T] = - amo.T_mask_T * spdiagm(0 => ones(Float64, amo.bmo.T_pts)) / cfg["τwk_SALT"]
        end

        return new(
            wksp,

            #mask_sT,
            #ev.topo.mask_T,

            amo_slab,
            amo,

            mtx,

            vd,

            cdatam,

        )
    end

end


text:
using SymEngine

function sabs(a)
    return abs.(a)
end

"""
    sum(i::Symbol,s::Integer ,e::Integer, expr::Expr)

Expands the expression `expr` by replacing the symol `i`
with the Integers between `s` and `e`
"""
macro sum(i::Symbol,s::Integer ,e::Integer, expr::Expr)
    expr = eval_macro(expr, Symbol("@sum"))
    q = []
    for j = s:e
        push!(q,replace_symbol(i, j, expr))
    end
    r = Expr(:call)
    r.args = [:+; q]
    return Meta.quot(r)
end

"""
    function replace_symbol(i::Symbol, v::Integer, expr::Expr)

Replaces the symbol with the number v
"""
function replace_symbol(i::Symbol, v::Integer, expr::Expr)
    temp = copy(expr)
    replace_symbol!(i, v, temp)
    return temp
end

"""
    function replace_symbol!(i::Symbol, v::Integer, expr::Expr)

Used by [`replace_symbol`](@ref).
"""
function replace_symbol!(i::Symbol, v::Integer, expr::Expr)
    #println(typeof(expr))
    for (j,a) in enumerate(expr.args)
        if a == i
            expr.args[j] = v
        elseif a isa Expr
            replace_symbol!(i,v,a)
        end
    end
    return  expr
end
"""
    function eval_macro(expr::Expr, mc::Vararg{Symbol})

Recursively searches the given macro(s) in a quoted block
and replaces it with the result of its evaluation.

!!! Use this, multiple times for each macro if the order of their evaluation
matters.
"""
function eval_macro(expr::Expr, mc::Vararg{Symbol})
    temp = copy(expr)
    eval_macro!(temp, mc...)
    return temp
end

"""
    function eval_macro!(expr::Any, mc::Vararg{Symbol})

Inplace function used by [`eval_macro`](@ref)
"""
function eval_macro!(expr::Any, mc::Vararg{Symbol})
    if !(expr isa Expr)
        return
    end

    if expr.head == :macrocall
        if findfirst(expr.args[1] .== mc) != nothing
            x = eval(expr)
            expr.args = x.args
            expr.head = x.head
            return
        end
    end
    if expr isa Expr
        for (i, e) = enumerate(expr.args)
            eval_macro!(e, mc...)
        end
    end
    return expr
end

"""
    function eval_func!(expr::Any, mc::Vararg{Symbol})

Recursively searches the given function(s) in a quoted block
and replaces it inplace with the result of its evaluation.

!!! Use this, multiple times for each function if the order of their evaluation
matters.
"""
function eval_func!(expr::Any, mc::Vararg{Symbol})
    if !(expr isa Expr)
        return expr
    end

    if expr.head == :call
        if findfirst(expr.args[1] .== mc) != nothing
            x = eval(expr)
            return x
        end
    end
    if reduce( ==, expr.head .== mc)
        x = eval(expr)
        return x
    end
    if expr isa Expr
        for (i, e) = enumerate(expr.args)
            res = eval_func!(e, mc...)
            if res isa Number
                expr.args[i] =res
            end
        end
    end
    return expr
end

"""
    function minmax_const(expr::Expr)

Recursively searches a quoted block, for the indexing expression (:ref)
and keeps track of the max/min values of indices offsets.
"""
function minmax_const(expr::Expr)
    minval = 100
    maxval = -100
    function rec_minmax(expr, mn, mx, inref)
        for e in expr.args
            if e isa Expr
                (mn1, mx1) = rec_minmax(e, mn, mx,e.head == :ref || inref)
                if mn1 < mn
                    mn = mn1
                end
                if mx1 > mx
                    mx = mx1
                end
            elseif e isa Number
                !inref && continue
                flag = false
                if expr.head == :(-)
                    flag = true
                elseif expr.head == :(call)
                    if expr.args[1] == :(-)
                        flag = true
                    end
                end
                sgn = flag   ? -1 : 1
                if e*sgn < mn
                    mn = e*sgn
                end
                if e*sgn > mx
                    mx = e*sgn
                end
            end
        end
        return (mn,mx)
    end
    return rec_minmax(expr, minval, maxval, false)
end

"""
    fix_indices(expr::Expr, offset::Integer)

Changes the indices of arrays (:ref) applying an offset
to transform them from zero centerd with positive and negative
indexing to one-based indexing.
"""
function fix_indices(expr::Expr, offset::Integer)
    ex = copy(expr)
    function fix!(expr, offset, inref)
        for (i,e) in enumerate(expr.args)
            if e isa Expr
                 #println("rec", e)
                 fix!(expr.args[i], offset, e.head == :ref || inref)
            elseif e isa Number
                !inref && continue
                #println(e, " ", i, " ", expr.head, " ## " ,expr.args)
                if ((expr.head == :(call) && expr.args[1] == :(-)) ||
                                                expr.head == :(-))
                    expr.args[i] = -e + offset
                    if (expr.head == :(-))
                        expr.args[1] = :(+)
                    end
                    break
                elseif ((expr.head == :(call) && expr.args[1] == :(+)) ||
                    expr.head == :(+))
                    expr.args[i] = e + offset
                    break
                end
                expr.args[i] = e + offset
            end
        end
    end
    fix!(ex,offset, false)
    return ex
    #expr.head = (:call)
end

# macro def_stencil(expr::Expr)#, data, current, t)
#     expr = Meta.quot(expr)
#     expr = eval_macro(expr, Symbol("@sum"))
#     # mn,mx = minmax_const(expr)
#     # println("MIN", mn, " MAX ", mx)
#     # offset = 0
#     # if mn <= 0
#     #     offset = -mn + 1
#     # end
#     # println(offset," ", mn)
#     #e = fix_indices(expr, offset)
#     return expr
# end
#
# macro def_stencil(expr::Expr, vars, radius)
#     expr = Meta.quot(expr)
#     expr = eval_macro(expr, Symbol("@sum"))
#     stencil = expr
#     for i in vars.args
#         println(i.value)
#         stencil = replace_symbol(i.value, 0, stencil)
#     end
#     stencil = fix_indices(stencil, radius)
#     return (s=stencil, v=expr)
# end

nothing


text:
cd(joinpath(dirname(@__FILE__), "src"))
include(joinpath(dirname(@__FILE__), "src", "ISP.jl"))

function doISP()
  include(joinpath(dirname(@__FILE__), "ex3.jl"))

  ed=ISP.ExData(a,b,c,s,u,ux,ut,sp,γf)
  sd=ISP.SolvData(ed,10)

  const tgrid = sd.grid.tgrid

  v_actual = convert(Vector{Float64}, [sd.grid.xbdy, [sd.g(t) for t in tgrid]])
  v0 = copy(v_actual)
  len = length(v0)
  v0[1:(len/2)] = mean(v0[1:(len/2)])
  v0[(len/2+1):len] = mean(v0[(len/2+1):len])

  mu_actual = Float64[u(s(t), t) for t in tgrid]
  nu_actual = Float64[u(zero(t), t) for t in tgrid]

  delta = 0.5
  l = 5.0
  R = 20.0

  func=ISP.SimpleTest(sd, mu_actual, nu_actual, v0, delta, l, R)

end
doISP()

max_stars_repo_name: jgoldfar/isp-julia-mwe

max_issues_repo_name: jgoldfar/isp-julia-mwe

max_forks_repo_name: jgoldfar/isp-julia-mwe


text:
using DelimitedFiles
using WriteVTK
using MeshCore: AbsShapeDesc, P1, L2, L3
using MeshCore: VecAttrib
using MeshCore: shapedesc, nshapes,  IncRel
using MeshCore: ShapeColl
using LinearAlgebra: norm
using StaticArrays

"""
    L2blockx(xs::Vector{T}; intbytes = 8) where {T}

Generate a graded mesh on an interval.

Keyword argument `intbytes` controls the size of the integer indexes.

# Return
By convention the function returns an incidence relation (`connectivity`)
between the elements (`connectivity.left`) and the vertices
(`connectivity.right`). The geometry is stored as the attribute "geom" of the
vertices.
"""
function L2blockx(xs::Vector{T}; intbytes = 8) where {T}
	inttype = Int64
    if intbytes == 4
        inttype = Int32
    end
    nL = length(xs)-1;
    nnodes = (nL+1);
    ncells = (nL);
    xys = zeros(T, nnodes, 1);
    conns = zeros(inttype, ncells, 2);

    f=1;
    for i in 1:(nL+1)
        xys[f, 1] = xs[i]
        f=f+1;
    end

    gc=1;
    for i in 1:nL
       conns[gc, 1] = i
       conns[gc, 2] = (i+1)
       gc=gc+1;
    end

    C = inttype.(conns[1:gc-1, :]);
    N, TC = size(xys, 2), eltype(xys)
    locs =  VecAttrib([SVector{N, TC}(xys[i, :]) for i in 1:size(xys, 1)])
    vertices = ShapeColl(P1, length(locs), "vertices")
    vertices.attributes["geom"] = locs
    elements = ShapeColl(L2, size(C, 1), "elements")
    return IncRel(elements, vertices, C)
end

"""
    L2block(Length, nL; kwargs...)

Generate a quadrilateral mesh  of the 2D block.

See also: L2blockx
"""
function L2block(Length, nL; kwargs...)
    return L2blockx(collect(linearspace(0.0, Length, nL+1)); kwargs...);
end


text:
# copied and tweaked from https://github.com/staticfloat/LLVMBuilder
# Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
using BinaryBuilder

# Collection of sources required to build LLVM
llvm_ver = "6.0.1"
sources = [
    "http://releases.llvm.org/$(llvm_ver)/llvm-$(llvm_ver).src.tar.xz" =>
    "b6d6c324f9c71494c0ccaf3dac1f16236d970002b42bb24a6c9e1634f7d0f4e2",
    "http://releases.llvm.org/$(llvm_ver)/cfe-$(llvm_ver).src.tar.xz" =>
    "7c243f1485bddfdfedada3cd402ff4792ea82362ff91fbdac2dae67c6026b667",
    "http://releases.llvm.org/$(llvm_ver)/compiler-rt-$(llvm_ver).src.tar.xz" =>
    "f4cd1e15e7d5cb708f9931d4844524e4904867240c306b06a4287b22ac1c99b9",
    "http://releases.llvm.org/$(llvm_ver)/libcxx-$(llvm_ver).src.tar.xz" =>
    "7654fbc810a03860e6f01a54c2297a0b9efb04c0b9aa0409251d9bdb3726fc67",
    "http://releases.llvm.org/$(llvm_ver)/libcxxabi-$(llvm_ver).src.tar.xz" =>
    "209f2ec244a8945c891f722e9eda7c54a5a7048401abd62c62199f3064db385f",
    "http://releases.llvm.org/$(llvm_ver)/polly-$(llvm_ver).src.tar.xz" =>
    "e7765fdf6c8c102b9996dbb46e8b3abc41396032ae2315550610cf5a1ecf4ecc",
    "http://releases.llvm.org/$(llvm_ver)/libunwind-$(llvm_ver).src.tar.xz" =>
    "a8186c76a16298a0b7b051004d0162032b9b111b857fbd939d71b0930fd91b96",
    "http://releases.llvm.org/$(llvm_ver)/lld-$(llvm_ver).src.tar.xz" =>
    "e706745806921cea5c45700e13ebe16d834b5e3c0b7ad83bf6da1f28b0634e11",
    "patches",
	# julia source
    "https://github.com/JuliaLang/julia/releases/download/v1.3.0/julia-1.3.0.tar.gz" =>
	"5e1bfa62b99921ec8f437a31e834802d8d357ab679611ed336bd6e214ea0ce75",
	# dlfcn-win32 source
    "https://github.com/dlfcn-win32/dlfcn-win32/archive/v1.2.0.tar.gz" =>
	"f18a412e84d8b701e61a78252411fe8c72587f52417c1ef21ca93604de1b9c55",
    # libcxxffi source
	"libcxxffi",
	"https://github.com/Gnimuc/JuliaBuilder/releases/download/v1.3.0/julia-1.3.0-x86_64-linux-gnu.tar.gz" =>
	"44099e27a3d9ebdaf9d67bfdaf745c3899654c24877c76cbeff9cade5ed79139",
    # LLVM binary
	"https://github.com/staticfloat/LLVMBuilder/releases/download/v6.0.1-7%2Bnowasm/LLVM.v6.0.1.x86_64-linux-gnu-gcc7.tar.gz" =>
	"f2c335eb912720a5b3318f909eda1650973043beb033dfb3adc0f7d150b91bf6",
]

# Since we kind of do this LLVM setup twice, this is the shared setup start:
script_setup = raw"""
# We want to exit the program if errors occur.
set -o errexit

cd $WORKSPACE/srcdir/

# First, move our other projects into llvm/projects
for f in *.src; do
    # Don't symlink llvm itself into llvm/projects...
    if [[ ${f} == llvm-*.src ]]; then
        continue
    fi

    # clang lives in tools/clang and not projects/cfe
    if [[ ${f} == cfe-*.src ]]; then
        mv $(pwd)/${f} $(echo llvm-*.src)/tools/clang
    elif [[ ${f} == polly-*.src ]]; then
        mv $(pwd)/${f} $(echo llvm-*.src)/tools/polly
    elif [[ ${f} == lld-*.src ]]; then
        mv $(pwd)/${f} $(echo llvm-*.src)/tools/lld
    else
        mv $(pwd)/${f} $(echo llvm-*.src)/projects/${f%-*}
    fi
done

# Next, boogie on down to llvm town
cd llvm-*.src

# Update config.guess/config.sub stuff
# update_configure_scripts

# Apply all our patches
for f in $WORKSPACE/srcdir/llvm_patches/*.patch; do
    echo "Applying patch ${f}"
    patch -p1 < ${f}
done
"""

# Next, build libcxxffi
script = script_setup * raw"""
# build libcxxffi
cd $WORKSPACE/srcdir/

if [[ ${target} == *mingw32* ]]; then
    cd $WORKSPACE/srcdir/dlfcn-win32-1.2.0
    mkdir build && cd build
    cmake .. -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=$CMAKE_TARGET_TOOLCHAIN
    make -j${nproc} VERBOSE=1
    make install VERBOSE=1
fi

# use CMake on all platforms except Windows
if [[ "${target}" == *linux* ]] || [[ "${target}" == *apple* ]] || [[ "${target}" == *freebsd* ]]; then
	mkdir build && cd build
	CMAKE_FLAGS="-DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=$CMAKE_TARGET_TOOLCHAIN"
	CMAKE_FLAGS="${CMAKE_FLAGS} -DJULIA_SOURCE_PREFIX=$WORKSPACE/srcdir/julia-1.3.0"
	CMAKE_FLAGS="${CMAKE_FLAGS} -DJULIA_BINARY_PREFIX=$WORKSPACE/srcdir/juliabin"
	CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVMBUILDER_PREFIX=$WORKSPACE/srcdir"
	cmake .. ${CMAKE_FLAGS}
	make -j${nproc} VERBOSE=1
	make install VERBOSE=1
fi

# use Make on Windows
if [[ ${target} == *mingw32* ]]; then
    cd $WORKSPACE/srcdir/
    make -j${nproc} VERBOSE=1 -f Windows.Makefile JULIA_SOURCE_PREFIX=$WORKSPACE/srcdir/julia-1.3.0 JULIA_BINARY_PREFIX=$WORKSPACE/srcdir/juliabin LLVMBUILDER_PREFIX=$WORKSPACE/srcdir
    cp -r ./build/* /workspace/destdir/bin/
fi

# copy LLVM/Clang assets
mkdir -p ${prefix}/src/clang-6.0.1/include
mkdir -p ${prefix}/src/llvm-6.0.1/include
mkdir -p ${prefix}/build/clang-6.0.1/lib/clang/6.0.1/include
mkdir -p ${prefix}/build/clang-6.0.1/include
mkdir -p ${prefix}/build/llvm-6.0.1/include
cp -r $WORKSPACE/srcdir/llvm-6.0.1.src/tools/clang/include/* ${prefix}/src/clang-6.0.1/include/
cp -r $WORKSPACE/srcdir/llvm-6.0.1.src/include/* ${prefix}/src/llvm-6.0.1/include/
cp -r $WORKSPACE/srcdir/include/* ${prefix}/build/llvm-6.0.1/include/
cp -r $WORKSPACE/srcdir/lib/* ${prefix}/lib/
cp -r $WORKSPACE/srcdir/lib/clang/6.0.1/include/* ${prefix}/build/clang-6.0.1/lib/clang/6.0.1/include/
cp -r ${prefix}/build/llvm-6.0.1/include/* ${prefix}/build/clang-6.0.1/include/

# generate clang_constants.jl
cd $WORKSPACE/srcdir
make -f GenerateConstants.Makefile BASE_LLVM_BIN=$WORKSPACE/srcdir BASE_JULIA_BIN=$WORKSPACE/srcdir/juliabin BASE_JULIA_SRC=$WORKSPACE/srcdir/julia LLVM_VERSION=6.0.1
cp $WORKSPACE/srcdir/clang_constants.jl ${prefix}/

# copy mingw headers
if [[ ${target} == *mingw32* ]] && [[ ${nbits} == 64 ]]; then
    cp -r $WORKSPACE/srcdir/bin/* ${prefix}/bin
    mkdir -p ${prefix}/mingw/include
    mkdir -p ${prefix}/mingw/sys-root/include
    cp -r /opt/x86_64-w64-mingw32/x86_64-w64-mingw32/include/* ${prefix}/mingw/include/
    cp -r /opt/x86_64-w64-mingw32/x86_64-w64-mingw32/sys-root/include/* ${prefix}/mingw/sys-root/include/
fi

if [[ ${target} == *mingw32* ]] && [[ ${nbits} == 32 ]]; then
    cp -r $WORKSPACE/srcdir/bin/* ${prefix}/bin
    mkdir -p ${prefix}/mingw/include
    mkdir -p ${prefix}/mingw/sys-root/include
    cp -r /opt/i686-w64-mingw32/i686-w64-mingw32/include/* ${prefix}/mingw/include/
    cp -r /opt/i686-w64-mingw32/i686-w64-mingw32/sys-root/include/* ${prefix}/mingw/sys-root/include/
fi

"""

platforms = supported_platforms()

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

build_tarballs(ARGS, "libcxxffi", v"0.0.0", sources, script, platforms, products, [], preferred_gcc_version=v"7")

max_stars_repo_name: Gnimuc/CxxBuilder

max_issues_repo_name: Gnimuc/CxxBuilder

max_forks_repo_name: Gnimuc/CxxBuilder


text:
"""
Get appropriate part of a body content by name

Returns found part or nothing

# Arguments
- `multipart::Vector{HTTP.Multipart}`: HTTP.jl provided vector of parts
- `name::String`: the name of part parameter
"""
get(multipart::Vector{HTTP.Multipart}, name::String) =
    findfirst(part -> part.name == name, multipart) |>
    idx -> isnothing(idx) ? nothing : multipart[idx]

"""
Process incomming query with image file
Generate output image in PNG format
"""
function process_image(c::WelcomeController)
    local multiparts = HTTP.parse_multipart_form(c.conn.request)

    # dump(multiparts)
    local mp_file = get(multiparts, "f_image")
    local mp_num = get(multiparts, "num")

    isnothing(mp_file) &&
        return render(Text, "Please upload some image file in f_image field")
    isnothing(mp_num) &&
        return render(Text, "Please specify number of sections in num field")

    local img = FileIO.load(Stream{format"PNG"}(mp_file.data))
    local out_img = JWebImageDemo.permute_image(
        img,
        tryparse(Int, String(take!(mp_num.data))) |> x -> something(x, 5),
    )
    local out_buf = IOBuffer()
    show(out_buf, MIME("image/png"), out_img)
    # return Render(mp_file.contenttype, identity, take!(out_buf))
    return Render("image/png", identity, take!(out_buf))
end


text:
using LoopFieldCalc
using Test

function test_many_loops()
	xs = 0.0:0.0
	ys = -0.4:0.005:0.4
	zs = -0.3:0.005:0.3

	nturns = 10

	inner_bobbin_radius = 0.1
	inner_bobbin_width = 0.15
	inner_bobbin_offset = -0.1
	inner_bobbin_current = 1.0

	outer_bobbin_radius = 0.25
	outer_bobbin_width = 0.15
	outer_bobbin_offset = -0.15
	outer_bobbin_current = 0.3 * inner_bobbin_current

	inner_bobbin = [
		LoopFieldCalc.CurrentLoop(
			radius = inner_bobbin_radius,
			current = inner_bobbin_current / nturns,
			center = LoopFieldCalc.CylindricalPoint(inner_bobbin_offset + z, 0.0, 0.0)
		)
		for z in LinRange(0.0, inner_bobbin_width, nturns)
	]

	outer_bobbin = [
		LoopFieldCalc.CurrentLoop(
			radius = outer_bobbin_radius,
			current = outer_bobbin_current / nturns,
			center = LoopFieldCalc.CylindricalPoint(outer_bobbin_offset + z, 0.0, 0.0)
		)
		for z in LinRange(0.0, outer_bobbin_width, nturns)
	]

	loops = [outer_bobbin; inner_bobbin]
	Bx, By, Bz = LoopFieldCalc.field_on_box(loops, xs, ys, zs)
	LoopFieldCalc.write_field("many_loops.dat", zs, ys, Bz[1, :, :]', By[1, :, :]')
end
test_many_loops()


text:
function cumulative_probability(obj::RealDistribution, arg0::jdouble)
    return jcall(obj, "cumulativeProbability", jdouble, (jdouble,), arg0)
end

function density(obj::RealDistribution, arg0::jdouble)
    return jcall(obj, "density", jdouble, (jdouble,), arg0)
end

function get_numerical_mean(obj::RealDistribution)
    return jcall(obj, "getNumericalMean", jdouble, ())
end

function get_numerical_variance(obj::RealDistribution)
    return jcall(obj, "getNumericalVariance", jdouble, ())
end

function get_support_lower_bound(obj::RealDistribution)
    return jcall(obj, "getSupportLowerBound", jdouble, ())
end

function get_support_upper_bound(obj::RealDistribution)
    return jcall(obj, "getSupportUpperBound", jdouble, ())
end

function inverse_cumulative_probability(obj::RealDistribution, arg0::jdouble)
    return jcall(obj, "inverseCumulativeProbability", jdouble, (jdouble,), arg0)
end

function is_support_connected(obj::RealDistribution)
    return jcall(obj, "isSupportConnected", jboolean, ())
end

function log_density(obj::RealDistribution, arg0::jdouble)
    return jcall(obj, "logDensity", jdouble, (jdouble,), arg0)
end

function probability(obj::RealDistribution, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "probability", jdouble, (jdouble, jdouble), arg0, arg1)
end


text:
using QuasiGeostrophy, BenchmarkTools

Ω = S¹(0, 2π) × S¹(0, 2π)
Nx = Ny = 2^10;
fourier_grid = FourierGrid((Nx, Ny), Ω)
fieldnames = ("u", "v", "ψ")
u, v, ψ = create_fields(names = fieldnames, grid = fourier_grid)
∂x, ∂y = create_operators(fourier_grid)
# initialize FourierField data
x, y = fourier_grid.grid
ψ(sin.(x) .* sin.(y))
u.data .= ∂y(ψ).data
v.data .= -∂x(ψ).data
Δ = FourierOperator((∂x * ∂x + ∂y * ∂y).op, FourierOperatorMetaData("Δ"))
Δ⁻¹ = inv(Δ)
Re = 0.1

function nsrhs(u,v, ∂x, ∂y, Δ, Δ⁻¹, Re)
    rhsu = -(∂x(u*u) + ∂y(u*v)) + (1/Re)*Δ(u)
    rhsv = -(∂x(u*v) + ∂y(v*v)) + (1/Re)*Δ(v)
    p = Δ⁻¹(∂x(rhsu) + ∂y(rhsv))
    rhsu -= ∂x(p)
    rhsv -= ∂y(p)
    return rhsu, rhsv
end

function nsrhs2(u,v, ∂x, ∂y, Δ, Δ⁻¹, Re)
    uu = u*u
    uv = u*v
    vu = v*u
    vv = v*v
    uux = ∂x(uu)
    uvy = ∂y(uv)
    vux = ∂x(vu)
    vvy = ∂y(vv)
    Δu = (1/Re)*Δ(u)
    Δv = (1/Re)*Δ(v)
    rhsu = -(uux + uvy) + Δu
    rhsv = -(vux + vvy) + Δv
    p = Δ⁻¹(∂x(rhsu) + ∂y(rhsv))
    rhsu -= ∂x(p)
    rhsv -= ∂y(p)
    return rhsu, rhsv
end

P = plan_fft(u.data, flags=FFTW.PATIENT)
iP = plan_ifft(u.data, flags=FFTW.PATIENT)
function nsrhs3(u,v, ∂x, ∂y, Δ, Δ⁻¹, Re, P, iP)
    ur = iP * u.data
    vr = iP * v.data
    fmd = FourierMetaData(" ", u.metadata.grid, u.metadata.transform)
    uu = FourierField( P * (ur .* ur), fmd)
    uv = FourierField( P * (ur .* vr),fmd)
    vv = FourierField( P * (vr .* vr),fmd)
    uux = ∂x(uu)
    uvy = ∂y(uv)
    vux = ∂x(uv)
    vvy = ∂y(vv)
    Δu = (1/Re)*Δ(u)
    Δv = (1/Re)*Δ(v)
    rhsu = -(uux + uvy) + Δu
    rhsv = -(vux + vvy) + Δv
    p = Δ⁻¹(∂x(rhsu) + ∂y(rhsv))
    rhsu -= ∂x(p)
    rhsv -= ∂y(p)
    return rhsu, rhsv
end

#nonfft operations
function srhs(u,v, ∂x, ∂y, Δ, Δ⁻¹, Re)
    uu = u
    Δu = (1/Re)*Δ(u)
    Δv = (1/Re)*Δ(v)
    rhsu = -(∂x(uu) + ∂y(uu)) + Δu
    rhsv = -(∂x(uu) + ∂y(uu)) + Δv
    p = Δ⁻¹(∂x(rhsu) + ∂y(rhsv))
    rhsu -= ∂x(p)
    rhsv -= ∂y(p)
    return rhsu, rhsv
end

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: sandreza/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
function FieldEquinoctialOrbit(arg0::FieldOrbit)
    return FieldEquinoctialOrbit((FieldOrbit,), arg0)
end

function FieldEquinoctialOrbit(arg0::FieldPVCoordinates, arg1::Frame, arg2::FieldAbsoluteDate, arg3::RealFieldElement)
    return FieldEquinoctialOrbit((FieldPVCoordinates, Frame, FieldAbsoluteDate, RealFieldElement), arg0, arg1, arg2, arg3)
end

function FieldEquinoctialOrbit(arg0::RealFieldElement, arg1::RealFieldElement, arg10::RealFieldElement, arg11::RealFieldElement, arg12::PositionAngle, arg13::Frame, arg14::FieldAbsoluteDate, arg15::RealFieldElement, arg2::RealFieldElement, arg3::RealFieldElement, arg4::RealFieldElement, arg5::RealFieldElement, arg6::RealFieldElement, arg7::RealFieldElement, arg8::RealFieldElement, arg9::RealFieldElement)
    return FieldEquinoctialOrbit((RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, PositionAngle, Frame, FieldAbsoluteDate, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement), arg0, arg1, arg10, arg11, arg12, arg13, arg14, arg15, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function FieldEquinoctialOrbit(arg0::RealFieldElement, arg1::RealFieldElement, arg2::RealFieldElement, arg3::RealFieldElement, arg4::RealFieldElement, arg5::RealFieldElement, arg6::PositionAngle, arg7::Frame, arg8::FieldAbsoluteDate, arg9::RealFieldElement)
    return FieldEquinoctialOrbit((RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, PositionAngle, Frame, FieldAbsoluteDate, RealFieldElement), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function FieldEquinoctialOrbit(arg0::TimeStampedFieldPVCoordinates, arg1::Frame, arg2::RealFieldElement)
    return FieldEquinoctialOrbit((TimeStampedFieldPVCoordinates, Frame, RealFieldElement), arg0, arg1, arg2)
end

function add_kepler_contribution(obj::FieldEquinoctialOrbit, arg0::PositionAngle, arg1::RealFieldElement, arg2::Vector{RealFieldElement})
    return jcall(obj, "addKeplerContribution", void, (PositionAngle, RealFieldElement, Vector{RealFieldElement}), arg0, arg1, arg2)
end

function eccentric_to_mean(::Type{FieldEquinoctialOrbit}, arg0::RealFieldElement, arg1::RealFieldElement, arg2::RealFieldElement)
    return jcall(FieldEquinoctialOrbit, "eccentricToMean", RealFieldElement, (RealFieldElement, RealFieldElement, RealFieldElement), arg0, arg1, arg2)
end

function eccentric_to_true(::Type{FieldEquinoctialOrbit}, arg0::RealFieldElement, arg1::RealFieldElement, arg2::RealFieldElement)
    return jcall(FieldEquinoctialOrbit, "eccentricToTrue", RealFieldElement, (RealFieldElement, RealFieldElement, RealFieldElement), arg0, arg1, arg2)
end

function get_a(obj::FieldEquinoctialOrbit)
    return jcall(obj, "getA", RealFieldElement, ())
end

function get_a_dot(obj::FieldEquinoctialOrbit)
    return jcall(obj, "getADot", RealFieldElement, ())
end

function get_e(obj::FieldEquinoctialOrbit)
    return jcall(obj, "getE", RealFieldElement, ())
end

function get_e_dot(obj::FieldEquinoctialOrbit)
    return jcall(obj, "getEDot", RealFieldElement, ())
end

function get_equinoctial_ex(obj::FieldEquinoctialOrbit)
    return jcall(obj, "getEquinoctialEx", RealFieldElement, ())
end

function get_equinoctial_ex_dot(obj::FieldEquinoctialOrbit)
    return jcall(obj, "getEquinoctialExDot", RealFieldElement, ())
end

function get_equinoctial_ey(obj::FieldEquinoctialOrbit)
    return jcall(obj, "getEquinoctialEy", RealFieldElement, ())
end

function get_equinoctial_ey_dot(obj::FieldEquinoctialOrbit)
    return jcall(obj, "getEquinoctialEyDot", RealFieldElement, ())
end

function get_hx(obj::FieldEquinoctialOrbit)
    return jcall(obj, "getHx", RealFieldElement, ())
end

function get_hx_dot(obj::FieldEquinoctialOrbit)
    return jcall(obj, "getHxDot", RealFieldElement, ())
end

function get_hy(obj::FieldEquinoctialOrbit)
    return jcall(obj, "getHy", RealFieldElement, ())
end

function get_hy_dot(obj::FieldEquinoctialOrbit)
    return jcall(obj, "getHyDot", RealFieldElement, ())
end

function get_i(obj::FieldEquinoctialOrbit)
    return jcall(obj, "getI", RealFieldElement, ())
end

function get_i_dot(obj::FieldEquinoctialOrbit)
    return jcall(obj, "getIDot", RealFieldElement, ())
end

function get_l(obj::FieldEquinoctialOrbit, arg0::PositionAngle)
    return jcall(obj, "getL", RealFieldElement, (PositionAngle,), arg0)
end

function get_l_dot(obj::FieldEquinoctialOrbit, arg0::PositionAngle)
    return jcall(obj, "getLDot", RealFieldElement, (PositionAngle,), arg0)
end

function get_le(obj::FieldEquinoctialOrbit)
    return jcall(obj, "getLE", RealFieldElement, ())
end

function get_le_dot(obj::FieldEquinoctialOrbit)
    return jcall(obj, "getLEDot", RealFieldElement, ())
end

function get_lm(obj::FieldEquinoctialOrbit)
    return jcall(obj, "getLM", RealFieldElement, ())
end

function get_lm_dot(obj::FieldEquinoctialOrbit)
    return jcall(obj, "getLMDot", RealFieldElement, ())
end

function get_lv(obj::FieldEquinoctialOrbit)
    return jcall(obj, "getLv", RealFieldElement, ())
end

function get_lv_dot(obj::FieldEquinoctialOrbit)
    return jcall(obj, "getLvDot", RealFieldElement, ())
end

function get_type(obj::FieldEquinoctialOrbit)
    return jcall(obj, "getType", OrbitType, ())
end

function has_derivatives(obj::FieldEquinoctialOrbit)
    return jcall(obj, "hasDerivatives", jboolean, ())
end

function interpolate(obj::FieldEquinoctialOrbit, arg0::FieldAbsoluteDate, arg1::Stream)
    return jcall(obj, "interpolate", FieldEquinoctialOrbit, (FieldAbsoluteDate, Stream), arg0, arg1)
end

function mean_to_eccentric(::Type{FieldEquinoctialOrbit}, arg0::RealFieldElement, arg1::RealFieldElement, arg2::RealFieldElement)
    return jcall(FieldEquinoctialOrbit, "meanToEccentric", RealFieldElement, (RealFieldElement, RealFieldElement, RealFieldElement), arg0, arg1, arg2)
end

function shifted_by(obj::FieldEquinoctialOrbit, arg0::RealFieldElement)
    return jcall(obj, "shiftedBy", FieldEquinoctialOrbit, (RealFieldElement,), arg0)
end

function shifted_by(obj::FieldEquinoctialOrbit, arg0::jdouble)
    return jcall(obj, "shiftedBy", FieldEquinoctialOrbit, (jdouble,), arg0)
end

function to_orbit(obj::FieldEquinoctialOrbit)
    return jcall(obj, "toOrbit", EquinoctialOrbit, ())
end

function to_string(obj::FieldEquinoctialOrbit)
    return jcall(obj, "toString", JString, ())
end

function true_to_eccentric(::Type{FieldEquinoctialOrbit}, arg0::RealFieldElement, arg1::RealFieldElement, arg2::RealFieldElement)
    return jcall(FieldEquinoctialOrbit, "trueToEccentric", RealFieldElement, (RealFieldElement, RealFieldElement, RealFieldElement), arg0, arg1, arg2)
end