text:
# Example (Node)
# -------
# <node id="61327001" visible="true" version="1" changeset="64956"
#       timestamp="2007-10-08T23:40:44Z" user="MassGIS Import"
#       uid="15750" lat="42.3627260" lon="-71.0936790">
#   <tag k="attribution"
#        v="Office of Geographic and Environmental Information (MassGIS)"/>
#   <tag k="created_by" v="JOSM"/>
#   <tag k="source" v="massgis_import_v0.1_20071008165629"/>
# </node>

type OSMData
    nodes::Vector{Node}
    ways::Vector{Way}
    relations::Vector{Relation}
    node_tags::Set{UTF8String}
    way_tags::Set{UTF8String}
    relation_tags::Set{UTF8String}
end
OSMData() = OSMData(Vector{Node}(), Vector{Way}(), Vector{Relation}(),
                    Set{UTF8String}(), Set{UTF8String}(), Set{UTF8String}())

type DataHandle
    element::Symbol
    osm::OSMData
    node::Node # initially undefined
    way::Way # initially undefined
    relation::Relation # initially undefined

    DataHandle() = new(:None, OSMData())
end

function parseElement(handler::LibExpat.XPStreamHandler,
                      name::AbstractString,
                      attr::Dict{AbstractString,AbstractString})
    data = handler.data::DataHandle
    if name == "node"
        data.element = :Node
        data.node = Node(parse(Int, attr["id"]),
                         (float(attr["lon"]), float(attr["lat"])))
    elseif name == "way"
        data.element = :Way
        data.way = Way(parse(Int, attr["id"]))
    elseif name == "relation"
        data.element = :Relation
        data.relation = Relation(parse(Int, attr["id"]))
    elseif name == "tag"
        k = attr["k"]; v = attr["v"]
        if data.element == :Node
            data_tags = tags(data.node)
            push!(data.osm.node_tags, k)
        elseif data.element == :Way
            data_tags = tags(data.way)
            push!(data.osm.way_tags, k)
        elseif data.element == :Relation
            data_tags = tags(data.relation)
            push!(data.osm.relation_tags, k)
        end
        data_tags[k] = v
    elseif name == "nd"
        push!(data.way.nodes, parse(Int, attr["ref"]))
    elseif name == "member"
        push!(data.relation.members, attr)
    end
end

function collectElement(handler::LibExpat.XPStreamHandler, name::AbstractString)
    if name == "node"
        push!(handler.data.osm.nodes, handler.data.node)
        handler.data.element = :None
    elseif name == "way"
        push!(handler.data.osm.ways, handler.data.way)
        handler.data.element = :None
    elseif name == "relation"
        push!(handler.data.osm.relations, handler.data.relation)
        handler.data.element = :None
    end
end

function parseOSM(filename::AbstractString; args...)
    callbacks = LibExpat.XPCallbacks()
    callbacks.start_element = parseElement
    callbacks.end_element = collectElement
    data = DataHandle()
    LibExpat.parsefile(filename, callbacks, data=data; args...)
    data.osm::OSMData
end

function parseNode(handler::LibExpat.XPStreamHandler,
                   name::AbstractString,
                   attr::Dict{AbstractString,AbstractString})
    data = handler.data
    if name == "node"
        data.element = :Node
        data.curr = Node(parse(Int, attr["id"]),
                         (float(attr["lon"]), float(attr["lat"])),
                         Dict{UTF8String,UTF8String}())
    elseif name == "tag" && data.element == :Node
        data.curr.tags[attr["k"]] = attr["v"]
    end
end

function collectNode(handler::LibExpat.XPStreamHandler, name::AbstractString)
    if name == "node"
        handler.data.element = :None
        push!(handler.data.nodes, handler.data.curr)
    end
end

function parseNodes(filename::AbstractString; args...)
    callbacks = LibExpat.XPCallbacks()
    callbacks.start_element = parseNode
    callbacks.end_element = collectNode
    nodehandle = NodeHandle()
    LibExpat.parsefile(filename, callbacks, data=nodehandle; args...)
    nodehandle.nodes::Vector{Node}
end

# Example (Way)
# -------
# <way id="5090250" visible="true" timestamp="2009-01-19T19:07:25Z" version="8"
#      changeset="816806" user="Blumpsy" uid="64226">
#     <nd ref="822403"/>
#     <nd ref="21533912"/>
#     ...
#     <nd ref="135791608"/>
#     <nd ref="823771"/>
#     <tag k="highway" v="residential"/>
#     <tag k="name" v="Clipstone Street"/>
#     <tag k="oneway" v="yes"/>
#   </way>

type WayHandle
    element::Symbol
    ways::Vector{Way}
    curr::Way # initially undefined
    WayHandle(element::Symbol, ways::Vector{Way}) = new(element, ways)
end

function parseWay(handler::LibExpat.XPStreamHandler,
                  name::AbstractString,
                  attr::Dict{AbstractString,AbstractString})
    if name == "way"
        handler.data.element = :Way
        handler.data.curr = Way(parse(Int, attr["id"]),
                                Vector{Int}(),
                                Dict{UTF8String,UTF8String}())
    elseif handler.data.element == :Way
        if name == "tag"
            handler.data.curr.tags[attr["k"]] = attr["v"]
        elseif name == "nd"
            push!(handler.data.curr.nodes, parse(Int, attr["ref"]))
        end
    end
end

function collectWay(handler::LibExpat.XPStreamHandler, name::AbstractString)
    if name == "way"
        push!(handler.data.ways, handler.data.curr)
        handler.data.element = :None
    end
end

function parseWays(filename::AbstractString; args...)
    wayhandle = WayHandle(:None, Vector{Way}())
    callbacks = LibExpat.XPCallbacks()
    callbacks.start_element = parseWay
    callbacks.end_element = collectWay
    LibExpat.parsefile(filename, callbacks, data=wayhandle; args...)
    wayhandle.ways::Vector{Way}
end

# Example (Relation)
# <relation id="1">
#   <tag k="type" v="multipolygon" />
#   <member type="way" id="1" role="outer" />
#   <member type="way" id="2" role="outer" />
#   ...
#   <member type="way" id="19" role="inner" />
#   <member type="way" id="20" role="outer" />
# </relation>

type RelationHandle
    element::Symbol
    relations::Vector{Relation}
    curr::Relation # initially undefined
    RelationHandle(e::Symbol, relations::Vector{Relation}) = new(e, relations)
end

function parseRelation(handler::LibExpat.XPStreamHandler,
                       name::AbstractString,
                       attr::Dict{AbstractString,AbstractString})
    if name == "relation"
        handler.data.element = :Relation
        handler.data.curr = Relation(parse(Int, attr["id"]),
                                     Vector{Dict{UTF8String,UTF8String}}(),
                                     Dict{UTF8String,UTF8String}())
    elseif handler.data.element == :Relation
        if name == "tag"
            handler.data.curr.tags[attr["k"]] = attr["v"]
        elseif name == "member"
            push!(handler.data.curr.members, attr)
        end
    end
end

function collectRelation(handler::LibExpat.XPStreamHandler, name::AbstractString)
    if name == "relation"
        push!(handler.data.relations, handler.data.curr)
        handler.data.element = :None
    end
end

function parseRelations(filename::AbstractString; args...)
    relationhandle = RelationHandle(:None, Vector{Relation}())
    callbacks = LibExpat.XPCallbacks()
    callbacks.start_element = parseRelation
    callbacks.end_element = collectRelation
    LibExpat.parsefile(filename, callbacks, data=relationhandle; args...)
    relationhandle.relations::Vector{Relation}
end


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

using DistancePlan
using StaticArrays
import Makie

scene = Makie.Scene()

DistancePlan.draw!(Sphere(SVector{2}([0. 0.]), 1.))

Makie.save("unit_circle.png", scene, resolution=(1000, 1000))

ext: jl

max_stars_repo_name: krzentner/DistancePlan

max_issues_repo_name: krzentner/DistancePlan

max_forks_repo_name: krzentner/DistancePlan


text:
################################################################################
#
#  Infinite places for number fields
#
#
#  (C) 2017 Tommy Hofmann
#
################################################################################

export is_complex, is_positive, is_totally_positive, signs, sign, real_places,
       complex_places, infinite_places, infinite_place

################################################################################
#
#  Equality and hashing
#
################################################################################

function Base.:(==)(P::InfPlc, Q::InfPlc)
  P.K != Q.K && error("Places are of different number fields")
  return P.i == Q.i
end

function Base.hash(P::InfPlc, h::UInt)
  return xor(Base.hash(P.K, h), Base.hash(P.i, h))
end

################################################################################
#
#  Printing
#
################################################################################

function Base.show(io::IO, P::InfPlc)
  if P.isreal
    print(io, "Real ")
  else
    print(io, "Complex ")
  end

  print(io, "place of \n$(P.K)\n")
  print(io, "corresponding to the root $(P.r)")
end

################################################################################
#
#  Predicates
#
################################################################################

function Base.isreal(P::InfPlc)
  return P.isreal
end

function is_complex(P::InfPlc)
  return !isreal(P)
end

################################################################################
#
#  Construction
#
################################################################################

@doc Markdown.doc"""
    infinite_place(K::AnticNumberField, i::Int) -> InfPlc

This function returns the infinite place of $K$ corresponding to the root
`conjugates_arb(a)[i]`, where `a` is the primitive element of $K$.
"""
function infinite_place(K::AnticNumberField, i::Int)
  !(1 <= i <= degree(K)) && error("Index must be between 1 and $(degree(K))")
  return InfPlc(K, i)
end

function infinite_places(K::AnticNumberField)
  _res = get_attribute(K, :infinite_places)
  if _res !== nothing
    return _res::Vector{InfPlc}
  end
  r1, r2 = signature(K)
  plcs = InfPlc[ InfPlc(K, i) for i in 1:(r1 + r2)]
  set_attribute!(K, :infinite_places => plcs)
  return plcs
end

@doc Markdown.doc"""
    real_places(K::AnticNumberField) -> Vector{InfPlc}

This function returns all infinite real places of $K$.
"""
function real_places(K::AnticNumberField)
  r1, r2 = signature(K)
  return [ InfPlc(K, i) for i in 1:r1]
end

@doc Markdown.doc"""
    complex_places(K::AnticNumberField) -> Vector{InfPlc}

This function returns all infinite complex places of $K$.
"""
function complex_places(K::AnticNumberField)
  r1, r2 = signature(K)
  return [ InfPlc(K, i) for i in (r1 + 1):r1 + r2]
end

################################################################################
#
#  Signs at places
#
################################################################################

@doc Markdown.doc"""
    signs(a::nf_elem)          -> Dict{InfPlc, Int}
    signs(a::FacElem{nf_elem}) -> Dict{InfPlc, Int}

This function returns a dictionary of the signs of $a$ at all infinite places
of the ambient number field. The keys are infinite places of the ambient
number field. The value is $1$ if the sign is positive and $-1$ if the sign
is negative.
"""
function signs(a::Union{nf_elem, FacElem{nf_elem, AnticNumberField}})
  K = _base_ring(a)
  r1, r2 = signature(K)
  D = Dict{InfPlc, Int}()
  s = _signs(a)
  for i in 1:r1
    P = InfPlc(K, i)
    D[P] = s[i]
  end
  return D
end

@doc Markdown.doc"""
    sign(a::nf_elem, P::InfPlc)          -> Int
    sign(a::FacElem{nf_elem}, P::InfPlc) -> Int

This function returns the sign of $a$ at the place $P$. The value is $1$ if
the sign is positive and $-1$ if the sign is negative.

"""
function sign(a::Union{nf_elem, FacElem{nf_elem, AnticNumberField}}, P::InfPlc)
  !isreal(P) && error("Place must be real")
  if a isa nf_elem && iszero(a)
    return 0
  end
  return signs(a, [P])[P]
end

@doc Markdown.doc"""
    signs(a::nf_elem, l::Vector{InfPlc})          -> Dict{InfPlc, Int}
    signs(a::FacElem{nf_elem}, l::Vector{InfPlc}) -> Dict{InfPlc, Int}

This function returns a dictionary of the signs of $a$ at places in $l$. The
keys are the elements of $l$. The value is $1$ if the sign is positive and
$-1$ if the sign is negative. The result will contain as many signs as there
are real places contained in $l$.
"""
function signs(a::Union{nf_elem, FacElem{nf_elem, AnticNumberField}}, l::Vector{InfPlc})
  K = _base_ring(a)
  r1, r2 = signature(K)
  D = Dict{InfPlc, Int}()
  s = _signs(a)
  for i in 1:r1
    P = InfPlc(K, i)
    if P in l
      D[P] = s[i]
    end
  end
  return D
end

# extend functionality to elements of orders

function signs(a::NfOrdElem, args...)
  return signs(a.elem_in_nf, args...)
end

function sign(a::NfOrdElem, args...)
  return sign(a.elem_in_nf, args...)
end

################################################################################
#
#  Positivity
#
################################################################################

@doc Markdown.doc"""
    is_positive(a::nf_elem, P::InfPlc)          -> Bool
    is_positive(a::FacElem{nf_elem}, P::InfPlc) -> Bool

Returns whether the element $a$ is positive at the embedding corresponding to
$P$. The place $P$ must be real.
"""
function is_positive(a::Union{nf_elem, FacElem{nf_elem, AnticNumberField}}, P::InfPlc)
  !isreal(P) && error("Place must be real")
  return sign(a, P) > 0
end

@doc Markdown.doc"""
    is_negative(a::nf_elem, P::InfPlc)          -> Bool
    is_negative(a::FacElem{nf_elem}, P::InfPlc) -> Bool

Returns whether the element $a$ is negative at the embedding corresponding to
$P$. The place $P$ must be real.
"""
function is_negative(a::Union{nf_elem, FacElem{nf_elem, AnticNumberField}}, P::InfPlc)
  !isreal(P) && error("Place must be real")
  return sign(a, P) < 0
end

@doc Markdown.doc"""
    is_positive(a::nf_elem, l::Vector{InfPlc})          -> Bool
    is_positive(a::FacElem{nf_elem}, l::Vector{InfPlc}) -> Bool

Returns whether the element $a$ is positive at the embeddings corresponding to
the real places of $l$.
"""
function is_positive(a::Union{nf_elem, FacElem{nf_elem, AnticNumberField}}, l::Vector{InfPlc})
  return all(x -> is_positive(a, x), (y for y in l if isreal(y)))
