text:
"""
    read_gmsh3d_mesh(filename) -> mesh::Mesh
"""
function read_gmsh3d_mesh(fn::AbstractString; kwargs...)
    open(fn, "r") do io
        read_gmsh3d_mesh(io; kwargs...)
    end
end

"""
    read_gmsh3d_mesh(iostream) -> mesh::Mesh

Reads the mesh nodes and elements stored in the input .msh file (`io`, output by GMSH)
    into arrays of node vectors and vertex vectors respectively.

Returns an object `mesh` of type `Mesh`, comprising both vector arrays.
"""
function read_gmsh3d_mesh(io; physical=nothing)

    entity_tag = 0
    if physical != nothing
        @assert physical isa String
        while true
            thisLine = io |> readline |> strip
            thisLine == "\$PhysicalNames" && break
            thisLine == "\$Nodes" && error("Mesh file does not define physical entities")
        end
        thisLine = io |> readline |> strip
        s = split(thisLine)
        NP = parse(Int,s[1])
        println("Mesh file defined $(NP) physical entities:")
        for i in 1:NP
            thisLine = io |> readline |> strip
            s = split(thisLine)
            entity_name = s[3]
            @show entity_name
            if entity_name == "\"$(physical)\""
                entity_tag = parse(Int, s[2])
                println("Target entity \"$(physical)\" is $i out of $(NP).")
                break
            end
        end
        entity_tag == 0 && error("Specified physical entitiy not found")
    end

    thisLine = io |> readline |> strip
    while thisLine != "\$Nodes"
        thisLine = io |> readline |> strip
    end

    thisLine = io |> readline |> strip
    s = split(thisLine)
    NV = parse(Int, s[1])

    P = SVector{3,Float64}
    v = Vector{P}(undef,NV)
    for i in 1:NV
        thisLine = io |> readline |>  strip
        d = readdlm(IOBuffer(thisLine), Float64)
        v[i] = P(d[2], d[3], d[4])
    end

    while thisLine != "\$Elements"
        thisLine = io |> readline |> strip
    end

    thisLine = io |> readline |> strip
    s = split(thisLine)
    NF = parse(Int, s[1])

    thisLine = io |> readline |> strip
    f = Vector{SVector{4,Int}}(undef,NF); i = 0
    while thisLine != "\$EndElements"
        d = readdlm(IOBuffer(thisLine), Int)
        thisLine = io |> readline |> strip
        d[2] == 4 || continue
        tag = d[4]
        if entity_tag == 0 || entity_tag == tag
            f[i +=1] = SVector(d[end-3], d[end-2], d[end-1], d[end-0])
        end
    end
    resize!(f,i)

    # gmsh meshes are not always consistently oriented
    for (i,fc) in enumerate(f)
        vs = v[fc]
        t1 = vs[1] - vs[4]
        t2 = vs[2] - vs[4]
        t3 = vs[3] - vs[4]
        orientation = dot(t1 × t2, t3)
        if orientation < 0
            f[i] = @SVector[fc[1],fc[2],fc[4],fc[3]]
        end
    end
    msh = Mesh(v, f)
    @assert isoriented(msh)
    return msh
end


text:
"""
    EvolvingGraph{V,T}(;is_directed=true; is_weighted=true)
    EvolvingGraph(;is_directed=true; is_weighted=true)

Construct an evolving graph with node type `V` and timestamp type `T`.
`EvolvingGraph()` constructs a simple evolving graph with integer nodes and timestamps.

# Example

```jldoctest
julia> using EvolvingGraphs

julia> g = EvolvingGraph{Node{String},Int}(is_weighted=false)
Directed EvolvingGraph 0 nodes, 0 static edges, 0 timestamps

julia> add_node!(g, "a")
Node(a)

julia> add_node!(g, "b")
Node(b)

julia> num_nodes(g)
2

julia> add_edge!(g, "a", "b", 2001)
Node(a)->Node(b) at time 2001

julia> add_edge!(g, "a", "c", 2002)
Node(a)->Node(c) at time 2002

julia> timestamps(g)
2-element Array{Int64,1}:
 2001
 2002

julia> active_nodes(g)
4-element Array{EvolvingGraphs.TimeNode{String,Int64},1}:
 TimeNode(a, 2001)
 TimeNode(b, 2001)
 TimeNode(a, 2002)
 TimeNode(c, 2002)

julia> g = EvolvingGraph()
Directed EvolvingGraph 0 nodes, 0 static edges, 0 timestamps

julia> add_edge!(g, 1, 2, 1)
Node(1)-1.0->Node(2) at time 1

julia> add_edge!(g, 2, 3, 2)
Node(2)-1.0->Node(3) at time 2

julia> add_edge!(g, 1, 3, 2)
Node(1)-1.0->Node(3) at time 2

julia> nodes(g)
3-element Array{EvolvingGraphs.Node{Int64},1}:
 Node(1)
 Node(2)
 Node(3)
```
"""
mutable struct EvolvingGraph{V, E, T, KV} <: AbstractEvolvingGraph{V, E, T}
    is_directed::Bool
    nodes::Vector{V}             # a vector of nodes
    node_indexof::Dict{KV, Int}   # map node keys to indices
    edges::Vector{E}           # a vector of edges
    timestamps::Vector{T}           # a vector of timestamps
    active_nodes::Vector{TimeNode{KV,T}}  # a vector of active nodes
    active_node_indexof::Dict{Tuple{KV,T},Int} # mapping active node keys to indices
end

function EvolvingGraph{V,T}(;is_directed::Bool=true, is_weighted::Bool=true) where {V,T}
    KV = eltype(V);
    E = is_weighted ? WeightedTimeEdge{V,T,Float64} : TimeEdge{V,T};
    return EvolvingGraph(is_directed, V[], Dict{KV, Int}(), E[], T[], TimeNode{KV, T}[], Dict{Tuple{KV,T},Int}())
end
EvolvingGraph(;is_directed::Bool=true, is_weighted::Bool=true) = EvolvingGraph{Node{Int},Int}(is_directed=is_directed, is_weighted=is_weighted)

"""
    evolving_graph_from_arrays(ils, jls, wls, timestamps; is_directed=true)
    evolving_graph_from_arrays(ils, jls, timestamps; is_directed=true)

Generate an EvolvingGraph type object from four input arrays: `ils`, `jls`, `wls` and `timestamps`, such that the ith entry `(ils[i], jls[i], wls[i], timestamps[i])` represents a WeightedTimeEdge from `ils[i]` to `jls[i]` with edge weight `wls[i]` at timestamp `timestamp[i]`. By default, `wls` is a vector of ones.

#  Example

```jldoctest
julia> using EvolvingGraphs

julia> g = evolving_graph_from_arrays([1,2,3], [4,5,2], [1,1,2])
Directed EvolvingGraph 5 nodes, 3 static edges, 2 timestamps

julia> nodes(g)
5-element Array{EvolvingGraphs.Node{Int64},1}:
 Node(1)
 Node(4)
 Node(2)
 Node(5)
 Node(3)

julia> edges(g)
3-element Array{EvolvingGraphs.WeightedTimeEdge{EvolvingGraphs.Node{Int64},Int64,Float64},1}:
 Node(1)-1.0->Node(4) at time 1
 Node(2)-1.0->Node(5) at time 1
 Node(3)-1.0->Node(2) at time 2

julia> g = evolving_graph_from_arrays([1,2], [2, 3], [2.5, 3.8], [1998,2001])
Directed EvolvingGraph 3 nodes, 2 static edges, 2 timestamps

julia> edges(g)
2-element Array{EvolvingGraphs.WeightedTimeEdge{EvolvingGraphs.Node{Int64},Int64,Float64},1}:
 Node(1)-2.5->Node(2) at time 1998
 Node(2)-3.8->Node(3) at time 2001
```
"""
function evolving_graph_from_arrays{V,T}(ils::Vector{V},
                            jls::Vector{V}, wls::Vector{<:Real},
                            timestamps::Vector{T};
                            is_directed::Bool = true)
    n = length(ils)
    n == length(jls) == length(timestamps) == length(wls)||
            error("4 input vectors must have the same length.")

    g = EvolvingGraph{Node{eltype(ils)}, eltype(timestamps)}(is_directed = is_directed)

    for i = 1:n
        v1 = add_node!(g, ils[i])
        v2 = add_node!(g, jls[i])
        w = wls[i]
        add_edge!(g, v1, v2, timestamps[i], weight = w)
    end
    g
end
evolving_graph_from_arrays{V,T}(ils::Vector{V}, jls::Vector{V}, timestamps::Vector{T}; is_directed::Bool = true) = evolving_graph_from_arrays(ils, jls, ones(Float64,length(ils)), timestamps, is_directed = is_directed)

deepcopy(g::EvolvingGraph) = EvolvingGraph(is_directed(g),
                                           deepcopy(g.nodes),
                                           deepcopy(g.node_indexof),
                                           deepcopy(g.edges),
                                           deepcopy(g.timestamps),
                                           deepcopy(g.active_nodes),
                                           deepcopy(g.active_node_indexof))

eltype{V,E,T,I}(g::EvolvingGraph{V,E,T,I}) = (V,E,T,I)

is_directed(g::EvolvingGraph) = g.is_directed

nodes(g::EvolvingGraph) = g.nodes
num_nodes(g::EvolvingGraph) = length(nodes(g))

has_node{V}(g::EvolvingGraph{V}, v::V) = v in g.nodes
has_node{V, E, T, KV}(g::EvolvingGraph{V, E, T, KV}, node_key::KV) = node_key in g.node_indexof

unique_timestamps(g::EvolvingGraph) = unique(g.timestamps)
timestamps(g::EvolvingGraph) = g.timestamps
num_timestamps(g::EvolvingGraph) = length(unique_timestamps(g))

active_nodes(g::EvolvingGraph) = g.active_nodes
num_active_nodes(g::EvolvingGraph) = length(g.active_nodes)
has_active_node{V,E,T,KV}(g::EvolvingGraph{V,E,T,KV}, v::TimeNode{KV,T}) =
v in g.active_nodes
has_active_node{V,E,T,KV}(g::EvolvingGraph{V,E,T,KV}, key::KV, t::T) =
get(g.active_node_indexof, (key,t), false) != false ? true : false

edges(g::EvolvingGraph) = g.edges
function edges(g::EvolvingGraph, t)
    inds = findin(g.timestamps, [t])
    if length(inds) == 0
        error("unknown timestamp $(t)")
    end
    return g.edges[inds]
end

num_edges(g::EvolvingGraph) = length(g.edges)

function add_node!{V}(g::EvolvingGraph{V}, v::V)
    push!(g.nodes, v)
    g.node_indexof[v.key] = node_index(g,v)
    v
end
function add_node!{V, E, T, KV}(g::EvolvingGraph{V, E, T, KV}, key::KV)

    id = get(g.node_indexof, key, 0)

    if id == 0
        v = V(g, key)
        return add_node!(g, v)
    else
        return V(id, key)
    end
end

function find_node{V, E, T, KV}(g::EvolvingGraph{V, E, T, KV}, key::KV)
    try
        id = g.node_indexof[v]
        return V(id, v)
    catch
        return false
    end
end

find_node{V}(g::EvolvingGraph{V}, v::V) = v in g.nodes ? v : false

function add_edge!{V,E,T,KV}(g::EvolvingGraph{V,E,T,KV}, v1::KV, v2::KV, t::T; weight::Real = 1.0)
    v1 = add_node!(g, v1)
    v2 = add_node!(g, v2)

    add_edge!(g, v1, v2, t, weight = weight)
end

function add_edge!{V,E,T,KV}(g::EvolvingGraph{V,E,T,KV}, v1::V, v2::V, t::T; weight::Real = 1.0)

    e1 = E <: WeightedTimeEdge? E(v1, v2, weight, t) : E(v1, v2, t)
    if !(e1 in edges(g))

        # we only add active nodes when we add edges to the graph.
        if !((node_key(v1),t) in keys(g.active_node_indexof))
            n1 = TimeNode{KV,T}(g, v1.key, t)
            push!(g.active_nodes, n1)
            g.active_node_indexof[(node_key(v1),t)] = node_index(n1)
        end

        if !((node_key(v2),t) in keys(g.active_node_indexof))
            n2 = TimeNode{KV,T}(g, v2.key, t)
            push!(g.active_nodes, n2)
            g.active_node_indexof[(node_key(v2),t)] = node_index(n2)
        end

        push!(g.edges, e1)
        push!(g.timestamps, t)

        if !(is_directed(g))
            push!(g.edges, edge_reverse(e1))
            push!(g.timestamps, t)
        end
    end
    return e1
end

function add_bunch_of_edges!(g::EvolvingGraph, ebunch)
    for e in ebunch
        if length(e) == 3
            i,j,t = e
            add_edge!(g, i, j, t)
        elseif length(e) == 4
            i, j, t, w = e
            add_edge!(g, i, j, t, weight = w)
        else
            error("Each edge in ebunch must have 3 or 4 elements")
        end

    end
    return g
end

"""
    adjacency_matrix(g, t)

Return an adjacency matrix representation of an evolving graph `g` at timestamp `t`.

# Example

```jldoctest
julia> using EvolvingGraphs

julia> g = evolving_graph_from_arrays([1,2,3], [4,5,2], [1,1,2])
Directed EvolvingGraph 5 nodes, 3 static edges, 2 timestamps

julia> adjacency_matrix(g, 1)
5×5 Array{Float64,2}:
 0.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
```
"""
function adjacency_matrix{V, E, T}(g::EvolvingGraph{V, E, T}, t::T)
    n = num_nodes(g)
    es = edges(g, t)
    A = zeros(Float64, n, n)
    for e in es
        i = node_index(source(e))
        j = node_index(target(e))
        v = E <: WeightedTimeEdge ? e.weight : 1.0
        A[i,j] = v
    end
    return A
end