end

@doc Markdown.doc"""
    is_totally_positive(a::nf_elem)          -> Bool
    is_totally_positive(a::FacElem{nf_elem}) -> Bool

Returns whether the element $a$ is totally positive, that is, whether it is
positive at all places of the ambient number field.
"""
function is_totally_positive(a::Union{nf_elem, FacElem{nf_elem, AnticNumberField}})
  K = _base_ring(a)
  return is_positive(a, real_places(K))
end

# extend functionality to elements of orders

function is_positive(a::NfOrdElem, args...)
  return is_positive(a.elem_in_nf, args...)
end

function is_totally_positive(a::NfOrdElem, args...)
  return is_totally_positive(a.elem_in_nf, args...)
end

is_totally_positive(x::fmpq) = x > 0

################################################################################
#
#  Action of a morphism on a infinite place
#
################################################################################

#The action of f on P is defined as f(P) = P\circ f^{-1} and not P\circ f
#In this way, (f\circ g)(P)= f(g(P)), otherwise it would fail.

@doc Markdown.doc"""
    induce_image(m::NfToNfMor, P::InfPlc) -> InfPlc

Find a place in the image of $P$ under $m$. If $m$ is an automorphism,
this is unique.
"""
function induce_image(m::NfToNfMor, P::InfPlc)
  k = number_field(P)
  @assert k == domain(m)
  Qx = parent(k.pol)
  h = Qx(preimage(m, gen(codomain(m))))
  #if isdefined(m, :inverse_data)
  #  h = Qx(m.prim_preimg)
  #else
  #  # I need to invert the map.
  #  image_primitive_element(inverse(m))
  #  inv_img = _compute_preimg(m)
  #  h = Qx(inv_img)
  #end
  im = h(P.r)
  lp = infinite_places(codomain(m))
  return lp[findfirst(x -> overlaps(lp[x].r, im), 1:length(lp))]
end

function number_field(P::InfPlc)
  return P.K
end

################################################################################
#
#  Uniformizers
#
################################################################################
@doc Markdown.doc"""
    uniformizer(P::InfPlc) -> nf_elem

Returns an element of the field which is positive at $P$ and negative at all the other real places.
Works only if $P$ is a real place.
"""
function uniformizer(P::InfPlc)
  if is_complex(P)
    error("The place is complex")
  end
  D = infinite_places_uniformizers(number_field(P))
  return deepcopy(D[P])
end

@doc Markdown.doc"""
    infinite_places_uniformizers(K::AnticNumberField)

Returns a dictionary having as keys the real places of $K$ and the values
are uniformizers for the corresponding real place.
A uniformizer of a real place $P$ is an element of the field which is negative
at $P$ and positive at all the other real places.
"""
function infinite_places_uniformizers(K::AnticNumberField)
  r, s = signature(K)
  if iszero(r)
    return Dict{InfPlc, nf_elem}()
  end

  return get_attribute!(K, :infinite_places_uniformizers) do
    return _infinite_places_uniformizers(K)
  end::Dict{InfPlc, nf_elem}
end

function _infinite_places_uniformizers(K::AnticNumberField)
  p = real_places(K) #Important: I assume these are ordered as the roots of the defining polynomial!
  S = abelian_group(Int[2 for i = 1:length(p)])

  s = Vector{GrpAbFinGenElem}(undef, length(p))
  g = Vector{nf_elem}(undef, length(p))
  ind = 1
  q, mq = quo(S, GrpAbFinGenElem[], false)
  b = 10
  cnt = 0
  pr = 16
  while b > 0
    a = rand(K, collect(-b:b))
    if iszero(a)
      continue
    end
    lc = conjugates(a, pr)
    pr = 16
    done = true
    while true
      for i = 1:length(p)
        if contains(reim(lc[i])[1], 0)
          pr *= 2
          done = false
          break
        end
      end
      if !done
        lc = conjugates(a, pr)
        done = true
      else
        break
      end
    end
    ar = zeros(Int, length(p))
    for i = 1:length(p)
      if is_positive(reim(lc[i])[1])
        ar[i] = 0
      else
        ar[i] = 1
      end
    end
    t = S(ar)
    if !iszero(mq(t))
      s[ind] = t
      g[ind] = a
      ind += 1
      q, mq = quo(S, s[1:ind-1], false)
      if ind > length(p)
        break
      end
    end
    cnt += 1
    if cnt > 1000
      b *= 2
      cnt = 0
    end
    for i = 1:length(p)
      good = true
      rep = reim(lc[i])[1]
      if is_positive(rep)
        y = -ceil(fmpz, BigFloat(rep))-1
      else
        y = -floor(fmpz, BigFloat(rep))+1
      end
      ar = zeros(Int, length(p))
      for j = 1:length(p)
        el = reim(lc[j])[1]+y
        if !is_positive(el) && !is_negative(el)
          good = false
          break
        end
        if is_positive(reim(lc[j])[1]+y)
          ar[j] = 0
        else
          ar[j] = 1
        end
      end
      if !good
        continue
      end
      t = S(ar)
      if !iszero(mq(t))
        s[ind] = t
        g[ind] = a + y
        ind += 1
        q, mq = quo(S, s[1:ind-1], false)
        if ind > length(p)
          break
        end
      end
    end
    if ind > length(p)
      break
    end
  end
  if b <= 0
    b = 10
    cnt = 0
    lllE = lll(EquationOrder(K))
    bas = basis(lllE, K, copy = false)
    while true
      @assert b>0
      a = rand(bas, 1:b)
      if iszero(a)
        continue
      end
      emb = signs(a, p)
      ar = Int[0 for i = 1:length(p)]
      for i=1:length(p)
        if emb[p[i]] == -1
          ar[i] = 1
        end
      end
      t = S(ar)
      if !iszero(mq(t))
        s[ind] = t
        g[ind] = a
        ind += 1
        q, mq = quo(S, s[1:ind-1], false)
        if ind > length(p)
          break
        end
      else
        cnt += 1
        if cnt > 1000
          b *= 2
          cnt = 0
        end
      end
    end
  end
  hS = hom(S, S, s)   # Change of coordinates so that the canonical basis elements are mapped to the elements found above
  r = Vector{nf_elem}(undef, length(p))
  hS1 = inv(hS)
  for i = 1:length(p)
    y = hS1(S[i])
    auxr = one(K)
    for w = 1:length(p)
      if !iszero(y[w])
        mul!(auxr, auxr, g[w])
      end
    end
    lc = conjugates(auxr, pr)
    while !is_negative(reim(lc[i])[1] + 1)
      auxr *= 2
      lc = conjugates(auxr, pr)
    end
    for j = 1:length(p)
      if j == i
        continue
      end
      while !is_positive(reim(lc[j])[1] - 1)
        auxr *= 2
        lc = conjugates(auxr, pr)
      end
    end
    r[i] = auxr
  end
  D = Dict{InfPlc, nf_elem}()
  for i = 1:length(p)
    D[p[i]] = r[i]
    @hassert :NfOrd 1 sign(r[i], p[i]) == -1
    #@hassert :NfOrd 1 all(x -> isone(x), values(signs(r[i], [P for P in p if P != p[i]])))
  end
  return D
end

function infinite_primes_map(O::NfOrd, p::Vector{InfPlc}, lying_in::NfOrdIdl)
  K = nf(O)

  if isempty(p)
    S = abelian_group(Int[])
    local exp1
    let S = S, lying_in = lying_in, O = O
      function exp1(A::GrpAbFinGenElem)
        return O(minimum(lying_in))
      end
    end

    local log1
    let S = S
      function log1(B::T) where T <: Union{nf_elem, FacElem{nf_elem, AnticNumberField}}
        return id(S)
      end
    end
    return S, exp1, log1
  end

  D = infinite_places_uniformizers(K)

  S = abelian_group(Int[2 for i = 1:length(p)])

  local exp
  let S = S, D = D, p = p, O = O, lying_in = lying_in, K = K
    function exp(A::GrpAbFinGenElem)
      s = one(K)
      if iszero(A)
        return minimum(lying_in)*O(s)
      end
      for i = 1:length(p)
        if isone(A[i])
          mul!(s, s, D[p[i]])
        end
      end
      return minimum(lying_in)*O(s)
    end
  end

  local log
  let S = S, D = D, p = p
    function log(B::T) where T <: Union{nf_elem, FacElem{nf_elem, AnticNumberField}}
      emb = signs(B, p)
      res = zeros(Int, length(p))
      for i=1:length(p)
        if emb[p[i]] == -1
          res[i] = 1
        end
      end
      return S(res)
    end
  end

  return S, exp, log
end


text:
struct Probe
    x::Float64
    y::Float64
    r::Float64
    angle::Float64
end

max_stars_repo_name: matsunagalab/BioAFM

max_issues_repo_name: matsunagalab/BioAFM

max_forks_repo_name: matsunagalab/BioAFM


text:
##### Extrapolations

export extrap!
export pre_post_colons

function extrap! end

#####
##### Self-extrapolation
#####

pre_post_colons(f, dim) = (;
    Ipre = ntuple(i->Colon(), Val(length(size(f)[1:dim-1]))),
    Ipost = ntuple(i->Colon(), Val(length(size(f)[dim+1:end])))
)

""" Simple linear extrapolation `f_g = 2f_b - f_i` """
extrap!(f::AbstractField, grid::Grid, dir::Int) =
    extrap!(f, grid, Val(dir))

function extrap!(f::AbstractField, grid::Grid)
    extrap!(f, grid, 1)
    extrap!(f, grid, 2)
    extrap!(f, grid, 3)
end

function extrap!(f::AbstractField, grid::Grid, ::Val{dim}) where {dim}
    @unpack Ipre, Ipost = pre_post_colons(f, dim)
    f[Ipre..., 1 , Ipost...] .= 2*f[Ipre...,  2  , Ipost...] .- f[Ipre...,  3  , Ipost...]
    f[Ipre...,end, Ipost...] .= 2*f[Ipre...,end-1, Ipost...] .- f[Ipre...,end-2, Ipost...]
    nothing
end

""" Simple linear extrapolation with two staggered fields: `f_g = 2g_b - f_i` """
extrap!(f::AbstractField, g::AbstractField, grid::Grid, dir::Int) =
    extrap!(f, g, grid, Val(dir))

extrap!(face::CellFace{dim}, cent::CellCenter, grid::Grid, ::Val{dim}) where {dim} =
    extrap!(face, cent, grid, Val(dim), Vertex1D(), Center1D())
extrap!(cent::CellCenter, face::CellFace{dim}, grid::Grid, ::Val{dim}) where {dim} =
    extrap!(cent, face, grid, Val(dim), Center1D(), Vertex1D())
extrap!(corn::CellCorner, edge::CellEdge{dim}, grid::Grid, ::Val{dim}) where {dim} =
    extrap!(corn, edge, grid, Val(dim), Vertex1D(), Center1D())
extrap!(edge::CellEdge{dim}, corn::CellCorner, grid::Grid, ::Val{dim}) where {dim} =
    extrap!(edge, corn, grid, Val(dim), Center1D(), Vertex1D())

extrap!(edge::CellEdge{dir}, face::CellFace{dir}, grid::Grid, ::Val{dir}) where {dir} =
    throw(MethodError("Cannot extrap CellEdge{dir} with CellFace{dir}"))
extrap!(face::CellFace{dir}, edge::CellEdge{dir}, grid::Grid, ::Val{dir}) where {dir} =
    throw(MethodError("Cannot extrap CellFace{dir} with CellEdge{dir}"))

extrap!(edge::CellEdge{edgedir}, face::CellFace{facedir}, grid::Grid, ::Val{dim}) where {dim,edgedir,facedir} =
    extrap!(edge, face, grid, Val(dim), Vertex1D(), Center1D())

extrap!(face::CellFace{facedir}, edge::CellEdge{edgedir}, grid::Grid, ::Val{dim}) where {dim,edgedir,facedir} =
    extrap!(face, edge, grid, Val(dim), Center1D(), Vertex1D())

function extrap!(f, g, grid::Grid, ::Val{dim}, ::Center1D, ::Vertex1D) where {dim}
    @unpack Ipre, Ipost = pre_post_colons(f, dim)
    f[Ipre..., 1 , Ipost...] .= 2*g[Ipre...,  2  , Ipost...] .- f[Ipre...,  2  , Ipost...]
    f[Ipre...,end, Ipost...] .= 2*g[Ipre...,end-1, Ipost...] .- f[Ipre...,end-1, Ipost...]
    nothing
end

function extrap!(f, g, grid::Grid, ::Val{dim}, ::Vertex1D, ::Center1D) where {dim}
    @unpack Ipre, Ipost = pre_post_colons(f, dim)
    f[Ipre..., 1 , Ipost...] .= 2*g[Ipre..., 1 , Ipost...] .- f[Ipre...,  2  , Ipost...]
    f[Ipre...,end, Ipost...] .= 2*g[Ipre...,end, Ipost...] .- f[Ipre...,end-1, Ipost...]
    nothing
end


text:
using ExploreFunction, Test

dataset = [ 0 0;
            0 1;
            0 3]

true_distance = [   0 1 3;
                    1 0 2;
                    3 2 0]

@testset "euclideandistance" begin
    @test euclideandistance(dataset) == true_distance
end

true_prune = [  0 3;
                3 0]
true_remaining = [1 3]

@testset "prunedataset" begin
    prune_dist, remaining = prunedataset(true_distance,1.5)
    @test prune_dist == true_prune
    @test true_remaining == true_remaining
end

dataset_ext = vcat(dataset, [0 6])
true_distance_ext = [   0 1 3 6;
                        1 0 2 5;
                        3 2 0 3;
                        6 5 3 0]
@testset "H0Barcode" begin
    barcode = H0Barcode(true_distance_ext)
    @test barcode[2,:] == [1, 2, 3]
end

@testset "lowerneighbors" begin
    adjacency_test = zeros(4,4)
    adjacency_test[1,2] = adjacency_test[2,1] = 1
    adjacency_test[2,3] = adjacency_test[3,2] = 1
    adjacency_test[3,4] = adjacency_test[4,3] = 1
    adjacency_test[2,4] = adjacency_test[4,2] = 1

    n1 = ExploreFunction.lowerneighbors(1,adjacency_test)
    @test n1 == []

    n2 = ExploreFunction.lowerneighbors(2,adjacency_test)
    @test n2 == [1]

    n3 = ExploreFunction.lowerneighbors(3,adjacency_test)
    @test n3 == [2]

    n4 = ExploreFunction.lowerneighbors(4,adjacency_test)
    @test n4 == [2,3]