"""
    sparse_adjacency_matrix(g, t)

Return a sparse adjacency matrix representation of an evolving graph
`g` at timestamp `t`.

# Example

```jldoctest
julia> using EvolvingGraphs

julia> g = evolving_graph_from_arrays([1,2,3], [4,5,2], [1,1,2])
Directed EvolvingGraph 5 nodes, 3 static edges, 2 timestamps

julia> sparse_adjacency_matrix(g,2)
5×5 SparseMatrixCSC{Float64,Int64} with 1 stored entry:
  [5, 3]  =  1.0

julia> sparse_adjacency_matrix(g,1)
5×5 SparseMatrixCSC{Float64,Int64} with 2 stored entries:
  [1, 2]  =  1.0
  [3, 4]  =  1.0
```
"""
function sparse_adjacency_matrix{V,E,T}(g::EvolvingGraph{V,E,T}, t::T, M::Type = Float64)
    n = num_nodes(g)
    is = Int[]
    js = Int[]
    vs = Float64[]
    es = edges(g, t)

    for e in es
        i = node_index(source(e))
        j = node_index(target(e))
        v = E <: WeightedTimeEdge ? e.weight : 1.0
        push!(is, i)
        push!(js, j)
        push!(vs, v)
    end

    return sparse(is, js, vs, n, n)
end

function neighbors{V,E,T,KV}(g::EvolvingGraph{V,E,T,KV}, v::TimeNode{KV,T}; mode::Symbol = :forward)
    r = TimeNode{KV,T}[]

    switch_mode = Dict(
        :forward => (source, target, >),
        :backward => (target, source, <)
    )

    this, that, compare = switch_mode[mode]

    if ! has_active_node(g,v)
        return r
    end

    v_t = node_timestamp(v)
    v_key = node_key(v)

    for e in edges(g)

        # forword neighbors at timestamp v_t
        if edge_timestamp(e) == v_t
            if node_key(this(e)) == v_key
                v_key_new = node_key(that(e))
                v_index = g.active_node_indexof[(v_key_new,v_t)]
                n = TimeNode(v_index, v_key_new, v_t)
                if !(n in r)
                    push!(r, n)
                end
            end
        end

        if compare(edge_timestamp(e), v_t)
            if node_key(source(e)) == v_key || node_key(target(e)) == v_key
                v_t_new = edge_timestamp(e)
                v_index = g.active_node_indexof[(v_key, v_t_new)]
                n = TimeNode(v_index, v_key, v_t_new)
                if !(n in r)
                    push!(r, n)
                end
            end
        end
    end
    return r
end

forward_neighbors{V,E,T,KV}(g::EvolvingGraph{V,E,T,KV}, v::TimeNode{KV,T}) = neighbors(g, v, mode = :forward)

backward_neighbors{V,E,T,KV}(g::EvolvingGraph{V,E,T,KV}, v::TimeNode{KV,T}) = neighbors(g, v, mode = :backward)

"""
    add_graph!(eg, g, t)

Add a static graph `g` at timestamp `t` to an evolving graph `eg`.
"""
function add_graph!{V, E, T}(eg::EvolvingGraph{V, E, T}, g::AbstractStaticGraph, t::T)
    es = edges(g)
    for e in es
        add_edge!(eg, node_key(source(e)), node_key(target(e)), t)
    end
    eg
end

ext: jl


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

name = "yaml_cpp"
version = v"0.7.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/jbeder/yaml-cpp/archive/refs/tags/yaml-cpp-$(version).zip",
                  "4d5e664a7fb2d7445fc548cc8c0e1aa7b1a496540eb382d137e2cc263e6d3ef5")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/yaml-cpp*/
mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN%.*}_gcc.cmake \
    -DCMAKE_BUILD_TYPE=Release \
    -DYAML_BUILD_SHARED_LIBS=ON \
    ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)
platforms = expand_cxxstring_abis(platforms; skip=p->false)

# The products that we will ensure are always built
products = [
    LibraryProduct("libyaml-cpp", :libyaml_cpp)
]

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

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"6")

max_stars_repo_name: singularitti/Yggdrasil

max_issues_repo_name: singularitti/Yggdrasil

max_forks_repo_name: singularitti/Yggdrasil


text:
# Autogenerated wrapper script for TracyProfiler_jll for x86_64-apple-darwin
export libtracy, tracy_capture_exe, tracy_csvexport_exe, tracy_import_chrome_exe, tracy_profiler_exe, tracy_update_exe

using GLFW_jll
using FreeType2_jll
using Capstone_jll
JLLWrappers.@generate_wrapper_header("TracyProfiler")
JLLWrappers.@declare_library_product(libtracy, "@rpath/libtracy-release.dylib")
JLLWrappers.@declare_executable_product(tracy_capture_exe)
JLLWrappers.@declare_executable_product(tracy_csvexport_exe)
JLLWrappers.@declare_executable_product(tracy_import_chrome_exe)
JLLWrappers.@declare_executable_product(tracy_profiler_exe)
JLLWrappers.@declare_executable_product(tracy_update_exe)
function __init__()
    JLLWrappers.@generate_init_header(GLFW_jll, FreeType2_jll, Capstone_jll)
    JLLWrappers.@init_library_product(
        libtracy,
        "lib/libtracy-release.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        tracy_capture_exe,
        "bin/capture-release",
    )

    JLLWrappers.@init_executable_product(
        tracy_csvexport_exe,
        "bin/csvexport-release",
    )

    JLLWrappers.@init_executable_product(
        tracy_import_chrome_exe,
        "bin/import-chrome-release",
    )

    JLLWrappers.@init_executable_product(
        tracy_profiler_exe,
        "bin/Tracy-release",
    )

    JLLWrappers.@init_executable_product(
        tracy_update_exe,
        "bin/update-release",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
load("partialarray.jl")
load("staged.jl")

quote_tuple(t) = expr(:tuple, {t...})
quote_expr(ex) = expr(:quote, ex)

# -- Test kernel construction -------------------------------------------------

function make_kernel_parts(dest, sources)
    body = expr(:block, dest.writes)
    for k = 1:ndims(dest)
        i = uinds[k]
        body = :(
            for ($i) = 1:size(($dest.name), $k)
                $body
            end
        )
    end

    args = append((dest.name,), map(arg->arg.name, sources))

    (args, body)
end

function make_kernel_function(dest, sources)
    args, body = make_kernel_parts(dest, sources)
    f = expr(:(->), quote_tuple(args), body)
end

function make_kernel(dest, sources)
    eval(make_kernel_function(dest, sources))
end
function make_staged_kernel(dest, sources)
    eval(make_staged_kernel_function(dest, sources))
end

function make_staged_kernel_function(dest, sources)
    args, body = make_kernel_parts(dest, sources)
    name = gensym("kernel")
    head = expr(:call, name, args...)

#    f = expr(:function, head, body)
    f = expr(:function, head, quote_expr(body))
    f = :(@staged $f)
    :($f; $name)
end

# -- Test code ----------------------------------------------------------------

VType = Float32

n = dimension(:n)
m = dimension(:m)

A = input_array(:A, VType, (n,m))
B = input_array(:B, VType, (n,m))
dest = output_array(:dest, VType, (n,m))

C = A.*B
dest[] = C

f = make_kernel(dest, (A, B))

(ni, mi) = (2, 3)
Ai = [1 2 3; 4 5 6]
Bi = [1 0 1; 0 1 0]
desti = Array(VType, (ni,mi))
f(desti, Ai, Bi)

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
#
# Function that reads atom information from PDB or mmCIF files
#
function read_atom(record::String;
                   mmCIF::Bool = false,
                   mmCIF_fields::Indexes_mmCIF_fields = Indexes_mmCIF_fields())

  atom = Atom()

  if length(record) < 6
    return nothing
  end

  if record[1:4] == "ATOM" || record[1:6] == "HETATM"

    if ! mmCIF

      atom.name = strip(record[13:16])

      atom.resname = strip(record[17:21])
      resname = alternate_conformation(atom)
      if resname == nothing
        return nothing
      else
        atom.resname = resname
      end

      atom.chain = strip(record[22:22])
      if atom.chain == " "
        atom.chain = "0"
      end

      atom.index = 1
      try
        atom.index_pdb = parse_int(record[7:11])
      catch
        atom.index_pdb = 0
      end
      atom.resnum = parse_int(record[23:26])
      atom.x = parse(Float64,record[31:38])
      atom.y = parse(Float64,record[39:46])
      atom.z = parse(Float64,record[47:54])
      try
        atom.beta = parse(Float64,record[61:66])
      catch
        atom.beta = 0.
      end
      try
        atom.occup = parse(Float64,record[56:60])
      catch
        atom.occup = 0.
      end
      atom.model = 1
      try
        atom.segname = strip(record[73:76])
        if  length(atom.segname) == 0
          atom.segname = "-"
        end
      catch
        atom.segname = "-"
      end

    else # if mmCIF

      mmcif_data = split(record)
      atom.index = 1
      try
        atom.index_pdb = parse(Int,mmcif_data[mmCIF_fields.index])
      catch
        atom.index_pdb = 0
      end
      atom.name = mmcif_data[mmCIF_fields.name]
      atom.resname = mmcif_data[mmCIF_fields.resname]
      resname = alternate_conformation(atom)
      if resname == nothing
        return nothing
      else
        atom.resname = resname
      end
      atom.chain = mmcif_data[mmCIF_fields.chain]
      try
        atom.resnum = parse(Int,mmcif_data[mmCIF_fields.resnum])
      catch
        atom.resnum = 0
      end
      try
        atom.segname = mmcif_data[mmCIF_fields.segname]
      catch
        atom.segname = ""
      end
      atom.x = parse(Float64,mmcif_data[mmCIF_fields.x])
      atom.y = parse(Float64,mmcif_data[mmCIF_fields.y])
      atom.z = parse(Float64,mmcif_data[mmCIF_fields.z])
      atom.beta = parse(Float64,mmcif_data[mmCIF_fields.beta])
      atom.occup = parse(Float64,mmcif_data[mmCIF_fields.occup])
      atom.model = 1

    end

    return atom

  else

    return nothing

  end

end


text:
# Run tests for the deflation

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra

@bench @testset "Deflation (combined with MSLP)" begin

    nep=nep_gallery("dep0");
    nep = DEP(nep.n, nep.A, [0, 0.8])

    n=size(nep,1);

    (λ,v)=newton(nep,v=ones(n))

    v=v/norm(v);
    dnep=deflate_eigpair(nep,λ,v)

    for k=1:3
        λ0=-0.1+0.1im;

        (λ2,v2)=mslp(dnep,maxit=1000,
                     λ=λ0,
                     logger=0
                     );
        @info "Computing eigenvalue k=$k"
        v2=v2/norm(v2);
        # Expand the partial Schur factorization with the computed solution
        dnep=deflate_eigpair(dnep,λ2,v2);
    end

    (λv,V)=get_deflated_eigpairs(dnep);

    # Test that λv and V are now eigpairs
    for i=1:size(λv,1)
        v=V[:,i]/norm(V[:,i]);
        λ=λv[i];
        @test norm(compute_Mlincomb(nep,λ,v))<sqrt(eps())
    end
end

@bench @testset "Deflation modes (new)" begin

    nep=nep_gallery("dep0_sparse",3);
    n=size(nep,1);
    local λ,v;
    (λ,v)=augnewton(nep,v=ones(n),λ=1.65,tol=1e-11)
    dnep=deflate_eigpair(nep,λ,v,mode=:Generic);
    dnep2=deflate_eigpair(nep,λ,v,mode=:SPMF);
    dnep3=deflate_eigpair(nep,λ,v,mode=:MM);
    λ=2+2im;
    @test compute_Mder(dnep,λ) ≈ compute_Mder(dnep2,λ)
    @test compute_Mder(dnep,λ) ≈ compute_Mder(dnep3,λ)

    X=[1 2; 3 4 ; 5 5.0; -1 -1];

    @test compute_Mlincomb(dnep,λ,X) ≈ compute_Mlincomb(dnep3,λ,X)
    @test compute_Mlincomb(dnep,λ,X) ≈ compute_Mlincomb(dnep2,λ,X)

    S=[2 4 ; 5 6.0]

    @test compute_MM(dnep,S,X) ≈ compute_MM(dnep3,S,X)
    @test compute_MM(dnep,S,X) ≈ compute_MM(dnep2,S,X)

    nep=nep_gallery("dep0_sparse",100);
    n = size(nep,1)
    dnep = nep; dnep2 = nep; dnep3 = nep;
    for i = 0:3
        (λ,v)=augnewton(dnep,v=ones(n+i),λ=-0.08,tol=1e-8,maxit=100)
        dnep=deflate_eigpair(dnep,λ,v,mode=:Generic);
        dnep2=deflate_eigpair(dnep2,λ,v,mode=:SPMF);
        dnep3=deflate_eigpair(dnep3,λ,v,mode=:MM);
        λ=2+2im;
        for der = 0:4
            @test compute_Mder(dnep,λ,der) ≈ compute_Mder(dnep2,λ,der)
            @test compute_Mder(dnep,λ,der) ≈ compute_Mder(dnep3,λ,der)
        end
    end

end


text:
using Tilde
using Tilde: ModelClosure
using ZigZagBoomerang
# using Makie
### Define the target distribution and its gradient
using ForwardDiff: gradient!
using LinearAlgebra
using SparseArrays
using NamedTupleTools: select
import TransformVariables as TV

"""
    bouncy(m, data, T = 1000.0; c=10.0, λ=0.1, ρ=0.0, adapt=false)
Draw  samples until time `T` from the posterior distribution of parameters defined in Tilde model `m`, conditional on `data`.
Samples are drawn using the Bouncy particle sampler.
Returns a `Trace` object.
"""
function bouncy(m::ModelClosure, T = 1000.0;  c=10.0, λref=0.1, ρ=0.0, adapt=false) where {A,B}

    ℓ(pars) = logdensity_def(m, pars)

    t = xform(m)

    function f(x)
        (θ, logjac) = transform_and_logjac(t, x)
        -ℓ(θ) - logjac
    end

    function ∇ϕ!(y, x)
        gradient!(y, f, x)
        y
    end

    # Draw a random starting points and velocity
    d = t.dimension

    tkeys = keys(t(zeros(d)))
    r = select(rand(m), tkeys)
    x0 = inverse(t, r)

    t0 = 0.0
    θ0 = randn(d)

    pdmp(∇ϕ!, t0, x0, θ0, T, c, BouncyParticle(sparse(I(d)), 0*x0, λref, ρ); adapt=adapt)

end

m = @model x begin
    α ~ Normal()
    β ~ Normal()
    yhat = α .+ β .* x
    y ~ For(eachindex(x)) do j
        Normal(yhat[j], 2.0)
    end
end

x = randn(3);
truth = [0.61, -0.34, -1.74];

post = m(x=x) | (y=truth,)

trace, final, (num, acc) = @time bouncy(post, c=10)

ts, xs = ZigZagBoomerang.sep(discretize(trace, 0.1))

p = lines(ts, getindex.(xs, 1))
lines!(ts, getindex.(xs, 2), color=:red)

p


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)

function isotropic_two(coef,strain,strainrate)
    isotropic_two_ = load_op_and_grad("./build/libIsotropicTwo","isotropic_two")
    coef,strain,strainrate = convert_to_tensor([coef,strain,strainrate], [Float64,Float64,Float64])
    isotropic_two_(coef,strain,strainrate)
end

N = 1
γ = rand(N,9)
ε = rand(N, 3)
dotε = rand(N,3)
out = zeros(N,3)
for i = 1:N
    A = [ε[i,1] ε[i,3]/2
        ε[i,3]/2 ε[i,2]]
    B = [dotε[i,1] dotε[i,3]/2
        dotε[i,3]/2 dotε[i,2]]
    T = γ[1] * I + γ[2] * A + γ[3] * B + γ[4] * A * A + γ[5] * (A*B+B*A) +
        γ[6] * B * B + γ[7] * (A * A * B + B * A * A) + γ[8] * (A * B * B + B * B * A) + γ[9] * (A * A * B * B + B * B * A * A)
    out[i, :] = [T[1,1];T[2,2];T[1,2]]
end
# TODO: specify your input parameters
u = isotropic_two(γ,ε,dotε)
sess = Session(); init(sess)
@show run(sess, u)-out

# uncomment it for testing gradients
# error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(m)
    # return sum(isotropic_two(γ,ε,dotε)^2)

    # return sum(isotropic_two(γ,ε,m)^2)

    return sum(isotropic_two(γ,m,dotε)^2)

    # return sum(isotropic_two(m,ε,dotε)^2)

end

# TODO: change `m_` and `v_` to appropriate values
# m_ = constant(rand(N, 9))
# v_ = rand(N, 9)

m_ = constant(rand(N,3))
v_ = [rand(N,1) rand(N,1) rand(N,1)]
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
using Test, DiscreteExteriorCalculus, DiscretePDEs
using UniqueVectors: UniqueVector

geo_code = GeoCode(gdspy_open(joinpath(@__DIR__, "two_pads.gds"), "main"),
    layer_translation=Dict(24 => -50),
    layer_extrusion=Dict(22 => -50, 23 => -50, 24 => 250),
    layer_volume_difference=[(23, 22, true), (24, 23, true), (24, 22, true)],
    layer_volume_groups=Dict("Substrate" => [22, 23], "Vacuum" => [24]))
file_name = joinpath(@__DIR__, "two_pads.geo")
geo_write!(file_name, geo_code, characteristic_length_factor=.35,
    footer="""
    Physical Surface("Ground") = {4};
    Physical Surface("Pad 1") = {5};
    Physical Surface("Pad 2") = {6};
    """)

initialize!()
gmsh_open!(file_name)
N, K = 3, 4
mesh!(K)

node_tags, points, tcomp = get_triangulated_complex(N, K, 1e-3)
group_dict = get_physical_groups(node_tags, points)
@test typeof(tcomp) <: TriangulatedComplex{N, K}
comp = tcomp.complex
orient!(comp)
sources = [get_charge_source(comp, group_dict[k])
    for k in ["Pad 1", "Pad 2"]]
exterior = boundary(comp)
m = Metric(N)
mesh = Mesh(tcomp, circumcenter(m))

ϵr = 2
ϵ_form = DiscretePDEs.ϵ₀ * (get_material(comp, 1, 2) +
    get_material(comp, group_dict["Substrate"], ϵr-1, 2))

bbox, null_basis = electrostatics_blackbox(m, mesh, sources,
    append!(subcomplex(comp, group_dict["Ground"]), boundary(comp)), ϵ_form)

capacitance = admittance_matrix(bbox, null_basis)

q3d_capacitance = [[80.995, -1.2219] [-1.2219, 80.994]] * 1e-12
@testset "capacitance two pads" begin
    @test isapprox(capacitance, q3d_capacitance, rtol=.25)
end

if false # set to true to display plots with Gmsh
    φ, q = solve_statics(bbox, null_basis, [1.0, -1.0])
    ρ = source_density(m, mesh, q, 1)
    vec_E = electric_field(m, mesh, φ)
    comp_points = UniqueVector([c.points[1] for c in comp.cells[1]])
    ordering = [findfirst(isequal(p), comp_points) for p in points]
    add_field!("Electric field", node_tags, vec_E[ordering])
    sub_points = UniqueVector([c.points[1] for c in subcomplex(comp,
        vcat(group_dict["Pad 1"], group_dict["Pad 2"], group_dict["Ground"])).cells[1]])
    sub_inds = findall([p in sub_points for p in points])
    add_field!("Charge density", node_tags[sub_inds], ρ[ordering[sub_inds]])
    add_field!("Voltage", node_tags[sub_inds], φ[ordering[sub_inds]])
    gui!()
end
rm(file_name)


text:
 using Flux

"""
`gate(x::T, y::T) where T <: AbstractArray`

* __input:__ arrays of equal length of Float32
* __output:__ array of Float32

The additional `y` vector is used to apply local and global conditioning.
`z = tanh(x[1:f÷2] + y[1:f÷2]) * σ(x[f÷2+1:end] + y[f÷2+1:end])`
"""
function gate(x)
    chan = size(x, 3) ÷ 2
    tanh.(x[:, :, 1:chan, :]) .* σ.(x[:, :, chan+1:end, :])
end

function gate(x, y)
    chan = size(x, 3) ÷ 2
    tanh.(x[:, :, 1:chan, :] .+ y[:, :, 1:chan, :]) .*
    σ.(x[:, :, chan+1:end, :] .+ y[:, :, chan+1:end, :])
end

"""
`WConv(dilate::Int, skpflt::Int, resflt::Int)`

* __input:__
    * `dilate` how wide the dilation should be for this conv layer
    * `resflt` how many filters to return for the residual layer
    * `skpflt` how many filters to return for the skip layer
* __output:__ a WConv struct

`WConv(x::AbstractArray)`

* __input:__ a (data, 1, channels, batch) Array like a 1D Conv
* __output:__ an Array processed by gated activation, a skip layer applied, and residual values stored in the `WConv` struct.

`WConv(x::Tuple)`

* __input:__ a Tuple of arrays
* __output:__ same as above but also with conditioning array applied in gated activation

"""
mutable struct WConv
    dilconv::Conv
    resconv::Conv
    skpconv::Conv
    skp::DenseArray
    condconv::Conv

    function WConv(dil::Int, res::Int, skp::Int, cond::Int)
        dc = gpu(Conv((2, 1), res => res*2, identity, dilation=(dil, 1), pad=(dil, 0)))
        rc = gpu(Conv((1, 1), res => res, identity))
        sc = gpu(Conv((1, 1), res => skp, identity))
        s = sc(gpu(zeros(Float32, 2, 1, res, 1)))
        cc = gpu(Conv((1, 1), cond => res*2, tanh))
        new(dc, rc, sc, s, cc)
    end

    function WConv(dc::Conv, rc::Conv, sc::Conv, s::DenseArray, cc::Conv)
        new(dc, rc, sc, s, cc)
    end
end

function (wc::WConv)(x)
    w,h,c,b = size(x)
    dat = gate(wc.dilconv(x))[1:w,:,:,:]
    wc.skp = wc.skpconv(dat)
    x .+ wc.resconv(dat)
end

function (wc::WConv)(input::Tuple)
    x, y = input
    w,h,c,b = size(x)
    dat = gate(wc.dilconv(x)[1:w,:,:,:], wc.condconv(y))
    wc.skp = wc.skpconv(dat)
    (x .+ wc.resconv(dat), y)
end

function Base.show(io::IO, wc::WConv)
    ln1 = ":: Wave.WConv ::\n"
    ln2 = "Dilated Conv:   " * string(wc.dilconv) *
    " dil: " * string(wc.dilconv.dilation) *
    " pad: " * string(wc.dilconv.pad) * "\n"
    ln3 = "Residual Conv:  " * string(wc.resconv) * "\n"
    ln4 = "Skip Conv:      " * string(wc.skpconv) * "\n"
    ln5 = "Cond Conv: " * string(wc.condconv) * "\n"
    print(io, ln1 * ln2 * ln3 * ln4 * ln5)
end

Flux.@treelike WConv

import Base: +
+(x::WConv, y::WConv) = x.skp + y.skp
+(x::DenseArray, y::WConv) = x + y.skp

"""
`Net(nlayers::Int; res=64, skip=256, cond=0)`

* __input:__
    * `nlayers` how many dilated layers to include
    * `res` how many channels to pass into residual layers
    * `skip` how many channels to pass into skip layers
    * `depth` reset dilation to 2 at every multiple of depth
    * `cond` how many variables are included in the conditioning tensor
* __output:__ a NetNet struct

`Net(x)`

* __input:__ array of shape (samples, 1, 1channel, batchsize)
* __output:__ probability distribution (softmax) of wave amplitudes

`Net(x::Tuple)`

* __input:__ Tuple of data and conditioning Array of shape (samples, condvars, batchsize)
* __output:__ probability distribution (softmax) of wave amplitudes

"""
mutable struct Net
    input::Conv
    layers::Chain
    output::Chain

    function Net(nlayers; res=64, skip=256, depth=9, cond=0)
        i = gpu(Conv((2, 1), 1=>res, tanh, pad=(1, 0)))
        l = Chain([WConv(2^(i%depth), res, skip, cond) for i in 1:nlayers]...)
        o = Chain(
            # x -> leakyrelu.(x),
            gpu(Conv((1, 1), skip=>256, leakyrelu)),
            gpu(Conv((1, 1), 256=>256, leakyrelu))
        )
        new(i, l, o)
    end

    function Net(i::Conv, l::Chain, o::Chain)
        new(i, l, o)
    end
end

Flux.@treelike Net

function (w::Net)(x)
    w.layers(w.input(x)[1:end-1,:,:,:])
    softmax(w.output(reduce(+, w.layers))[end, 1, :, :])
end

function (w::Net)(input::Tuple)
    x, y = input
    dat = (w.input(x)[1:end-1,:,:,:], y)
    w.layers(dat)
    softmax(w.output(reduce(+, w.layers))[end, 1, :, :])
end

function Base.show(io::IO, wave::Net)
    ln1 = ":: Wave.Net ::\n"
    ln2 = "Hidden layers:     " * string(length(wave.layers)) * "\n"
    ln3 = "Receptive field:   " * string(2^length(wave.layers)) * "\n"
    ln4 = "Residual channels: " * string(size(wave.layers[end].resconv.weight, 4)) * "\n"
    ln5 = "Skip channels:     " * string(size(wave.layers[end].skpconv.weight, 4)) * "\n"
    ln6 = "Conditioning vars: " * string(size(wave.layers[end].condconv.weight, 3)) * "\n"
    ln7 = "GPU enabled:       " * ifelse(typeof(wave.input.weight) <: Array, "no", "yes")
    print(io, ln1 * ln2 * ln3 * ln4 * ln5 * ln6 * ln7)
end


text:
#
println("An) Test of a mean-field basis and configuration-interaction (CI) expansion.")
#

name        = "Oxygen 1s^2 2s^2 2p^4 ground configuration"
refConfigs  = [Configuration("[He] 2s^2 2p^4")]
mfSettings  = MeanFieldSettings()
#
wa          = Representation(name, Nuclear.Model(8.), Radial.Grid(true), refConfigs, MeanFieldBasis(mfSettings) )
println("wa = $wa")

wb = generate(wa, output=true)

orbitals    = wb["mean-field basis"].orbitals
ciSettings  = CiSettings(CoulombInteraction(), LevelSelection())
from        = [Shell("2s")]
#
frozen      = [Shell("1s")]
to          = [Shell("2s"), Shell("2p")]
excitations = RasStep()
#             RasStep(RasStep(), seFrom=from, seTo=deepcopy(to), deFrom=from, deTo=deepcopy(to), frozen=deepcopy(frozen))
#
wc          = Representation(name, Nuclear.Model(8.), Radial.Grid(true), refConfigs,
                             CiExpansion(orbitals, excitations, ciSettings) )
println("wc = $wc")

wd = generate(wc, output=true)


text:
module PlotlyBaseTest
using TestSetExtensions

using Base.Test

using PlotlyBase
const M = PlotlyBase

try
    @testset ExtendedTestSet "PlotlyJS Tests" begin
        @includetests ARGS
    end
catch
    exit(-1)
end

end


text:
module Scryfall

export getOracle, getImgurl

import HTTP
import JSON