end

@testset "edgelist" begin
    adjacency_test = zeros(4,4)
    adjacency_test[1,2] = adjacency_test[2,1] = 1
    adjacency_test[2,3] = adjacency_test[3,2] = 1
    adjacency_test[3,4] = adjacency_test[4,3] = 1
    adjacency_test[2,4] = adjacency_test[4,2] = 2

    edges, weights = ExploreFunction.edgelist(adjacency_test)
    true_edges = [  1 2 2 3;
                    2 3 4 4]
    @test true_edges == edges
    @test weights == [1, 1, 2, 1]
end

@testset "neighborintersect" begin
    adjacency_test = zeros(4,4)
    adjacency_test[1,2] = adjacency_test[2,1] = 1
    adjacency_test[2,3] = adjacency_test[3,2] = 1
    adjacency_test[3,4] = adjacency_test[4,3] = 1
    adjacency_test[2,4] = adjacency_test[4,2] = 1

    s1 = [3, 4]
    i1 = ExploreFunction.neighborintersect(s1,adjacency_test)
    @test [2] == i1

    s2 = [2, 3]
    i2 = ExploreFunction.neighborintersect(s2,adjacency_test)
    @test [] == i2
end

@testset "inductiveVR" begin
    adjacency_test = zeros(4,4)
    adjacency_test[1,2] = adjacency_test[2,1] = 1
    adjacency_test[2,3] = adjacency_test[3,2] = 1
    adjacency_test[3,4] = adjacency_test[4,3] = 2
    adjacency_test[2,4] = adjacency_test[4,2] = 1

    VR, weights = inductiveVR(adjacency_test,2)
    # Test the 1-simplices
    true_s1 = [ 1 2 2 3;
                2 3 4 4]
    @test VR[1] == true_s1

    # Test the 2-simplices
    true_s2 = [2, 3, 4]
    @test VR[2][:] == true_s2

    # Test the 3-simplices
    @test VR[3][:] == []

    # Test Weights
    @test weights[1] == [1, 1, 1, 2]
    @test weights[2] == [2]
end


text:
using Calculus

################################################################################
# Function 1
################################################################################
f = function(x)
    n = length(x)
    total = 0
    for i in 1:n-1
        total += (x[i+1]-x[i])^2 + (exp(x[i]) - 1)^2
    end
    return(total + (exp(x[n]) - 1)^2)
end

grad_f = function(x)
    n = length(x)
    J = zeros(n)
    for i in 1:n
        if i == 1
            J[1] = -2 * (x[2]-x[1]) + 2*(exp(2*x[1]) - exp(x[1]))
        elseif i == n
            J[n] = 2 * (x[n]-x[n-1]) + 2*(exp(2*x[n]) - exp(x[n]))
        else
            J[i] = -2 * (x[i+1]-x[i]) + 2 * (x[i]-x[i-1]) + 2*(exp(2*x[i])-exp(x[i]))
        end
    end
    return(J)
end

hess_f = function(x)
    n = length(x)
    H = zeros(n,n)
    for i in 1:n
        for j in 1:n
            if i == j
                H[i,j] = 2 + 4*exp(2*x[i]) - 2*exp(x[i])
                if (i != n) && (i != 1)
                    H[i,j] += 2
                end

            elseif abs(i-j) == 1
                H[i,j] = -2
            else
                H[i,j] = 0
            end
        end
    end
    return(H)
end

# test that the numerical gradients equal the analytically computed gradients
#x = [1.;2.;3.]
#println(sum(grad_f(x)-Calculus.gradient(f,x)))
#println(sum(hess_f(x)-Calculus.hessian(f,x)))

f_g = Calculus.gradient(f)
f_h = Calculus.hessian(f)

################################################################################
# Rosenbrock
################################################################################

function rosenbrock(x)
    """
    Assumes x has an even number of elemnts.
    """
    n = length(x)
    @assert n%2 == 0
    total = 0
    for i in 1:Integer(n/2)
        total += (1-x[2*i-1])^2 + 10*(x[2*i] - x[2*i-1]^2)^2
    end
    return total
end

rosenbrock_g = Calculus.gradient(rosenbrock)
rosenbrock_h = Calculus.hessian(rosenbrock)

################################################################################
# Cute functions
################################################################################

function cute(x)
    n = length(x)
    total = 0
    for i in 1:n-4
        total += (-4x[i]+3)^2 + (x[i]^2 + 2*x[i+1]^2 + 3*x[i+2]^2 + 4*x[i+3]^2 + 5*x[n]^2)^2
    end
	return total
end

cute_g = Calculus.gradient(cute)
cute_h = Calculus.hessian(cute)

function cute2(x)
    """
    the cute2 function
    """
    n = length(x)
    inds = 2:n-3
    return sum((-4*x[inds]+3.0).^2) +
        sum((x[inds].^2 + 2*x[inds+1].^2 +
		3*x[inds+2].^2 + 4*x[inds+3].^2 + 5*x[1].^2).^2)
end

cute2_g = Calculus.gradient(cute2)
cute2_h = Calculus.hessian(cute2)

################################################################################
# Fenton
################################################################################

function fenton(x)
    """
    Using Newton's method, the starting point [3;2] converges but [3;4] diverges
    """
    x1 = x[1]; x2 = x[2];
    return((12 + x1*x1 + (1+x2*x2)/(x1*x1) + (x1*x1*x2*x2+100)/((x1*x2)^4)) / 10)
end

fenton_g = Calculus.gradient(fenton)
fenton_h = Calculus.hessian(fenton)


text:
# Global variable referring to the standard out stream.
stdout

# Global variable referring to the standard error stream.
stderr

# Global variable referring to the standard input stream.
stdin

# open(filename::AbstractString, [mode::AbstractString]) -> IOStream
io = open("myfile.txt", "w");
write(io, "Hello world!");
close(io);

io = open("myfile.txt", "r");
println(read(io, String))

# IOBuffer
# Create an in-memory I/O stream, which may optionally operate on a pre-existing array.

# flush
# Commit all currently buffered writes to the given stream.

# seek
# Seek a stream to the given position.
io = IOBuffer("JuliaLang is a GitHub organization.");
seek(io, 5);
println(read(io, Char))

# truncate
# Resize the file or buffer given by the first argument to exactly n bytes, filling previously unallocated space with '\0' if the file or buffer is grown.
io = IOBuffer();
write(io, "JuliaLang is a GitHub organization.")
truncate(io, 15);
println(String(take!(io)))

max_stars_repo_name: startprogress/Hello-Julia

max_issues_repo_name: startprogress/Hello-Julia

max_forks_repo_name: startprogress/Hello-Julia


text:
uniform(key::JaxRNGKey, dims...; kwargs...) =
  uniform(Float64, key, dims...; kwargs...)
uniform(T::Type{<:Number}, key::JaxRNGKey, dims...; kwargs...) =
  uniform(T, key, Dims(dims); kwargs...)
function uniform(
  T::Type{<:Number},
  key::JaxRNGKey,
  dims::Dims;
  minval = 0.0,
  maxval = 1.0,
)
  _random.uniform(key, reverse(dims), jl_to_np_type(T); minval = minval, maxval = maxval)
end

function uniform(
  T::Type{<:Integer},
  key::JaxRNGKey,
  dims::Dims;
  minval::Integer = 0,
  maxval::Integer = 1,
)
  _random.randint(
    key,
    dims,
    dtype = jl_to_np_type(T),
    minval = minval,
    maxval = maxval,
  )
end

normal(key::JaxRNGKey, dims...; kwargs...) = uniform(Float64, key, dims...)
normal(T::Type{<:Number}, key::JaxRNGKey, dims...) = normal(T, key, Dims(dims))
function normal(T::Type{<:Number}, key::JaxRNGKey, dims::Dims)
  _random.normal(key, reverse(dims), dtype = jl_to_np_type(T))
end


text:
##
using Plots
plotlyjs()

theme(:juno)

##
using IJulia
display(1)

print("hello doug")

##
plot([1,2,1])

function f(x)
    return x
end

f(3)

3+3

using IJulia
notebook()

## Simple testing suite

## Preloads
using Statistics
using FFTW
using Plots
using BenchmarkTools
using Profile
using LinearAlgebra

## BLAS.dot

test_img = rand(256,256)

function variance_BLAS(A)
    (Nx, Ny)  = size(A)
    norm_im_vec = reshape(A, Nx*Ny)
    BLAS.dot(norm_im_vec,norm_im_vec)/(Nx*Ny)
end

function variance(A)
    (Nx, Ny)  = size(A)
    sum(A .* A)/(Nx*Ny)
end

function variance_BLAS_sq(A)
    (Nx, Ny)  = size(A)
    BLAS.dot(A,A)/(Nx*Ny)
end

@time variance(test_img)

@benchmark variance(test_img)

@time variance_BLAS(test_img)

@benchmark variance_BLAS(test_img)

#BLAS dot 30x faster; Was the unflattened version really wrong

@time variance_BLAS_sq(test_img)

@benchmark variance_BLAS_sq(test_img)

function variance_lc(A)
    (Nx, Ny)  = size(A)
    sum([x^2 for x in A])/(Nx*Ny)
end

@time variance_lc(test_img)

@benchmark variance_lc(test_img)

function variance_lc_abs2(A)
    (Nx, Ny)  = size(A)
    sum([abs2(x) for x in A])/(Nx*Ny)
end

@time variance_lc_abs2(test_img)

@benchmark variance_lc_abs2(test_img)

## Test the * versus mul!

function Doug_mul(A)
    B_test = zeros(256,256)
    B_test     .= A' * A
end

@time Doug_mul(test_img)

@benchmark Doug_mul(test_img)