#=
"oracle_text"      => "{T}, Sacrifice Black Lotus: Add three mana of any one color."
"scryfall_set_uri" => "https://scryfall.com/sets/vma?utm_source=api"
"set"              => "vma"
"lang"             => "en"
"rarity"           => "mythic"
"colors"           => Any[]
"legalities"       => Dict{String,Any}("frontier"=>"not_legal","1v1"=>"banned","modern"=>"not_legal","penny"=>"not_legal","commander"=>"banned","standard"=>"not_legal","duel"=>"…
"name"             => "Black Lotus"
"oracle_id"        => "5089ec1a-f881-4d55-af14-5d996171203b"
"reserved"         => true
"uri"              => "https://api.scryfall.com/cards/bd8fa327-dd41-4737-8f19-2cf5eb1f7cdd"
"type_line"        => "Artifact"
"mana_cost"        => "{0}"
"id"               => "bd8fa327-dd41-4737-8f19-2cf5eb1f7cdd"
"scryfall_uri"     => "https://scryfall.com/card/vma/4/black-lotus?utm_source=api"
"rulings_uri"      => "https://api.scryfall.com/cards/bd8fa327-dd41-4737-8f19-2cf5eb1f7cdd/rulings"
"image_uris"       => Dict{String,Any}("normal"=>"https://img.scryfall.com/cards/normal/en/vma/4.jpg?1517813031","border_crop"=>"https://img.scryfall.com/cards/border_crop/en/vm…
"frame"            => "2015"
"full_art"         => false
"color_identity"   => Any[]
"reprint"          => true
"colorshifted"     => false
"digital"          => true
"multiverse_ids"   => Any[382866]
"illustration_id"  => "da62ded1-bedd-44c6-8950-ca56e691a899"
"highres_image"    => true
"cmc"              => 0.0
"set_name"         => "Vintage Masters"
"set_uri"          => "https://api.scryfall.com/sets/vma"
=#

function getRaw(fuzzyName::AbstractString;kwargs...)
    fuzzyName = HTTP.URIs.escapeuri(fuzzyName) #encode the fuzzyName for GET purpose
    kwargs = Dict(kwargs) #make keywords arguments a dictonary
    if haskey(kwargs, :setCode) #if 3-letter keyset is specified
        setCode = kwargs[:setCode]
        if length(setCode) == 3
            fuzzyRequest = HTTP.request("GET","https://api.scryfall.com/cards/named?set=$setCode&fuzzy=$fuzzyName")
        else
            fuzzyRequest = HTTP.request("GET","https://api.scryfall.com/cards/named?set=&fuzzy=$fuzzyName")
            @warn ("Length of setCode is incorrect, 3-letter only")
        end
    else #fallback to using anyset, let scryfall decide
        fuzzyRequest = HTTP.request("GET","https://api.scryfall.com/cards/named?set=&fuzzy=$fuzzyName")
    end
    cardDict = JSON.parse(String(fuzzyRequest.body)) #return a dictonary for wrapper functions to use
    #= println(cardDict["name"]) =#
    return cardDict
end

# wrapper function for getting oracle
function getOracle(fuzzyName::AbstractString; kwargs...)
    try
        result = getRaw(fuzzyName; kwargs...)["oracle_text"]
        sleep(0.1)
        return result
    catch
        @warn ("error in fetching from scryfall")
        return "No matching on scryfall"
    end
end

# wrapper function for getting imgurl, the "normal" one
function getImgurl(fuzzyName::AbstractString; kwargs...)
    try
        result = getRaw(fuzzyName; kwargs...)["image_uris"]["normal"]
        sleep(0.1)
        return result
    catch
        @warn ("error in fetching from scryfall, may due to too many results")
        return "No matching on scryfall"
    end
end
end


text:
"""
	SeisCopy(in,out)

Copy a seis input file to a seis output file. In and out should be of type AbstractString

*Credits: AS, 2015*
"""
function SeisCopy(in::AbstractString,out::AbstractString)

	filename_data_in = ParseDataName(in)
	filename_headers_in = ParseHeaderName(in)
	DATAPATH = get(ENV,"DATAPATH",join([pwd(),"/"]))
	filename_data_out = join([DATAPATH out "@data@"])
	if filename_headers_in != "NULL"
		filename_headers_out = join([DATAPATH out "@headers@"])
	else
		filename_headers_out = "NULL"
	end
	extent = ReadTextHeader(in)
	WriteTextHeader(out,extent,"native_float",4,filename_data_out,filename_headers_out)
	run(`cp $filename_data_in $filename_data_out`)
	if filename_headers_in != "NULL"
		run(`cp $filename_headers_in $filename_headers_out`)
	end

end

function SeisCopy(in::Array{AbstractString,1},out::Array{AbstractString,1})

	for j = 1 : length(in)
		SeisCopy(in[j],out[j])
	end

end


text:
using Pkg, Test

using CLIMA.TurbulenceConvection.Grids
using CLIMA.TurbulenceConvection.StateVecs
using CLIMA.TurbulenceConvection.StateVecFuncs

n_subdomains = 3 # number of sub-domains
n_elems_real = 10 # number of elements

grid = Grid(0.0, 1.0, n_elems_real)
vars = ( (:ρ_0, 1), (:a, n_subdomains), (:w, n_subdomains),
         (:ϕ, n_subdomains), (:ψ, n_subdomains) )
state_vec = StateVec(vars, grid)
vars = ((:w_ave, 1), (:cv_ϕ_ψ, n_subdomains), (:TCV_ϕ_ψ, 1))
tmp = StateVec(vars, grid)

@testset "Export fields" begin
  for k in over_elems(grid)
    state_vec[:ρ_0, k] = 1*k
    for i in over_sub_domains(state_vec)
      state_vec[:w, k, i] = 2*k
      state_vec[:a, k, i] = 3*k
    end
  end
  export_state(state_vec, "./", "state_vec")
  export_state(state_vec, grid, "./", "state_vec")
  export_state(state_vec, grid, "./", "state_vec", UseVTK())
end

@testset "Assign ghost, extrapolate, surface funcs" begin
  for k in over_elems(grid)
    state_vec[:a, k] = 2
  end

  assign_ghost!(state_vec, :a, 0.0, grid)
  @test all(state_vec[:a, k, 1] ≈ 0.0 for k in over_elems_ghost(grid))

  @test surface_val(state_vec, :a, grid) ≈ 1
  @test first_elem_above_surface_val(state_vec, :a, grid) ≈ 2

  extrap!(state_vec, :a, grid, 1)
  k = 1+grid.n_ghost
  @test state_vec[:a, k-1] ≈ 2*state_vec[:a, k] - state_vec[:a, k+1]
  k = grid.n_elem-grid.n_ghost
  @test state_vec[:a, k+1] ≈ 2*state_vec[:a, k] - state_vec[:a, k-1]
end

@testset "Domain average" begin
  state_vec[:a, 1, 1] = 0.25
  state_vec[:a, 1, 2] = 0.75
  state_vec[:a, 1, 3] = 0
  state_vec[:w, 1, 1] = 2
  state_vec[:w, 1, 2] = 2
  state_vec[:w, 1, 3] = 2
  domain_average!(tmp, state_vec, state_vec, (:w_ave,), (:w,), :a, grid)
  @test tmp[:w_ave, 1] ≈ 2
end

@testset "Distribute" begin
  tmp[:w_ave, 1] = 2

  distribute!(state_vec, tmp, (:w,), (:w_ave,), grid)
  @test state_vec[:w, 1, 1] ≈ 2
  @test state_vec[:w, 1, 2] ≈ 2
  @test state_vec[:w, 1, 3] ≈ 2
end

@testset "Total covariance" begin
  state_vec[:a, 1, 1] = 0.1
  state_vec[:a, 1, 2] = 0.2
  state_vec[:a, 1, 3] = 0.3
  state_vec[:ϕ, 1, 1] = 1
  state_vec[:ϕ, 1, 2] = 2
  state_vec[:ϕ, 1, 3] = 3
  state_vec[:ψ, 1, 1] = 2
  state_vec[:ψ, 1, 2] = 3
  state_vec[:ψ, 1, 3] = 4
  tmp[:cv_ϕ_ψ, 1, 1] = 1.0
  tmp[:cv_ϕ_ψ, 1, 2] = 1.0
  tmp[:cv_ϕ_ψ, 1, 3] = 1.0
  decompose_ϕ_ψ(cv) = cv == :cv_ϕ_ψ ?  (:ϕ , :ψ) : error("Bad init")
  total_covariance!(tmp, state_vec, tmp, state_vec,
                    (:TCV_ϕ_ψ,), (:cv_ϕ_ψ,), :a, grid, decompose_ϕ_ψ)
  @test tmp[:TCV_ϕ_ψ, 1] ≈ 0.8
end

@testset "Integrate ODE" begin
  c_0 = 3.0
  y_0 = 2.0
  func(z, args) = y_0 + z
  sol_correct = c_0 .+ [y_0*get_z(grid, k) + 1/2*get_z(grid, k)^2 for k in over_elems_real(grid)]
  integrate_ode!(state_vec, :w, grid, func, c_0, NamedTuple(), 1)
  sol_approx = [state_vec[:w, k] for k in over_elems_real(grid)]
  sol_percent_error = [(sol_approx[k]-sol_correct[k])/sol_correct[k]*100 for k in 1:length(sol_correct)]
  @test all(sol_percent_error .< 100*grid.dz)
end

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.3
#     language: julia
#     name: julia-1.6
# ---

# %%
using Distributions
using StatsPlots

function plot_tstat(;
        n = 10,
        L = 10^5,
        dist = Normal(1, 2),
        xlim = (-5, 5),
        xtick = -5:5,
        bin = range(extrema(xtick)...; length=51),
    )
    tstat(sample) = (mean(sample) - mean(dist))/√(var(sample)/n)
    T = [tstat(rand(dist, n)) for _ in 1:L]
    histogram(T; norm=true, alpha=0.3, label="(x - μ)/√(u²/n)", xlim, xtick, bin)
    plot!(TDist(n-1); lw=1.5, label="TDist(n - 1)")
    sk = skewness(dist)
    diststr = replace("$dist", r"\{.*\}"=>"")
    title!("$diststr,  n = $n\nskewness = $(round(sk; digits=3))"; titlefontsize=10)
end

function plot_chisqstat(;
        n = 10,
        L = 10^5,
        dist = Normal(1, 2),
    )
    X = (n-1)*[var(rand(dist, n)) for _ in 1:L]/var(dist)
    m, s = mean(X), std(X)
    xlim = (max(-1, m - 4s), n + 4s)
    bin = range(0, xlim[2]; length=round(Int, 51/(1 - xlim[1]/xlim[2])))
    histogram(X; norm=true, alpha=0.3, label="(n - 1)u²/σ²", xlim, bin)
    plot!(Chisq(n-1); lw=1.5, label="Chisq(n - 1)")
    kurt = kurtosis(dist)
    if isfinite(kurt)
        M, S = n - 1, (n - 1)*√((kurt + 3)/n - (n - 3)/(n*(n - 1)))
        plot!(Normal(M, S); lw=1.5, ls=:dash, label="normal dist.")
    end
    diststr = replace("$dist", r"\{.*\}"=>"")
    title!("$diststr,  n = $n\nkurtosis = $(round(kurt; digits=3))"; titlefontsize=10)
end

function plot_both(;
        n = 10,
        L = 10^5,
        dist = Normal(1, 2),
        xlim1 = (-5, 5),
        xtick1 = -5:5,
        bin1 = range(extrema(xtick1)...; length=51),
    )
    P = plot_tstat(; n, L, dist, xlim=xlim1, xtick=xtick1, bin=bin1)
    Q = plot_chisqstat(; n, L, dist)
    plot(P, Q; size=(800, 300))
end

# %%
for n in (10, 20, 50, 100, 200, 400)
    plot_both(; n) |> display
end

# %%
for n in (10, 20, 50, 100, 200, 400)
    plot_both(; dist = Uniform(0, 1), n) |> display
end

# %%
for n in (10, 20, 50, 100, 200, 400)
    plot_both(; dist = TDist(4), n) |> display
end

# %%
for n in (10, 20, 50, 100, 200, 400)
    plot_both(; dist = TDist(5), n) |> display
end

# %%
for n in (10, 20, 50, 100, 200, 400)
    plot_both(; dist = Gamma(2, 1), n) |> display
end

# %%
for n in (10, 20, 50, 100, 200, 400)
    plot_both(; dist = Gamma(4, 1), n) |> display
end

# %%
for n in (10, 20, 50, 100, 200, 400)
    plot_both(; dist = Gamma(10, 1), n) |> display
end

# %%
for n in (10, 20, 50, 100, 200, 400)
    plot_both(; dist = Exponential(), n) |> display
end

# %%

max_stars_repo_path: 0021/plot tstat and chisqatat.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0021/plot tstat and chisqatat.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0021/plot tstat and chisqatat.jl

max_forks_repo_name: genkuroki/public


text:
@safetestset "Router tests" begin

  @safetestset "Basic routing" begin
    using Genie, Genie.Router

    route("/hello") do
      "Hello"
    end

  end;

  @safetestset "router_delete" begin
    using Genie, Genie.Router

    x = route("/caballo") do
      "caballo"
    end

    @test (x in routes()) == true
    Router.delete!(:get_caballo)
    @test (x in routes()) == false
  end;

end;


text:
using OpenEphysLoader, DSP

type AnalogData
  #data, with each row being a different channel and each column a time value
  x_all::Array{Float64,2}
  t::Vector{Float64} #time
  original_fs::Int64
  fs::Int64
  channel_nums::Vector{Int64} #desired channel numbers
end

"""
    AnalogData(x_all::Array{Float64,2}, t::Vector{Float64};
    original_fs::Int64=30000, channel_nums::Vector{Int64}=[0,0])

Create analogdata object given only data and time- assumes all channels are
desired.
"""
function AnalogData(x_all::Array{Float64,2}, t::Vector{Float64};
  original_fs::Int64=30000, channel_nums::Vector{Int64}=[0,0])
  if channel_nums == [0,0]
    channel_nums = collect(1:size(x_all)[1])
  end
  fs = original_fs
  AD = AnalogData(x_all::Array{Float64,2}, t::Vector{Float64},
  original_fs::Int64, fs::Int64, channel_nums::Vector{Int64})
end

"""
    ad_equals(adone::AnaloglData, adtwo::AnalogData)
Check whether two analogdata objects are equal.
"""
function ad_equals(adone::AnalogData, adtwo::AnalogData)
  if ((adone.x_all == adtwo.x_all) && (adone.t == adtwo.t) &&
    (adone.original_fs == adtwo.original_fs)
    && (adone.fs == adtwo.fs) && (adone.channel_nums == adtwo.channel_nums))
    return true
  else
    return false
  end
end

#Base.==(x::AnalogData, y::AnalogData) =  ((x.x_all == y.x_all) && (x.t == y.t)
#&& (x.original_fs == y.original_fs) && (x.fs == y.fs) && (x.channel_nums == y.channel_nums))

"""
    load_continuous(path::String, fs::Int64)

Return an array of the data and an array of the time values given a file path
and sampling rate
"""
function load_continuous(path::String, fs::Int64)
  A = nothing
  open("$path", "r") do io
    A = Array(SampleArray(Float64, io))
    #sleep(0.1)
    seekstart(io)
    T = Array(TimeArray(Float64, io))
    return(A,T)
  end
end

"""
    load_continuous_channels(prefix::String, data_directory::String,fs::Int64,
    channel_nums::Vector{Int64}, recording_num::Int64=1)

Creates an analogdata object given the desired channel numbers, data directory, and data prefix.
Recording num is to account for if there are multiple recordings,
in which case the appropriate number is added to the filename.
"""
function load_continuous_channels(prefix::String, data_directory::String,
  fs::Int64, channel_nums::Vector{Int64}, recording_num::Int64=1)
  x_all = Nullable{Array{Float64,2}}()
  t = nothing
  if recording_num==1
    recording_num_str = ""
  elseif recording_num>1
    recording_num_str = "_$recording_num"
  else
    error("load_continuous_channels: invalid recording number")
  end
  for chan_name in channel_nums
    #TODO how to make this not have just my username?
    filename =
    "/Users/mj2/data/eeg_tests/$data_directory/$prefix$chan_name$recording_num_str.continuous"
    (x,t) = load_continuous("$filename", fs)
    if isnull(x_all)
      x_all = x'
    else
      x_all = vcat(x_all, x')
    end
  end
  data = AnalogData(x_all, t; original_fs=fs, channel_nums=channel_nums)
  return data
end

type DigitalData
  #data, with each row being a different channel and each column a time value
  x_all::Array{Bool,2}
  t::Vector{Float64} #time
  original_fs::Int64
  fs::Int64
  channel_nums::Vector{Int64} #desired channel numbers
end

"""
    DigitalData(x_all::Array{Bool,2}, t::Vector{Float64};
    original_fs::Int64=30000, channel_nums::Vector{Int64}=[0,0])

Create digitaldata object given only data and time- assumes all channels are
desired.
"""
function DigitalData(x_all::Array{Bool,2}, t::Vector{Float64};
  original_fs::Int64=30000, channel_nums::Vector{Int64}=[0,0])
  if channel_nums == [0,0]
    channel_nums = collect(1:size(x_all)[1])
  end
  fs = original_fs
  DD = DigitalData(x_all::Array{Bool,2}, t::Vector{Float64},
  original_fs::Int64, fs::Int64, channel_nums::Vector{Int64})
end

"""
    dd_equals(ddone::DigitalData, ddtwo::DigitalData)
Check whether two digitaldata objects are equal.
"""
function dd_equals(ddone::DigitalData, ddtwo::DigitalData)
  if ((ddone.x_all == ddtwo.x_all) && (ddone.t == ddtwo.t) &&
    (ddone.original_fs == ddtwo.original_fs)
    && (ddone.fs == ddtwo.fs) && (ddone.channel_nums == ddtwo.channel_nums))
    return true
  else
    return false
  end
end

type Spectrogram
  analog_data::AnalogData
  power_all
  freq_bins::DSP.Util.Frequencies
  time_bins::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
  time::Array{Float64,1}
  n::Int64
end

"""
    Spectrogram(analog_data::AnalogData, n=1024)
Create a Spectrogram object from an AnalogData object. N decides the size of time bins
"""
function Spectrogram(analog_data::AnalogData, n::Int64=256)
  power_all = Array{Float64}[]
  freq_bins = nothing
  time_bins = nothing
  for chan_name = 1:(size(analog_data.channel_nums)[1])
    s = spectrogram(analog_data.x_all[chan_name, :], n)
    freq_bins = freq(s)
    time_bins = time(s)
    #TODO think about clamp range below!
    power_x = 20.*log10(clamp!(power(s), .00001 , 100000000))
    push!(power_all, power_x)
  end
  ti = analog_data.t
  S = Spectrogram(analog_data, power_all, freq_bins, time_bins, ti, n)
end

type Session
  name::String
  fs_openephys::Int64
  directory::String
  eeg_data::Nullable{AnalogData}
  spectrum::Nullable{Spectrogram}
  ica::Nullable{AnalogData}
  ica_spectrum::Nullable{Spectrogram}
end

"""
    Session(name::String, directory::String, eeg_data=Nullable{AnalogData}(), n::Int64=256))
Create a Session object. Providing eeg_data will automatically create a spectrogram, so change n
here if necessary.
"""
function Session(name::String, directory::String, eeg_data=Nullable{AnalogData}(), n::Int64=256)
  if isnull(eeg_data)
    S = Session(name, 30000, directory, eeg_data, Nullable{Spectrogram}(),
    Nullable{AnalogData}(), Nullable{Spectrogram}())
  else S = Session(name, 30000, directory, eeg_data, Spectrogram(eeg_data, n),
    Nullable{AnalogData}(), Nullable{Spectrogram}())
  end
end

"""
    function session_equals(s1::Session, s2::Session)
Check whether two session objects are equal.
"""
function session_equals(s1::Session, s2::Session)
  if ((s1.name==s2.name) && (s1.fs_openephys==s2.fs_openephys) && (s1.directory==s2.directory) &&
    ad_equals(get(s1.eeg_data), get(s2.eeg_data)) && (get(s1.spectrum)==get(s2.spectrum)) &&
    ad_equals(get(s1.ica), get(s2.ica)) && (get(s1.ica_spectrum)==get(s2.ica_spectrum)))
    return true
  else
    return false
  end
end

"""
    load_eeg(session::Session, channel_nums::Vector{Int64})
Load EEG data into a session
"""
function load_eeg(session::Session, channel_nums::Vector{Int64}, prefix::String="100_CH")
  if !isnull(session.eeg_data)
    error("load_eeg: eeg_data already loaded")
  else
    session.eeg_data = load_continuous_channels(prefix, session.directory,
    session.fs_openephys, channel_nums)
    session.spectrum = Spectrogram(get(session.eeg_data))
  end
end

#you can get the eeg_data in the form of a non-nullable analog data type by using get(session.eeg_data)

"""
    lowpass_session(session::Session, lowpass_cutoff::Float64, lowpass_fs::Int64,
    down_sample_factor::Int64, lowpass_order::Int64=5, prefix::String="100_CH")
Filter the eeg_data of a session with both lowpass and downsampling and return
the new filtered session. If no downsampling factor is provided, it is assumed to be 1
and no downsampling occurs. Change n here if you want a different value for the
spectrogram of the filtered data.
"""
function lowpass_session(session::Session, lowpass_cutoff::Float64,
   down_sample_factor::Int64=1, lowpass_order::Int64=5,
  prefix::String="100_CH"; n::Int64=256)
  if isnull(session.eeg_data)
    error("lowpass_session: no eeg_data loaded")
  else
    new_eeg_data = down_sample(lowpass(get(session.eeg_data), lowpass_cutoff,
     lowpass_order), down_sample_factor)
    filt_session = Session("Filtered $(session.name) (fL $(lowpass_cutoff) Hz, $(new_eeg_data.fs) fs)",
    session.directory, new_eeg_data, n)
    return filt_session
  end
end

"""
    highpass_session(session::Session, highpass_cutoff::Float64, highpass_order::Int64=5,
    prefix::String="100_CH"; n::Int64=256)
Highpass filter the eeg_data of a session and return the new filtered session. Change n here if you
want a different value for the spectrogram of the filtered data.
"""
function highpass_session(session::Session, highpass_cutoff::Float64, highpass_order::Int64=5,
  prefix::String="100_CH"; n::Int64=256)
  if isnull(session.eeg_data)
    error("highpass_session: no eeg_data loaded")
  else
    new_eeg_data = highpass(get(session.eeg_data), highpass_cutoff, get(session.eeg_data).fs, highpass_order)
    filt_session = Session("Filtered $(session.name) (fH $(highpass_cutoff) Hz, $(new_eeg_data.fs) fs)",
    session.directory, new_eeg_data, n)
    return filt_session
  end
end


text:
module LaleLibOps

using PyCall
import Pandas

using DataFrames
using Random
using ..AbsTypes
using ..Utils
using ..LaleAbsTypes

import ..AbsTypes: fit!, transform!
import ..LaleAbsTypes: fit, transform, predict

export fit!, transform!
export fit, transform, predict

import Base: >>,|,+,|>,&
export  >>,|,+,|>,&

export LalePipeOptimizer, lalepipeoptimizers
export LalePipe, visualize

const optim_dict    = Dict{String, PyObject}()
const LALELIBS      = PyNULL()
const LALEOPS       = PyNULL()

function __init__()
   copy!(LALELIBS, pyimport("lale.lib.lale"))
   copy!(LALEOPS,  pyimport("lale.operators"))

   optim_dict["Hyperopt"]     = LALELIBS.Hyperopt
   optim_dict["GridSearchCV"] = LALELIBS.GridSearchCV
   global _ConcatFeatures     = LALELIBS.ConcatFeatures
   global _make_pipeline      = LALEOPS.make_pipeline
   global _make_choice        = LALEOPS.make_choice
   global _make_union         = LALEOPS.make_union
   global _make_union_nc      = LALEOPS.make_union_no_concat
end

# hide pyobject
struct LalePipe <: LaleOperator
   model::Dict{Symbol, Any}
   LalePipe(x::PyObject) = new(Dict{Symbol,Any}(:laleobj=>x))
end

⊖(a::PyObject,b::PyObject)          = _make_pipeline(a,b)
|>(a::LaleOperator,b::LaleOperator) = a.model[:laleobj] ⊖ b.model[:laleobj] |> LalePipe
|>(a::PyObject,b::LaleOperator)     = a ⊖ b.model[:laleobj] |> LalePipe
|>(a::LaleOperator,b::PyObject)     = a.model[:laleobj] ⊖ b |> LalePipe
>>(a::LaleOperator,b::LaleOperator) = a.model[:laleobj] ⊖ b.model[:laleobj] |> LalePipe
>>(a::PyObject,b::LaleOperator)     = a ⊖ b.model[:laleobj] |> LalePipe
>>(a::LaleOperator,b::PyObject)     = a.model[:laleobj] ⊖ b |> LalePipe

⊕(a::PyObject,b::PyObject)          = _make_union(a,b)
+(a::LaleOperator,b::LaleOperator)  = a.model[:laleobj] ⊕ b.model[:laleobj] |> LalePipe
+(a::PyObject,b::LaleOperator)      = a ⊕ b.model[:laleobj] |> LalePipe
+(a::LaleOperator,b::PyObject)      = a.model[:laleobj] ⊕ b |> LalePipe

⨸(a::PyObject,b::PyObject)           = _make_union_nc(a,b)
(&)(a::LaleOperator,b::LaleOperator) = a.model[:laleobj] ⨸ b.model[:laleobj] |> LalePipe
(&)(a::PyObject,b::LaleOperator)     = a ⨸ b.model[:laleobj] |> LalePipe
(&)(a::LaleOperator,b::PyObject)     = a.model[:laleobj] ⨸ b |> LalePipe

⊗(a::PyObject,b::PyObject)          = _make_choice(a,b)
|(a::LaleOperator,b::LaleOperator)  = a.model[:laleobj] ⊗ b.model[:laleobj] |> LalePipe
|(a::PyObject,b::LaleOperator)      = a ⊗ b.model[:laleobj] |> LalePipe
|(a::LaleOperator,b::PyObject)      = a.model[:laleobj] ⊗ b |> LalePipe

mutable struct LalePipeOptimizer <: LaleOperator
   name::String
   model:: Dict{Symbol,Any}

   function LalePipeOptimizer(args=Dict{Symbol,Any}())
      default_args=Dict{Symbol,Any}(
            :lalepipe => PyNULL(),
            :name => "LalePipeOptimizer",
            :impl_args => Dict{Symbol,Any}(
                :cv        => 3,
                :max_evals => 10,
            )
      )
      cargs = nested_dict_merge(default_args, args)
      cargs[:name] = cargs[:name]*"_"*randstring(3)
      @assert cargs[:lalepipe] != PyNULL()
      new(cargs[:name],cargs)
   end
end

function LalePipeOptimizer(pipe::LalePipe, args::Dict)
   LalePipeOptimizer(Dict(:lalepipe => pipe.model[:laleobj], args...))
end

function LalePipeOptimizer(pipe::LalePipe; args...)
   LalePipeOptimizer(Dict(:lalepipe => pipe.model[:laleobj], :impl_args => Dict(pairs(args))))
end

function lalepipeoptimizers()
  opts = keys(optim_dict) |> collect |> x -> sort(x,lt=(x,y)->lowercase(x)<lowercase(y))
  println("syntax: LalePipeOptimizer(name::String, args::Dict=Dict())")
  println("where 'name can be one of:")
  println()
  [print(opt," ") for opt in opts]
  println()
  println()
  println("and 'args' are the corresponding learner's initial parameters.")
  println("Note: Consult Lale online help for more details about the auto_configure arguments.")
end

function fit!(lopt::LalePipeOptimizer, xx::DataFrame, y::Vector=Vector())
   Xpd = Pandas.DataFrame(xx).pyo
   Ypd = Pandas.DataFrame(y).pyo
   margs = lopt.model[:impl_args]
   pipe = lopt.model[:lalepipe]
   hyperopt = LALELIBS.Hyperopt(estimator=pipe; margs...)
   trained = hyperopt.fit(Xpd,Ypd)
   lopt.model[:trained] = trained
end

function transform!(lopt::LalePipeOptimizer, xx::DataFrame)
   Xpd = Pandas.DataFrame(xx).pyo
   trainedmodel = lopt.model[:trained]
   trainedmodel.predict(Xpd) |> Pandas.DataFrame |> DataFrame |> x -> x[:,1]
end

function fit(lopt::LalePipeOptimizer, xx::DataFrame, y::Vector=Vector())
   fit!(lopt,xx,y)
   loptcopy = deepcopy(lopt)
   return loptcopy
end

transform(lopt::LalePipeOptimizer, xx::DataFrame) = transform!(lopt,xx)
predict(lopt::LalePipeOptimizer, xx::DataFrame) = transform!(lopt,xx)

function visualize(lopt::LalePipeOptimizer)
   auto_trained = lopt.model[:trained]
   best_pipeline = auto_trained.get_pipeline()
   best_pipeline.visualize(ipython_display=false)
end

function visualize(lpipe::LalePipe)
   lpipe.model[:laleobj].visualize(ipython_display=false)
end

end


text:
# This file autogenerated through https://github.com/timholy/SnoopCompile.jl.  Precompile every day!
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Nettle.Decryptor, (String, Vector{UInt8},))
    precompile(Nettle.HMACState, (String, Vector{UInt8},))
    precompile(Nettle.CipherType, (Nettle.NettleCipher,))
    precompile(Nettle.HashType, (Nettle.NettleHash, Ptr{Cvoid},))
    precompile(Nettle.Encryptor, (String, Vector{UInt8},))
    precompile(Nettle.Hasher, (String,))
    precompile(Nettle.HMACState, (String, String,))
    precompile(Nettle.HashType, (String, UInt32, UInt32, UInt32, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},))
    precompile(Nettle.CipherType, (String, UInt32, UInt32, UInt32, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},))
    precompile(Nettle.encrypt!, (Nettle.Encryptor, Vector{UInt8}, Vector{UInt8},))
    precompile(Nettle.decrypt!, (Nettle.Decryptor, Vector{UInt8}, Vector{UInt8},))
    precompile(Nettle.digest!, (Nettle.HMACState,))
    precompile(Nettle.get_cipher_types, ())
    precompile(Nettle.get_hash_types, ())
    precompile(Nettle.digest!, (Nettle.Hasher,))
    precompile(Nettle.update!, (Nettle.HMACState, Vector{UInt8},))
    precompile(Nettle.update!, (Nettle.Hasher, String,))
    precompile(Nettle.decrypt, (String, Vector{UInt8}, Vector{UInt8},))
    precompile(Nettle.update!, (Nettle.HMACState, String,))
    precompile(Nettle.encrypt, (String, Vector{UInt8}, Vector{UInt8},))
    precompile(Nettle.hexdigest!, (Nettle.Hasher,))
end


text:
# Note: These methods need a more helpfull error message than a `NoMethodError`,
#       when the deprecation is removed
import Base.@deprecate
import Base.Operators: /
@deprecate (/)(x::(@compat Union{NAtype,Number}),A::AbstractDataArray)    x ./ A

#' @description
#'
#' Turn a DataArray into an Array. Raises an error if NA's are encountered.
#'
#' @param da::DataArray{T} DataArray that will be converted to an Array.
#'
#' @returns a::Array{T} Array containing values of `da`.
#'
#' @examples
#'
#' dv = @data [1, 2, 3, 4]
#' v = convert(Vector, dv)
#'
#' dm = @data [1 2; 3 4]
#' m = convert(Matrix, dm)
function array{T}(da::DataArray{T}) # -> Array{T}
    Base.depwarn(
        """
        array(da::DataArray{T}) is deprecated.
        Use convert(Array, da).
        """,
        :array
    )
    res = Array(T, size(da))
    for i in 1:length(da)
        if da.na[i]
            throw(NAException())
        else
            res[i] = da.data[i]
        end
    end
    return res
end

#' @description
#'
#' Turn a DataArray into an Array. Replace any NA's with the value
#' of second argument, `replacement`.
#'
#' @param da::DataArray{T} DataArray that will be converted to an Array.
#' @param replacement::T Value that will replace NA's in `da`.
#'
#' @returns a::Array{T} Array containing values of `da` plus replacements.
#'
#' @examples
#'
#' dv = @data [1, 2, NA, 4]
#' v = convert(Vector, dv, 3)
#'
#' dm = @data [1 2; NA 4]
#' m = convert(Matrix, dm, 3)
function array{T}(da::DataArray{T}, replacement::Any) # -> Array{T}
    Base.depwarn(
        """
        array(da::DataArray{T}, replacement::Any) is deprecated.
        Use convert(Array, da, replacement) instead.
        """,
        :array
    )
    res = Array(T, size(da))
    replacementT = convert(T, replacement)
    for i in 1:length(da)
        if da.na[i]
            res[i] = replacementT
        else
            res[i] = da.data[i]
        end
    end
    return res