function noraml_mul(A)
    B_test = zeros(256,256)
    mul!(B_test,A',A)
end

@time noraml_mul(test_img)

@benchmark noraml_mul(test_img)

BLAS.syrk!(test_img,B_test)

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2020
#
# This is a Julia package for discrete event simulation
#

println("... doctests  ...")

using Documenter

resetClock!(𝐶)
doctest(DiscreteEvents; manual = false)


text:
using HiQGA.transD_GP
fname_dfn = "L103004.dfn"
zipsaveprefix = splitpath(fname_dfn)[end]
# electronics file
fname_specs_halt = "electronics_halt.jl"
# column numbers from hdr file
X, Y, Z = "E", "N", "DTM_AHD"
fid = "Fid"
linenum = "Line"
frame_height ="Height"
frame_dx = "dx"
frame_dy = "dy"
frame_dz = "dz"
HM_Z = "HM_Z_"
LM_Z = "LM_Z_"  # first LM gate dropped, also noise and times in electronics_halt
HM_σ = "RUNC_HM_Z_"
LM_σ = "RUNC_LM_Z_"
relerror = false
units = 1e-12
LM_drop,HM_drop = 9,15
##
soundings = transD_GP.SkyTEM1DInversion.read_survey_files(fname_dfn,
						             fname_specs_halt = fname_specs_halt,
						             LM_Z             = LM_Z,
									 HM_Z             = HM_Z,
									 LM_σ             = LM_σ,
									 HM_σ             = HM_σ,
									 LM_drop          = LM_drop,
									 HM_drop          = HM_drop,
									 frame_height     = frame_height,
									 frame_dz         = frame_dz,
									 frame_dy         = frame_dy,
									 frame_dx         = frame_dx,
									 X                = X,
									 Y                = Y,
                                     Z                = Z,
									 fid              = fid,
									 units            = units,
									 relerror         = relerror,
									 linenum          = linenum,
                                     startfrom        = 1,
									 skipevery        = 20,
									 dotillsounding   = nothing,
									 makesounding     = true)

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA

max_forks_repo_name: GeoscienceAustralia/HiQGA


text:
# implement QZ-40 algorithm
# Dellnitz, M., Schütze, O., & Zheng, Q. (2002). Locating all the zeros of an analytic function in one complex variable. Journal of Computational and Applied Mathematics, 138(2), 325–333. https://doi.org/10.1016/S0377-0427(01)00371-5
# currently assumes all roots have multiplicity one

using QuadGK

function winding_number(dlogfdz::Function, rectangle::Array{Float64,1}, maxevals::Int = 10^4)::Int
    (x_min, x_max, y_min, y_max) = rectangle
    (bl, br, tr, tl) = (x_min + y_min*im, x_max + y_min*im, x_max + y_max*im, x_min + y_max*im)
    round(Int,real(1 ./(2*pi*im)*quadgk(dlogfdz,bl,br,tr,tl,bl,rtol=0, atol=0.45*2*pi,maxevals=maxevals)[1]))
end

function qz40(dlogfdz::Function, rectangle::Array{Float64,1}, verbose::Bool=false, maxevals::Int=10^4)

    rectangles = Array{Float64,1}[rectangle]
    zeroes = Complex[]
    wn = -1
    qz40_recursive(dlogfdz,rectangles,zeroes,wn,verbose,maxevals)
end

function locate_root(dlogfdz::Function, rectangle::Array{Float64,1}, maxevals::Int = 10^4)
    (x_min, x_max, y_min, y_max) = rectangle
    (bl, br, tr, tl) = (x_min + y_min*im, x_max + y_min*im, x_max + y_max*im, x_min + y_max*im)
    integrand(z::Complex) = z*dlogfdz(z)
    location = 1 ./(2*pi*im)*quadgk(integrand,bl,br,tr,tl,bl,maxevals=maxevals)[1]
    if (real(location) < x_min || real(location) > x_max || imag(location) < y_min || imag(location) > y_max)
        error("Root location $(location) outside of rectangle (x_min, x_max, y_min, y_max) $(x_min),$(x_max),$(y_min),$(y_max).")
    end
    return location
end

function random_subdivide(rectangle::Array{Float64,1})

    (x_min, x_max, y_min, y_max) = rectangle
    x_mid = x_min + rand()*(x_max - x_min)
    y_mid = y_min + rand()*(y_max - y_min)
    # rectangles = Array{Array{Float64,1},1}[]
    rectangles = Array{Array{Float64,1},1}(undef,0)
    push!(rectangles,[x_min, x_mid, y_min, y_mid]);
    push!(rectangles,[x_mid, x_max, y_min, y_mid]);
    push!(rectangles,[x_mid, x_max, y_mid, y_max]);
    push!(rectangles,[x_min, x_mid, y_mid, y_max]);
    return rectangles
end

function qz40_recursive(dlogfdz::Function, rectangles::Array{Array{Float64,1},1},zeroes::Vector{Complex},nroots::Int,verbose::Bool,maxevals::Int)
    # the rows of rectangles are the delimiters (x_min, x_max, y_min, y_max) of the rectangles
    new_rectangles = Array{Float64,1}[]
    n_rectangles = size(rectangles)[1]
    if (verbose)
        println("Number of rectangles: $(n_rectangles)")
    end
    nroots_new = length(zeroes)
    for i in 1:n_rectangles
        rect = rectangles[i]
        wn = winding_number(dlogfdz,rect,maxevals)
        nroots_new += wn
        if (wn == 1)
            push!(zeroes,locate_root(dlogfdz,rect,maxevals))
        elseif (wn > 1)
            append!(new_rectangles,random_subdivide(rect))
        elseif (wn < 0)
            error("Negative winding numbers")
        end
    end
    if (nroots > 0 && nroots_new != nroots)
        error("Expected number of roots: $(nroots), new number of roots: $(nroots_new).")
    end
    if (nroots == -1)
        println("Number of roots: $(nroots_new)")
    end
    if size(new_rectangles)[1] == 0
        return zeroes
    else
        return qz40_recursive(dlogfdz, new_rectangles, zeroes, nroots_new, verbose,maxevals)
    end
end

function newton(f, dfdz, z_init, n_iter)

    z_current = z_init
    for i in 1:n_iter
        z_current = z_current - f(z_current)/dfdz(z_current)
    end
    z_current
end

max_stars_repo_name: jbierkens/spectral

max_issues_repo_name: jbierkens/spectral

max_forks_repo_name: jbierkens/spectral


text:
using Mmap

########################
# fractions of formula #
########################

factorial_div(i::Integer, j::Integer) = factorial(big(i)) / factorial(big(j))

z(x::Real, p::Real) = sqrt(2.)*(x + p*im)
z(x::Vector{<:Real}, p::Vector{<:Real}) = z.(x, p')

gaussian_function(x::Real, p::Real) = exp(-0.5 * abs2(z(x,p))) / π
gaussian_function(x::Vector{<:Real}, p::Vector{<:Real}) = gaussian_function.(x, p')

neg_one_to_power_of(i::Integer) = (i % 2 == 0) ? 1 : -1

function coefficient_of_wave_function(m::Integer, n::Integer)
    if n ≥ m
        # adjust index bases for number state on `m`
        return neg_one_to_power_of(m-1) * sqrt(factorial_div(m, n))
    else
        # adjust index bases for number state on `n`
        return neg_one_to_power_of(n-1) * sqrt(factorial_div(n, m))
    end
end

function z_to_power(m::Integer, n::Integer, x::Real, p::Real)
    if n ≥ m
        return conj(z(x, p'))^(n - m)
    else
        return z(x, p')^(m - n)
    end
end

function z_to_power(m::Integer, n::Integer, x::Vector{<:Real}, p::Vector{<:Real})
    z_to_power.(m, n, x, p')
end

function laguerre(m::Integer, n::Integer, x::Real, p::Real)
    if n ≥ m
        # adjust index bases for number state on `m`
        return laguerre(m-1, n - m, abs2(z(x, p)))
    else
        # adjust index bases for number state on `n`
        return laguerre(n-1, m - n, abs2(z(x, p)))
    end
end

function laguerre(m::Integer, n::Integer, x::Vector{<:Real}, p::Vector{<:Real})
    return laguerre.(m, n, x, p')
end

#########
# utils #
#########

function save_𝐰(bin_path::String, 𝐰::Array{ComplexF64,4})
    @info "Save W_{m,n,x,p} to $bin_path"
    mem = open(bin_path, "w+")
    write(mem, 𝐰)
    close(mem)
end

function load_𝐰(
    m_dim::Integer,
    n_dim::Integer,
    x_range::AbstractRange,
    p_range::AbstractRange,
    bin_path::String
)
    @info "Load W_{m,n,x,p} from $bin_path"
    mem = open(bin_path)
    𝐰 = Mmap.mmap(
        mem,
        Array{ComplexF64,4},
        (m_dim, n_dim, length(x_range), length(p_range))
    )
    close(mem)

    return 𝐰
end

function gen_wigner_bin_path(
    m_dim::Integer,
    n_dim::Integer,
    x_range::AbstractRange,
    p_range::AbstractRange,
)
    bin_path = joinpath(
        mkpath(joinpath(datadep_root(), "wigner_function")),
        "W " *
        "m=$(m_dim) n=$(n_dim) " *
        "x=$(range2str(x_range)) p=$(range2str(p_range)).bin"
    )

    return bin_path
end

range2str(range::AbstractRange) = replace(string(range), r":|," => "_")

check_zero(m_dim, n_dim) = !iszero(m_dim) && !iszero(n_dim)

check_empty(x_range, p_range) = !isempty(x_range) && !isempty(p_range)

function check_argv(m_dim, n_dim, x_range, p_range)
    return check_zero(m_dim, n_dim) && check_empty(x_range, p_range)
end


text:
using Printf, Distributions

newgrid(p::Float64, M::Int=15, N::Int=15) = rand(Bernoulli(p), M, N)

function walkmaze!(grid::Matrix{Int}, r::Int, c::Int, indx::Int)
    NOT_VISITED = 1 # const
    N, M = size(grid)
    dirs = [[1, 0], [-1, 0], [0, 1], [1, 0]]
    # fill cell
    grid[r, c] = indx

    # is the bottom line?
    rst = r == N

    # for each direction, if has not reached the bottom yet and can continue go to that direction
    for d in dirs
        rr, cc = (r, c) .+ d
        if !rst && checkbounds(Bool, grid, rr, cc) && grid[rr, cc] == NOT_VISITED
            rst = walkmaze!(grid, rr, cc, indx)
        end
    end
    return rst
end

function checkpath!(grid::Matrix{Int})
    NOT_VISITED = 1 # const
    N, M = size(grid)
    walkind = 1
    for m in 1:M
        if grid[1, m] == NOT_VISITED
            walkind += 1
            if walkmaze!(grid, 1, m, walkind)
                return true
            end
        end
    end
    return false
end

function printgrid(G::Matrix{Int})
    LETTERS = vcat(' ', '#', 'A':'Z')
    for r in 1:size(G, 1)
        println(r % 10, ") ", join(LETTERS[G[r, :] .+ 1], ' '))
    end
    if any(G[end, :] .> 1)
        println("!) ", join((ifelse(c > 1, LETTERS[c+1], ' ') for c in G[end, :]), ' '))
    end
end

const nrep = 1000 # const
sampleprinted = false

p = collect(0.0:0.1:1.0)
f = similar(p)
for i in linearindices(f)
    c = 0
    for _ in 1:nrep
        G = newgrid(p[i])
        perc = checkpath!(G)
        if perc
            c += 1
            if !sampleprinted
                @printf("Sample percolation, %i×%i grid, p = %.2f\n\n", size(G, 1), size(G, 2), p[i])
                printgrid(G)
                sampleprinted = true
            end
        end
    end
    f[i] = c / nrep
end

println("\nFrequencies for $nrep tries that percolate through\n")
for (pi, fi) in zip(p, f)
    @printf("p = %.1f ⇛ f = %.3f\n", pi, fi)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@testset "verifykeys" begin
    @testset "JSON3.Object" begin
        json = JSON3.read("""
        {
            "a": 1,
            "b": 2,
            "c": 3
        }
        """)

        @test verifykeys(json, (:a, :b, :c)) == true
        @test verifykeys(json, ["a", "b", "c"]) == true

        @test verifykeys(json, [:a, :d, :c]) == false
        @test verifykeys(json, ("a", "d", "c")) == false
    end

    @testset "WebAPI.Params" begin
        params = WebAPI.Params(:a => "1", :b => "2", :c => "3")

        @test verifykeys(params, (:a, :b, :c)) == true
        @test verifykeys(params, ["a", "b", "c"]) == true

        @test verifykeys(params, [:a, :d, :c]) == false
        @test verifykeys(params, ("a", "d", "c")) == false
    end

    @testset "WebAPI.Query" begin
        query = WebAPI.Query(:a => "1", :b => "2", :c => "3")

        @test verifykeys(query, (:a, :b, :c)) == true
        @test verifykeys(query, ["a", "b", "c"]) == true

        @test verifykeys(query, [:a, :d, :c]) == false
        @test verifykeys(query, ("a", "d", "c")) == false
    end
end


text:
type ExactIterator
    inneritr::Base.Combinations{Array{Int64,1}}
    n::Int
end

exact(n::Int) = ExactIterator(combinations(collect(1:(2n-1)), n), n)

eltype(itr::ExactIterator) = Array{Int64, 1}

start(itr::ExactIterator) = start(itr.inneritr)

function next(itr::ExactIterator, s)
    c, s = next(itr.inneritr, s)
    v = Array(Int, itr.n)
    j = 1
    p = 1
    for k = 1:itr.n
        while !(j in c)
            j += 1
            p += 1
        end
        v[k] = p
        j += 1
    end
    return v, s
end

done(itr::ExactIterator, s) = done(itr.inneritr, s)

length(itr::ExactIterator) = length(itr.inneritr)

type ExactIterator2
    inneritr::Base.Combinations{Array{Int64,1}}
    n::Int
end

exact2(n::Int) = ExactIterator2(combinations(collect(1:(2n-1)), n-1), n)

function next(itr::ExactIterator2, s)
    (c, s) = next(itr.inneritr, s)
    c3 = [c, 2*itr.n; ] .- [0, c; ]
    (vcat([fill(i, c3[itr.n-i+1]-1) for i = 1:itr.n]...), s)
end

done(itr::ExactIterator2, s) = done(itr.inneritr, s)

length(itr::ExactIterator2) = length(itr.inneritr)

eltype(itr::ExactIterator2) = Array{Int64, 1}

start(itr::ExactIterator2) = start(itr.inneritr)

using Iterators

function sample_set(c::Array) ## Integer,1
    v = Array(Int, length(c))
    j = 1
    p = 1
    for k = 1:length(c)
        while !(j in c)
            j += 1
            p += 1
        end
        v[k] = p
        j += 1
    end
    return v
end

function exact3(n::Integer)
    return imap(sample_set, combinations([1:(2n-1);], n))
end

@time s1 = collect(exact(12));
# n = 12: 2.58s, 6.76 M allocations: 711.773 MB

@time s2 = collect(exact2(12)); ## slow and much memory
# n = 12: 30.23s, 104.11 M allocations: 6.358 GB

@time s3 = collect(exact3(12));
# n = 12: 11.22s, 16.22 M allocations: 1.145 GB

sort([hash(x) for x in s1]) == sort([hash(x) for x in s2])

sort([hash(x) for x in s1]) == sort([hash(x) for x in s3])


text:
struct SystemMatrix
    rows::Any
    cols::Any
    vals::Any
    function SystemMatrix(rows, cols, vals)
        @assert length(rows) == length(cols) == length(vals)
        new(rows, cols, vals)
    end
end

function SystemMatrix()
    rows = Int[]
    cols = Int[]
    vals = zeros(0)
    SystemMatrix(rows, cols, vals)
end

function Base.show(io::IO, sysmatrix::SystemMatrix)
    numvals = length(sysmatrix.rows)
    str = "SystemMatrix with $numvals entries"
    print(io, str)
end

function assemble!(matrix, rows, cols, vals)
    @assert length(rows) == length(cols) == length(vals)
    append!(matrix.rows, rows)
    append!(matrix.cols, cols)
    append!(matrix.vals, vals)
end

struct SystemRHS
    rows::Any
    vals::Any
    function SystemRHS(rows, vals)
        @assert length(rows) == length(vals)
        new(rows, vals)
    end
end

function SystemRHS()
    SystemRHS(Int[], zeros(0))
end

function Base.show(io::IO, sysrhs::SystemRHS)
    numvals = length(sysrhs.rows)
    str = "SystemRHS with $numvals entries"
    print(io, str)
end

function assemble!(systemrhs, rows, vals)
    @assert length(rows) == length(vals)
    append!(systemrhs.rows, rows)
    append!(systemrhs.vals, vals)
end

function node_to_dof_id(nodeid, dofs, dofspernode)
    return (nodeid - 1) * dofspernode + dofs
end

function element_dofs(nodeids, dofs, dofspernode)
    numnodes = length(nodeids)
    numdofs = length(dofs)

    extnodeids = repeat(nodeids, inner = numdofs)
    extdofs = repeat(dofs, outer = numnodes)

    edofs = [
        node_to_dof_id(n, d, dofspernode) for
        (n, d) in zip(extnodeids, extdofs)
    ]
    return edofs
end

function element_dofs(nodeids, dofspernode)
    dofs = 1:dofspernode
    return element_dofs(nodeids, dofs, dofspernode)
end

function element_dofs_to_operator_dofs(rowdofs, coldofs)
    nr = length(rowdofs)
    nc = length(coldofs)
    rows = repeat(rowdofs, outer = nc)
    cols = repeat(coldofs, inner = nr)
    return rows, cols
end

function assemble_couple_cell_matrix!(
    sysmatrix,
    nodeids1,
    dofs1,
    nodeids2,
    dofs2,
    dofspernode,
    vals,
)

    edofs1 = element_dofs(nodeids1, dofs1, dofspernode)
    edofs2 = element_dofs(nodeids2, dofs2, dofspernode)

    rows, cols = element_dofs_to_operator_dofs(edofs1, edofs2)
    assemble!(sysmatrix, rows, cols, vals)
end

function assemble_couple_cell_matrix!(
    sysmatrix,
    nodeids1,
    nodeids2,
    dofspernode,
    vals,
)
    edofs1 = element_dofs(nodeids1, dofspernode)
    edofs2 = element_dofs(nodeids2, dofspernode)

    rows, cols = element_dofs_to_operator_dofs(edofs1, edofs2)
    assemble!(sysmatrix, rows, cols, vals)
end

function assemble_cell_matrix!(sysmatrix, nodeids, dofspernode, vals)
    assemble_couple_cell_matrix!(sysmatrix, nodeids, nodeids, dofspernode, vals)
end

function assemble_cell_rhs!(sysrhs, nodeids, dofs, dofspernode, vals)
    rows = element_dofs(nodeids, dofs, dofspernode)
    assemble!(sysrhs, rows, vals)
end

function assemble_cell_rhs!(sysrhs, nodeids, dofspernode, vals)
    dofs = 1:dofspernode
    assemble_cell_rhs!(sysrhs,nodeids,dofs,dofspernode,vals)
end

function sparse_operator(sysmatrix, ndofs)
    return dropzeros!(
        sparse(sysmatrix.rows, sysmatrix.cols, sysmatrix.vals, ndofs, ndofs),
    )
end

function sparse_operator(sysmatrix, mesh, dofspernode)
    numnodes = number_of_nodes(mesh)
    totaldofs = dofspernode * numnodes
    return sparse_operator(sysmatrix, totaldofs)
end

function sparse_displacement_operator(sysmatrix, mesh)
    dim = dimension(mesh)
    return sparse_operator(sysmatrix, mesh, dim)
end

function rhs_vector(sysrhs, ndofs)
    return Array(sparsevec(sysrhs.rows, sysrhs.vals, ndofs))
end

function rhs_vector(sysrhs, mesh, dofspernode)
    numnodes = number_of_nodes(mesh)
    totaldofs = dofspernode * numnodes
    return rhs_vector(sysrhs, totaldofs)
end

function displacement_rhs_vector(sysrhs, mesh)
    dim = dimension(mesh)
    return rhs_vector(sysrhs, mesh, dim)
end


text:
using LinearAlgebra: svd, svdvals, cond, pinv, I, diagm, norm
using Statistics
using Distributions, TransformVariables
using Random
using JLD

include("../../../src/PotentialUQ.jl")

######################################################
##################### Setup ##########################
######################################################

ϵ = 1.0
σ = 1.0
Temp = 0.25
seed = randperm(100)
ljtrue = PotentialUQ.Potentials.LennardJones(ϵ, σ)
num_configs = 40
r_train = Vector{PotentialUQ.Potentials.Configuration}(undef, num_configs)
radii = 3.4
for j = 1:num_configs
    file_path = "examples/Argon7/DATA_training/TEMP_$(Int(Temp*100))/seed_1/DATA.$(Int(10000*j*2.5))"
    c_temp = PotentialUQ.Potentials.Configuration(file_path; atom_names = [:Ar],
                    radii = [radii], weights = [1.0])
    r_train[j] = c_temp
end

num_configs_test = 25
r_test = Vector{PotentialUQ.Potentials.Configuration}(undef, num_configs_test)
for j = 1:num_configs_test
    file_path = "examples/Argon7/DATA_training/TEMP_$(Int(Temp*100))/seed_1/DATA.$(Int(10000*j*3))"
    c_temp = PotentialUQ.Potentials.Configuration(file_path; atom_names = [:Ar],
                    radii = [radii], weights = [1.0])
    r_test[j] = c_temp
end

etrain = PotentialUQ.Potentials.potential_energy(r_train, ljtrue)
etest = PotentialUQ.Potentials.potential_energy(r_test, ljtrue)

ftrain = PotentialUQ.Potentials.force(r_train, ljtrue)
ftest = PotentialUQ.Potentials.force(r_test, ljtrue)

vtrain = PotentialUQ.Potentials.virial(r_train, ljtrue)
vtest = PotentialUQ.Potentials.virial(r_test, ljtrue)

vttrain = PotentialUQ.Potentials.virial_stress(r_train, ljtrue)
vttest = PotentialUQ.Potentials.virial_stress(r_train, ljtrue)

###################### Setup SNAP ####################
# SNAP
b = vcat(etrain, reduce(vcat, reduce(vcat, ftrain)), reduce(vcat, vttrain))
# b = etrain
twojmax = 6
rcutfac = 0.75/radii
snap = PotentialUQ.Potentials.SNAP(rcutfac, twojmax, r_train[1].num_atom_types)
A = PotentialUQ.Potentials.get_snap(r_train, snap)

# Fit
snap.β = A\b
# show(stdout, "text/plain", snap.β)
# println(" ")
# Test
e_snap = PotentialUQ.Potentials.potential_energy(r_test, snap)
error = mean(abs.(etest - e_snap)./abs.(etest))

println("Relative Error = ", error)

f_snap = PotentialUQ.Potentials.force(r_test, snap)
ff = vcat( (ftest - f_snap)... )
fff = vcat( ftest... )
error = mean( norm.(ff)./norm.(fff)  )
println("Force Relative Error = ",  error)

v_snap = PotentialUQ.Potentials.virial(r_test, snap)
error = mean(abs.(vtest - v_snap)./abs.(vtest))
println("Virial Error = ", error)
println(" ")

######################################################
##################### UQ #############################
######################################################

# Unknowns
n, m = size(A)

# β = pinv(A' * A, 1e-6)*(A'*b)
β = snap.β
nt = (β = β, )
# println("v = ", nt)
t = PotentialUQ.namedtp_to_vec(nt)
x = TransformVariables.inverse(t, nt)

# # Prior
F = svd(A)
Σ = diagm(F.S)
scale_A = vec(std(A, dims = 2))
Q =  F.U' * ( diagm(scale_A) ) * (F.U)
Q = I(m) + 0.5*Q + 0.5*Q'
y =  copy(F.U' * b)

# sampler(rng::Random.AbstractRNG, d::ContinuousMultivariateDistribution) = β + 1e-1*randn(length(d))
# logpdf(x::Vector) = loglikelihood(MvNormal(β, 1e-1*I(m)), x)
# snap_prior = PotentialUQ.SNAP_Prior_Distribution(A, b, x, t, sampler, logpdf)
snap_prior = PotentialUQ.SNAP_Prior_Distribution(A, b, nt)
println("logpdf = ", Distributions.logpdf(snap_prior, x))
snap_likelihood(x::Vector) = loglikelihood( MvNormal( Σ*F.V'*x - y, Q ),  zeros(m) )
println("Likelihood = ", snap_likelihood(x))
dsnap = PotentialUQ.PotentialDistribution(snap, nt, snap_prior, snap_likelihood)

@time chain, samp = PotentialUQ.Sample(dsnap; return_chain = true, num_samples = 1000, num_adapts = 500, verbose = true)

snap.β = dsnap.x.β
println("MAP Fitted β = ")
show(stdout, "text/plain", snap.β)
println(" ")
e_snap = PotentialUQ.Potentials.potential_energy(r_test, snap)
println("MAP SNAP potential_energy = ", e_snap)
println("MAP Relative Error = ", mean(abs.(etest - e_snap)./abs.(etest)))
v_snap = PotentialUQ.Potentials.virial(r_test, snap)
f_snap = PotentialUQ.Potentials.force(r_test, snap)
ff = vcat( (ftest - f_snap)... )
fff = vcat( ftest... )
println("Force Relative Error = ", mean( norm.(ff)./norm.(fff)  ) )
println("MAP SNAP Virial = ", v_snap)
println("MAP Virial Error = ", mean(abs.(vtest - v_snap)./abs.(vtest)))

save("examples/Argon7/snap_estimation/TEMP_25/samples_snap.jld", "samples", samp, "map", snap.β, "A", A, "mle", A\b)


text:
push!(LOAD_PATH, "./src")
using DelimitedFiles

using WriteVTK

using CommonUtils
using Basis1D
using Basis2DQuad
using UniformQuadMesh
using UnPack

using SetupDG

N = 2
K1D = 20
# t=0.027652618923773267
# t=0.0008610198724528106
#t=0.004593014339084991
#t=0.009755907564170784
# t=0.021238249430540698
# t=0.01107131919723909
# t=0.016685099657334425
#t=0.038174565986393724
# t=0.09433024491945503
# t=0.09655400705355822
# t=0.19575285632508105
# t=0.07550828108140312
# t=0.07590867790313385
# t=0.024393258642572562
# t=0.025518477789014574
# t=0.024393255476498787
# t=0.025518475001162794
# t=0.10674597963531857
# t = 0.2
# t =0.02549578894035154
# t=0.013226118481316584
# t=0.013037857428159683
t=0.014477691777086832
t=0.009145127272670283
t=0.013037781513546847
t=0.19060893512806573
t=0.2
t=0.06000064065592031
t=0.0783551163019341
t=0.09945653897728977
t=0.2
t=0.17009267106494133
t=0.0
flname = "dmr"
CFL = 0.9
XLENGTH = 2.0
Np = (N+1)*(N+1)

# Initial condition 2D shocktube
const γ = 1.4
const rhoL  = 8.0
const uL    = 8.25*cos(pi/6)
const vL    = -8.25*sin(pi/6)
const pL    = 116.5
const rhouL = rhoL*uL
const rhovL = rhoL*vL
const EL    = pL/(γ-1)+.5*rhoL*(uL^2+vL^2)
const betaL = rhoL/(2*pL)
const rhoR  = 1.4
const uR    = 0.0
const vR    = 0.0
const pR    = 1.0
const rhouR = rhoR*uR
const rhovR = rhoR*vR
const ER    = pR/(γ-1)+.5*rhoR*(uR^2+vR^2)
const betaR = rhoR/(2*pR)

"Mesh related variables"
# VX, VY, EToV = uniform_quad_mesh(Int(round(XLENGTH*K1D)),K1D)
# @. VX = (VX+1)/2*XLENGTH
# @. VY = (VY+1)/2
VX, VY, EToV = uniform_quad_mesh(Int(round(XLENGTH*K1D)),K1D)
# @. VX = (VX+1)/2
# @. VY = (VY+1)/4
@. VX = (VX+1)/2*XLENGTH*5
@. VY = (VY+1)/2*5

rd = init_reference_quad(N,gauss_lobatto_quad(0,0,N))
md = init_mesh((VX,VY),EToV,rd)
@unpack x,y = md;

@unpack VDM = rd
rp,sp = equi_nodes_2D(10)
Vp = vandermonde_2D(N,rp,sp)/VDM

r_vtk = x[1:4]
x_vtk = []
y_vtk = []
for k = 1:Int(round(XLENGTH*K1D))
    push!(x_vtk,r_vtk[1]+(k-1)*5/K1D)
    push!(x_vtk,r_vtk[2]+(k-1)*5/K1D)
    push!(x_vtk,r_vtk[3]+(k-1)*5/K1D)
    # push!(x_vtk,r_vtk[4]+(k-1)*1/K1D)
end
for k = 1:K1D
    push!(y_vtk,r_vtk[1]+(k-1)*5/K1D)
    push!(y_vtk,r_vtk[2]+(k-1)*5/K1D)
    push!(y_vtk,r_vtk[3]+(k-1)*5/K1D)
    # push!(y_vtk,r_vtk[4]+(k-1)*1/K1D)
end

x_vtk = x_vtk .* ones((N+1)*K1D,1)'
y_vtk = y_vtk' .* ones((N+1)*Int(round(XLENGTH*K1D)),1)
# y_vtk = y_vtk' .* ones((N+1)*Int(round(2*K1D)),1)
rho_vtk,rhou_vtk,rhov_vtk,E_vtk,sigmax_vtk,sigmay_vtk = (zeros(size(x_vtk)) for i = 1:6)

# vtk_grid("N=$N,K1D=$K1D,t=$t,CFL=$CFL,x=$XLENGTH",x,y) do vtk
#     vtk["rho"] = readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,xt=0.014317166411179642=$XLENGTH,rho,dmr.txt")
#     vtk["rhou"] = readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,x=$XLENGTH,rhou,dmr.txt")
#     vtk["rhov"] = readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,x=$XLENGTH,rhov,dmr.txt")
#     vtk["E"] = readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,x=$XLENGTH,E,dmr.txt")
# end

vtk_grid("N=$N,K1D=$K1D,t=$t,CFL=$CFL,x=$XLENGTH,$flname",x_vtk,y_vtk) do vtk
    rho = readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,x=$XLENGTH,rho,$flname.txt")
    rhou = readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,x=$XLENGTH,rhou,$flname.txt")
    rhov = readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,x=$XLENGTH,rhov,$flname.txt")
    E = readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,x=$XLENGTH,E,$flname.txt")
    # rho = readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,rho,dmr.txt")
    # rhou = readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,rhou,dmr.txt")
    # rhov = readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,rhov,dmr.txt")
    # E = readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,E,dmr.txt")
    # sigmax = readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,x=$XLENGTH,sigmax,dmr.txt")
    # sigmay = readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,x=$XLENGTH,sigmay,dmr.txt")

    for k = 1:Int(round(XLENGTH*K1D))*K1D
        kx = mod1(k,Int(round(XLENGTH*K1D)))
        ky = div(k-1,Int(round(XLENGTH*K1D)))+1

        # if kx > 3.3*K1D
        #     rho[:,k] .= rhoR
        #     rhou[:,k] .= rhouR
        #     rhov[:,k] .= rhovR
        #     E[:,k] .= ER
        # end

        rho_vtk[(kx-1)*(N+1)+1:(kx-1)*(N+1)+N+1, (ky-1)*(N+1)+1:(ky-1)*(N+1)+N+1] = rho[:,k]
        rhou_vtk[(kx-1)*(N+1)+1:(kx-1)*(N+1)+N+1, (ky-1)*(N+1)+1:(ky-1)*(N+1)+N+1] = rhou[:,k]
        rhov_vtk[(kx-1)*(N+1)+1:(kx-1)*(N+1)+N+1, (ky-1)*(N+1)+1:(ky-1)*(N+1)+N+1] = rhov[:,k]
        E_vtk[(kx-1)*(N+1)+1:(kx-1)*(N+1)+N+1, (ky-1)*(N+1)+1:(ky-1)*(N+1)+N+1] = E[:,k]
        # sigmax_vtk[(kx-1)*(N+1)+1:(kx-1)*(N+1)+N+1, (ky-1)*(N+1)+1:(ky-1)*(N+1)+N+1] = sigmax[:,k]
        # sigmay_vtk[(kx-1)*(N+1)+1:(kx-1)*(N+1)+N+1, (ky-1)*(N+1)+1:(ky-1)*(N+1)+N+1] = sigmay[:,k]
    end

    vtk["rho"]  = rho_vtk
    vtk["rhou"] = rhou_vtk
    vtk["rhov"] = rhov_vtk
    vtk["E"]    = E_vtk
    # vtk["sigmax"]    = sigmax_vtk
    # vtk["sigmay"]    = sigmay_vtk
end

# t_arr = []
# for file in readdir()
#     if length(file) >= length("N=$N,K1D=$K1D,") && file[1:length("N=$N,K1D=$K1D,")] == "N=$N,K1D=$K1D," && file[end-3:end] == ".txt"
#         idx = findnext(",",file,length("N=$N,K1D=$K1D,t="))
#         t = parse(Float64,file[length("N=$N,K1D=$K1D,t=")+1:idx[1]-1])
#         if !(t in t_arr)
#             push!(t_arr,t)
#         end
#     end
# end

# pvd = paraview_collection("N=$N,K1D=$K1D,CFL=$CFL,x=$XLENGTH")

# for t in t_arr
#     vtk_grid("N=$N,K1D=$K1D,t=$t,CFL=$CFL,x=$XLENGTH",Vp*x,Vp*y) do vtk
#         vtk["rho"] = Vp*readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,x=$XLENGTH,rho,dmr.txt")
#         vtk["rhou"] = Vp*readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,x=$XLENGTH,rhou,dmr.txt")
#         vtk["rhov"] = Vp*readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,x=$XLENGTH,rhov,dmr.txt")
#         vtk["E"] = Vp*readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,x=$XLENGTH,E,dmr.txt")
#         pvd[t] = vtk
#     end
# end

# vtk_save(pvd)

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
# load required tensor flow functionality

ENV["PYTHON"] = "/usr/bin/python3.6"

using Pkg
Pkg.build("PyCall")

using PyCall

# tf = pyimport("tensorflow")
# tf.keras.models.load_model

py"""
import sys
sys.path.insert(0, "/home/singhk/learning-odometry/")
"""

# @show tfmodeldir = joinpath(ENV["HOME"],"learning-odometry/")
# Base.cd(tfmodeldir)

PyTFOdoPredictorPoint2 = pyimport("queryAllModels")."predictSinglePose" # was ["predictSinglePose"]

# test that its working
A = [rand(4) for i in 1:25]
PyTFOdoPredictorPoint2(A)
# returns  10×3 Array{Real,2}:

# mmodel = tf.keras.models.load_model(joinpath(tfmodeldir, "original_network0.h5"))
# mmodel.predict
# mmodel.predict([[A]])
# 1×2 Array{Float32,2}:
#  0.0169379  1.58055


text:
##### Eddy-diffusivity models

abstract type EddyDiffusivityModel end

struct SCAMPyEddyDiffusivity{FT} <: EddyDiffusivityModel
    tke_ed_coeff::FT
end

"""
    compute_eddy_diffusivities_tke!

Define eddy-diffusivity fields

 - `aux[:K_m, k, i]`
 - `aux[:K_h, k, i]`

 for all `k` and all `i`
"""
function compute_eddy_diffusivities_tke! end

function compute_eddy_diffusivities_tke!(
    grid::Grid{FT},
    q,
    aux,
    params,
    model::SCAMPyEddyDiffusivity,
) where {FT}
    gm, en, ud, sd, al = allcombinations(q)
    @inbounds for k in over_elems_real(grid)
        l_mix = aux[:l_mix, k, gm]
        K_m_k = model.tke_ed_coeff * l_mix * sqrt(max(q[:tke, k, en], FT(0)))
        aux[:K_m, k, gm] = K_m_k
        aux[:K_h, k, gm] = K_m_k / params[:prandtl_number]
    end
end


text:
module SplineHazard

using StatsBase, Statistics, Dierckx, Distributions, DataFrames, QuadGK

using LinearAlgebra: rank, cholesky!, diagm

import Base.summary

export Sampler, Prior, Param, Tuner

export update!,
    create_sampler,
    setup_sampler,
    set_initial_state!,
    sample!,
    sample_prior,
    set_tuner!,
    extract_spline,
    extract_hazard,
    extract_cumhaz,
    extract_survival,
    summary

include("Spline.jl")

using .Spline

include("tuner.jl")
include("sampler.jl")
include("update.jl")

"""
A Julia package for fitting spline based non- and semi-parametric hazard models.
"""
SplineHazard

end


text:
# Dictionary with resnames for ionic liquids

export dict
dict = Dict( "DCA" => "NC" ,
             "BF4" => "BF4",
             "NO3" => "NO3",
 	         "SCN" => "SCN",
	         "TFO" => "TFO",
             "Cl " => "Cl ",
	     	 "MS " => "MS ",
			  "Cl" => "Cl",
	     	 "EMIM" => "EMI",
 	     	 "BMIM" => "BMI",
	     	 "NTF2" => "NSC",
 	     	 "Br"   => "Br",
	     	 "MS"   => "MS",
	     	 "PF6"  => "PF6",
	     	 "TCM"  => "TCM")


text:
module GeoInterfaceRFC

include("types.jl")
include("interface.jl")
include("defaults.jl")
# include("primitives.jl")  # needs rethinking
include("utils.jl")

export test_interface_for_geom

export geomtype
export ncoord
export getcoord
export ngeom
export getgeom

end # module


text:
@testset "extrema_polynomial" begin
    # Define a number of problems for which we expect the method to
    # work without issues. The problems are defined as (p, (a, b), x,
    # y). Where the minimum is attained at `x` and the maximum at `y`.
    problems = [
        # Min and max at interior points
        (
            Arblib.integral(Arblib.fromroots(ArbPoly, Arb[0.25, 1.75])),
            (Arf(0), Arf(2)),
            Arb(1.75),
            Arb(0.25),
        ),

        # Min at interior point, max at endpoint
        (ArbPoly([0, 0, 1]), (Arf(-1), Arf(2)), Arb(0), Arb(2)),
        (
            Arblib.integral(Arblib.fromroots(ArbPoly, [-1, 0, 1])),
            (Arf(-1.5), Arf(2)),
            Arb(1),
            Arb(2),
        ),
        (
            Arblib.integral(Arblib.fromroots(ArbPoly, [0, 1, 4])),
            (Arf(-1), Arf(6)),
            Arb(4),
            Arb(6),
        ),

        # Min at endpoint, max at interior point
        (ArbPoly([0, 0, -1]), (Arf(-1), Arf(2)), Arb(2), Arb(0)),
        (
            -Arblib.integral(Arblib.fromroots(ArbPoly, [-1, 0, 1])),
            (Arf(-1.5), Arf(2)),
            Arb(2),
            Arb(1),
        ),
        (
            -Arblib.integral(Arblib.fromroots(ArbPoly, [0, 1, 4])),
            (Arf(-1), Arf(6)),
            Arb(6),
            Arb(4),
        ),

        # Min and max at endpoints
        (ArbPoly([0, 1]), (Arf(-1), Arf(1)), Arb(-1), Arb(1)),
        (ArbPoly([1, -1, 3, -1]), (Arf(-3), Arf(3)), Arb(3), Arb(-3)),

        # Derivative has double root so cant isolate roots
        (
            Arblib.integral(Arblib.fromroots(ArbPoly, [0.25, 0.25, 1.75, 1.75])),
            (Arf(0), Arf(2)),
            Arb(0),
            Arb(2),
        ),

        # Special Cases

        # Thin interval
        (ArbPoly([1, 1, 1]), (Arf(1), Arf(1)), Arb(1), Arb(1)),

        # Constant polynomial
        (ArbPoly(5), (Arf(0), Arf(1)), Arb(0.5), Arb(0.5)),

        # Linear polynomial
        (ArbPoly([5, 1]), (Arf(0), Arf(1)), Arb(0), Arb(1)),
        (ArbPoly([5, -1]), (Arf(0), Arf(1)), Arb(1), Arb(0)),
    ]

    i = 0
    for (p, (a, b), x, y) in problems
        min, max = p(x), p(y)
        min1, max1 = ArbExtras.extrema_polynomial(p, a, b)
        min2 = ArbExtras.minimum_polynomial(p, a, b)
        max2 = ArbExtras.maximum_polynomial(p, a, b)

        @test isequal(min1, min2)
        @test isequal(max1, max2)

        @test Arblib.overlaps(min, min1)
        @test Arblib.overlaps(max, max1)

        @test Arblib.radius(min1) < sqrt(eps(Arb))
        @test Arblib.radius(max1) < sqrt(eps(Arb))
    end

    # Polynomial with NaN
    @test all(isnan, ArbExtras.extrema_polynomial(ArbPoly(NaN), Arf(0), Arf(1)))
    @test isnan(ArbExtras.minimum_polynomial(ArbPoly(NaN), Arf(0), Arf(1)))
    @test isnan(ArbExtras.maximum_polynomial(ArbPoly(NaN), Arf(0), Arf(1)))

    # TODO: Isolate roots fails
    #(Arblib.integral(Arblib.fromroots(ArbPoly, [0, 1, 2, 3, 4])), (Arf(-1), Arf(6)), Arb(0), Arb(6)),

    @testset "abs_value = true" begin
        # Define a number of problems for which we expect the method to
        # work without issues. The problems are defined as (p, (a, b), x,
        # y). Where the minimum is attained at `x` and the maximum at `y`.
        problems = [
            # Min and max at interior points
            (
                Arblib.integral(Arblib.fromroots(ArbPoly, [-1, 1])) + 2,
                (Arf(-1.5), Arf(1.5)),
                Arf(1),
                Arf(-1),
            ),
            (ArbPoly([-1, 0, 1]), (Arf(-1.1), Arf(1.2)), Arf(1), Arf(0)),

            # Min is non-zero at interior point, max at endpoint
            (ArbPoly([1, 0, 1]), (Arf(-1), Arf(2)), Arf(0), Arf(2)),
            (
                Arblib.integral(Arblib.fromroots(ArbPoly, [-0.5, 0, 0.5])) - 9 // 1024,
                (Arf(-0.9), Arf(1)),
                Arf(0.75),
                Arf(1),
            ),

            # Min is zero at interior point, max at endpoint
            (ArbPoly([-1, 0, 1]), (Arf(0), Arf(2)), Arf(1), Arf(2)),
            (ArbPoly([-1, 0, 1]), (Arf(-2), Arf(2)), Arf(1), Arf(2)),

            # Min at endpoint, max at interior point
            (ArbPoly([5, 0, -1]), (Arf(-1), Arf(2)), Arf(2), Arf(0)),

            # Min and max at endpoints
            (
                Arblib.integral(Arblib.fromroots(ArbPoly, [-1, 1])) + 2,
                (Arf(-2.2), Arf(2.2)),
                Arf(-2.2),
                Arf(2.2),
            ),

            # Derivative has double root so cant isolate roots
            (ArbPoly([5, 0, 0, 1]), (Arf(-1), Arf(2)), Arb(-1), Arb(2)),

            # Special Cases

            # Thin interval
            (ArbPoly([1, 1, 1]), (Arf(1), Arf(1)), Arb(1), Arb(1)),

            # Constant polynomial
            (ArbPoly(5), (Arf(0), Arf(1)), Arb(0.5), Arb(0.5)),

            # Linear polynomial
            (ArbPoly([0, 1]), (Arf(-2), Arf(1)), Arb(0), Arb(-2)),
            (ArbPoly([0, 1]), (Arf(-1), Arf(2)), Arb(0), Arb(2)),
            (ArbPoly([5, -1]), (Arf(0), Arf(1)), Arb(1), Arb(0)),
            (ArbPoly([-5, -1]), (Arf(0), Arf(1)), Arb(0), Arb(1)),

            # Zero at endpoint
            (ArbPoly([0, 1, 1]), (Arf(0), Arf(2)), Arb(0), Arb(2)),
        ]

        i = 0
        for (p, (a, b), x, y) in problems
            min, max = abs(p(x)), abs(p(y))
            # Run the problems with both p and -p, they should give
            # identical results
            min11, max11 = ArbExtras.extrema_polynomial(p, a, b, abs_value = true)
            min12, max12 = ArbExtras.extrema_polynomial(-p, a, b, abs_value = true)
            min21 = ArbExtras.minimum_polynomial(p, a, b, abs_value = true)
            min22 = ArbExtras.minimum_polynomial(-p, a, b, abs_value = true)
            max21 = ArbExtras.maximum_polynomial(p, a, b, abs_value = true)
            max22 = ArbExtras.maximum_polynomial(-p, a, b, abs_value = true)

            @test isequal(min11, min12)
            @test isequal(min11, min21)
            @test isequal(min11, min22)
            @test isequal(max11, max12)
            @test isequal(max11, max21)
            @test isequal(max11, max22)

            @test Arblib.overlaps(min, min11)
            @test Arblib.overlaps(max, max11)

            @test Arblib.radius(min11) < sqrt(eps(Arb))
            @test Arblib.radius(max11) < sqrt(eps(Arb))
        end

        # Polynomial with NaN
        @test all(
            isnan,
            ArbExtras.extrema_polynomial(ArbPoly(NaN), Arf(0), Arf(1), abs_value = true),
        )
        @test isnan(
            ArbExtras.minimum_polynomial(ArbPoly(NaN), Arf(0), Arf(1), abs_value = true),
        )
        @test isnan(
            ArbExtras.maximum_polynomial(ArbPoly(NaN), Arf(0), Arf(1), abs_value = true),
        )
    end
end


text:
#---------------------------------------------------------------------#
#This code computes the element mass and differentiation matrices
#Written by F.X. Giraldo on April 24, 2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function element_matrices(ψ,dψ,ξ_x,ξ_y,η_x,η_y,jac,ωq,Ne,Np,Nq,DFloat)

    #Initialize Matrices
    M=zeros(DFloat,Np,Np,Np,Np,Ne)
    D_x=zeros(DFloat,Np,Np,Np,Np,Ne)
    D_y=zeros(DFloat,Np,Np,Np,Np,Ne)

    #Construct Mass and Differentiation Matrices
    for e=1:Ne
        for l=1:Nq, k=1:Nq
            wq=ωq[k]*ωq[l]*jac[k,l,e]
            for j=1:Np, i=1:Np
                Ψ_JK=ψ[i,k]*ψ[j,l] #h_ik*h_jl
                dΨdx_JK=dψ[i,k]*ψ[j,l]*ξ_x[k,l,e] + ψ[i,k]*dψ[j,l]*η_x[k,l,e]
                dΨdy_JK=dψ[i,k]*ψ[j,l]*ξ_y[k,l,e] + ψ[i,k]*dψ[j,l]*η_y[k,l,e]
                for n=1:Np, m=1:Np
                    Ψ_IK=ψ[m,k]*ψ[n,l] #h_ik*h_jl
                    dΨdx_IK=dψ[m,k]*ψ[n,l]*ξ_x[k,l,e] + ψ[m,k]*dψ[n,l]*η_x[k,l,e]
                    dΨdy_IK=dψ[m,k]*ψ[n,l]*ξ_y[k,l,e] + ψ[m,k]*dψ[n,l]*η_y[k,l,e]
                    M[m,n,i,j,e]+=wq*Ψ_IK*Ψ_JK
                    D_x[m,n,i,j,e]+=wq*dΨdx_IK*Ψ_JK
                    D_y[m,n,i,j,e]+=wq*dΨdy_IK*Ψ_JK
                end #m,n
            end #j,i
        end #k,l
    end #e

    return (M,D_x,D_y)
end


text:
"""
$(TYPEDEF)

Grid component key type for storing node in parent array
"""
abstract type NodeInParent <: AbstractGridIntegerArray1D end

"""
$(TYPEDEF)

Grid component key type for storing parent grid
"""
abstract type ParentGrid <: AbstractGridComponent end

"""
$(TYPEDSIGNATURES)

Default transform for subgrid creation
"""
function _copytransform!(a::AbstractArray,b::AbstractArray)
    for i=1:length(a)
        a[i]=b[i]
    end
end

"""
$(TYPEDSIGNATURES)

Create subgrid from list of regions.

- `parent`: parent grid
- `subregions`:  Array of subregions
- `transform` (kw parameter): transformation function between
   grid and subgrid coordinates acting on one point.
   Default: `copytransform`
- `boundary`: if true, create codimension 1 subgrid from boundary region.

A subgrid is of type `ExtendableGrid` and stores two additional components:
[`ParentGrid`](@ref) and [`NodeInParent`](@ref)

"""
function subgrid(parent,
                 subregions::AbstractArray;
                 transform::Function=_copytransform!,
                 boundary=false,
                 project=true)

    Tc=coord_type(parent)
    Ti=index_type(parent)

    @inline function insubregions(xreg)
        for i in eachindex(subregions)
            if subregions[i]==xreg
                return true
            end
        end
        return false
    end

    if boundary
        xregions=parent[BFaceRegions]
        xnodes=parent[BFaceNodes]
        sub_gdim=dim_grid(parent)-1
        xct=parent[BFaceGeometries]
        sub_gdim=dim_grid(parent)-1
    else
        xregions=parent[CellRegions]
        xnodes=parent[CellNodes]
        xct=parent[CellGeometries]
        sub_gdim=dim_grid(parent)
    end

    nodemark=zeros(Ti,num_nodes(parent))

    nsubcells=0
    nsubnodes=0
    for icell in eachindex(xregions)
        if insubregions(xregions[icell])
            nsubcells+=1
            for inode=1:num_targets(xnodes,icell)
                ipnode=xnodes[inode,icell]
                if nodemark[ipnode]==0
                    nsubnodes+=1
                    nodemark[ipnode]=nsubnodes
                end
            end
        end
    end

    sub_xnodes=VariableTargetAdjacency(Ti)
    sub_nip=zeros(Ti,nsubnodes)
    sub_ct=Vector{DataType}(undef,0)
    sub_cr=Vector{Ti}(undef,0)
    for inode in eachindex(nodemark)
        if nodemark[inode]>0
            sub_nip[nodemark[inode]]=inode
        end
    end

    isubcell=0
    for icell in eachindex(xregions)
        if insubregions(xregions[icell])
            ncn=num_targets(xnodes,icell)
            col=zeros(Ti,0)
            for inode=1:ncn
                push!(col,nodemark[xnodes[inode,icell]])
            end
            append!(sub_xnodes,col)
            push!(sub_ct,xct[icell])
            push!(sub_cr,xregions[icell])
        end
    end

    if project
        sub_coord=zeros(Tc,sub_gdim,nsubnodes)
    else
        sub_coord=zeros(Tc,dim_space(parent),nsubnodes)
    end
    coord=parent[Coordinates]
    @views for inode=1:nsubnodes
        transform(sub_coord[:,inode],coord[:,sub_nip[inode]])
    end
    subgrid=ExtendableGrid{Tc,Ti}()
    subgrid[Coordinates]=sub_coord
    subgrid[CellRegions]=sub_cr
    subgrid[CellGeometries]=sub_ct
    subgrid[CellNodes]=tryfix(sub_xnodes)
    subgrid[ParentGrid]=parent
    subgrid[NodeInParent]=sub_nip

    if boundary
        subgrid[NumBFaceRegions]=0
    else
        bfacenodes=parent[BFaceNodes]
        bfaceregions=parent[BFaceRegions]
        bfacetypes=parent[BFaceGeometries]

        sub_bfacenodes=VariableTargetAdjacency(Ti)
        sub_bfaceregions=Vector{Ti}(undef,0)
        sub_bfacetypes=Vector{DataType}(undef,0)

        for ibface in eachindex(bfaceregions)
            nbn=num_targets(bfacenodes,ibface)
            insubgrid=true
            for inode=1:nbn
                if nodemark[bfacenodes[inode,ibface]]==0
                    insubgrid=false
                    continue
                end
            end
            if insubgrid
                col=zeros(Ti,0)
                for inode=1:nbn
                    push!(col,nodemark[bfacenodes[inode,ibface]])
                end
                append!(sub_bfacenodes,col)
                push!(sub_bfacetypes,bfacetypes[ibface])
                push!(sub_bfaceregions,bfaceregions[ibface])
            end
        end

        subgrid[BFaceRegions]=sub_bfaceregions
        subgrid[BFaceGeometries]=sub_bfacetypes
        subgrid[BFaceNodes]=tryfix(sub_bfacenodes)
        subgrid[NumBFaceRegions]=maximum(sub_bfaceregions)
    end

    subgrid
end

"""
$(TYPEDEF)
Vector view on subgrid

$(TYPEDFIELDS)
"""
struct SubgridVectorView{Tv,Ti} <: AbstractVector{Tv}
    sysarray::AbstractVector{Tv}
    node_in_parent::Vector{Ti}
end

##################################################################
"""
$(TYPEDSIGNATURES)

Create a view of the vector on a subgrid.
"""
Base.view(a::AbstractVector,subgrid::ExtendableGrid)  = SubgridVectorView(a,subgrid[NodeInParent])

##############################################################################
"""
$(TYPEDSIGNATURES)

Accessor method for subgrid vector view.
"""
Base.getindex(aview::SubgridVectorView,inode::Integer) = aview.sysarray[aview.node_in_parent[inode]]

##############################################################################
"""
$(TYPEDSIGNATURES)

Accessor method for subgrid vector view.
"""
@inline function Base.setindex!(aview::SubgridVectorView,v,inode::Integer)
    aview.sysarray[aview.node_in_parent[inode]]=v
    return aview
end

##################################################################
"""
$(TYPEDSIGNATURES)

Return size of vector view.
"""
Base.size(a::SubgridVectorView)=(size(a.node_in_parent,1),)


text:
function generate_cartesian_orbit(obj::OPMFile)
    return jcall(obj, "generateCartesianOrbit", CartesianOrbit, ())
end

function generate_keplerian_orbit(obj::OPMFile)
    return jcall(obj, "generateKeplerianOrbit", KeplerianOrbit, ())
end

function generate_spacecraft_state(obj::OPMFile)
    return jcall(obj, "generateSpacecraftState", SpacecraftState, ())
end

function get_has_maneuver(obj::OPMFile)
    return jcall(obj, "getHasManeuver", jboolean, ())
end

function get_maneuver(obj::OPMFile, arg0::jint)
    return jcall(obj, "getManeuver", OPMFile_Maneuver, (jint,), arg0)
end

function get_maneuvers(obj::OPMFile)
    return jcall(obj, "getManeuvers", List, ())
end

function get_meta_data(obj::OPMFile)
    return jcall(obj, "getMetaData", ODMMetaData, ())
end

function get_meta_data_comment(obj::OPMFile)
    return jcall(obj, "getMetaDataComment", List, ())
end

function get_nb_maneuvers(obj::OPMFile)
    return jcall(obj, "getNbManeuvers", jint, ())
end

function get_position(obj::OPMFile)
    return jcall(obj, "getPosition", Vector3D, ())
end

function get_pv_coordinates(obj::OPMFile)
    return jcall(obj, "getPVCoordinates", PVCoordinates, ())
end

function get_velocity(obj::OPMFile)
    return jcall(obj, "getVelocity", Vector3D, ())
end


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

mutable struct ConnectionStateSnapshot <: SwaggerModel
    connectionState::Any # spec type: Union{ Nothing, String } # spec name: connectionState
    startTime::Any # spec type: Union{ Nothing, DateTime } # spec name: startTime
    endTime::Any # spec type: Union{ Nothing, DateTime } # spec name: endTime
    evaluationState::Any # spec type: Union{ Nothing, String } # spec name: evaluationState
    avgLatencyInMs::Any # spec type: Union{ Nothing, Int32 } # spec name: avgLatencyInMs
    minLatencyInMs::Any # spec type: Union{ Nothing, Int32 } # spec name: minLatencyInMs
    maxLatencyInMs::Any # spec type: Union{ Nothing, Int32 } # spec name: maxLatencyInMs
    probesSent::Any # spec type: Union{ Nothing, Int32 } # spec name: probesSent
    probesFailed::Any # spec type: Union{ Nothing, Int32 } # spec name: probesFailed
    hops::Any # spec type: Union{ Nothing, Vector{ConnectivityHop2} } # spec name: hops

    function ConnectionStateSnapshot(;connectionState=nothing, startTime=nothing, endTime=nothing, evaluationState=nothing, avgLatencyInMs=nothing, minLatencyInMs=nothing, maxLatencyInMs=nothing, probesSent=nothing, probesFailed=nothing, hops=nothing)
        o = new()
        validate_property(ConnectionStateSnapshot, Symbol("connectionState"), connectionState)
        setfield!(o, Symbol("connectionState"), connectionState)
        validate_property(ConnectionStateSnapshot, Symbol("startTime"), startTime)
        setfield!(o, Symbol("startTime"), startTime)
        validate_property(ConnectionStateSnapshot, Symbol("endTime"), endTime)
        setfield!(o, Symbol("endTime"), endTime)
        validate_property(ConnectionStateSnapshot, Symbol("evaluationState"), evaluationState)
        setfield!(o, Symbol("evaluationState"), evaluationState)
        validate_property(ConnectionStateSnapshot, Symbol("avgLatencyInMs"), avgLatencyInMs)
        setfield!(o, Symbol("avgLatencyInMs"), avgLatencyInMs)
        validate_property(ConnectionStateSnapshot, Symbol("minLatencyInMs"), minLatencyInMs)
        setfield!(o, Symbol("minLatencyInMs"), minLatencyInMs)
        validate_property(ConnectionStateSnapshot, Symbol("maxLatencyInMs"), maxLatencyInMs)
        setfield!(o, Symbol("maxLatencyInMs"), maxLatencyInMs)
        validate_property(ConnectionStateSnapshot, Symbol("probesSent"), probesSent)
        setfield!(o, Symbol("probesSent"), probesSent)
        validate_property(ConnectionStateSnapshot, Symbol("probesFailed"), probesFailed)
        setfield!(o, Symbol("probesFailed"), probesFailed)
        validate_property(ConnectionStateSnapshot, Symbol("hops"), hops)
        setfield!(o, Symbol("hops"), hops)
        o
    end
end # type ConnectionStateSnapshot

const _property_map_ConnectionStateSnapshot = Dict{Symbol,Symbol}(Symbol("connectionState")=>Symbol("connectionState"), Symbol("startTime")=>Symbol("startTime"), Symbol("endTime")=>Symbol("endTime"), Symbol("evaluationState")=>Symbol("evaluationState"), Symbol("avgLatencyInMs")=>Symbol("avgLatencyInMs"), Symbol("minLatencyInMs")=>Symbol("minLatencyInMs"), Symbol("maxLatencyInMs")=>Symbol("maxLatencyInMs"), Symbol("probesSent")=>Symbol("probesSent"), Symbol("probesFailed")=>Symbol("probesFailed"), Symbol("hops")=>Symbol("hops"))
const _property_types_ConnectionStateSnapshot = Dict{Symbol,String}(Symbol("connectionState")=>"String", Symbol("startTime")=>"DateTime", Symbol("endTime")=>"DateTime", Symbol("evaluationState")=>"String", Symbol("avgLatencyInMs")=>"Int32", Symbol("minLatencyInMs")=>"Int32", Symbol("maxLatencyInMs")=>"Int32", Symbol("probesSent")=>"Int32", Symbol("probesFailed")=>"Int32", Symbol("hops")=>"Vector{ConnectivityHop2}")
Base.propertynames(::Type{ ConnectionStateSnapshot }) = collect(keys(_property_map_ConnectionStateSnapshot))
Swagger.property_type(::Type{ ConnectionStateSnapshot }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ConnectionStateSnapshot[name]))}
Swagger.field_name(::Type{ ConnectionStateSnapshot }, property_name::Symbol) =  _property_map_ConnectionStateSnapshot[property_name]

const _allowed_ConnectionStateSnapshot_connectionState = ["Reachable", "Unreachable", "Unknown"]

const _allowed_ConnectionStateSnapshot_evaluationState = ["NotStarted", "InProgress", "Completed"]

function check_required(o::ConnectionStateSnapshot)
    true
end

function validate_property(::Type{ ConnectionStateSnapshot }, name::Symbol, val)
    if name === Symbol("connectionState")
        Swagger.validate_param(name, "ConnectionStateSnapshot", :enum, val, _allowed_ConnectionStateSnapshot_connectionState)
    end
    if name === Symbol("evaluationState")
        Swagger.validate_param(name, "ConnectionStateSnapshot", :enum, val, _allowed_ConnectionStateSnapshot_evaluationState)
    end
end


text:
module log

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
function hello2()
    println("hello")
end

function cpa()
    println("cpa")
end

function add2()
    println("")
end


text:
@doc doc"""
    Circle{F} <: Manifold

The circle $\mathbb S^1$ as a manifold ere manifold represented by
real-valued data in $[-\pi,\pi)$ or complex-valued data $z\in \mathbb C$ of absolute value
$\lvert z\rvert = 1$.
# Constructor

    Circle(f=ℝ)

Generate the `ℝ`-valued Circle represented by angles, which
alternatively can be set to use the [`AbstractNumbers`](@ref) `f=ℂ` to obtain the `Circle`
represented by `ℂ`-valued `Circle` of unit numbers.
"""
struct Circle{F} <: Manifold where {F<:AbstractNumbers} end

Circle(f::AbstractNumbers = ℝ) = Circle{f}()

@doc doc"""
    check_manifold_point(M::Circle, x)

Check whether `x` is a point on the [`Circle`](@ref) `M`.
For the real-valued case, `x` is an angle and hence it checks that $x \in [-\pi,\pi)$.
for the complex-valued case its a unit number, $x \in \mathbb C$ with $\lvert x \rvert = 1$.
"""
check_manifold_point(::Circle, ::Any...)

function check_manifold_point(M::Circle{ℝ}, x; kwargs...)
    if !isapprox(sym_rem(x), x; kwargs...)
        return DomainError(
            x,
            "The point $(x) does not lie on $(M), since its is not in [-π,π).",
        )
    end
    return nothing
end
function check_manifold_point(M::Circle{ℂ}, x; kwargs...)
    if !isapprox(sum(abs.(x)), 1.0; kwargs...)
        return DomainError(
            abs(x),
            "The point $(x) does not lie on the $(M) since its norm is not 1.",
        )
    end
    return nothing
end

"""
    check_tangent_vector(M::Circle, x, v)

Check whether `v` is a tangent vector in the tangent space of `x` on the
[`Circle`](@ref) `M`.
For the real-valued case represented by angles all `v` are valid, since the tangent space is
the whole real line.
For the complex-valued case `v` has to lie on the line parallel to the tangent line at `x`
in the complex plane, i.e. the inner product is zero.
"""
check_tangent_vector(::Circle{ℝ}, ::Any...; ::Any...)

function check_tangent_vector(M::Circle{ℝ}, x, v; kwargs...)
    perr = check_manifold_point(M, x)
    return perr # if x is valid all v that are real numbers are valid
end
function check_tangent_vector(M::Circle{ℂ}, x, v; kwargs...)
    perr = check_manifold_point(M, x)
    perr === nothing || return perr
    if !isapprox(abs(complex_dot(x, v)), 0.0; kwargs...)
        return DomainError(
            abs(complex_dot(x, v)),
            "The value $(v) is not a tangent vector to $(x) on $(M), since it is not orthogonal in the embedding.",
        )
    end
    return nothing
end

@doc doc"""
    complex_dot(a, b)

Compute the inner product of two (complex) numbers with in the complex plane.
"""
complex_dot(a, b) = dot(map(real, a), map(real, b)) + dot(map(imag, a), map(imag, b))
complex_dot(a::Number, b::Number) = (real(a) * real(b) + imag(a) * imag(b))

@doc doc"""
    distance(M::Circle, x, y)

Compute the distance on the [`Circle`](@ref) `M`, which is
the absolute value of the symmetric remainder of `x` and `y` for the real-valued
case and the angle between both complex numbers in the Gaussian plane for the
complex-valued case.
"""
distance(::Circle, ::Any...)
distance(::Circle{ℝ}, x::Real, y::Real) = abs(sym_rem(x - y))
distance(::Circle{ℝ}, x, y) = abs(sum(sym_rem.(x - y)))
distance(::Circle{ℂ}, x, y) = acos(clamp(complex_dot(x, y), -1, 1))

@doc doc"""
    exp(M::Circle, x, v)

Compute the exponential map on the [`Circle`](@ref).
````math
\exp_xv = (x+v)_{2\pi},
````
where $(\cdot)$ is the (symmetric) remainder with respect to division by $2\pi$,
i.e. in $[-\pi,\pi)$.

For the complex-valued case the formula is the same as for the [`Sphere`](@ref)
applied to valuedin the complex plane.
"""
exp(::Circle, ::Any...)
exp(::Circle{ℝ}, x::Real, v::Real) = sym_rem(x + v)
function exp(M::Circle{ℂ}, x::Number, v::Number)
    θ = norm(M, x, v)
    return cos(θ) * x + usinc(θ) * v
end

exp!(::Circle{ℝ}, y, x, v) = (y .= sym_rem(x + v))
function exp!(M::Circle{ℂ}, y, x, v)
    θ = norm(M, x, v)
    y .= cos(θ) * x + usinc(θ) * v
    return y
end

flat(M::Circle, x::Number, w::TFVector) = FVector(CotangentSpace, w.data)

flat!(::Circle, v::CoTFVector, x, w::TFVector) = copyto!(v, w)

function get_basis(M::Circle{ℝ}, x, B::DiagonalizingOrthonormalBasis)
    sbv = sign(B.v[1])
    vs = @SVector [@SVector [sbv == 0 ? one(sbv) : sbv]]
    return PrecomputedDiagonalizingOrthonormalBasis(vs, @SVector [0])
end

get_coordinates(M::Circle{ℝ}, x, v, B::ArbitraryOrthonormalBasis) = v
function get_coordinates(M::Circle{ℝ}, x, v, B::DiagonalizingOrthonormalBasis)
    sbv = sign(B.v[1])
    return v .* (sbv == 0 ? 1 : sbv)
end
"""
    get_coordinates(M::Circle{ℂ}, x, v, B::ArbitraryOrthonormalBasis)

Return tangent vector coordinates in the Lie algebra of the circle.
"""
function get_coordinates(M::Circle{ℂ}, x, v, B::ArbitraryOrthonormalBasis)
    v, x = v[1], x[1]
    w = imag(v) * real(x) - real(v) * imag(x)
    return @SVector [w]
end

get_vector(M::Circle{ℝ}, x, v, B::ArbitraryOrthonormalBasis) = v
function get_vector(M::Circle{ℝ}, x, v, B::DiagonalizingOrthonormalBasis)
    sbv = sign(B.v[1])
    return v .* (sbv == 0 ? 1 : sbv)
end
"""
    get_vector(M::Circle{ℂ}, x, v, B::ArbitraryOrthonormalBasis)

Return tangent vector from the coordinates in the Lie algebra of the circle.
"""
get_vector(M::Circle{ℂ}, x, v, B::ArbitraryOrthonormalBasis) = @SVector [1im * v[1] * x[1]]

@doc doc"""
    injectivity_radius(M::Circle[, x])

Return the injectivity radius on the [`Circle`](@ref) `M`, i.e. $\pi$.
"""
injectivity_radius(::Circle, args...) = π

@doc doc"""
    inner(M::Circle, x, w, v)

Compute the inner product of the two tangent vectors `w,v` from the tangent plane at `x` on
the [`Circle`](@ref) `M` using the restriction of the metric from the embedding,
i.e.

````math
g_x(v,w) = w*v
````

for the real case and

    ````math
g_x(v,w) = v^\mathrm{T}w
````

for the complex case interpreting complex numbers in the Gaussian plane.
"""
inner(::Circle, ::Any...)
@inline inner(::Circle{ℝ}, x, w, v) = dot(v, w)
@inline inner(::Circle{ℝ}, x::Real, w::Real, v::Real) = v * w
@inline inner(::Circle{ℂ}, x, w, v) = complex_dot(w, v)

function inverse_retract(M::Circle, x::Number, y::Number)
    return inverse_retract(M, x, y, LogarithmicInverseRetraction())
end
function inverse_retract(M::Circle, x::Number, y::Number, ::LogarithmicInverseRetraction)
    return log(M, x, y)
end

@doc doc"""
    log(M::Circle, x, y)

Compute the logarithmic map on the [`Circle`](@ref) `M`.
````math
\exp_xv = (y,x)_{2\pi},
````
where $(\cdot)$ is the (symmetric) remainder with respect to division by $2\pi$,
i.e. in $[-\pi,\pi)$.

For the complex-valued case the formula is the same as for the [`Sphere`](@ref)
applied to valuedin the complex plane.
"""
log(::Circle, ::Any...)
log(::Circle{ℝ}, x::Real, y::Real) = sym_rem(y - x)
function log(M::Circle{ℂ}, x::Number, y::Number)
    cosθ = complex_dot(x, y)
    if cosθ ≈ -1  # appr. opposing points, return deterministic choice from set-valued log
        v = real(x) ≈ 1 ? 1im : 1 + 0im
        v = v - complex_dot(x, v) * x
        v *= π / norm(v)
    else
        cosθ = cosθ > 1 ? one(cosθ) : cosθ
        θ = acos(cosθ)
        v = (y - cosθ * x) / usinc(θ)
    end
    return project_tangent(M, x, v)
end

log!(::Circle{ℝ}, v, x, y) = (v .= sym_rem(y - x))
function log!(M::Circle{ℂ}, v, x, y)
    cosθ = complex_dot(x, y)
    if cosθ ≈ -1
        v .= sum(real.(x)) ≈ 1 ? 1.0im : 1.0 + 0.0im
        v .= v - complex_dot(x, v) * x
        v .*= π / norm(v)
    else
        cosθ = cosθ > 1 ? one(cosθ) : cosθ
        θ = acos(cosθ)
        v .= (y - cosθ * x) / usinc(θ)
    end
    return project_tangent!(M, v, x, v)
end

@doc doc"""
    manifold_dimension(M::Circle)

Return the dimension of the [`Circle`](@ref) `M`,
i.e. $\operatorname{dim}(\mathbb S^1) = 1$.
"""
manifold_dimension(::Circle) = 1

@doc doc"""
    mean(M::Circle, x::AbstractVector[, w::AbstractWeights])

Compute the Riemannian [`mean`](@ref mean(M::Manifold, args...)) of `x` on the
[`Circle`](@ref) $\mathbb S^1$ by the wrapped mean, i.e. the remainder of the
mean modulo 2π.
"""
mean(::Circle, ::Any)
mean(::Circle, x::Array{<:Real}; kwargs...) = sym_rem(sum(x))
mean(::Circle, x::Array{<:Real}, w::AbstractVector; kwargs...) = sym_rem(sum(w .* x))

@inline norm(::Circle, x, v) = sum(abs, v)

@doc doc"""
    project_point(M::Circle, x)

Project a point `x` onto the [`Circle`](@ref) `M`.
For the real-valued case this is the remainder with respect to modulus $2\pi$.
For the complex-valued case the result is the projection of `x` onto the unit circle in the
complex plane.
"""
project_point(::Circle, ::Any)
project_point(::Circle{ℝ}, x::Real) = sym_rem(x)
project_point(::Circle{ℂ}, x::Number) = x / abs(x)

project_point!(::Circle{ℝ}, x) = (x .= sym_rem(x))
project_point!(::Circle{ℂ}, x) = (x .= x / sum(abs.(x)))

@doc doc"""
    project_tangent(M::Circle, x, v)

Project a value `v` onto the tangent space of the point `x` on the [`Circle`](@ref) `M`.

For the real-valued case this is just the identity.
For the complex valued case `v` is projected onto the line in the complex plane
that is parallel to the tangent to `x` on the unit circle and contains `0`.
"""
project_tangent(::Circle{ℝ}, x::Real, v::Real) = v
project_tangent(::Circle{ℂ}, x::Number, v::Number) = v - complex_dot(x, v) * x

project_tangent!(::Circle{ℝ}, w, x, v) = (w .= v)
project_tangent!(::Circle{ℂ}, w, x, v) = (w .= v - complex_dot(x, v) * x)

retract(M::Circle, x, y) = retract(M, x, y, ExponentialRetraction())
retract(M::Circle, x, y, m::ExponentialRetraction) = exp(M, x, y)

representation_size(::Circle) = ()

sharp(M::Circle, x::Number, w::CoTFVector) = FVector(TangentSpace, w.data)

sharp!(M::Circle, v::TFVector, x, w::CoTFVector) = copyto!(v, w)

@doc doc"""
    sym_rem(x,[T=π])

Compute symmetric remainder of `x` with respect to the interall 2*`T`, i.e.
`(x+T)%2T`, where the default for `T` is $\pi$
"""
function sym_rem(x::N, T = π) where {N<:Number}
    return (x ≈ T ? convert(N, -T) : rem(x, convert(N, 2 * T), RoundNearest))
end
sym_rem(x, T=π) where N = map(sym_rem, x, Ref(T))

@doc doc"""
    vector_transport_to(M::Circle, x, v, y, ::ParallelTransport)

Compute the parallel transport of `v` from the tangent space at `x` to the tangent space at
`y` on the [`Circle`](@ref) `M`.
For the real-valued case this results in the identity.
For the complex-valud case, the formula is the same as for the [`Sphere`](@ref)`(1)` in the
complex plane.
````math
P_{y\gets x}(v) = v - \frac{\langle \log_xy,v\rangle_x}{d^2_{\mathbb C}(x,y)}
\bigl(\log_xy + \log_yx \bigr),
````
where [`log`](@ref) denotes the logarithmic map on `M`.
"""
vector_transport_to(::Circle, ::Any, ::Any, ::Any, ::ParallelTransport)
vector_transport_to(::Circle{ℝ}, x::Real, v::Real, y::Real, ::ParallelTransport) = v
function vector_transport_to(
    M::Circle{ℂ},
    x::Number,
    v::Number,
    y::Number,
    ::ParallelTransport,
)
    v_xy = log(M, x, y)
    vl = norm(M, x, v_xy)
    vto = v
    if vl > 0
        factor = 2 * complex_dot(v, y) / (abs(x + y)^2)
        vto -= factor .* (x + y)
    end
    return vto
end

vector_transport_to!(::Circle{ℝ}, w, x, v, y, ::ParallelTransport) = (w .= v)
function vector_transport_to!(M::Circle{ℂ}, vto, x, v, y, ::ParallelTransport)
    v_xy = log(M, x, y)
    vl = norm(M, x, v_xy)
    vto .= v
    if vl > 0
        factor = 2 * complex_dot(v, y) / (sum(abs.(x + y) .^ 2))
        vto .-= factor .* (x + y)
    end
    return vto
end

function vector_transport_direction(
    M::Circle,
    x::Number,
    v::Number,
    vdir::Number,
    m::AbstractVectorTransportMethod,
)
    y = exp(M, x, vdir)
    return vector_transport_to(M, x, v, y, m)
end

zero_tangent_vector(::Circle, x::Number) = zero(x)

zero_tangent_vector!(::Circle, v, x) = fill!(v, 0)