end

# Turn a PooledDataArray into an Array. Fail on NA
function array{T, R}(da::PooledDataArray{T, R})
    Base.depwarn(
        """
        array(pda::PooledDataArray{T, R}) is deprecated.
        Use convert(Array, pda) instead.
        """,
        :array
    )
    n = length(da)
    res = Array(T, size(da))
    for i in 1:n
        if da.refs[i] == zero(R)
            throw(NAException())
        else
            res[i] = da.pool[da.refs[i]]
        end
    end
    return res
end

function array{T, R}(da::PooledDataArray{T, R}, replacement::T)
    Base.depwarn(
        """
        array(pda::PooledDataArray{T, R}, replacement::T) is deprecated.
        Use convert(Array, pda, replacement) instead.
        """,
        :array
    )
    n = length(da)
    res = Array(T, size(da))
    for i in 1:n
        if da.refs[i] == zero(R)
            res[i] = replacement
        else
            res[i] = da.pool[da.refs[i]]
        end
    end
    return res
end

@deprecate head(dv::AbstractDataVector) dv[1:min(6, end)]
@deprecate tail(dv::AbstractDataVector) dv[max(end-6, 1):end]

function rep{T <: Integer}(x::AbstractVector, lengths::AbstractVector{T})
    Base.depwarn(
        """
        rep{T <: Integer}(x::AbstractVector, lengths::AbstractVector{T}) is deprecated.
        """,
        :rep
    )
    if length(x) != length(lengths)
        throw(DimensionMismatch("vector lengths must match"))
    end
    res = similar(x, sum(lengths))
    i = 1
    for idx in 1:length(x)
        tmp = x[idx]
        for kdx in 1:lengths[idx]
            res[i] = tmp
            i += 1
        end
    end
    return res
end

@deprecate rep(x::AbstractVector, times::Integer, each::Integer = 1) Compat.repeat(x; inner=each, outer=times)
@deprecate rep(x::AbstractVector; times::Integer = 1, each::Integer = 1) Compat.repeat(x; inner=each, outer=times)
@deprecate rep(x::Any, times::Integer) Compat.repeat(x; inner=times)


text:
using NCDatasets
using JSON

function write_to_json(data, target_file)
	d = JSON.json(data)
	open(target_file, "w") do x
		write(x,d)
	end
end

function main()
    ds = Dataset("geolocation.nc")
    # geoJSON takes points in a [lon,lat] format
    sites = Dict(site => (ds["lon"][site],ds["lat"][site]) for site in ds["site"][:])
    write_to_json(sites,"../src/assets/geolocations.json")
end

main()

max_stars_repo_name: CliMA/Nimbus

max_issues_repo_name: CliMA/Nimbus

max_forks_repo_name: CliMA/Nimbus


text:
function convert(::Type{Float32}, val::Float16)
    ival::Uint32 = reinterpret(Uint16, val)
    sign::Uint32 = (ival & 0x8000) >> 15
    exp::Uint32  = (ival & 0x7c00) >> 10
    sig::Uint32  = (ival & 0x3ff) >> 0
    ret::Uint32

    if exp == 0
        if sig == 0
            sign = sign << 31
            ret = sign | exp | sig
        else
            n_bit = 1
            bit = 0x0200
            while (bit & sig) == 0
                n_bit = n_bit + 1
                bit = bit >> 1
            end
            sign = sign << 31
            exp = (-14 - n_bit + 127) << 23
            sig = ((sig & (~bit)) << n_bit) << (23 - 10)
            ret = sign | exp | sig
        end
    elseif exp == 0x1f
        if sig == 0  # Inf
            if sign == 0
                ret = 0x7f800000
            else
                ret = 0xff800000
            end
        else  # NaN
            ret = 0x7fc00000 | (sign<<31)
        end
    else
        sign = sign << 31
        exp  = (exp - 15 + 127) << 23
        sig  = sig << (23 - 10)
        ret = sign | exp | sig
    end
    return reinterpret(Float32, ret)
end

# Float32 -> Float16 algorithm from:
#   "Fast Half Float Conversion" by Jeroen van der Zijp
#   ftp://ftp.fox-toolkit.org/pub/fasthalffloatconversion.pdf

const basetable = Array(Uint16, 512)
const shifttable = Array(Uint8, 512)

for i = 0:255
    e = i - 127
    if e < -24  # Very small numbers map to zero
        basetable[i|0x000+1] = 0x0000
        basetable[i|0x100+1] = 0x8000
        shifttable[i|0x000+1] = 24
        shifttable[i|0x100+1] = 24
    elseif e < -14  # Small numbers map to denorms
        basetable[i|0x000+1] = (0x0400>>(-e-14))
        basetable[i|0x100+1] = (0x0400>>(-e-14)) | 0x8000
        shifttable[i|0x000+1] = -e-1
        shifttable[i|0x100+1] = -e-1
    elseif e <= 15  # Normal numbers just lose precision
        basetable[i|0x000+1] = ((e+15)<<10)
        basetable[i|0x100+1] = ((e+15)<<10) | 0x8000
        shifttable[i|0x000+1] = 13
        shifttable[i|0x100+1] = 13
    elseif e < 128  # Large numbers map to Infinity
        basetable[i|0x000+1] = 0x7C00
        basetable[i|0x100+1] = 0xFC00
        shifttable[i|0x000+1] = 24
        shifttable[i|0x100+1] = 24
    else  # Infinity and NaN's stay Infinity and NaN's
        basetable[i|0x000+1] = 0x7C00
        basetable[i|0x100+1] = 0xFC00
        shifttable[i|0x000+1] = 13
        shifttable[i|0x100+1] = 13
    end
end

function convert(::Type{Float16}, val::Float32)
    f = reinterpret(Uint32, val)
    i = (f >> 23) & 0x1ff + 1
    sh = shifttable[i]
    f &= 0x007fffff
    h::Uint16 = uint16(basetable[i] + (f >> sh))
    # round
    # NOTE: we maybe should ignore NaNs here, but the payload is
    # getting truncated anyway so "rounding" it might not matter
    nextbit = (f >> (sh-1)) & 1
    if nextbit != 0
        if h&1 == 1 ||  # round halfway to even
            (f & ((1<<(sh-1))-1)) != 0  # check lower bits
            h += 1
        end
    end
    reinterpret(Float16, h)
end

convert(::Type{Integer}, x::Float16) = convert(Int,float32(x))

round(x::Float16) = float16(round(float32(x)))
trunc(x::Float16) = float16(trunc(float32(x)))
floor(x::Float16) = float16(floor(float32(x)))
 ceil(x::Float16) = float16( ceil(float32(x)))

isnan(x::Float16)    = reinterpret(Uint16,x)&0x7fff  > 0x7c00
isinf(x::Float16)    = reinterpret(Uint16,x)&0x7fff == 0x7c00
isfinite(x::Float16) = reinterpret(Uint16,x)&0x7c00 != 0x7c00

function ==(x::Float16, y::Float16)
    ix = reinterpret(Uint16,x)
    iy = reinterpret(Uint16,y)
    if (ix|iy)&0x7fff > 0x7c00 #isnan(x) || isnan(y)
        return false
    end
    if (ix|iy)&0x7fff == 0x0000
        return true
    end
    return ix == iy
end

-(x::Float16) = reinterpret(Float16, reinterpret(Uint16,x) $ 0x8000)
abs(x::Float16) = reinterpret(Float16, reinterpret(Uint16,x) & 0x7fff)
for op in (:+,:-,:*,:/,:\,:^)
    @eval ($op)(a::Float16, b::Float16) = float16(($op)(float32(a), float32(b)))
end
for op in (:<,:<=,:isless)
    @eval ($op)(a::Float16, b::Float16) = ($op)(float32(a), float32(b))
end
for func in (:sin,:cos,:tan,:asin,:acos,:atan,:sinh,:cosh,:tanh,:asinh,:acosh,
             :atanh,:exp,:log,:sqrt)
    @eval begin
        $func(a::Float16) = float16($func(float32(a)))
        $func(a::Complex32) = complex32($func(complex64(a)))
    end
end
atan2(a::Float16, b::Float16) = float16(atan2(float32(a), float32(b)))
hypot(a::Float16, b::Float16) = float16(hypot(float32(a), float32(b)))
ldexp(a::Float16, b::Integer) = float16(ldexp(float32(a), b))
exponent(x::Float16) = exponent(float32(x))
^(x::Float16, y::Integer) = x^float16(y)

hash(x::Float16) = hash(reinterpret(Uint16, isnan(x) ? NaN16 : x))

max_stars_repo_name: tobig/julia

max_issues_repo_name: tobig/julia

max_forks_repo_name: tobig/julia


text:
module ModelicaPipe_GPSSM

using GP, Plots, CSV, DataFrames

function read_csv(file_suffix="1", x_i=4:6, y_i=2)
    # read csv and create datasets
    data_df = CSV.read("data/gp_data_$file_suffix.csv", DataFrame)
    input_df = CSV.read("data/input_data_$file_suffix.csv", DataFrame)
    validation_df = CSV.read("data/validation_data_$file_suffix.csv", DataFrame)

    # read csv and create datasets
    data_df_p = CSV.read("data/gp_data_$file_suffix.csv", DataFrame)
    y = data_df_p[:,y_i]
    X = convert(Array, data_df_p[:,x_i])

    # indicies
    #ix = [(data_df.t_samples[i], data_df.x_samples[i]) for i in 1:size(data_df,1)]
    ix = [(data_df.x_samples[i], data_df.t_samples[i]) for i in 1:size(data_df,1)]

    return X, y, ix, data_df, input_df, validation_df
end

function update_inputs(T_sim_t, V_sim_t, T_in, m_flow_in, Σs_in, nx)
    # initialize
    Xs_t = zeros(nx, 3)
    Σs_t = zeros(nx, 3)

    # update inputs for temperature
    Xs_t[1, 1] = T_in
    Xs_t[2:end, 1] = T_sim_t[2:end]
    # update inputs for Δ temperature
    Xs_t[1, 2] = T_sim_t[1] - T_in
    Xs_t[2:end, 2] = T_sim_t[2:end] - T_sim_t[1:end-1]
    # update inputs for m_flow
    Xs_t[:, 3] = m_flow_in*ones(size(Xs_t,1))

    # propagate forward Σs, set nx=1 from input
    Σs_t[1, :] = Σs_in
    Σs_t[2:end, 1] = V_sim_t[2:end]
    Σs_t[2:end, 2] = (V_sim_t[1:end-1] + V_sim_t[2:end])
    Σs_t[2:end, 3] = Σs_in[3]*ones(size(Σs_t,1)-1)

    return Xs_t, Σs_t
end

function update_states!(T_sim_t, V_sim_t, T_in, Δ_sim, σ_sim)
    # update T and V
    T_sim_t[1] = T_in
    T_sim_t[2:end] = T_sim_t[2:end] + Δ_sim[2:end]
    V_sim_t[1:end] = σ_sim
end

function forward_simulation(t_range, input_df, GP_pipe::GP.GP_obj, k, nx, calc_σ, nmc=0)
    # time series of inputs
    Σs_in = [5.0,5.0,1]
    #Σs_in = [0,0,0]

    T_in = Array(input_df[1:length(t_range)+1,2])
    T_in = hcat(T_in, T_in .+ Σs_in[2]*randn(length(t_range)+1,nmc))
    m_flow_in = Array(input_df[1:length(t_range)+1,3])
    m_flow_in = hcat(m_flow_in, m_flow_in .+ Σs_in[3]*randn(length(t_range)+1,nmc))

    # inputs at time t
    #Xs_t = zeros(Float64, nx, 3)
    #Σs_t = zeros(Float64, nx, 3)

    # new values
    Δ_sim = zeros(nx)
    σ_sim = zeros(nx)

    print("Simulating pipe [")
    # initialize
    T_sim_mc, V_sim_mc = zeros((nmc+1), nx, length(t_range)), zeros((nmc+1), nx, length(t_range))

    for n in 1:(nmc+1)
        print("[n = $n]")
        T_sim, V_sim = (273.15+40)*ones(nx), zeros(nx)
        T_sim_t, V_sim_t = (273.15+40.0)*ones(nx), zeros(nx)
        # run in t_range
        for t in t_range
            print(".")
            simulation_step!(T_sim_t, V_sim_t, GP_pipe, T_in[t,n], m_flow_in[t,n], Σs_in, nx)
            T_sim, V_sim = hcat(T_sim, T_sim_t), hcat(V_sim, V_sim_t)
        end
        # store results
        T_sim, V_sim = T_sim[:, 2:end] .- 273.15, V_sim[:, 2:end]
        T_sim_mc[n, :, :], V_sim_mc[n, :, :] = T_sim, V_sim
        # re-initialize
        T_sim, V_sim = (273.15+40)*ones(nx), zeros(nx)
        T_sim_t, V_sim_t = (273.15+40.0)*ones(nx), zeros(nx)
    end

    print("]")

    return T_sim_mc, V_sim_mc
end

function simulation_step!(T_sim_t, V_sim_t, GP_pipe, T_in, m_flow_in, Σs_in, nx)
    Xs_t, Σs_t = update_inputs(T_sim_t, V_sim_t, T_in, m_flow_in, Σs_in, nx)
    Δ_sim, σ_sim = GP.posterior(GP_pipe, Xs_t, Σs_t, k_pipe, true)
    update_states!(T_sim_t, V_sim_t, T_in, Δ_sim, σ_sim)
end

function plot_results_t(T_sim_mc, V_sim_mc, T_val, Δt, x, t_range, xlabel, ylabel, filename, save_file, plot_samples)
    t_plot = 1:Δt:Δt*(length(t_range))

    p = plot(xlabel=xlabel, ylabel=ylabel, legend=:topright)

    if plot_samples
        for i in 2:size(T_sim_mc, 1)
            T_sim = T_sim_mc[i, :, :]
            V_sim = V_sim_mc[i, :, :]
            μs_x = T_sim[x, t_range];
            σs_x = V_sim[x, t_range];
            ref = T_val[x, t_range];

            # plot mean of GP
            plot!(p, t_plot, μs_x, linewidth=0.1, linecolor=:blue, alpha=0.3, label="");
        end
    end

    T_sim = T_sim_mc[1, :, :]
    V_sim = V_sim_mc[1, :, :]
    μs_x = T_sim[x, t_range];
    σs_x = V_sim[x, t_range];
    ref = T_val[x, t_range];

    # plot mean of GP
    plot!(p, t_plot, μs_x, label="f(x)", linecolor=:green);

    # plot validation measurement
    plot!(p, t_plot, ref, label="y", linecolor=:red, linestyle=:dash);

    # plot posterior marginal std. dev.
    plot!(p, t_plot, [μs_x μs_x];
        linewidth=0.0,
        fillrange=[μs_x .- 3  .* σs_x μs_x .+ 3 .* σs_x],
        fillalpha=0.1,
        fillcolor=:blue,
        label="");

    if save_file
        savefig(p, filename)
    else
        display(p)
    end
end

function plot_results_x(T_sim_mc, V_sim_mc, T_val, Δx, x_range, t, xlabel, ylabel, filename, save_file)
    T_sim = T_sim_mc[1, :, :]
    V_sim = V_sim_mc[1, :, :]

    x_plot = 1:Δx:Δx*(length(x_range))
    μs_x = T_sim[x_range, t];
    σs_x = V_sim[x_range, t];
    ref = T_val[x_range, t];

    p = plot(xlabel=xlabel, ylabel=ylabel, legend=:topright)

    # plot mean of GP
    plot!(p, x_plot, μs_x, label="f(x)", linecolor=:green);

    # plot validation measurement
    plot!(p, x_plot, ref, label="y", linecolor=:red, linestyle=:dash);

    # plot posterior marginal std. dev.
    plot!(p, x_plot, [μs_x μs_x];
        linewidth=0.0,
        fillrange=[μs_x .- 3  .* σs_x μs_x .+ 3 .* σs_x],
        fillalpha=0.1,
        fillcolor=:blue,
        label="");

    if save_file
        savefig(p, filename)
    else
        display(p)
    end
end

# pipe specific kernel
function k_pipe(X₁, X₂, Θ::GP.hyperparameters, d=1:3)
    # kernels to multiply
    mul_1 = [2,3]; mul_2 = [1]

    # distance matrixes
    r_1 = GP.r_SE_mul(X₁[:, mul_1], X₂[:, mul_1], Θ.ℓ[Θ.ℓx[mul_1]])
    r_2 = GP.r_SE_mul(X₁[:, mul_2], X₂[:, mul_2], Θ.ℓ[Θ.ℓx[mul_2]])
    # squared exponential kernel
    k_1 =  GP.k_SE(r_1, Θ.σ[1])
    k_2 = GP.k_SE(r_2, Θ.σ[2])
    if d in mul_1
        return k_1
    elseif d in mul_2
        return k_2
    else
        return k_1 + k_2
    end
end

function POD(X,y,ix)
    return X, X, 0
end

function generate_GP(n_pre, n_large, m_sparse, Θ₀)

    # read csv files
    Xp, yp, ixp = read_csv("pre", 4:6, 2)
    Xl, yl, ixl = read_csv("large", 4:6, 2)

    # pre-optimize hyperparameters on the n dataset
    println("Pre-optimizing hyperparameters...")
    Θp = GP.optimize_hyperparameters(Xp, yp, k_pipe, Θ₀)

    # find a sparse approximation with m_sparse points
    println("Finding sparse approximation...")
    Xm, ym, ixm = GP.sparse_subset(Xl, yl, ixl, k_pipe, Θp, m_sparse)

    # re-optimize hyperparameters on the sparse dataset
    println("Post-optimizing hyperparameters...")
    Θm = GP.optimize_hyperparameters(Xm, ym, k_pipe, Θp)

    # cholesky factorization and inverse of K
    Lxxm, αm = GP.calcLxx(Xm, ym, k_pipe, Θm)

    return GP.GP_obj(Xm, Lxxm, αm, Θm), ixl, ixm
end

end

max_stars_repo_name: smnssn/GPDH

max_issues_repo_name: smnssn/GPDH

max_forks_repo_name: smnssn/GPDH


text:
export fixed_field, subfields

# TODO !!!: Claus says this is wrong :(
function _subfield_basis(K, elt)
  n = degree(K)
  bas = elem_type(K)[one(K)]
  phase = 1
  local B::FakeFmpqMat

  for e = elt
    if phase == 2
      C = basis_mat([e])
      fl, _ = can_solve(matrix(FlintQQ, B.num), matrix(FlintQQ, C.num), side = :left)
      fl && continue
    end
    df = n-1
    f = one(K)
    for i=1:df
      mul!(f, f, e)
      if phase == 2
        C = matrix(FlintQQ, basis_mat(elem_type(K)[f]).num)
        reduce_mod!(C, matrix(FlintQQ, B.num))
        fl = iszero(C)
        fl && break
      end
      b = elem_type(K)[e*x for x in bas]
      append!(bas, b)
      if length(bas) >= n
        B = basis_mat(bas)
        hnf!(B)
        rk = nrows(B) - n + 1
        while iszero_row(B, rk)
          rk += 1
        end
        B = sub(B, rk:nrows(B), 1:n)
        phase = 2
        bas = elem_type(K)[ elem_from_mat_row(K, B.num, i, B.den) for i = 1:nrows(B) ]
      end
    end
  end

  if length(bas) >= n
    B = basis_mat(bas)
    hnf!(B)
    rk = nrows(B) - n + 1
    # This is wrong. But I need to see an actual error
    if iszero_row(B.num, rk)
      error("data does not define an order: dimension to small")
    end
    B = sub(B, rk:nrows(B), 1:n)
    bas = elem_type(K)[ elem_from_mat_row(K, B.num, i, B.den) for i = 1:nrows(B) ]
  end

  return bas
end

function _improve_subfield_basis(K, bas)
  # First compute the maximal order of <bas> by intersecting and saturating
  # Then B_Ok = N * B_LLL_OK
  # Then B' defined as lllN * B_LLL_OK will hopefully be small
  OK = maximal_order(K)
  OKbmatinv = basis_mat_inv(OK, copy = false)
  basinOK = bas * matrix(FlintQQ, OKbmatinv.num) * fmpq(1, OKbmatinv.den)
  deno = fmpz(1)
  for i in 1:nrows(basinOK)
    for j in 1:ncols(basinOK)
      deno = lcm(deno, denominator(basinOK[i, j]))
    end
  end
  S = saturate(matrix(FlintZZ, basinOK * deno))
  SS = S * basis_mat(OK, copy = false)
  lllOK = lll(OK)
  N = (SS * basis_mat_inv(lllOK)).num
  lllN = lll(N)
  maybesmaller = lllN * basis_mat(lllOK)
  return maybesmaller
end

function _subfield_primitive_element_from_basis(K, elt)
  if length(elt) == 0
    return gen(K)
  end

  s = zero(K)

  # First check basis elements
  for i in 1:length(elt)
    if degree(minpoly(elt[i])) == length(elt)
      return elt[i]
    end
  end

  while true
    rand!(s, elt, 0:1)
    if (degree(minpoly(s)) == length(elt))
      return s
    end
  end
end

#returns minimal subfield of K containing A
function subfield(K::S, elt::Array{T, 1}; isbasis::Bool = false) where {S <: Union{AnticNumberField, Hecke.NfRel}, T <: Union{nf_elem, Hecke.NfRelElem}}
  if length(elt) == 1
    return _subfield_from_primitive_element(K, elt[1])
  end

  if isbasis
    s = _subfield_primitive_element_from_basis(K, elt)
  else
    bas = _subfield_basis(K, elt)
    s = _subfield_primitive_element_from_basis(K, bas)
  end

  return _subfield_from_primitive_element(K, s)
end

function _subfield_from_primitive_element(K, s)
  f = minpoly(s)
  L,_ = NumberField(f)
  return L, hom(L, K, s)
end

################################################################################
#
#  Fixed field
#
################################################################################

@doc Markdown.doc"""
    fixed_field(K::NumberField, sigma::NfToNfMor;
                                simplify::Bool = true) -> NumberField, NfToNfMor

Given a number field $K$ and an automorphisms $\sigma$ of $K$, this function
returns the fixed field of $\sigma$ as a pair $(L, i)$ consisting of a number
field $L$ and an embedding of $L$ into $K$.

By default, the function tries to find a small defining polynomial of $L$. This
can be disabled by setting `simplify = false`.
"""
function fixed_field(K::S, sigma::T; simplify::Bool = true) where {S <: Union{AnticNumberField, NfRel}, T <: Union{NfToNfMor, NfRelToNfRelMor}}
  return fixed_field(K, T[sigma], simplify = simplify)
end

@doc Markdown.doc"""
    fixed_field(K::NumberField, A::Vector{NfToNfMor}) -> NumberField, NfToNfMor

Given a number field $K$ and a set $A$ of automorphisms of $K$, this function
returns the fixed field of $A$ as a pair $(L, i)$ consisting of a number
field $L$ and an embedding of $L$ into $K$.

By default, the function tries to find a small defining polynomial of $L$. This
can be disabled by setting `simplify = false`.
"""
function fixed_field(K::S, A::Array{T, 1}; simplify::Bool = true) where {S <: Union{AnticNumberField, NfRel}, T <: Union{NfToNfMor, NfRelToNfRelMor}}

  autos = A

  # Everything is fixed by nothing :)
  if length(autos) == 0
    return K, id_hom(K)
  end

  F = base_field(K)
  a = gen(K)
  n = degree(K)
  ar_mat = Vector{dense_matrix_type(elem_type(F))}()
  v = Vector{elem_type(K)}(undef, n)
  for i in 1:length(autos)
    domain(autos[i]) !== codomain(autos[i]) && throw(error("Maps must be automorphisms"))
    domain(autos[i]) !== K && throw(error("Maps must be automorphisms of K"))
    o = one(K)
    # Compute the image of the basis 1,a,...,a^(n - 1) under autos[i] and write
    # the coordinates in a matrix. This is the matrix of autos[i] with respect
    # to 1,a,...a^(n - 1).
    as = autos[i](a)
    if a == as
      continue
    end
    v[1] = o
    for j in 2:n
      o = o * as
      v[j] = o
    end

    bm = basis_mat(v)

    if S === AnticNumberField
      # We have to be a bit careful (clever) since in the absolute case the
      # basis matrix is a FakeFmpqMat

      m = matrix(FlintQQ, bm.num)
      for j in 1:n
        m[j, j] = m[j, j] - bm.den # This is autos[i] - identity
      end
    else
      # In the generic case just subtract the identity
      m = bm - identity_matrix(F, degree(K))
    end

    push!(ar_mat, m)
  end

  if length(ar_mat) == 0
    return K, id_hom(K)
  else
    bigmatrix = hcat(ar_mat)
    k, Ker = kernel(bigmatrix, side = :left)
    bas = Vector{elem_type(K)}(undef, k)
    if S === AnticNumberField
      # Try to find a small basis for absolute simple number fields
      if simplify
        KasFMat = _improve_subfield_basis(K, Ker)
        for i in 1:k
          bas[i] = elem_from_mat_row(K, KasFMat.num, i, KasFMat.den)
        end
      else
        KasFMat = FakeFmpqMat(Ker)
        Ksat = saturate(KasFMat.num)
        Ksat = lll(Ksat)
        onee = one(fmpz)
        for i in 1:k
          bas[i] = elem_from_mat_row(K, Ksat, i, onee)
        end
      end
    else
      for i in 1:k
        bas[i] = elem_from_mat_row(K, Ker, i)
      end
    end

    return subfield(K, bas, isbasis = true)
  end
end

################################################################################
#
#  Principal subfields
#
################################################################################

function _principal_subfields_basis(K::AnticNumberField)
  f = K.pol
  Kx, x = PolynomialRing(K, "x", cached = false)
  n = degree(K)
  #f in Kx
  #fk = Kx([coeff(f,i) for i in 0:n])
  #determine roots
  rts = roots(f, K)
  ar_lin_fac = elem_type(Kx)[x - root for root in rts]
  fK = change_ring(f, Kx)
  ##divide by roots
  for lin_fac in ar_lin_fac
    fK = div(fK, lin_fac)#divaxa
  end

  fac = factor(fK)

  for (g, e) in fac
    push!(ar_lin_fac, g)
  end

  factor_ar = ar_lin_fac
  k = base_field(K)
  principal_subfields_ar = dense_matrix_type(elem_type(k))[]

  # #compute kernel of (phi - id)
  for fi in factor_ar
    M = zero_matrix(k, n, n * degree(fi))
    im_ar = elem_type(Kx)[(mod(x^l,fi)-gen(K)^l) for l in 0:n-1]
    for j in 1:n
      for l in 0:degree(fi)-1
        t = basis_mat([coeff(im_ar[j], l)])
        for m in 1:n
          M[j, l * n + m] = (t.num[1, m])//t.den
        end
      end
    end
    nu, ker = kernel(M, side = :left)

    # This might be expensive for bigger fields?
    ker_rref = matrix(FlintQQ, lll(saturate(FakeFmpqMat(rref(ker)[2]).num)))

    # TODO: Remove this madness once we switched to Nemo >=0.13.2
    if ker_rref in [ b for b in principal_subfields_ar if nrows(b) == nrows(ker_rref)]
      continue
    else
      push!(principal_subfields_ar, ker_rref)
    end
  end
  return principal_subfields_ar
end

function principal_subfields(K::T) where {T}
  ba = _principal_subfields_basis(K)
  elts = Vector{Vector{nf_elem}}(undef, length(ba))
  for i in 1:length(ba)
    baf = FakeFmpqMat(ba[i])
    elts[i] = elem_type(K)[elem_from_mat_row(K, baf.num, j, baf.den) for j in 1:nrows(ba[i])]
  end
  return Tuple{T, morphism_type(T)}[ subfield(K, elts[i], isbasis = true) for i in 1:length(elts)]
end

# Computes the minpoly of a over M if k(a)=K/M/k
# TODO: Improve this.
function _get_sfpoly(Kx, M)
  K = base_ring(Kx)
  k = base_field(K)
  n = degree(K)
  @assert k === base_ring(M)
  M = transpose(M)
  rk = rank(M)
  if rk == 1
      return Kx([K(coeff(minpoly(gen(K)),i)) for i in 0:degree(minpoly(gen(K)))])
  end
  my = div(n, rk)
  ar_basis = Vector{elem_type(K)}(undef, rk)
  for i in 1:rk
    elem_basis = zero(K)
    for j in 1:n
      elem_basis += M[j,i] * gen(K)^(j-1)
    end
    ar_basis[i] = elem_basis
  end
  null_mat = zero_matrix(k, n, 0)
  for i in 1:length(ar_basis)
    null_mat = hcat(null_mat,representation_matrix(ar_basis[i])'[:,1:my])
  end
  temp = zero_matrix(k,n,1)
  temp[my+1,1] = one(k)
  null_mat = hcat(null_mat,temp)
  mat_poly = nullspace(null_mat)[2]
  ar_coeff = Array{elem_type(K),1}(undef,my)
  for i = 1:my
    indx = 1
    temp_coeff = K(0)
    for j = i:my:n
      temp_coeff += mat_poly[j,1] * ar_basis[indx]
      indx += 1
    end
    ar_coeff[i] = temp_coeff
  end
  #w.l. we cannot assume basis ordered in any way
  push!(ar_coeff,K(mat_poly[n+1,1]))
  if !isone(mat_poly[n+1,1])
      error("LC of minpoly not 1")
  end
  return Kx(ar_coeff)
end

# TODO:
# The algorithm works with subfields represented as vector spaces and it needs
# to compute intersections and test containment
# - Improve this by exploiting that everything should be in rref (use reduce_mod)
# - Maybe also cache the pivots

# Computes the intersection of subfields A,B of K/k, represented as k-VS
# TODO (easy): Get rid of the transpose :)
function intersect_spaces(A::Hecke.AbstractAlgebra.Generic.MatElem, B::Hecke.AbstractAlgebra.Generic.MatElem)
  A = transpose(A)
  B = transpose(B)
  M = nullspace(hcat(A,-B))[2][1:size(A)[2],:]
  intersect_mat = zero_matrix(base_ring(A), nrows(A), ncols(M))
  for i in 1:ncols(intersect_mat)
    intersect_mat[:,i] = A * M[:,i]
  end
  return transpose(intersect_mat)
end

# Computes the intersection of subfields A = [a1,...,an] of K/k, represented as k-VS
function intersect_spaces(A::Vector{T}) where T
  if length(A) < 1
    throw(error("Number of spaces must be non-zero"))
  elseif length(A) == 1
    return A[1]
  elseif length(A) == 2
    return intersect_spaces(A[1],A[2])
  else
    intersection_temp = intersect_spaces(A[1],A[2])
    for i in 3:length(A)
      intersection_temp = intersect_spaces(intersection_temp,A[i])
    end
    return intersection_temp
  end
end

# Returns true if A is subspace of B,otherwise false, for A,B k-VS
# TODO (easy): Make this great again :)
function issubspace(A::MatElem, B::MatElem, proper_subspace::Bool = false)     #or cmpr remark7
  intersectAB = intersect_spaces(A, B)
  Bol = rank(intersectAB) == rank(A)
  if proper_subspace
    return Bol && (rank(intersectAB) < rank(B))
  else
    return Bol
  end
end

function generating_subfields(S, len::Int64 = -1)
  # I don't know (yet) why it is necessary, but we have to do it
  if length(S) == 1
    return S
  end
  ar_2delete = Bool[false for i in 1:length(S)]
  for i in 1:length(S)
    if nrows(S[i]) >= len
      ar_intersection = typeof(S)()
      for j in 1:length(S)
        if issubspace(S[i],S[j],true)
          push!(ar_intersection, S[j])
        end
      end
      #K always principal_subfield
      if length(ar_intersection) == 0
        ar_2delete[i] = true
      else
        intersection = intersect_spaces(ar_intersection)
        if issubspace(intersection, S[i]) && issubspace(S[i], intersection)
          ar_2delete[i] == true
        end
      end
    else
      ar_2delete[i] = true
    end
  end
  generating_subfields_ar = S
  indx = length(ar_2delete)
  while indx > 0
    if ar_2delete[indx]
      deleteat!(generating_subfields_ar, indx)
    end
    indx -= 1
  end
  return generating_subfields_ar
end

function allSubfields(K, S::Vector{T}, len::Int64 = -1) where {T}
    Kx, _  = PolynomialRing(K, "x", cached = false)
    if length(S) == 0
        return S
    end
    #compute associated tuple of K
    e = Int[0 for i in 1:length(S)]
    K_mat = identity_matrix(base_field(K), degree(K))
    sf_ar = elem_type(Kx)[_get_sfpoly(Kx, M) for M in S]
    ListSubfields = T[K_mat]
    nextSubfields(ListSubfields, Kx, S ,K_mat,e,0,sf_ar, len)
    return ListSubfields
end

#Hilfsfunktion for allSubfields
function nextSubfields(ListSubfields, Kx, S::Vector{T}, L::T, e::Array{Int64,1}, s::Int64, sf_ar, len::Int64) where {T}
  i = s + 1
  while i <= length(S)
    if e[i] == 0
      M = intersect_spaces(L,S[i])
      #check len
      if rank(M) >= len
        #constr associated tuple of M
        ee = Int[0 for j in 1:length(S)]
        for j in 1:length(S)
          #    if issubspace(M,S[j])
          if iszero(mod(_get_sfpoly(Kx, M),sf_ar[j]))
            ee[j] = 1
          end
        end
        #check whether to append M
        # This should be checked when creating the e'
        appnd = true
        for j in 1:i-1
          if ee[j] > e[j]
            appnd = false
            break
          end
        end
        if appnd
          push!(ListSubfields,M)
          hi = i
          nextSubfields(ListSubfields, Kx, S,M,ee,hi,sf_ar, len)
        end
      end
    end
    i += 1
  end
  #check length
  if len != -1
    indx = length(ListSubfields)
    while indx > 0
      if rank(ListSubfields[indx]) != len
        deleteat!(ListSubfields,indx)
      end
      indx -= 1
    end
  end
end

function subfields(K::T; degree::Int64 = -1) where {T <: Union{AnticNumberField, Hecke.NfRel}}
  n = Hecke.degree(K) # I want to keep the degree keyword
  k = base_field(K)
  #K = k[x]/f
  # K no generating subfield

  # TODO (medium)
  # I don't know why we have to do this.
  # This needs to be fixed properly
  if n == 1
    return Tuple{T, morphism_type(T)}[(K, id_hom(K))]
  end

  if degree == n
    return Tuple{T, morphism_type(T)}[(K, id_hom(K))]
  end
  princ_subfields = _principal_subfields_basis(K)
  gg = generating_subfields(princ_subfields)
  sf_asmat_ar = allSubfields(K, gg, degree)
  #compute embedding
  Res = Vector{Tuple{typeof(K), morphism_type(typeof(K))}}()
  #get minimal polynomial of primitive elem k(pe) = M, over k
  for sf_mat in sf_asmat_ar
    #interpret column vectors as field elems
    #sf_mat = transpose(sf_mat)
    sf_mat_f = FakeFmpqMat(sf_mat)
    basis_ar = Array{elem_type(K),1}(undef, nrows(sf_mat_f))
    for i in 1:nrows(sf_mat_f)
      basis_ar[i] = elem_from_mat_row(K, sf_mat_f.num, i, sf_mat_f.den)
      #field_elem = K(0)
      #for j in 1:n
      #  field_elem += sf_mat[j,i] * gen(K)^(j-1)
      #end
      #basis_ar[i] = field_elem
    end
    push!(Res, subfield(K, basis_ar, isbasis = true))
  end
  return Res
end

# TODO: Write a dedicated function for the normal case and use the subgroup functions