text:
const tokenfile = joinpath(tempdir(), "24cb59ae-a5a4-11e9-0ec9-534fb6590556")
const SERVER = "http://3.9.208.37/graphql/"
const query = """
    {
      experiments {
        objectId
        description
        name
        runs {
          objectId
          factorLevels {
            factor {
              name
            }
            level {
              value
            }
            createdDate
          }
          runDate
        }
      }
    }
"""
tokenquery(username, password) = """
    mutation {
      tokenAuth(username: "$username", password: "$password") {
        token
      }
    }
"""
validityquery(token) = """
mutation {
  verifyToken(token: "$token") {
    payload
  }
}
"""
const username_dialog = Dialog()

function getcredentials()
    username = requestᵐ("Username: ", username_dialog)
    io = Base.getpass("Password")
    password = read(io, String)
    Base.shred!(io)
    username, password
end

function gettoken()
    if isfile(tokenfile)
        open(tokenfile) do io
            token = readline(io)
        end
        r = Queryclient(SERVER, validityquery(token))
        r.Info.status == 200 || error("couldn't verify token")
        j = JSON3.read(r.Data)
        if unix2datetime(j[:data][:verifyToken][:payload][:exp]) > now() + Minute(1)
            return token
        end
    end
    username, password = getcredentials()
    r = Queryclient(SERVER, tokenquery(username, password))
    r.Info.status == 200 || error("couldn't authenticate with server")
    j = JSON3.read(r.Data)
    token = j[:data][:tokenAuth][:token]
    open(tokenfile, "w") do io
        println(io, token)
    end
    return token
end

function getexperimentsruns(token)
    r = Queryclient(SERVER, query, auth = "JWT $token")
    r.Info.status == 200 || error("couldn't query server")
    r.Data
end

ext: jl


text:
module RuntimeGeneratedFunctions

using ExprTools, Serialization, SHA

export @RuntimeGeneratedFunction

"""
    RuntimeGeneratedFunction

This type should be constructed via the macro @RuntimeGeneratedFunction.
"""
struct RuntimeGeneratedFunction{argnames, cache_tag, context_tag, id} <: Function
    body::Expr
    function RuntimeGeneratedFunction(cache_tag, context_tag, ex)
        def = splitdef(ex)
        args, body = normalize_args(def[:args]), def[:body]
        id = expr_to_id(body)
        cached_body = _cache_body(cache_tag, id, body)
        new{Tuple(args), cache_tag, context_tag, id}(cached_body)
    end
end

"""
    @RuntimeGeneratedFunction(function_expression)
    @RuntimeGeneratedFunction(context_module, function_expression)

Construct a function from `function_expression` which can be called immediately
without world age problems. Somewhat like using `eval(function_expression)` and
then calling the resulting function. The differences are:

* The result can be called immediately (immune to world age errors)
* The result is not a named generic function, and doesn't participate in
  generic function dispatch; it's more like a callable method.

You need to use `RuntimeGeneratedFunctions.init(your_module)` a single time at
the top level of `your_module` before any other uses of the macro.

If provided, `context_module` is module in which symbols within
`function_expression` will be looked up. By default this is module in which
`@RuntimeGeneratedFunction` is expanded.

# Examples
```
RuntimeGeneratedFunctions.init(@__MODULE__) # Required at module top-level

function foo()
    expression = :((x,y)->x+y+1) # May be generated dynamically
    f = @RuntimeGeneratedFunction(expression)
    f(1,2) # May be called immediately
end
```
"""
macro RuntimeGeneratedFunction(code)
    _RGF_constructor_code(:(@__MODULE__), esc(code))
end
macro RuntimeGeneratedFunction(context_module, code)
    _RGF_constructor_code(esc(context_module), esc(code))
end

function _RGF_constructor_code(context_module, code)
    quote
        code = $code
        cache_module = @__MODULE__
        context_module = $context_module
        if #==# !isdefined(cache_module,   $(QuoteNode(_tagname))) ||
                !isdefined(context_module, $(QuoteNode(_tagname)))
            init_mods = unique([context_module, cache_module])
            error("""You must use `RuntimeGeneratedFunctions.init(@__MODULE__)` at module
                     top level before using runtime generated functions in $init_mods""")
        end
        RuntimeGeneratedFunction(cache_module.$_tagname, context_module.$_tagname, $code)
    end
end

function Base.show(io::IO, ::MIME"text/plain", f::RuntimeGeneratedFunction{argnames, cache_tag, context_tag, id}) where {argnames,cache_tag,context_tag,id}
    cache_mod = parentmodule(cache_tag)
    context_mod = parentmodule(context_tag)
    func_expr = Expr(:->, Expr(:tuple, argnames...), f.body)
    print(io, "RuntimeGeneratedFunction(#=in $cache_mod=#, #=using $context_mod=#, ", repr(func_expr), ")")
end

(f::RuntimeGeneratedFunction)(args::Vararg{Any,N}) where N = generated_callfunc(f, args...)

# We'll generate a method of this function in every module which wants to use
# @RuntimeGeneratedFunction
function generated_callfunc end

function generated_callfunc_body(argnames, cache_tag, id, __args)
    setup = (:($(argnames[i]) = @inbounds __args[$i]) for i in 1:length(argnames))
    body = _lookup_body(cache_tag, id)
    @assert body !== nothing
    quote
        $(setup...)
        $(body)
    end
end

### Body caching and lookup
#
# Looking up the body of a RuntimeGeneratedFunction based on the id is a little
# complicated because we want the `id=>body` mapping to survive precompilation.
# This means we need to store the mapping created by a module in that module
# itself.
#
# For that, we need a way to lookup the correct module from an instance of
# RuntimeGeneratedFunction. Modules can't be type parameters, but we can use
# any type which belongs to the module as a proxy "tag" for the module.
#
# (We could even abuse `typeof(__module__.eval)` for the tag, though this is a
# little non-robust to weird special cases like Main.eval being
# Base.MainInclude.eval.)

# It appears we can't use a ReentrantLock here, as contention seems to lead to
# deadlock. Perhaps because it triggers a task switch while compiling the
# @generated function.
_cache_lock = Threads.SpinLock()
_cachename = Symbol("#_RuntimeGeneratedFunctions_cache")
_tagname = Symbol("#_RGF_ModTag")

function _cache_body(cache_tag, id, body)
    lock(_cache_lock) do
        cache = getfield(parentmodule(cache_tag), _cachename)
        # Caching is tricky when `id` is the same for different AST instances:
        #
        # Tricky case #1: If a function body with the same `id` was cached
        # previously, we need to use that older instance of the body AST as the
        # canonical one rather than `body`. This ensures the lifetime of the
        # body in the cache will always cover the lifetime of the parent
        # `RuntimeGeneratedFunction`s when they share the same `id`.
        #
        # Tricky case #2: Unless we hold a separate reference to
        # `cache[id].value`, the GC can collect it (causing it to become
        # `nothing`). So root it in a local variable first.
        #
        cached_body = haskey(cache, id) ? cache[id].value : nothing
        cached_body = cached_body !== nothing ? cached_body : body
        # Use a WeakRef to allow `body` to be garbage collected. (After GC, the
        # cache will still contain an empty entry with key `id`.)
        cache[id] = WeakRef(cached_body)
        return cached_body
    end
end

function _lookup_body(cache_tag, id)
    lock(_cache_lock) do
        cache = getfield(parentmodule(cache_tag), _cachename)
        cache[id].value
    end
end

"""
    RuntimeGeneratedFunctions.init(mod)

Use this at top level to set up your module `mod` before using
`@RuntimeGeneratedFunction`.
"""
function init(mod)
    lock(_cache_lock) do
        if !isdefined(mod, _cachename)
            mod.eval(quote
                const $_cachename = Dict()
                struct $_tagname
                end

                # We create method of `generated_callfunc` in the user's module
                # so that any global symbols within the body will be looked up
                # in the user's module scope.
                #
                # This is straightforward but clunky.  A neater solution should
                # be to explicitly expand in the user's module and return a
                # CodeInfo from `generated_callfunc`, but it seems we'd need
                # `jl_expand_and_resolve` which doesn't exist until Julia 1.3
                # or so. See:
                #   https://github.com/JuliaLang/julia/pull/32902
                #   https://github.com/NHDaly/StagedFunctions.jl/blob/master/src/StagedFunctions.jl#L30
                @inline @generated function $RuntimeGeneratedFunctions.generated_callfunc(
                        f::$RuntimeGeneratedFunctions.RuntimeGeneratedFunction{argnames, cache_tag, $_tagname, id}, __args...) where {argnames, cache_tag, id}
                    $RuntimeGeneratedFunctions.generated_callfunc_body(argnames, cache_tag, id, __args)
                end
            end)
        end
    end
end

###
### Utilities
###
normalize_args(args::Vector) = map(normalize_args, args)
normalize_args(arg::Symbol) = arg
function normalize_args(arg::Expr)
    arg.head === :(::) || error("argument malformed. Got $arg")
    arg.args[1]
end

function expr_to_id(ex)
    io = IOBuffer()
    Serialization.serialize(io, ex)
    return Tuple(reinterpret(UInt32, sha1(take!(io))))
end

end


text:
```
what makes an idea Value-Able?
<needs time to ponder upon>

it's not valuable on its own (per se)

but it's the value, in the Relation it has
to other ideas

 - an idea , stems from ... Nothing (or is it a `Complex Plane` ?....)
Q. is an idea fixed (or changing?)

there must be a meta-idea : explaining an idea

an idea is indifferent on its own; it makes no Difference ....
so...*Ponders....*

what makes an idea, an idea?
```

struct idea <: Nothing # is it nothing, or something.. of value??
    ```
    is it nothing, in the real plane, or in the complex plane
    (before that) where does the real plane start, & the complex plane meet with it?

        that is the question

        an idea is a thing you can't grasp...

    ```

    #space , is it/

    ϕ::nothing
    ϵ=10^-6 # is it a smallest thing - of value

end

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
using ColorTypes

# creating HdrImage struct
"""
HdrImage(w, h; array=[RGB() for i in 1:h*w])

It creates a **HdrImage**.

## Arguments
- *width* -> integer;
- *height* -> integer;
- *pixels* -> a monodimensional array of RGB colors;

If not specified, the pixels array is this: ``[RGB() for i in 1:h*w]``.
"""
mutable struct HdrImage
    width::Int
    height::Int
    pixels::Array{RGB, 1}

    HdrImage(w, h, array=[RGB() for i in 1:h*w]) = new(w, h, array )
end

# Check if the coordinates passed are valid
valid_coordinates(img::HdrImage, x, y) = return((x > 0) && (x <= img.width) && (y >= 0) && (y <= img.height))

# Check that the (x,y) pixel of the image is in the right place in the linear array pixels of the struct
pixel_offset(img::HdrImage, x, y) = (y-1) * img.width + x

# Get and set methods

"""
    get_pixel(img, x, y)

It returns the positional index of the (x, y) pixel of the image (img).
If the coordinates are not valid it returns a Boolean.
"""
get_pixel(img::HdrImage, x, y) =  valid_coordinates(img, x, y) && return( pixel_offset(img, x, y))

set_pixel(img::HdrImage, x, y, new_color::RGB) =  valid_coordinates(img, x, y) && (img.pixels[ get_pixel(img, x, y)] = new_color)

## SAVING an HdrImage on a stream or an output file in PFM format ################################################################

function Base.write(io::IO, img::HdrImage)
    header = transcode(UInt8, "PF\n$(img.width) $(img.height)\n$(-1.0)\n")
    write(io, header)

    for y in img.height:-1:1
        for x in 1:img.width

            color = img.pixels[ get_pixel(img, x, y)]
            write(io, convert(Vector{Float32}, [color.r,color.g,color.b] ) )

        end
    end
end

function Base.write(file_output::String, img::HdrImage)
    io = open(file_output, "w")
    write(io, img)
end

## READ an HdrImage in PMF format from a file #############################################################################

# new error message InvalidPfmFileFormat
struct InvalidPfmFileFormat <: Exception
    msg::String

    InvalidPfmFileFormat(msg::String) = new(msg)
end

# Supporting function for read_pfm_image
function _read_line(stream::IO)
    result = ""
    while true
        cur_byte = read(stream, 1)
        if cur_byte == transcode(UInt8, "") || cur_byte == transcode(UInt8, "\n")
            return result
        end
        result *= transcode(String, cur_byte)
    end
end

function _read_float(io::IO, endianness::String)
    if endianness == "LE"
        try
            return ltoh(reinterpret(Float32, read(io, 4))[1])
        catch e
            throw(InvalidPfmFileFormat("impossible to read binary data from the file"))
        end
    else
        a = read(io,4)
        try
            return ntoh(reinterpret(Float32, a)[1])
        catch e
            throw(InvalidPfmFileFormat("impossible to read binary data from the file"))
        end
    end
end

function _parse_img_size(line::String)
    elements = split(line, " ")
    if length(elements) != 2
        throw(InvalidPfmFileFormat("invalid image size specification!"))
    end

    (width, height) = (parse(Int, elements[1]), parse(Int, elements[2]))

    if (width < 0) || (height < 0)
        throw(InvalidPfmFileFormat("invalid width/size!"))
    end

    return (width, height)
end

function _parse_endianness(line::String)
    value = 0
    try
        value = parse(Float32, line)
    catch e
        throw(InvalidPfmFileFormat("missing endianness specification"))
    end

    if value == 1.0
        return "BE"
    elseif value == -1.0
        return "LE"
    else
        throw(InvalidPfmFileFormat("invalid endianness specification"))
    end
end

# finally, the real READING method:
function read_pfm_image(io::IO)

    magic =  _read_line(io)
    if magic != "PF"
    throw(InvalidPfmFileFormat("invalid magic in PFM file"))
    end

    img_size =  _read_line(io)
    (width, height) =  _parse_img_size(img_size)

    endianness_line =  _read_line(io)
    endianness =  _parse_endianness(endianness_line)

    result =  HdrImage(width, height)
    for y in height:-1:1
        for x in 1:width
            (r, g, b) = [ _read_float(io, endianness) for i in 1:3]
            set_pixel(result, x, y, RGB(r, g, b))
        end
    end

    return result
end

function read_pfm_image(filein::String)
    io = open(filein, "r")
    read_pfm_image(io)
end

## SAVE an LdrImage on an output file #############################################################################à
function average_luminosity(img::HdrImage, delta)
    sum = 0.0f0
    for pixel in img.pixels
        sum += log10(delta + luminosity(pixel))
    end
    return Float32(10^(sum/length(img.pixels)))
end

function average_luminosity(img::HdrImage)
    return average_luminosity(img, Float32(10^(-10)))
end

# Normalizing the luminosity of a image
function normalize_image(img::HdrImage, a_factor, luminosity)

    for i in 1:length(img.pixels)
        img.pixels[i] = img.pixels[i]::RGB * (a_factor/luminosity)::Float64
    end

end

function normalize_image(img::HdrImage, a_factor)

    lum = average_luminosity(img)
    normalize_image(img, a_factor, lum)

end

# clamping method
function _clamp(x)
    return x/(x+1)
end

function clamp_image(img)
    for i in 1:length(img.pixels)
        (r, g, b) = (_clamp(img.pixels[i].r), _clamp(img.pixels[i].g) , _clamp(img.pixels[i].b) )
        img.pixels[i] = RGB(r, g, b)
    end
end

function tone_mapping(img, γ)
    for i in 1:length(img.pixels)
        (r, g, b) = (floor(255*(img.pixels[i].r)^(1/γ)),
                     floor(255*(img.pixels[i].g)^(1/γ)),
                     floor(255*(img.pixels[i].b)^(1/γ)))
        img.pixels[i] = RGB(r/255, g/255, b/255)
    end
end


text:
export RegexDict

#Miscellanious utility's for users

"""
    RegexDict(::Tuple{Union{Regex, String}, T}...)

Creates a dictionary that uses Regexes as keys and tests against those
when looking up keys.
Symbols can be used as lookup keys, by using their name.
If multiple Regexes match a key, the value associated with any of them may be
returned.
"""
struct RegexDict{T}
    entries::Vector{Tuple{Regex, T}}
end

function RegexDict{T}() where {T}
    RegexDict(Vector{Tuple{Regex, T}}(undef, 0))
end

function RegexDict(pairs::Tuple{Union{Regex, String}, T}...) where {T}
    entries = Vector{Tuple{Regex, T}}(undef, length(pairs))
    for i in 1:length(pairs)
        (key, val) = pairs[i]
        if key isa String
            key = Regex(key)
        end
        entries[i] = (key, val)
    end
    RegexDict(entries)
end

function Base.get(dict::RegexDict, key::Symbol, default)
    get(dict, String(key), default)
end

function Base.get(dict::RegexDict, key::String, default)
    for (entryKey, entryVal) in dict.entries
        if occursin(entryKey, key)
            return entryVal
        end
    end
    return default
end

function Base.getindex(dict::RegexDict, index::Symbol)
    getindex(dict, String(index))
end

function Base.getindex(dict::RegexDict, index::String)
    for (entryKey, entryVal) in dict.entries
        if occursin(entryKey, index)
            return entryVal
        end
    end
    throw(BoundsError(dict, index))
end


text:
"""The cost to build a single ship."""
const SHIP_COST = constants['NEW_ENTITY_ENERGY_COST']

"""The cost to build a dropoff."""
const DROPOFF_COST = constants['DROPOFF_COST']

"""The maximum amount of halite a ship can carry."""
const MAX_HALITE = constants['MAX_ENERGY']

"""
The maximum number of turns a game can last. This reflects the fact
that smaller maps play for fewer turns.
"""
const MAX_TURNS = constants['MAX_TURNS']

"""1/EXTRACT_RATIO halite (truncated) is collected from a square per turn."""
const EXTRACT_RATIO = constants['EXTRACT_RATIO']

"""1/MOVE_COST_RATIO halite (truncated) is needed to move off a cell."""
const MOVE_COST_RATIO = constants['MOVE_COST_RATIO']

"""Whether inspiration is enabled."""
const INSPIRATION_ENABLED = constants['INSPIRATION_ENABLED']

"""
A ship is inspired if at least INSPIRATION_SHIP_COUNT opponent
ships are within this Manhattan distance.
"""
const INSPIRATION_RADIUS = constants['INSPIRATION_RADIUS']

"""
A ship is inspired if at least this many opponent ships are within
INSPIRATION_RADIUS distance.
"""
const INSPIRATION_SHIP_COUNT = constants['INSPIRATION_SHIP_COUNT']

"""An inspired ship mines 1/X halite from a cell per turn instead."""
const INSPIRED_EXTRACT_RATIO = constants['INSPIRED_EXTRACT_RATIO']

"""An inspired ship that removes Y halite from a cell collects X*Y additional halite."""
const INSPIRED_BONUS_MULTIPLIER = constants['INSPIRED_BONUS_MULTIPLIER']

"""An inspired ship instead spends 1/X% halite to move."""
const INSPIRED_MOVE_COST_RATIO = constants['INSPIRED_MOVE_COST_RATIO']

max_stars_repo_name: Affos/Halite-III

max_issues_repo_name: Affos/Halite-III

max_forks_repo_name: Affos/Halite-III


text:
# Basic types for scRNA-seq data

"""
ScRNAseqData with metadata
Fields:
* `X`: gene-by-cell expression matrix
* `g`: gene list
* `s`: coordinates of embedding of cells
* `c`: class ID of cells
"""
struct ScRNAseqData
    X::AbstractMatrix
    g::Array{String,1}
    s::Array{Float64,2}
    c::Array{Integer,1}
end

# Constructs a ScRNAseqData
function ScRNAseqData(X::AbstractMatrix, g::String, s::Array{Float64,2}, c::Array{Integer,1})
    # a=ScRNAseqData(rand(3,4),["g1","g2","g3"],[4.8 3.2 1.2 2.3; 1.0 2.0 3.0 4.0],[3, 1, 2, 1])
    @assert size(X,1)==length(g)
    @assert size(X,2)==size(s,1)
    @assert size(X,2)==length(c)
    new(X, g, s, c)
end

# using TableReader

rescale(A; dims=1) = (A .- mean(A, dims=dims)) ./ max.(std(A, dims=dims), eps())

# vecnorm(x) = x./norm.(x[:,i] for i in 1:size(x,2))'
vecnorm(x::AbstractMatrix) = norm.(x[:,i] for i in 1:size(x,2))
function normc!(x)
    for i in 1:size(x,2)
        x[:,i]=x[:,i]./norm(x[:,i])
    end
end

#=
function [A,G,H] = pythagoreanMeans(list)
    A = mean(list);           % arithmetic mean
    G = exp(mean(log(list))); % geometric mean
    H = 1./mean(1./list);     % harmonic mean
end
=#


text:
# interface functions for mapping between xyz and CAD parametric coordinates

"""
  This function gets the vector of CAD parametric coordinates from the
  mesh database.  See [`coords_xyzToXi`](@ref) for a description of the
  parametric coordinate system.

  **Inputs**

   * mesh

  **Inputs/Outputs**

   * xivec: vector to be overwritten with parametric coordinates.  Length
            `mesh.geoNums.numXiDofs`.
"""
function getXiCoords(mesh::PumiMeshDG, xivec::AbstractVector)
# retrieve Xi coordinates

  g = apf.getModel(mesh.m_ptr)
  if !gmi.can_eval(g)
    error("geometric model does not support evaluating parametric coordinates")
  end

  xi_j = zeros(Float64, 3)
  xyz_j = zeros(Float64, 3)

  for (e, dim) in apf.FieldEntityIt(mesh.m_ptr, mesh.coordshape_ptr)
    me = apf.toModel(mesh.m_ptr, e)
    me_dim = apf.getModelType(mesh.m_ptr, me)

    for j=1:mesh.coord_numNodesPerType[dim+1]
      firstidx = apf.getNumberJ(mesh.geoNums.xiNums, e, j-1, 0)

      # put the new coordinates in xyz_j
      if me_dim == mesh.dim  # xyz coordinates == xi coordinates
        getCoords(mesh, e, j-1, xi_j)
        ndof = mesh.dim
      elseif firstidx == mesh.geoNums.numXiDofs + 1  # no geometric dofs
        fill!(xi_j, 0)
        ndof = 0
      else  # constrained entity, convert parametric to xyz
        ndof = me_dim
        getCoordsXi(mesh, e, j-1, xi_j)
      end

      # write into output vector
      for k=1:ndof
        idx = apf.getNumberJ(mesh.geoNums.xiNums, e, j-1, k-1)
        xivec[idx] = xi_j[k]
      end
    end  # end j
  end  # end iterator

  return nothing
end

"""
  Similar to the other method, but allocates a new vector and returns it.

  **Inputs**

   * mesh

  **Outputs**

   * vector containing xi coordinates.  The vector has element type Float64
"""
function getXiCoords(mesh::PumiMeshDG)

  xivec = zeros(Float64, mesh.geoNums.numXiDofs)
  getXiCoords(mesh, xivec)

  return xivec
end

"""
  Gets the vector of xyz coordinates for the entire mesh.

  **Inputs**

   * mesh

  **Inputs/Outputs**

   * xvec: vector to be overwritten, length `mesh.dim*mesh.coord_numNodes`.
"""
function getXCoords(mesh::PumiMeshDG, xvec::AbstractVector)

  xyz_j = zeros(Float64, 3)
  for (e, dim) in apf.FieldEntityIt(mesh.m_ptr, mesh.coordshape_ptr)
    for j=1:mesh.coord_numNodesPerType[dim+1]

      getCoords(mesh, e, j-1, xyz_j)
      for k=1:mesh.dim
        idx = apf.getNumberJ(mesh.geoNums.coordNums, e, j-1, k-1)
        xvec[idx] = xyz_j[k]
      end
    end  # end j
  end # end iterator

  return nothing
end

"""
  Similar to the other method, but allocates a new vector and returns it.

  **Inputs**

   * mesh

  **Outputs**

   * vector containing the xyz coordinates.  Element type Float64.
"""
function getXCoords(mesh::PumiMeshDG)

  xvec = zeros(Float64, mesh.dim*mesh.coord_numNodes)
  getXCoords(mesh, xvec)

  return xvec
end

"""
  This function takes a vector of the mesh xyz coordinate field and turns it
  into a parametric representation using the CAD system.  The transformation
  has the following rules:

   1. MeshEntities that do not have any degrees of freedom in the CAD
      parametric system do not appear in the output (verticies classified on
      geometric vertices)
   2. MeshEntities that are on the interior of the domain retain their xyz
      coordinates (not all CAD system have a parametric representation of
      3D regions)
   3. All other entities (those that lie on a curve/surface that bounds a
      higher-dimension geometric entity) have their xyz coordinates transformed
      to the CAD parametric coordinates (which are fewer number than the xyz
      coordinates).

  It is assumed that the input coordinates have not been modified in such a way
  that MeshEntities no longer lie on the geometric entity they are classified
  on.

  This function cannot be complex-stepped.  The complex part of the input
  vector will be ignored, and the output vector will have zero complex part.

  This function has limited accuracy for entities classified on boundaries.
  TODO: what function to use instead.

  **Inputs**

   * mesh: DG mesh
   * xvec: a vector of mesh coordinates, perhaps obtained by
           [`coords3DTo1D`](@ref), length `mesh.dim*mesh.coord_numNodes`

  **Inputs/Outputs*

   * xivec: vector to be overwritten with parametric coordinates, length
            `mesh.geoNums.numXiDofs`
"""
function coords_xyzToXi(mesh::PumiMeshDG, xvec::AbstractVector,
                        xivec::AbstractVector)

  g = apf.getModel(mesh.m_ptr)
  if !gmi.can_eval(g)
    error("geometric model does not support evaluating parametric coordinates")
  end

  geonums = mesh.geoNums

  fill!(xivec, 0)
  xyz_j = zeros(Float64, 3)  # Pumi uses Float64 only
  xi_j = zeros(Float64, 3)

  newxyz_j = zeros(Float64, 3)

  for (e, dim) in apf.FieldEntityIt(mesh.m_ptr, mesh.coordshape_ptr)
    me = apf.toModel(mesh.m_ptr, e)
    me_dim = apf.getModelType(mesh.m_ptr, me)

    for j=1:mesh.coord_numNodesPerType[dim+1]
      firstidx = apf.getNumberJ(geonums.xiNums, e, j-1, 0)

      # get xyz coordinates
      for k=1:mesh.dim
        idx = apf.getNumberJ(geonums.coordNums, e, j-1, k-1)
        xyz_j[k] = real(xvec[idx])
      end

      # decide parametric representation
      if me_dim == mesh.dim  # keep xyz coordinates
        for k=1:mesh.dim
          xi_j[k] = xyz_j[k]
        end

      elseif firstidx == geonums.numXiDofs + 1  # no geometric dofs
        for k=1:mesh.dim
          xi_j[k] = 0
        end
      else  # this is a constrained entity, find parametric representation
        getSnappedCoords(mesh, e, true, xyz_j, newxyz_j, xi_j)
#        gmi.closest_point(g, me, xyz_j, newxyz_j, sview(xi_j, 1:2))
#          @assert norm(newxyz_j - xyz_j) < 1e-8  # xyz_j should lie on this
#                                                 # geometric entity
      end

      # put in output vector
      for k=1:me_dim
        idx = apf.getNumberJ(geonums.xiNums, e, j-1, k-1)
        xivec[idx] = xi_j[k]
      end

    end  # end j
  end  # end iterator

  return nothing
end

"""
  Methods that takes a `mesh.vert_coords` shaped array as input, rather
  than a vector.

  **Inputs**

   * mesh
   * xarr: 3D array of xyz coordinates

  **Inputs/Outputs**

   * xivec: vector, length `mesh.geoNums.numXiDofs`
"""
function coords_xyzToXi(mesh::PumiMeshDG, xarr::Abstract3DArray{T},
                        xivec::AbstractVector) where {T}

  xvec = zeros(T, mesh.geoNums.numCoordDofs)
  op = AssignReduction{T}()
  coords3DTo1D(mesh, xarr, xvec, op, parallel=false)

  coords_xyzToXi(mesh, xvec, xivec)
end

"""
  For non-periodic geometric entities, this function asserts the xi coordinates
  are in the allowed range.  If periodic geometric entities, this function wraps
  the coordinates to ensure they are in the correct range.  This function
  correctly preserves complex values if `xivec` is complex.

  This function is called by the other geometric interface functions when
  needed, users generally do not need to call it.

  **Inputs**

   * mesh

  **Inputs/Outputs**

   * xivec: vector of xi coordinates, length `mesh.geoNums.numXiDofs`.
"""
function wrapXiCoords(mesh::PumiMeshDG, xivec::AbstractVector)

  g = apf.getModel(mesh.m_ptr)
  if !gmi.can_eval(g)
    error("geometric model does not support evaluating parametric coordinates")
  end

  xi_idx = zeros(Int, 2)
  rng = zeros(Float64, 2)
  geonums = mesh.geoNums

  for (e, dim) in apf.FieldEntityIt(mesh.m_ptr, mesh.coordshape_ptr)
    me = apf.toModel(mesh.m_ptr, e)
    me_dim = apf.getModelType(mesh.m_ptr, me)

    for j=1:mesh.coord_numNodesPerType[dim+1]

      for k=1:2
        xi_idx[k] = apf.getNumberJ(geonums.xiNums, e, j-1, k-1)
      end

      # if not in geometric interior and has geometric degrees of freedom
      if !(me_dim == mesh.dim) && !(xi_idx[1] == geonums.numXiDofs + 1)
        for k=1:me_dim
          idx = xi_idx[k]
          gmi.range(g, me, k-1, rng)

          if gmi.periodic(g, me, k-1)
            xivec[idx] = wrapPeriodicXi(rng, xivec[idx])
          else
            @assert real(xivec[idx]) >= rng[1]
            @assert real(xivec[idx]) <= rng[2]
          end
        end  # end for
      end  # end if
    end  # end j
  end  # end iterator

  return nothing
end

"""
  Wraps a periodic xi value such that it represents the same location, but has
  value within the given range.  Methods are available for real and complex
  numbers.

  **Inputs**

   * rng: vector of length 2, containing the range of parametric values
   * xival: the xi value

  **Outputs**

   * xival: a new xi value
"""
function wrapPeriodicXi(rng::AbstractVector, xival::Real)

  # the algorithm is to shift to the range (0, rng[2] - rng[1]),
  # compute the modulus with the length of the parameter value range, then
  # shift back to the original range.
  delta = rng[2] - rng[1]
  t1 = (xival - rng[1]) % delta

  t1 += rng[1]

  if t1 < rng[1]
    t1 += delta
  end

  return t1
end

function wrapPeriodicXi(rng::AbstractVector, xival::T) where {T <: Complex}

  xival_imag = imag(xival)
  xival_real = wrapPeriodicXi(rng, real(xival))

  # because the geometric entity is periodic, wrapping by a multiple of the
  # period shouldn't change the derivative
  return T(xival_real, xival_imag)
end

"""
  This function is the inverse of [`coords_XiToXYZ`](@ref).  It recomputes
  the xyz coordinates from the parametric coordinates.  For MeshEntities which
  do not have any geometric degrees of freedom (and therefore no entries in
  `xivec`, the coordinates are retrieved from Pumi.

  **Inputs**

   * mesh
   * xivec: vector of CAD parametric coordinates, length
            `mesh.geoNums.numXiDofs`.  This vector may be modified if any
            periodic parametric values are out of range.  see
            [`wrapXiCoords`](@ref)

  **Inputs/Outputs**

   * xvec: vector to be overwritten with xyz coordinates, length `mesh.dim*mesh.coord_numNodes`.
"""
function coords_XiToXYZ(mesh::PumiMeshDG, xivec::AbstractVector,
                        xvec::AbstractVector)
  g = apf.getModel(mesh.m_ptr)
  if !gmi.can_eval(g)
    error("geometric model does not support evaluating parametric coordinates")
  end

  # make sure all values are valid
  wrapXiCoords(mesh, xivec)

  geonums = mesh.geoNums

  xyz_j = zeros(Float64, 3)
  xi_j = zeros(Float64, 2)
  xi_idx = zeros(Cint, mesh.dim)

  for (e, dim) in apf.FieldEntityIt(mesh.m_ptr, mesh.coordshape_ptr)
    me = apf.toModel(mesh.m_ptr, e)
    me_dim = apf.getModelType(mesh.m_ptr, me)

    for j=1:mesh.coord_numNodesPerType[dim+1]

      for k=1:mesh.dim
        xi_idx[k] = apf.getNumberJ(geonums.xiNums, e, j-1, k-1)
      end

      # put the new coordinates in xyz_j
      if me_dim == mesh.dim  # xivec has xyz coordinates
        for k=1:mesh.dim
          xyz_j[k] = xivec[xi_idx[k]]
        end
      elseif xi_idx[1] == geonums.numXiDofs + 1  # no geometric dofs
        getCoords(mesh, e, j-1, xyz_j)
      else  # constrained entity, convert parametric to xyz
        for k=1:me_dim
          xi_j[k] = xivec[xi_idx[k]]
        end
        gmi.geval(g, me, xi_j, xyz_j)
      end

      # write xyz_j to xvec
      for k=1:mesh.dim
        idx = apf.getNumberJ(geonums.coordNums, e, j-1, k-1)
        xvec[idx] = xyz_j[k]
      end

    end   # end j
  end  # end iterator

  return nothing
end

"""
  Method that takes a 3D array of xyz coordinates as the output.

  Note that this method does not update the Pumi coordinate database nor
  does it recalculate the metrics.
  Users should call [`update_coords`](@ref) if this is required.

  **Inputs**

   * mesh
   * xivec: vector if CAD parametric coordinates, length `mesh.geoNums.numXiDofs`.

  **Inputs/Outputs**

   * xarr: array, same shape as `mesh.vert_coords`, to be overwritten with
           new coordinates.
"""
function coords_XiToXYZ(mesh::PumiMeshDG, xivec::AbstractVector,
                        xarr::Abstract3DArray{T}) where {T}

  xvec = zeros(T, mesh.geoNums.numCoordDofs)
  coords_XiToXYZ(mesh, xivec, xvec)

  coords1DTo3D(mesh, xvec, xarr)
end

"""
  This function takes a vector that contains derivatives wrt to the xyz
  coordinates and transforms it to derivative wrt geometric coordinates xi.
  The derivative dx/dxi is evaluated at the x coordinates stored in the Pumi
  database.

  **Inputs**

   * mesh
   * xvec: vector containing dJ/dx

  **Inputs/Outputs**

   * xivec: vector to be overwritten with dJ/dxi
"""
function coords_dXTodXi(mesh::PumiMeshDG, xvec::AbstractVector,
                        xivec::AbstractVector)

  g = apf.getModel(mesh.m_ptr)
  if !gmi.can_eval(g)
    error("geometric model does not support evaluating parametric coordinates")
  end

  geonums = mesh.geoNums

  fill!(xivec, 0)
  dx_j = zeros(Float64, 3)
  dxi_j = zeros(Float64, 3)
  xidx = zeros(Cint, 3)

#  x_j = zeros(Float64, 3)
#  xnew_j = zeros(Float64, 3)
  xi_j = zeros(Float64, 2)
  dx_dxi = zeros(Float64, 3, 2)
  dx_dxi1 = sview(dx_dxi, :, 1)
  dx_dxi2 = sview(dx_dxi, :, 2)

  for (e, dim) in apf.FieldEntityIt(mesh.m_ptr, mesh.coordshape_ptr)
    me = apf.toModel(mesh.m_ptr, e)
    me_dim = apf.getModelType(mesh.m_ptr, me)

    for j=1:mesh.coord_numNodesPerType[dim+1]
      firstidx = apf.getNumberJ(geonums.xiNums, e, j-1, 0)

       # get xyz derivative
      for k=1:mesh.dim
        idx = apf.getNumberJ(geonums.coordNums, e, j-1, k-1)
        dx_j[k] = real(xvec[idx])
      end

      # decide parametric representation
      if me_dim == mesh.dim  # keep xyz coordinates
        for k=1:mesh.dim
          dxi_j[k] = dx_j[k]
        end
      elseif firstidx == geonums.numXiDofs + 1  # no geometric dofs
        for k=1:mesh.dim
          dxi_j[k] = 0
        end
      else  # this is a constrained entity, compute dx/dxi
        getCoordsXi(mesh, e, j-1, xi_j)
        gmi.first_derivative(g, me, xi_j, dx_dxi1, dx_dxi2)

        # apply chain rule to compute dJ/dxi
        fill!(dxi_j, 0)
        for k=1:me_dim
          for d=1:mesh.dim
            dxi_j[k] += dx_j[d]*dx_dxi[d, k]
          end
        end

      end  # end if

      # put in output vector
      for k=1:me_dim
        idx = apf.getNumberJ(geonums.xiNums, e, j-1, k-1)
        xivec[idx] = dxi_j[k]
      end

    end  # end j
  end  # end iterator

  return nothing
end

"""
  Takes a vector of the xyz coordinates for the entire mesh and updates the
  mesh with them.  Prefer the per-element method of this function when possible,
  this method is only preferable when the vector of all coordinates is required
  for some other reason.

  Unlike the per-element method, this function calls [`commit_coords`](@ref).

  Prefer the method that takes both `xvec` and `xivec` as arguments.  This
  method may lead to somewhat less accurate derivatives computed by
  [`coords_dXTodXi`](@ref).  This method is useful only when the caller
  has no interest in the parametric coordinates.

  **Inputs**

   * mesh
   * sbp
   * opts

  **Inputs/Outputs**

   * xvec: vector of new coordinates, length `mesh.dim*mesh.coord_numNodes`

  **Keyword Arguments**

   * snap: if true, the coordinates in `xvec` will be snapped to the
           nearest point on the geometry.  The updates coordinates will
           be saved to the mesh database and written back into `xvec`.
           If false, the coordinates in `xvec` will be entered into the
           mesh database without snapping.  This can lead to inconsistent
           parametric and Cartesian coordinates, as well as inaccurate
           derivatives.  Use only if the caller absolutely requires the
           coordinates be entered exactly as they appear in `xvec`
   * verify: run Pumi's verifier on the new mesh, default true
   * write_vis: write visulzation files after updating the coordinates but
                before recalculating the metrics (which errors if an element
                is turned inside out), default false.
"""
function update_coords(mesh::PumiMesh, sbp, opts, xvec::AbstractVector;
                       snap::Bool=true, verify=true, write_vis=false)

  @assert length(xvec) == mesh.dim*mesh.coord_numNodes

  xyz_j = zeros(Float64, 3)

  g = apf.getModel(mesh.m_ptr)
  _snap::Bool = snap && gmi.can_eval(g)   # make type-stable
  can_eval = gmi.can_eval(g)

  for (e, dim) in apf.FieldEntityIt(mesh.m_ptr, mesh.coordshape_ptr)
    me = apf.toModel(mesh.m_ptr, e)
    for j=1:mesh.coord_numNodesPerType[dim+1]

      for k=1:mesh.dim
        idx = apf.getNumberJ(mesh.coord_nodenums_Nptr, e, j-1, k-1)
        xyz_j[k] = xvec[idx]
      end

      setCoords(mesh, e, j-1, xyz_j, _snap)
      # write back into xvec
      for k=1:mesh.dim
        idx = apf.getNumberJ(mesh.coord_nodenums_Nptr, e, j-1, k-1)
        xvec[idx] = xyz_j[k]
      end
    end  # end j
  end  # end iterator

  commit_coords(mesh, sbp, opts, verify=verify, write_vis=write_vis)

  return nothing
end

#TODO: single elemet version of update_coords

# It looks like the MeshAdapt Snapper updates the parametric coords for
# vertices, but not mid edge nodes (I also recall that MeshAdapt doesn't
# support quadratic triangles either).

"""
  Convenience method that computes the xyz coordinates from the parametric
  coordinates and updates the mesh database.  This method is more accurate than
  the function that updates the mesh database using the xyz coordinates

  **Inputs**

   * mesh
   * sbp
   * opts
   * xivec: vector of xi coordinates, length `mesh.geoNums.numXiDofs`

  **Keyword Arguments**

   * verify: run Pumi's verifier on the new mesh, default true
   * write_vis: write visulzation files after updating the coordinates but
                before recalculating the metrics (which errors if an element
                is turned inside out), default false.
"""
function update_coordsXi(mesh::PumiMesh, sbp, opts, xivec::AbstractVector;
                       verify=true, write_vis=false)
  @assert length(xivec) == mesh.geoNums.numXiDofs

  # wrap periodic coordinates
  wrapXiCoords(mesh, xivec)

  xi_j = zeros(Float64, 3)
  geonums = mesh.geoNums
  for (e, dim) in apf.FieldEntityIt(mesh.m_ptr, mesh.coordshape_ptr)

    me = apf.toModel(mesh.m_ptr, e)
    me_dim = apf.getModelType(mesh.m_ptr, me)

    for j=1:mesh.coord_numNodesPerType[dim+1]

      # decide parametric representation
      firstidx = apf.getNumberJ(geonums.xiNums, e, j-1, 0)
      if me_dim == mesh.dim
        # parametric coords == xyz coordinates
        fill!(xi_j, 0)
        for k=1:me_dim
          xi_j[k] = xivec[apf.getNumberJ(geonums.xiNums, e, j-1, k-1)]
        end
        setCoords(mesh, e, j-1, xi_j)
      elseif firstidx == geonums.numXiDofs + 1  # no geometric dofs
        # this must be classified on a vertex, which always has parametric
        # coordinates 0
        fill!(xi_j, 0)
      else  # this is a constrained entity, get the values from xivec
        fill!(xi_j, 0)
        for k=1:me_dim
          xi_j[k] = xivec[apf.getNumberJ(geonums.xiNums, e, j-1, k-1)]
        end
        setCoordsXi(mesh, e, j-1, xi_j)
      end  # end if
    end  # end j
  end  # end iterator

  commit_coords(mesh, sbp, opts, verify=verify, write_vis=write_vis)

  return nothing
end


text:
__precompile__(false)

module SCRIPT

using DataStructures #OrderedDict

using GLFW, ModernGL, SharedArrays, StaticArrays, Quaternions

using RessourceManager
#using WindowManager
using GraphicsManager
using DefaultModelData
using CameraManager
using FrustumManager
using ChunkManager
using MeshManager
using TextureManager
using ShaderManager
#using ..ScriptManager
using TimeManager
using LogManager
using MathManager

const GPU = GraphicsManager

mutable struct Script
  args::Dict{Symbol,Any}
  vars::Dict{Symbol,Any}
  shaderManager::JShaderManager
  Script() = new(Dict(),Dict(),JShaderManager())
end

this = Script()

function reset()
  this.vars = typeof(this.vars)()
  delPrograms(this.shaderManager)
end

""" TODO """
function main(args::Dict{Symbol,Any})
  this.args = args
  println("Script: $(basename(@__FILE__)), $(this.args), time: $(mtime(@__FILE__))")
  reset()
end

""" TODO """
function OnDestroy()
  println("OnDestroy")
  reset()
  GraphicsManager.cleanUp() #remove all
  GC.gc()
end

function resizeWindow(window, width, height)
  this.vars[:RESOLUTION] = (
    width = width,
    height = height,
    ratio = width / (height*1f0),
    size = width * height
  )

  GLFW.SetWindowSize(window, width, height)
  glViewport(0, 0, width, height)
end

""" TODO """
function OnInit()
  SS = this.vars[:SECTION_SHADER] = Dict{Symbol,Any}()
  SM = this.vars[:SECTION_MESH] = Dict{Symbol,Any}()
  ST = this.vars[:SECTION_TEXTURES] = Dict{Symbol,Any}()

  UF = SS[:UNIFORMS] = OrderedDict{Symbol,Any}()
  #UF = SS[:UNIFORMS] = OrderedDict{Symbol,Any}()
  #TS = ST[:TEXTURES] = OrderedDict{Symbol,Any}()

  resizeWindow(this.args[:WINDOW], 800, 600)

  glEnable(GL_DEPTH_TEST)
  glEnable(GL_BLEND)
  glEnable(GL_CULL_FACE)
  #glEnable(GL_MULTISAMPLE)

  glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)
  #glBlendEquation(GL_FUNC_ADD)
  #glFrontFace(GL_CCW)

  glCullFace(GL_BACK) #default: GL_BACK
  glDepthMask(GL_TRUE) #default: GL_TRUE
  glDepthFunc(GL_LEQUAL) #default: GL_LESS

  #glPolygonMode(GL_FRONT_AND_BACK, GL_LINE)#GL_FILL,GL_LINE
  glClearColor(0.5, 0.5, 0.5, 0.0)

  resolution = values(this.vars[:RESOLUTION])[1:2]

  UF[:iResolution] = Uniform(value=Float32[resolution...,1])
  UF[:iTime] = Uniform(value=Float32(0),default=0)
  UF[:iFrame] = Uniform(value=Float32(0),default=0)
  UF[:iTimeDelta] = Uniform(value=Float32(0),default=0)
  UF[:iFrameRate] = Uniform(value=Float32(0),default=0)
  UF[:iChannelTime] = Uniform(value=zeros(Float32,4))
  UF[:iMouse] = Uniform(value=zeros(Float32,4))
  UF[:iDate] = Uniform(value=zeros(Float32,4))
  UF[:iSampleRate] = Uniform(value=Float32(0))
  UF[:iChannelResolution] = Uniform(value=[zeros(Float32,3),zeros(Float32,3),zeros(Float32,3),zeros(Float32,3)])
  UF[:iChannel0] = Uniform(value=Int32(0),typ=GL_SAMPLER_2D,layout=Dict(:binding=>0))
  UF[:iChannel1] = Uniform(value=Int32(1),typ=GL_SAMPLER_2D,layout=Dict(:binding=>1))
  UF[:iChannel2] = Uniform(value=Int32(2),typ=GL_SAMPLER_2D,layout=Dict(:binding=>2))
  UF[:iChannel3] = Uniform(value=Int32(3),typ=GL_SAMPLER_2D,layout=Dict(:binding=>3))

  this.vars[:renderTime] = this.vars[:startTime] = time()

  SM[:SCREEN] = MeshData()
  linkData(SM[:SCREEN], :vertices=>(DATA_PLANE2D_VERTEX_STRIP,2))

  size = (UInt32(resolution[1]), UInt32(resolution[2]))

  txp(path::String) = RessourceManager.getPath(:TEXTURES, path)

  global iChannel0=createTexture(:iChannel0, size) #txp("A.png"))
  global iChannel1=createTexture(:iChannel1, size) #txp("B.png"))
  global iChannel2=createTexture(:iChannel2, size) #txp("C.png"))
  global iChannel3=createTexture(:iChannel3, size) #txp("D.png"))

  glBindImageTexture(0, iChannel0, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RGBA32F)
  glBindImageTexture(1, iChannel1, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RGBA32F)
  glBindImageTexture(2, iChannel2, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RGBA32F)
  glBindImageTexture(3, iChannel3, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RGBA32F)

  SS[:PROGRAMS] = OrderedDict(
  :SCREEN => (shaders=[:SCREEN_VSH, :MC_FSH],attributes=[:SCREEN],),
  :MC_BA => (shaders=[:MC_BA_CSH],),
  :MC_BB => (shaders=[:MC_BB_CSH],),
  :MC_BC => (shaders=[:MC_BC_CSH],),
  :MC_BD => (shaders=[:MC_BD_CSH],),
  )

  dispatch = (x=round(UInt32, resolution[1]/16f0), y=round(UInt32, resolution[2]/16f0), z=UInt32(1))
  dispatchInFile = (x=16,y=16,z=1)

  SS[:DISPATCH] = [
  (:MC_BA, dispatch),
  (:MC_BB, dispatch),
  (:MC_BC, dispatch),
  (:MC_BD, dispatch),
  ]

  SS[:GLOBALS] = Dict(
  :iDispatchX=>dispatchInFile.x,
  :iDispatchY=>dispatchInFile.y,
  :iDispatchZ=>dispatchInFile.z,
  )

  ST[:TEXTURES] = OrderedDict(
  :iChannel0 => (location=0,),
  :iChannel1 => (location=1,),
  :iChannel2 => (location=2,),
  :iChannel3 => (location=3,),
  )

  reloadShaderPrograms()
end

""" TODO """
function OnReload()
end

""" TODO """
function OnUpdate()
  checkForUpdate()

  UF = this.vars[:SECTION_SHADER][:UNIFORMS]

  deltaTime = Float32(this.vars[:renderTime] - this.vars[:startTime]) #now timeIntervalSinceDate:time()
  if deltaTime > 1f0/20f0
      deltaTime = 1f0/20f0
  end

  UF[:iTimeDelta].value = Float32(deltaTime)
  UF[:iFrameRate].value = Float32(1f0/deltaTime)
end

function setShaderProperties()
  for (k,v) in this.vars[:SECTION_SHADER][:UNIFORMS]
    setShaderProperty(string(k), v.value)
  end
end

function bindTextures()
  for (k,v) in this.vars[:SECTION_TEXTURES][:TEXTURES]
    bindTexture(k, v.location)
  end
end

function dispatchPrograms(list)
  for (p,d) in list
    useShaderProgram(p)
    setShaderProperties()
    glDispatchCompute(d.x, d.y, d.z)
    glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT)
  end
end

""" TODO """
function OnRender()
  SS = this.vars[:SECTION_SHADER]
  SM = this.vars[:SECTION_MESH]
  UF = SS[:UNIFORMS]

  iframe = UF[:iFrame]
  UF[:iTime].value = Float32(programTime())

  bindTextures()

  dispatchPrograms(SS[:DISPATCH])

  mesh = SM[:SCREEN]

  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT | GL_STENCIL_BUFFER_BIT)
  useShaderProgram(:SCREEN)
  glBindVertexArray(mesh.vao)
  setShaderProperties()
  glDrawArrays(GL_TRIANGLE_STRIP, 0, mesh.draw.count)

  UF[:iFrame].value += 1
  this.vars[:renderTime] = time()
end

""" TODO """
function checkForUpdate()
  keyValue, keyPressed = getKey()

  if keyPressed
    if keyValue == 45 #ß
    elseif keyValue == 48 #0
    elseif keyValue >= 49 && keyValue <= 57 #1-9

    elseif keyValue == 65 #a
    elseif keyValue == 66 #b
    elseif keyValue == 67 #c
    elseif keyValue == 68 #d
    elseif keyValue == 69 #e
    elseif keyValue == 70 #f
    elseif keyValue == 71 #g
    elseif keyValue == 72 #h
    elseif keyValue == 73 #i
    elseif keyValue == 74 #j
    elseif keyValue == 75 #k
    elseif keyValue == 76 #l
    elseif keyValue == 77 #m
    elseif keyValue == 78 #n
    elseif keyValue == 79 #o
    elseif keyValue == 80 #p
    elseif keyValue == 81 #q
    elseif keyValue == 82 #r
      reloadShaderPrograms()
    elseif keyValue == 83 #s
    elseif keyValue == 84 #t
    elseif keyValue == 85 #u
    elseif keyValue == 86 #v
    elseif keyValue == 86 #w
    elseif keyValue == 87 #x
    elseif keyValue == 88 #y
    elseif keyValue == 89 #z

    elseif (keyValue >= 290 && keyValue <= 301) # F1 - F12
    end
  end
end

useShaderProgram(program::Union{Nothing, Symbol}) = ShaderManager.useProgram(this.shaderManager, program)
getShaderProgram(program::Symbol) = ShaderManager.getProgram(this.shaderManager, program).id
setShaderProperty(name::String, value::Any; program::Union{Nothing,Symbol}=getActiveProgram(this.shaderManager)) = ShaderManager.setProperty(this.shaderManager, name::String, value::Any; program=program)
#if mode != "" println("MODE(",stringColor(mode;color=:yellow),"): ",stringColor(value;color=:yellow)) end

function reloadShaderPrograms()
  println("Load shader files...")
  SS = this.vars[:SECTION_SHADER]
  SM = this.vars[:SECTION_MESH]

  shaders = loadShaders(SS)

  for (p,t) in SS[:PROGRAMS]
    ShaderManager.reloadProgram(this.shaderManager, p, (x->shaders[x]).(t.shaders))
    programID = getShaderProgram(p)
    if haskey(t,:attributes)
      for a in t.attributes
        setAttributes(SM[a], programID)
      end
    end
  end

end

end # SCRIPT

max_stars_repo_name: Gilga/JLGLAssets

max_issues_repo_name: Gilga/JLGLAssets

max_forks_repo_name: Gilga/JLGLAssets


text:
isdefined(Base, :__precompile__) && __precompile__()

module MichiganAutonomousVehicles
using NLOptControl
import YAML.load

include("VehicleModels/VehicleModels.jl")
using .VehicleModels
using Parameters # needed to export @unpack and @pack

include("CaseModule.jl")
using .CaseModule

include("AutonomousControl.jl")
using .AutonomousControl

export
    # CaseModule.jl
    setConfig,
    case2dfs,

    # AutonomousControl.jl
    initializeAutonomousControl,
    updateAutoParams!,
    avMpc,
    solverConfig,
    fixYAML,
    load, # from YAML
    configProb!,
    obstacleAvoidanceConstraints!,
    lidarConstraints!,
    objFunc!,
    goalRange,
    simMPC!,

    # VehicleModels.jl
    # Objects
    ###########
    Vpara,

    # Functions
    ###########
    checkCrash,
    FNames,

    # Three DOF
    ThreeDOFv1,
    ThreeDOFv2,
    ThreeDOFv2_expr,
    ThreeDOFv3,
    ThreeDOFv3_expr,

    # KinematicBicycle
    KinematicBicycle,
    KinematicBicycle_expr,
    KinematicBicycle2,
    KinematicBicycle_expr2,

    # Macros and support functions
    ###############################
    # Three DOF
    @F_YF,
    @F_YR,
    @FZ_RL,
    @FZ_RR,
    @FZ_FL,
    @FZ_FR,
    @Ax_min,
    @Ax_max,

    # parameters
    @unpack_Vpara,
    @pack_Vpara,

    # Parameters.jl
    @unpack,
    @pack,

    # Plots
    obstaclePlot,
    trackPlot,
    mainSim,
    posterP,
    posPlot,
    vtPlot,
    pSimGR,
    axLimsPlot
end


text:
"""
   Add the virtual source to pressure (just add to pz), the source-side wavefield
of all time-steps are saved in memory.
"""
function add_virtual_source!(spt::Snapshot, m::Vector{Tv}, wfd1::Wavefield, wfd2::Wavefield,
                             params::TdParams) where {Tv<:AbstractFloat}

    # total number of element
    N = params.nz * params.nx

    for i = 1 : N
        j = params.spt2wfd[i]
        dpdt = 2.0 * (wfd2.p[i] - wfd1.p[i]) / params.vel[i]
        spt.px[j] = spt.px[j] +  dpdt * m[i]
    end

    return nothing
end

"""
   reconstruct source-side wavefield for born approximation and FWI
"""
function sourceside_reconstruct_forward(path_bnd::Tp, src::Ts,
         params::TdParams) where {Tp<:String, Ts<:Union{Source,Vector{Source}}}

    # make the source term as a vector has one element
    if Ts <: Source
       srcs    = Vector{Source}(undef,1)
       srcs[1] = src
    else
       srcs    = src
    end
    ns = length(srcs)

    # length of one-step pressure field
    N = params.nz * params.nx

    # initialize intermediate variables
    wfd1 = Wavefield(params)
    wfd2 = Wavefield(params)
    tmp_z1 = zeros(params.data_format, params.nz)
    tmp_z2 = zeros(params.data_format, params.nz)
    tmp_x1 = zeros(params.data_format, params.nx)
    tmp_x2 = zeros(params.data_format, params.nx)

    # save the sourceside wavefield at one time step
    dpdt   = zeros(params.data_format, N)

    # initialize boundary value
    fid_bnd = open(path_bnd, "r")
    bnd     = WavefieldBound(params)

    # allocate memory for saving pressure field
    pre = zeros(params.data_format, N * (params.nt-1))

    # add source to the first wavefield
    add_source!(wfd1, srcs, 1)

    # the start index for saving the current source-side wavefield
    idx_o = 1

    # loop over time stepping
    for it = 2 : params.nt

        # read the boundary value
        read_one_boundary!(bnd, fid_bnd, it, params)

        # forward time steping and correcting boundaries
        one_step_forward!(wfd2, wfd1, bnd, params,
                          tmp_z1, tmp_z2, tmp_x1, tmp_x2)

        # subtract the source term if source located in the saved boundary area
        for i = 1 : ns

            if (srcs[i].isz <= params.order || srcs[i].isz > params.nz-params.order ||
                srcs[i].isx <= params.order || srcs[i].isx > params.nx-params.order)

                subtract_source!(wfd2, srcs[i], it)
            end
        end

        # compute the source-side wave field
        for i = 1 : N
            dpdt[i] = 2.0 * (wfd2.p[i] - wfd1.p[i]) / params.vel[i]
        end

        # save the source-side field
        copyto!(pre, idx_o, dpdt, 1, N)
        idx_o = idx_o + N

        # add source to wavefield and prepare for next time step
        add_source!(wfd2, srcs, it)

        # prepare for next step
        copy_wavefield!(wfd1, wfd2)
    end

    close(fid_bnd)
    return reshape(pre, params.nz, params.nx, params.nt-1)

end

"""
   Born approximation for sourceside wavefield is generated by simultaneouse source
"""
function born_approximation_forward(irz::Ti, irx::Ti, m::Vector{Tv}, path_bnd::Tp,
         src::Ts, params::TdParams; location_flag="index") where {Ti<:Vector, Tv<:AbstractFloat, Tp<:String, Ts<:Union{Source,Vector{Source}}}

    # create recordigns
    rec = Recordings(irz, irx, params; location_flag=location_flag)

    # make the source term as a vector has one element
    if Ts <: Source
       srcs    = Vector{Source}(undef,1)
       srcs[1] = src
    else
       srcs    = src
    end
    ns = length(srcs)

    # model length
    N = params.nz * params.nx

    # allocate memory for computing adjoint wavefield
    spt1 = Snapshot(params)
    spt2 = Snapshot(params)
    tmp_z1 = zeros(params.data_format, params.Nz)
    tmp_z2 = zeros(params.data_format, params.Nz)
    tmp_x1 = zeros(params.data_format, params.Nx)
    tmp_x2 = zeros(params.data_format, params.Nx)

    # allocate memory for reconstructing source-side wavefield backward
    wfd1 = Wavefield(params); add_source!(wfd1, srcs, 1); # get the wavefield at first time step
    wfd2 = Wavefield(params)
    wfd_z1 = zeros(params.data_format, params.nz)
    wfd_z2 = zeros(params.data_format, params.nz)
    wfd_x1 = zeros(params.data_format, params.nx)
    wfd_x2 = zeros(params.data_format, params.nx)

    # initialize the boundary value as zero
    bnd = WavefieldBound(params)
    fid_bnd = open(path_bnd, "r")

    # forward time stepping
    for it = 2 : params.nt

        # time stepping of scattered wavefield
        one_step_forward!(spt2, spt1, params, tmp_z1, tmp_z2, tmp_x1, tmp_x2)

        # reconstruct source-side wavefield
        read_one_boundary!(bnd, fid_bnd, it, params)
        one_step_forward!(wfd2, wfd1, bnd, params, wfd_z1, wfd_z2, wfd_x1, wfd_x2)

        # subtract source term if it being added during update boundary part
        # loop over all the sources
        for i = 1 : ns

            if (srcs[i].isz <= params.order || srcs[i].isz > params.nz-params.order ||
                srcs[i].isx <= params.order || srcs[i].isx > params.nx-params.order)

                subtract_source!(wfd2, srcs[i], it)
            end
        end

        # add the virtual source
        add_virtual_source!(spt2, m, wfd1, wfd2, params)

        # sampling scatter wavefield
        sample_spt2rec!(rec, spt2, it)

        # prepare for the next step backward reconstruction
        add_source!(wfd2, srcs, it)
        copy_wavefield!(wfd1, wfd2)
        copy_snapshot!(spt1, spt2)

    end

    # close the boundary value file
    close(fid_bnd)

    return rec
end

"""
   parallel forward born approximation for multiple shots.
"""
function born_approximation_forward(dir_born::Ts, path_m::Ts, irz::Ti, irx::Ti, dir_sourceside::Ts,
         fidiff::TdParams; location_flag="index", normalization_flag=true, mute_index::Int64=0) where {Ts<:String, Ti<:Vector, T<:Union{Source, Vector{Source}}}

    # wrapper of forward born approximation
    function wrap_born_forward(params::NamedTuple)

        # read source wavelet
        src = read_source(params.path_src)

        rec = born_approximation_forward(params.irz, params.irx, params.m, params.path_bnd, src, params.fidiff; location_flag=params.location_flag)
        write_recordings(params.path_born, rec)

        return nothing
    end

    # create folder to save the result
    rm(dir_born, force=true, recursive=true)
    mkdir(dir_born)
    if !isdir(dir_born) # check the directory is created
       error("can't create directory for forward born approximation")
    end

    dir_bnd = joinpath(dir_sourceside, "boundary")
    dir_src = joinpath(dir_sourceside, "source")
    isdir(dir_bnd) || error("boundary file doesn't exist")
    isdir(dir_src) || error("source wavelet file doesn't exist")

    # determine the number of shot
    file_src = readdir(dir_src)
    ns       = length(file_src)
    if file_src[1] == ".DS_Store" # in case of MacOS
       ns = ns - 1
    end

    # apply preconditioner to model parameter
    (hdr, m) = read_RSdata(path_m)
    if normalization_flag
       path_normalization = joinpath(dir_sourceside, "normalization.rsf")
       (hdr, scale) = read_RSdata(path_normalization)
       m .= vec(scale) .* m
    end

    if mute_index > 0
       m  = reshape(m, fidiff.nz, fidiff.nx)
       m[1:mute_index,:] .= 0.0
       m  = vec(m)
    end

    # prepare argument
    argument_collection = Vector{NamedTuple}(undef, ns)
    for i = 1 : ns

        file_name = join(["recordings_" "$i" ".bin"])
        path_born = joinpath(dir_born, file_name)

        file_name = join(["boundary_" "$i" ".rsf"])
        path_bnd  = joinpath(dir_bnd, file_name)

        file_name = join(["source_" "$i" ".bin"])
        path_src  = joinpath(dir_src, file_name)

        # OBN acquisition geometry
        if eltype(irz) <: Real
           argument_collection[i] = (path_born=path_born, irz=irz, irx=irx, location_flag=location_flag, m=m,
                                     path_bnd=path_bnd, path_src=path_src, fidiff=fidiff)

        # towed streamer
        elseif eltype(irz) <: Vector
           argument_collection[i] = (path_born=path_born, irz=irz[i], irx=irx[i], location_flag=location_flag, m=m,
                                     path_bnd=path_bnd, path_src=path_src, fidiff=fidiff)
        else
           error("wrong type receiver locations")
        end
    end

    # do simulation parallel
    if nprocs() == 1
       for i = 1 : ns
           wrap_born_forward(argument_collection[i])
       end

    else

       pmap(wrap_born_forward, argument_collection)
    end

    return nothing
end


text:
module FastSplat

export @fastsplat

macro fastsplat(call)
    @assert call.head == :call
    f = call.args[1]
    args = call.args[2:end]

    # figure out which arguments are splats
    splats = map(arg->isa(arg,Expr) && arg.head==:(...), args)

    gen = define_generator(splats)
    return esc(call_generator(gen, f, args, splats))
end

# define a generated function of signature `gen(f, args...)`
# where every arg corresponds with the original arg,
# or the arg and its length if it was splatted
function define_generator(splats::Vector{Bool})
    # build an arglist of form `(arg1, ... argN, ::Type{Val{LenN}}, ...)`
    arglist = Union{Expr,Symbol}[]
    typevars = Symbol[]
    for (i,splat) in enumerate(splats)
        push!(arglist, Symbol(:arg, i))
        if splat
            typevar = Symbol(:Len, i)
            push!(arglist, :(::Type{Val{$typevar}}))
            push!(typevars, typevar)
        end
    end

    # build an invocation of the original function, with every splat expanded
    args = Any[]
    for (i,splat) in enumerate(splats)
        arg = Symbol(:arg, i)
        if splat
            len = Symbol(:Len, i)
            push!(args, Expr(:(...), :((Expr(:ref, $(QuoteNode(arg)), i) for i in 1:$len))))
        else
            push!(args, QuoteNode(arg))
        end
    end

    @gensym gen
    @eval begin
        @generated function $gen(f, $(arglist...)) where {$(typevars...)}
            args = ($(args...),)
            quote
                Base.@_inline_meta
                @inbounds ret = f($(args...))
                ret
            end
        end
    end
    return gen
end

# emit an expression to call a generator as defined by `define_generator`
function call_generator(gen, f, args, splats)
    call = Expr(:call, :($FastSplat.$gen), f)
    for (i,splat) in enumerate(splats)
        if splat
            arg = args[i].args[1]
            push!(call.args, arg)
            push!(call.args, :(Val{length($arg)}))
        else
            arg = args[i]
            push!(call.args, arg)
        end
    end
    return call
end

end


text:
using Base
using Test
using LinearAlgebra
using Surrogates

@testset "Kriging" begin
    @testset "ND" begin

        #WITHOUT ADD POINT
        x = [(1.0,2.0,3.0),(4.0,5.0,6.0),(7.0,8.0,9.0)]
        s = sample(3,[0,0],[10,10],UniformSample())
        x = Tuple.(s)
        y = [1.0,2.0,3.0]
        p = [1.0 1.0 1.0]
        theta = [2.0 2.0 2.0]
        my_k = Kriging(x,y,p,theta)
        est = my_k((1.0,2.0,3.0))
        std_err = std_error_at_point(my_k,(1.0,2.0,3.0))

        #WITH ADD POINT adding singleton
        x = [(1.0,2.0,3.0),(4.0,5.0,6.0),(7.0,8.0,9.0)]
        y = [1.0,2.0,3.0]
        p = [1.0 1.0 1.0]
        theta = [2.0 2.0 2.0]
        my_k = Kriging(x,y,p,theta)
        add_point!(my_k,(10.0,11.0,12.0),4.0)
        est = my_k((10.0,11.0,12.0))
        std_err = std_error_at_point(my_k,(10.0,11.0,12.0))
        @test std_err < 10^(-6)

        #WITH ADD POINT ADDING MORE
        x = [(1.0, 2.0, 3.0),(4.0,5.0,6.0),(7.0, 8.0, 9.0)]
        y = [1.0,2.0,3.0]
        p = [1.0 1.0 1.0]
        theta = [2.0 2.0 2.0]
        my_k = Kriging(x,y,p,theta)
        add_point!(my_k,[(10.0, 11.0, 12.0),(13.0,14.0,15.0)],[4.0,5.0])
        est = my_k((10.0,11.0,12.0))
        std_err = std_error_at_point(my_k,(10.0,11.0,12.0))
        @test std_err < 10^(-6)

    end
    @testset "1D" begin
        #WITHOUT ADD POINT
        x = [1.0,2.0,3.0]
        y = [4.0,5.0,6.0]
        p = 1.3
        my_k = Kriging(x,y,p)
        est = my_k(1.0)
        std_err = std_error_at_point(my_k,1.0)
        @test std_err < 10^(-6)

        #WITH ADD POINT adding singleton
        x = [1.0,2.0,3.0]
        y = [4.0,5.0,6.0]
        p = 1.3
        my_k = Kriging(x,y,p)
        add_point!(my_k,4.0,9.0)
        est = my_k(4.0)
        std_err = std_error_at_point(my_k,4.0)
        @test std_err < 10^(-6)

        #WITH ADD POINT adding more
        x = [1.0,2.0,3.0]
        y = [4.0,5.0,6.0]
        p = 1.3
        my_k = Kriging(x,y,p)
        add_point!(my_k,[4.0,5.0,6.0],[9.0,13.0,15.0])
        est = my_k(4.0)
        std_err = std_error_at_point(my_k,4.0)
        @test std_err < 10^(-6)
    end

end


text:
"""
======================== VORTEX SHEET BOUNDARY =================================
HEADER FILE

Vortex Sheet Boundary method for coupling with Vortex Particle Method. Written
Julia 1.0

See README.md for details on how to install and docs/Theory.pdf for derivation.

@author: Damyn Chipman - 2019
================================================================================
"""
module VSB

# ===== Imports =====
using SciTools
using QuadGK
using LinearAlgebra
using SimpleVPM
SVPM = SimpleVPM

# ===== Exports =====
export Boundary, Parameters
export CalcRhoCoefs, CalcRho
export Circle, NACA4
export RBF_gauss
export CalcVSCoefs, CalcVS
export CalcVSDiffusionCoefs, CalcDiffusion

# ===== GLOBAL VARIABLES =====
const NUMB_MATRIX = Array{Array{T,1},1} where {T<:Real} # Numerical matrix type
#const module_path = splitdir(@__file__)[1]             # File path to module VSB.jl

# ===== Files =====
file_names = ["RBF",
              "Boundary",
              "CalcVortexSheet",
              "VortexDiffusion",
              "Parameters"]
for header_name in file_names
    include("VSB_"*header_name*".jl")
end

# ========================= END OF VSB.jl MODULE ===============================
end


text:
# Advent of Code, Day 9
#
# Find lowest points in a 2-dimensional space, and the areas around these low
# points, bounded by edges or 9s.
#
# AK, 9 Dec 2021

# Read input, convert to list of vectors, one per row
f = "test_input.txt"
f = "input.txt"
d = split(strip(read(f, String)), "\n")
m = map(r -> [x - '0' for x in r], d)
#display(m)
#println("\n")

# Initialize lists of basins: every basin is centred around a low point, and
# consists of list of points around that point (a point is [row, col] tuple).
basins = Vector{Vector{Vector{Int64}}}()

# Part 1: Find the low points - the locations that are lower than any of its adjacent
# locations -- and add together their "risk", i.e., the value in the cell plus 1.
# Also initialize a list of basins around each low point, and set the value of the low
# point to -1 to set up for part 2.
function part1()

    totalRisk = 0
    for ri in 1:length(m)
        r = m[ri]
        for ci in 1:length(r)

            # If this is this a low point, add up total risk,
            # and set cell to -1 and initialize basin to set up for part 2
            x = r[ci]
            lowX = (ci == 1 || x < r[ci-1]) && (ci == length(r) || x < r[ci+1])
            lowY = (ri == 1 || x < m[ri-1][ci]) && (ri == length(m) || x < m[ri+1][ci])
            if lowX && lowY
                totalRisk += x + 1
                m[ri][ci] = -1   # mark the low point, first in this basin
                #println("Found low point at $ri,$ci")
                push!(basins, [[ri,ci]])  # create start of a new basin
            end
        end
    end

    println("Part 1: total risk = $totalRisk")

end

# Part 2: find the size of the three largest "basins" around the low points, i.e.,
# the areas around low points bounded by edges or 9s, and multiply these three
# sizes together. Since a basin can include points adjacent to a cell found by
# search left/right/up/down from the original low point, need to iterate until
# no more points are found. Then, just add up the lengths of each basin and find
# the highest values.
function part2()

    # Find all basins
    println("Part 2: finding adjacent points")
    iter = 1
    while findBasins()
        println("  iteration ", iter)
        iter += 1
    end

     # Sort the basin lengths, and multiply the top 3 together
     #println("Ending basins: ")
     #for b in basins
     #   println(b)
     #end
     blens = [length(b) for b in basins]
     blens = reverse(sort(blens))
     println("Sorted basin lengths: ", blens)
     result = blens[1] * blens[2] * blens[3]
     println("Part 2 result = $result")
end

# For part 2, find basins iteratively, i.e., areas around -1 bounded by
# 9s or edges. Looks for every -1, and expands in all four directions,
# assigning each point to the basin of the starting point.
# Returns true if there were changes, so we know when to stop iterating.
function findBasins()

    changed = false  # if something changed, will need to iterate further
    for ri in 1:length(m)
        r = m[ri]
        for ci in 1:length(r)

            # Skip if this point is not within a basin
            if m[ri][ci] != -1
                continue
            end

            # Find the basin to which this point belongs
            p = [ri, ci]
            basin = nothing
            for b in basins
                if in(p, b)
                    basin = b
                    break
                end
            end

            # If no basin was found, something is wrong,
            # so stop here (does not happen with data provided)
            if basin == nothing
                println("No basin found for", p)
                return false  # to stop iterations
            end

            # Check points to the left/right/up/down of this point,
            # and mark/save any points until we reach a boundary,
            # i.e., edge of matrix or a 9

            # Check to right of this point
            i = ci + 1
            while i <= length(r) && r[i] != 9 && r[i] != -1
                r[i] = -1  # mark point as -1
                push!(basin, [ri, i])
                i += 1
                changed = true
            end

            # Check to left of this point
            i = ci - 1
            while i >= 1 && r[i] != 9 && r[i] != -1
                r[i] = -1
                push!(basin, [ri, i])
                i -= 1
                changed = true
            end

            # Check above this point
            i = ri - 1
            while i >= 1 && m[i][ci] != 9 && m[i][ci] != -1
                m[i][ci] = -1
                push!(basin, [i, ci])
                i -= 1
                changed = true
            end

            # Check below this point
            i = ri + 1
            while i <= length(m) && m[i][ci] != 9 && m[i][ci] != -1
                m[i][ci] = -1
                push!(basin, [i,  ci])
                i += 1
                changed = true
            end
        end
    end

    # Return status indicates whether something changed, so will need
    # to iterate further until there are no more changes.
    return changed

end

# Execute both parts (need to do part 1 in order to do part 2)
part1()
part2()


text:
#=
In Julia, as a developer you can define your own types to structure data used in
applications. For Example, if you need to represent points in a three-dimensional
space, you can define a type Point, as follows:
=#

mutable struct Point
    x::Float64
    y::Float64
    z::Float64
end

struct Vector3D
    x::Float64
    y::Float64
    z::Float64
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
@testset "Special Functions" begin
    @testset "Gamma Function" begin
        @test gamma(Arb(2)) ≈ gamma(2)
        @test gamma(Arb(3)) ≈ gamma(3)
        @test gamma(Acb(2 + 2im)) ≈ gamma(2 + 2im)
        @test gamma(Acb(3 + 3im)) ≈ gamma(3 + 3im)
        @test gamma(ArbSeries([2, 1]))[0] ≈ gamma(2)
        @test gamma(AcbSeries([2 + 2im, 1]))[0] ≈ gamma(2 + 2im)

        @test loggamma(Arb(2)) ≈ loggamma(2)
        @test loggamma(Arb(3)) ≈ loggamma(3)
        @test loggamma(Acb(2 + 2im)) ≈ loggamma(2 + 2im)
        @test loggamma(Acb(3 + 3im)) ≈ loggamma(3 + 3im)
        @test loggamma(ArbSeries([2, 1]))[0] ≈ loggamma(2)
        @test loggamma(AcbSeries([2 + 2im, 1]))[0] ≈ loggamma(2 + 2im)

        @test digamma(Arb(2)) ≈ digamma(2)
        @test digamma(Arb(3)) ≈ digamma(3)
        @test digamma(Acb(2 + 2im)) ≈ digamma(2 + 2im)
        @test digamma(Acb(3 + 3im)) ≈ digamma(3 + 3im)
        @test digamma(ArbSeries([2, 1]))[0] ≈ digamma(2)
        @test digamma(AcbSeries([2 + 2im, 1]))[0] ≈ digamma(2 + 2im)

        @test trigamma(Acb(2 + 2im)) ≈ trigamma(2 + 2im)
        @test trigamma(Acb(3 + 3im)) ≈ trigamma(3 + 3im)

        @test polygamma(Acb(2), Acb(3 + 3im)) ≈ polygamma(2, 3 + 3im)
        @test polygamma(Acb(3), Acb(4 + 4im)) ≈ polygamma(3, 4 + 4im)

        @test gamma(Arb(2), Arb(3)) ≈ gamma(2, 3)
        @test gamma(Arb(3), Arb(4)) ≈ gamma(3, 4)
        @test gamma(Acb(2 + 2im), Acb(3 + 3im)) ≈ gamma(2 + 2im, 3 + 3im)
        @test gamma(Acb(3 + 3im), Acb(4 + 4im)) ≈ gamma(3 + 3im, 4 + 4im)
        @test gamma(Arb(2), ArbSeries([3, 1]))[0] ≈ gamma(2, 3)
        @test gamma(Acb(2 + 2im), AcbSeries([3 + 3im, 1]))[0] ≈ gamma(2 + 2im, 3 + 3im)

        @test all(gamma_inc(Arb(2), Arb(3)) .≈ gamma_inc(2, 3))
        @test all(gamma_inc(Arb(3), Arb(4)) .≈ gamma_inc(3, 4))
        @test all(gamma_inc(Acb(2), Acb(3)) .≈ gamma_inc(2, 3))
        @test all(gamma_inc(Acb(3), Acb(4)) .≈ gamma_inc(3, 4))
        @test all(getindex.(gamma_inc(Arb(2), ArbSeries([3, 1])), 0) .≈ gamma_inc(2, 3))
        @test all(getindex.(gamma_inc(Acb(2), AcbSeries([3, 1])), 0) .≈ gamma_inc(2, 3))

        @test all(beta_inc(Arb(2), Arb(3), Arb(1 // 2)) .≈ beta_inc(2.0, 3.0, 0.5))
        @test all(beta_inc(Arb(3), Arb(4), Arb(1 // 4)) .≈ beta_inc(3.0, 4.0, 0.25))
        @test all(beta_inc(Acb(2), Acb(3), Acb(1 // 2)) .≈ beta_inc(2.0, 3.0, 0.5))
        @test all(beta_inc(Acb(3), Acb(4), Acb(1 // 4)) .≈ beta_inc(3.0, 4.0, 0.25))
        @test all(
            getindex.(beta_inc(Arb(2), Arb(3), ArbSeries([1 // 2, 1])), 0) .≈
            beta_inc(2.0, 3.0, 0.5),
        )
        @test all(
            getindex.(beta_inc(Acb(2), Acb(3), AcbSeries([1 // 2, 1])), 0) .≈
            beta_inc(2.0, 3.0, 0.5),
        )
    end

    @testset "Trigonometric Integrals" begin
        @test expint(Arb(2), Arb(3)) ≈ expint(2, 3)
        @test expint(Arb(3), Arb(4)) ≈ expint(3, 4)
        @test expint(Acb(2 + 2im), Acb(3 + 3im)) ≈ expint(2 + 2im, 3 + 3im)
        @test expint(Acb(3 + 3im), Acb(4 + 4im)) ≈ expint(3 + 3im, 4 + 4im)
        @test expint(Arb(2), ArbSeries([3, 1]))[0] ≈ expint(2, 3)
        @test expint(Acb(2 + 2im), AcbSeries([3 + 3im, 1]))[0] ≈ expint(2 + 2im, 3 + 3im)

        @test expint(Arb(2)) ≈ expint(2)
        @test expint(Arb(3)) ≈ expint(3)
        @test expint(Acb(2 + 2im)) ≈ expint(2 + 2im)
        @test expint(Acb(3 + 3im)) ≈ expint(3 + 3im)
        @test expint(ArbSeries([2, 1]))[0] ≈ expint(2)
        @test expint(AcbSeries([2 + 2im, 1]))[0] ≈ expint(2 + 2im)

        @test expinti(Arb(2)) ≈ expinti(2)
        @test expinti(Arb(3)) ≈ expinti(3)
        @test expinti(Acb(2)) ≈ expinti(2)
        @test expinti(Acb(3)) ≈ expinti(3)
        @test expinti(ArbSeries([2, 1]))[0] ≈ expinti(2)
        @test expinti(AcbSeries([2, 1]))[0] ≈ expinti(2)

        @test sinint(Arb(2)) ≈ sinint(2)
        @test sinint(Arb(3)) ≈ sinint(3)
        @test sinint(Acb(2)) ≈ sinint(2)
        @test sinint(Acb(3)) ≈ sinint(3)
        @test sinint(ArbSeries([2, 1]))[0] ≈ sinint(2)
        @test sinint(AcbSeries([2, 1]))[0] ≈ sinint(2)

        @test cosint(Arb(2)) ≈ cosint(2)
        @test cosint(Arb(3)) ≈ cosint(3)
        @test cosint(Acb(2)) ≈ cosint(2)
        @test cosint(Acb(3)) ≈ cosint(3)
        @test cosint(ArbSeries([2, 1]))[0] ≈ cosint(2)
        @test cosint(AcbSeries([2, 1]))[0] ≈ cosint(2)
    end

    @testset "Error Functions" begin
        @test erf(Arb(2)) ≈ erf(2)
        @test erf(Arb(3)) ≈ erf(3)
        @test erf(Acb(2 + 2im)) ≈ erf(2 + 2im)
        @test erf(Acb(3 + 3im)) ≈ erf(3 + 3im)
        @test erf(ArbSeries([2, 1]))[0] ≈ erf(2)
        @test erf(AcbSeries([2 + 2im, 1]))[0] ≈ erf(2 + 2im)

        @test erfc(Arb(2)) ≈ erfc(2)
        @test erfc(Arb(3)) ≈ erfc(3)
        @test erfc(Acb(2 + 2im)) ≈ erfc(2 + 2im)
        @test erfc(Acb(3 + 3im)) ≈ erfc(3 + 3im)
        @test erfc(ArbSeries([2, 1]))[0] ≈ erfc(2)
        @test erfc(AcbSeries([2 + 2im, 1]))[0] ≈ erfc(2 + 2im)

        @test erfi(Arb(2)) ≈ erfi(2)
        @test erfi(Arb(3)) ≈ erfi(3)
        @test erfi(Acb(2 + 2im)) ≈ erfi(2 + 2im)
        @test erfi(Acb(3 + 3im)) ≈ erfi(3 + 3im)
        @test erfi(ArbSeries([2, 1]))[0] ≈ erfi(2)
        @test erfi(AcbSeries([2 + 2im, 1]))[0] ≈ erfi(2 + 2im)
    end

    @testset "Airy Functions" begin
        @test airyai(Arb(2)) ≈ airyai(2)
        @test airyai(Arb(3)) ≈ airyai(3)
        @test airyai(Acb(2 + 2im)) ≈ airyai(2 + 2im)
        @test airyai(Acb(3 + 3im)) ≈ airyai(3 + 3im)
        @test airyai(ArbSeries([2, 1]))[0] ≈ airyai(2)
        @test airyai(AcbSeries([2 + 2im, 1]))[0] ≈ airyai(2 + 2im)

        @test airyaiprime(Arb(2)) ≈ airyaiprime(2)
        @test airyaiprime(Arb(3)) ≈ airyaiprime(3)
        @test airyaiprime(Acb(2 + 2im)) ≈ airyaiprime(2 + 2im)
        @test airyaiprime(Acb(3 + 3im)) ≈ airyaiprime(3 + 3im)
        @test airyaiprime(ArbSeries([2, 1]))[0] ≈ airyaiprime(2)
        @test airyaiprime(AcbSeries([2 + 2im, 1]))[0] ≈ airyaiprime(2 + 2im)

        @test airybi(Arb(2)) ≈ airybi(2)
        @test airybi(Arb(3)) ≈ airybi(3)
        @test airybi(Acb(2 + 2im)) ≈ airybi(2 + 2im)
        @test airybi(Acb(3 + 3im)) ≈ airybi(3 + 3im)
        @test airybi(ArbSeries([2, 1]))[0] ≈ airybi(2)
        @test airybi(AcbSeries([2 + 2im, 1]))[0] ≈ airybi(2 + 2im)

        @test airybiprime(Arb(2)) ≈ airybiprime(2)
        @test airybiprime(Arb(3)) ≈ airybiprime(3)
        @test airybiprime(Acb(2 + 2im)) ≈ airybiprime(2 + 2im)
        @test airybiprime(Acb(3 + 3im)) ≈ airybiprime(3 + 3im)
        @test airybiprime(ArbSeries([2, 1]))[0] ≈ airybiprime(2)
        @test airybiprime(AcbSeries([2 + 2im, 1]))[0] ≈ airybiprime(2 + 2im)
    end

    @testset "Bessel Functions" begin
        @test besselj(Arb(2), Arb(3)) ≈ besselj(2, 3)
        @test besselj(Arb(3), Arb(4)) ≈ besselj(3, 4)
        @test besselj(Acb(2), Acb(3 + 3im)) ≈ besselj(2, 3 + 3im)
        @test besselj(Acb(3), Acb(4 + 4im)) ≈ besselj(3, 4 + 4im)

        @test besselj0(Arb(2)) ≈ besselj0(2)
        @test besselj0(Arb(3)) ≈ besselj0(3)
        @test besselj0(Acb(2 + 2im)) ≈ besselj0(2 + 2im)
        @test besselj0(Acb(3 + 3im)) ≈ besselj0(3 + 3im)

        @test besselj1(Arb(2)) ≈ besselj1(2)
        @test besselj1(Arb(3)) ≈ besselj1(3)
        @test besselj1(Acb(2 + 2im)) ≈ besselj1(2 + 2im)
        @test besselj1(Acb(3 + 3im)) ≈ besselj1(3 + 3im)

        @test sphericalbesselj(Arb(2), Arb(3)) ≈ sphericalbesselj(2, 3)
        @test sphericalbesselj(Arb(3), Arb(4)) ≈ sphericalbesselj(3, 4)

        @test bessely(Arb(2), Arb(3)) ≈ bessely(2, 3)
        @test bessely(Arb(3), Arb(4)) ≈ bessely(3, 4)
        @test bessely(Acb(2), Acb(3 + 3im)) ≈ bessely(2, 3 + 3im)
        @test bessely(Acb(3), Acb(4 + 4im)) ≈ bessely(3, 4 + 4im)

        @test bessely0(Arb(2)) ≈ bessely0(2)
        @test bessely0(Arb(3)) ≈ bessely0(3)
        @test bessely0(Acb(2 + 2im)) ≈ bessely0(2 + 2im)
        @test bessely0(Acb(3 + 3im)) ≈ bessely0(3 + 3im)

        @test bessely1(Arb(2)) ≈ bessely1(2)
        @test bessely1(Arb(3)) ≈ bessely1(3)
        @test bessely1(Acb(2 + 2im)) ≈ bessely1(2 + 2im)
        @test bessely1(Acb(3 + 3im)) ≈ bessely1(3 + 3im)

        @test besselh(Arb(2), Arb(3)) ≈ besselh(2, 3)
        @test besselh(Arb(3), Arb(4)) ≈ besselh(3, 4)
        @test besselh(Acb(2), Acb(3)) ≈ besselh(2, 3)
        @test besselh(Acb(3), Acb(4)) ≈ besselh(3, 4)
        @test besselh(Arb(2), 1, Arb(3)) ≈ besselh(2, 1, 3)
        @test besselh(Arb(3), 1, Arb(4)) ≈ besselh(3, 1, 4)
        @test besselh(Acb(2), 1, Acb(3)) ≈ besselh(2, 1, 3)
        @test besselh(Acb(3), 1, Acb(4)) ≈ besselh(3, 1, 4)
        @test besselh(Arb(2), 2, Arb(3)) ≈ besselh(2, 2, 3)
        @test besselh(Arb(3), 2, Arb(4)) ≈ besselh(3, 2, 4)
        @test besselh(Acb(2), 2, Acb(3)) ≈ besselh(2, 2, 3)
        @test besselh(Acb(3), 2, Acb(4)) ≈ besselh(3, 2, 4)
        @test_throws SpecialFunctions.AmosException(1) besselh(Arb(2), 3, Arb(3))
        @test_throws SpecialFunctions.AmosException(1) besselh(Acb(2), 3, Acb(3))

        @test besseli(Arb(2), Arb(3)) ≈ besseli(2, 3)
        @test besseli(Arb(3), Arb(4)) ≈ besseli(3, 4)
        @test besseli(Acb(2), Acb(3 + 3im)) ≈ besseli(2, 3 + 3im)
        @test besseli(Acb(3), Acb(4 + 4im)) ≈ besseli(3, 4 + 4im)

        @test besselk(Arb(2), Arb(3)) ≈ besselk(2, 3)
        @test besselk(Arb(3), Arb(4)) ≈ besselk(3, 4)
        @test besselk(Acb(2), Acb(3 + 3im)) ≈ besselk(2, 3 + 3im)
        @test besselk(Acb(3), Acb(4 + 4im)) ≈ besselk(3, 4 + 4im)

        @test besselkx(Arb(2), Arb(3)) ≈ besselkx(2, 3)
        @test besselkx(Arb(3), Arb(4)) ≈ besselkx(3, 4)
        @test besselkx(Acb(2), Acb(3 + 3im)) ≈ besselkx(2, 3 + 3im)
        @test besselkx(Acb(3), Acb(4 + 4im)) ≈ besselkx(3, 4 + 4im)
    end

    @testset "Elliptic Integrals" begin
        @test ellipk(Acb(-1)) ≈ ellipk(-1)
        @test ellipk(Acb(-2)) ≈ ellipk(-2)
        @test ellipk(AcbSeries([-1, 1]))[0] ≈ ellipk(-1)

        @test ellipe(Acb(-1)) ≈ ellipe(-1)
        @test ellipe(Acb(-2)) ≈ ellipe(-2)
    end

    @testset "Zeta Functions" begin
        @test eta(Acb(2 + 2im)) ≈ eta(2 + 2im)
        @test eta(Acb(3 + 3im)) ≈ eta(3 + 3im)

        @test zeta(Arb(2)) ≈ zeta(2)
        @test zeta(Arb(3)) ≈ zeta(3)
        @test zeta(Acb(2 + 2im)) ≈ zeta(2 + 2im)
        @test zeta(Acb(3 + 3im)) ≈ zeta(3 + 3im)
        @test zeta(ArbSeries([2, 1]))[0] ≈ zeta(2)
        @test zeta(AcbSeries([2 + 2im, 1]))[0] ≈ zeta(2 + 2im)
        @test zeta(Arb(2), Arb(3)) ≈ zeta(2, 3)
        @test zeta(Arb(3), Arb(4)) ≈ zeta(3, 4)
        @test zeta(Acb(2 + 2im), Acb(3 + 3im)) ≈ zeta(2 + 2im, 3 + 3im)
        @test zeta(Acb(3 + 3im), Acb(4 + 4im)) ≈ zeta(3 + 3im, 4 + 4im)
        @test zeta(ArbSeries([2, 1]), Arb(3))[0] ≈ zeta(2, 3)
        @test zeta(AcbSeries([2 + 2im, 1]), Acb(3 + 3im))[0] ≈ zeta(2 + 2im, 3 + 3im)
    end
end


text:
# Bogoliubov-de Gennes problems

# Idea: find eigenvalues of B^2, with :SR.  The positive eigenvalues of B are the square roots.  Find eigenvectors of B belonging to those eigenvalues, then flip the modes with negative norms.
# TODO investigate why Arpack fails to find the eigenvalues.

"""
    ωs, us, vs = bdg_modes(s, d, ψ, Ω, nmodes; nev=nmodes, raw=false)

Return vectors of the lowest energy BdG modes

TODO Normalised with ∥u∥² - ∥v∥² = 1.
"""
function bdg_modes(s, d, ψ, Ω, nmodes; nev = nmodes, raw = false)
    # TODO how many hartree modes to include?
    # write test to double nev, compare first half of modes
    @info "starting"
    _, us = hartree_modes(s, d, ψ, nev, Ω)
    U = hcat((u[:] for u in us)...)
    @info "expanding over modes"
    B = expand(bdg_operator(s, d, ψ, Ω), us, conj.(us))
    @info "diagonizing"
    ew, ev = eigs(B; nev = 2nmodes, which = :SM)[1:2]
    @info "finishing up"
    ev = [U zero(U); zero(U) conj(U)] * ev
    if raw
        ew, ev
    else
        bdg_output(d, ew, ev)
    end
end

"""
    bdg_output(d, ew, ev; safe=true)

Select the positive norm modes, and reshape evs to 2D grid
"""
function bdg_output(d, ew, ev; safe=true, a_tol=sqrt(eps()))
    us = [reshape(ev[1:d.n^2, j], d.n, d.n) for j in eachindex(ew)]
    vs = [reshape(ev[d.n^2+1:end, j], d.n, d.n) |> conj for j in eachindex(ew)]
    safe || return ew, us, vs
    ixs = @. norm(us) > norm(vs)
    # There are 3 cases:
    # ω ≈ 0, norm(u) ≈ norm(v)
    # ω real, norm(u) > norm(v)
    # ω imaginary, v = exp(iφ)*u
    # TODO think about which degenerate pair to take in the imaginary case
    if count(ixs) != length(ew) ÷ 2
        for j = 1:2:length(ixs)
            ixs[j] ⊻ ixs[j+1] && continue
            @assert all(@. abs(real(ew[j:j+1])) < a_tol)
            ixs[j] = !ixs[j]
        end
    end
    ew[ixs], us[ixs], vs[ixs]
end

"expand f(u) as a matrix over us"
function expand(f, us)
    [dot(us[j], f(us[k])) for j in eachindex(us), k in eachindex(us)]
end

"expand u', v' = f(u,v) as a matrix over us, vs"
function expand(f, us, vs)
    x(g, as, bs) = [dot(as[j], g(bs[k])) for j in eachindex(as), k in eachindex(bs)]
    [
        x(u -> f(u, zero(u))[1], us, us) x(v -> f(zero(v), v)[1], us, vs)
        x(u -> f(u, zero(u))[2], vs, us) x(v -> f(zero(v), v)[2], vs, vs)
    ]
end

"""
    bdg_operator(s,d,ψ,Ω)

v is conjugated
"""
function bdg_operator(s, d, ψ, Ω)
    T, V, U, J, L = operators(s, d, :T, :V, :U, :J, :L)

    μ = dot(L(ψ; Ω), ψ) |> real
    (u, v) -> (
        T(u) + V(u) + 2U(u, abs2.(ψ)) - μ * u - Ω * J(u) - U(v, ψ .^ 2),
        U(u, conj(ψ) .^ 2) - T(v) - V(v) - 2U(v, abs2.(ψ)) + μ * v - Ω * J(v),
    )
end

@defaults BdGmatrix(s::Superfluid, d::Discretisation, Ω, ψ)

"""
    BdGmatrix([s], [d], Ω, ψ)

Return the matrix of the BdG operator.
"""
BdGmatrix

function BdGmatrix(s::Superfluid, d::Discretisation, Ω, ψ)
    # TODO make this a BlockBandedMatrix
    T, V, U, J, L = matrices(s, d, Ω, ψ, :T, :V, :U, :J, :L)
    Q = diagm(0 => ψ[:])
    μ = sum(conj(ψ[:]) .* (L * ψ[:])) |> real |> UniformScaling
    [
        T+V+2U-μ-Ω*J -s.C/d.h^2*Q .^ 2
        s.C/d.h^2*conj.(Q) .^ 2 -T-V-2U+μ-Ω*J
    ]
end

"Expand u over us, v* over us*"
function BdGmatrix(s, d, Ω, ψ, us)
    T!, V!, U!, J!, L = operators(s, d, :T!, :V!, :U!, :J!, :L)
    nev = length(us)
    M = Array{Complex{Float64}}(undef, 2nev, 2nev)
    μ = dot(ψ, L(ψ; Ω)) |> real
    norm(L(ψ; Ω) - μ * ψ) > 1e-3 &&
        @warn "Non-stationary order parameter" residual = norm(L(ψ; Ω) - μ * ψ)
    w = similar(us[1])
    for j = 1:nev, k = 1:j
        # Diagonal blocks
        u = us[k]
        w .= 0
        T!(w, u)
        V!(w, u)
        U!(w, 2, ψ, u)
        a = dot(us[j], w)
        j == k && (a -= μ)

        w .= 0
        b = dot(us[j], J!(w, Ω, u))
        M[j, k] = a - b
        M[k, j] = conj(a - b)
        # conj(v) ⋅L(conj(u))) = conj(u ⋅L(v))
        M[k+nev, j+nev] = -a - b
        M[j+nev, k+nev] = conj(-a - b)
    end

    for j = 1:nev, k = 1:nev
        # Off-diagonal blocks
        M[j, k+nev] = -s.C / d.h^2 * dot(us[j], @. ψ^2 * us[k])
        M[j+nev, k] = s.C / d.h^2 * dot(conj(us[j]), @. conj(ψ)^2 * conj(us[k]))
    end

    M
end

"""
    ωs, us = hartree_modes(s, d, ψ, n, Ω; raw=false)

Return the lowest `n` hartree eigenfunctions

Modes are returned as a vector of N-dimensional arrays.  If `raw` is set true, a unitary matrix of [u1[:] u2[:] ...] is returned instead.
"""
function hartree_modes end

function hartree_modes(s, d, ψ, nev, Ω; raw = false)
    L = operators(s, d, :L) |> only
    Lop = LinearMap{Complex{Float64}}(d.n^2) do u
        u = reshape(u, d.n, d.n)
        L(u, abs2.(ψ); Ω)[:]
    end
    ws, us = eigs(Lop; nev, which = :SR)[1:2]
    @info "Maximum imag eigenvalue" im = norm(imag(ws), Inf)
    ws = real(ws)
    raw || (us = [reshape(u, d.n, d.n) for u in eachcol(us)])
    ws, us
end

"Matrix whose columns are a basis for the space orthogonal to ψ"
function householder1(ψ)
    ψ = ψ[:]
    n = length(ψ)
    u = ψ[1]/abs(ψ[1])
    v = ψ + [u*norm(ψ); zeros(n-1)]
    normalize!(v)
    R = Matrix(I, n, n) - 2*v*v'
    R[:, 2:end]
end

function householder2(ψ)
    R = householder1(ψ)
    S = householder1(conj(ψ))
    Z = zero(R)
    [R Z; Z S]
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 4fb1c47c-8d70-11eb-3d5c-01c5bb7190fe
using Pkg, DrWatson

# ╔═╡ 5838d806-8d70-11eb-2a78-fb586f056012
begin
	using Distributions
	using StatsPlots
	using StatsBase
	using LaTeXStrings
	using CSV
	using DataFrames
	using LinearAlgebra
	using Random
	using StatisticalRethinking
end

# ╔═╡ 3562c0e2-8d70-11eb-0f89-9d039fe996b9
md" ## Clip-11-01-05s.jl"

# ╔═╡ 07c85e9e-8d95-11eb-295d-63a99f7aa921
md" ##### Fig 10.7."

# ╔═╡ 9bfd33fc-8d70-11eb-3a9c-bb9a990cbcf9
x = collect(-1.0:0.01:1.0);

# ╔═╡ ee25cbe6-8d8c-11eb-3323-c521a67b10f0
begin
	α = 0.0
	β = 2.0

	# p = exp.(α .+ β .* x) ./ (1 .+ exp.(α .+ β .* x))

	p = logistic.(α .+ β .* x)
end

# ╔═╡ 1cba8970-8d91-11eb-01d8-09d8d076dc73
begin
	odds = p ./ (1 .- p)
	logodds = log.(odds)
end

# ╔═╡ 6b96d95e-8d91-11eb-25ac-75f21ee9c98e
logistic.(x)

# ╔═╡ 15dcde1e-8d8c-11eb-0964-d9c5f616e59c
begin
	fig1 = plot(x, logodds, leg=false)
	fig2 = plot(x, p, ylims=(0,1), leg=false)
	plot(fig1, fig2, layout=(1, 2))
end

# ╔═╡ 1404a820-8d95-11eb-27d2-531cea578fda
md" ##### Fig 10.8."

# ╔═╡ c859986c-8d93-11eb-1211-db9c435de034
logsig = exp.(α .+ β .* x)

# ╔═╡ f3bac6a4-8d93-11eb-2614-1572443e3415
begin
	fig3 = plot(x, log.(logsig), leg=false)
	fig4 = plot(x, logsig, ylims=(0,10), leg=false)
	plot(fig3, fig4, layout=(1, 2))
end

# ╔═╡ 7e069492-8d70-11eb-2093-8375b07f4099
md" ## End of clip-11-01-05s.jl"

# ╔═╡ Cell order:
# ╠═3562c0e2-8d70-11eb-0f89-9d039fe996b9
# ╠═4fb1c47c-8d70-11eb-3d5c-01c5bb7190fe
# ╠═5838d806-8d70-11eb-2a78-fb586f056012
# ╟─07c85e9e-8d95-11eb-295d-63a99f7aa921
# ╠═9bfd33fc-8d70-11eb-3a9c-bb9a990cbcf9
# ╠═ee25cbe6-8d8c-11eb-3323-c521a67b10f0
# ╠═1cba8970-8d91-11eb-01d8-09d8d076dc73
# ╠═6b96d95e-8d91-11eb-25ac-75f21ee9c98e
# ╠═15dcde1e-8d8c-11eb-0964-d9c5f616e59c
# ╟─1404a820-8d95-11eb-27d2-531cea578fda
# ╠═c859986c-8d93-11eb-1211-db9c435de034
# ╠═f3bac6a4-8d93-11eb-2614-1572443e3415
# ╟─7e069492-8d70-11eb-2093-8375b07f4099


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

name = "SharedMimeInfo"
version = v"1.12"

sources = [
    "https://gitlab.freedesktop.org/xdg/shared-mime-info/uploads/80c7f1afbcad2769f38aeb9ba6317a51/shared-mime-info-$(version.major).$(version.minor).tar.xz" =>
    "18b2f0fe07ed0d6f81951a5fd5ece44de9c8aeb4dc5bb20d4f595f6cc6bd403e",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/shared-mime-info-*/
apk add intltool

./configure --prefix=$prefix --host=$target
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    ExecutableProduct("update-mime-database", :update_mime_database),
    FileProduct("share/locale", :locale_dir),
]

# Dependencies that must be installed before this package can be built
# Based on http://www.linuxfromscratch.org/blfs/view/8.3/general/shared-mime-info.html
dependencies = [
    "Glib_jll",
    "XML2_jll",
]

# Build the tarballs, and possibly a `build.jl` as well.  We use GCC 8 because it is the only GCC version that links
# properly on powerpc64le.  Shocking, I know, but this is the world we live in.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"8")


text:
using RegistryTools
using Test

@testset "RegistryTools" begin

include("regedit.jl")

end


text:
using DataFrames
import Unitful
using Dates: Date

@testset "unit" begin
    @testset "unit" begin
        @system SUnit(Controller) begin
            a => 2 ~ track(u"m")
            b => 1 ~ track(u"s")
            c(a, b) => a / b ~ track(u"m/s")
        end
        s = instance(SUnit)
        @test s.a' == 2u"m" && s.b' == 1u"s" && s.c' == 2u"m/s"
    end

    @testset "unitless" begin
        @system SUnitless(Controller) begin
            a => 200 ~ track(u"cm^2")
            b => 10 ~ track(u"m^2")
            c(a, b) => a / b ~ track(u"cm^2/m^2")
            d(a, b) => a / b ~ track(u"m^2/m^2")
            e(a, b) => a / b ~ track
        end
        s = instance(SUnitless)
        @test s.a' == 200u"cm^2"
        @test s.b' == 10u"m^2"
        @test s.c' == 20u"cm^2/m^2"
        @test s.d' == 0.002
        @test s.e' == 0.002
    end

    @testset "nounit" begin
        @system SNounit(Controller) begin
            a => 1 ~ track(u"m")
            b(nounit(a)) ~ track
        end
        s = instance(SNounit)
        @test s.a' == u"1m"
        @test s.b' == 1
    end

    @testset "nounit with alias" begin
        @system SNounitAlias(Controller) begin
            a: aa => 1 ~ track(u"m")
            b(nounit(aa)): bb ~ track
        end
        s = instance(SNounitAlias)
        @test s.aa' == u"1m"
        @test s.bb' == 1
    end

    @testset "nounit with call" begin
        @system SNounitCall(Controller) begin
            a => 1 ~ track(u"m")
            b(nounit(a); x) => (a + x) ~ call
            c(b) => b(1) ~ track
        end
        s = instance(SNounitCall)
        @test s.a' == u"1m"
        @test s.b'(1) == 2
        @test s.c' == 2
    end

    @testset "nounit nested" begin
        @system SNounitNestedComponent begin
            a => 1 ~ track(u"m")
        end
        @system SNounitNested(Controller) begin
            n(context) ~ ::SNounitNestedComponent
            a(n.a) ~ track(u"m")
            b(nounit(n.a)) ~ track
            c(x=nounit(n.a)) => 2x ~ track
        end
        s = instance(SNounitNested)
        @test s.n.a' == s.a' == u"1m"
        @test s.b' == 1
        @test s.c' == 2
    end

    @testset "nothing" begin
        @test Cropbox.unitfy(nothing, u"m") === nothing
        @test Cropbox.unitfy(nothing, nothing) === nothing
    end

    @testset "missing" begin
        @test Cropbox.unitfy(missing, u"m") === missing
        @test Cropbox.unitfy(missing, nothing) === missing
    end

    @testset "nothing units" begin
        @test Cropbox.unitfy(1, nothing) === 1
        @test Cropbox.unitfy(1u"cm", nothing) === 1u"cm"
    end

    @testset "single" begin
        @test Cropbox.unitfy(1, u"m") === 1u"m"
        @test Cropbox.unitfy(1.0, u"m") === 1.0u"m"
        @test Cropbox.unitfy(1u"m", u"cm") === 100u"cm"
        @test Cropbox.unitfy(1.0u"m", u"cm") === 100.0u"cm"
    end

    @testset "percent" begin
        @test Cropbox.unitfy(1, u"percent") === 1u"percent"
        @test 1 |> u"percent" === 100u"percent"
        @test Cropbox.unitfy(1u"percent", nothing) === 1u"percent"
        @test Cropbox.unitfy(1u"percent", u"NoUnits") === 1//100
    end

    @testset "array" begin
        @test Cropbox.unitfy([1, 2, 3], u"m") == [1, 2, 3]u"m"
        @test Cropbox.unitfy([1, 2, 3]u"m", u"cm") == [100, 200, 300]u"cm"
        @test Cropbox.unitfy([1u"cm", 0.02u"m", 30u"mm"], u"cm") == [1, 2, 3]u"cm"
    end

    @testset "tuple" begin
        @test Cropbox.unitfy((1, 2, 3), u"m") === (1u"m", 2u"m", 3u"m")
        @test Cropbox.unitfy((1u"m", 2u"m", 3u"m"), u"cm") === (100u"cm", 200u"cm", 300u"cm")
        @test Cropbox.unitfy((1u"cm", 0.02u"m", 30u"mm"), u"cm") === (1u"cm", 2.0u"cm", 3//1*u"cm")
    end

    @testset "array with missing" begin
        @test isequal(Cropbox.unitfy([1, 2, missing], u"m"), [1, 2, missing]u"m")
        @test isequal(Cropbox.unitfy([1, 2, missing]u"m", u"cm"), [100, 200, missing]u"cm")
        @test isequal(Cropbox.unitfy([1u"cm", 0.02u"m", missing], u"cm"), [1, 2, missing]u"cm")
    end

    @testset "tuple with missing" begin
        @test isequal(Cropbox.unitfy((1, 2, missing), u"m"), (1u"m", 2u"m", missing))
        @test isequal(Cropbox.unitfy((1u"m", 2u"m", missing), u"cm"), (100u"cm", 200u"cm", missing))
        @test isequal(Cropbox.unitfy((1u"cm", 0.02u"m", missing), u"cm"), (1u"cm", 2.0u"cm", missing))
    end

    @testset "mixed units" begin
        #TODO: reduce subtlety
        @test_throws Unitful.DimensionError Cropbox.unitfy([1u"m", 2u"cm", 3], u"mm")
        @test Cropbox.unitfy(Any[1u"m", 2u"cm", 3], u"mm") == [1000u"mm", 20u"mm", 3u"mm"]
        @test Cropbox.unitfy((1u"m", 2u"cm", 3), u"mm") == (1000u"mm", 20u"mm", 3u"mm")
    end

    @testset "deunitfy" begin
        @test Cropbox.deunitfy(1) == 1
        @test Cropbox.deunitfy(1u"m") == 1
        @test Cropbox.deunitfy([1, 2, 3]u"m") == [1, 2, 3]
        @test Cropbox.deunitfy((1u"m", 2u"cm", 3)) === (1, 2, 3)
    end

    @testset "deunitfy with units" begin
        @test Cropbox.deunitfy(1, u"m") == 1
        @test Cropbox.deunitfy(1u"m", u"cm") == 100
        @test Cropbox.deunitfy([1, 2, 3]u"m", u"cm") == [100, 200, 300]
        @test Cropbox.deunitfy([1u"m", 2u"cm", 3u"mm"], u"mm") == [1000, 20, 3]
        @test Cropbox.deunitfy((1u"m", 2u"cm", 3u"mm"), u"mm") === (1000, 20, 3)
    end

    @testset "dataframe" begin
        df = DataFrame(a=[0], b=[0])
        U = [u"s", nothing]
        r = Cropbox.unitfy(df, U)
        @test Cropbox.unittype(eltype(r[!, 1])) == u"s"
        @test Cropbox.unittype(eltype(r[!, 2])) == u"NoUnits"
    end

    @testset "dataframe auto unit" begin
        df = DataFrame()
        z = [0]
        df."a (g/cm^2)" = z
        df."c (a)(b)(s)" = z
        df."b ()" = z
        df."d" = z
        r = Cropbox.unitfy(df)
        @test Cropbox.unittype(eltype(r[!, 1])) == u"g/cm^2"
        @test Cropbox.unittype(eltype(r[!, 2])) == u"s"
        @test Cropbox.unittype(eltype(r[!, 3])) == u"NoUnits"
        @test Cropbox.unittype(eltype(r[!, 4])) == u"NoUnits"
        N = names(r)
        @test N[1] == "a"
        @test N[2] == "c (a)(b)"
        @test N[3] == "b ()"
        @test N[4] == "d"
    end

    @testset "dataframe auto type" begin
        df = DataFrame()
        z = ["2020-10-08"]
        df."a (:Date)" = z
        df."b (:String)" = z
        df."c (:Symbol)" = z
        df."d ()" = z
        df."e" = z
        r = Cropbox.unitfy(df)
        @test eltype(r[!, 1]) == Date
        @test eltype(r[!, 2]) == String
        @test eltype(r[!, 3]) == Symbol
        @test eltype(r[!, 4]) == String
        @test eltype(r[!, 5]) == String
        N = names(r)
        @test N[1] == "a"
        @test N[2] == "b"
        @test N[3] == "c"
        @test N[4] == "d ()"
        @test N[5] == "e"
    end

    @testset "dataframe auto datetime" begin
        df = DataFrame()
        d = Date("2020-11-10")
        t = Time("10:12")
        df."d1 (:Date)" = [d]
        df."d2" = [d]
        df."t1 (:Time)" = [t]
        df."t2" = [t]
        r = Cropbox.unitfy(df)
        @test r[1, 1] == d
        @test r[1, 2] == d
        @test r[1, 3] == t
        @test r[1, 4] == t
        @test eltype(r[!, 1]) == Date
        @test eltype(r[!, 2]) == Date
        @test eltype(r[!, 3]) == Time
        @test eltype(r[!, 4]) == Time
        N = names(r)
        @test N[1] == "d1"
        @test N[2] == "d2"
        @test N[3] == "t1"
        @test N[4] == "t2"
    end
end


text:
using PyCall
using LinearAlgebra
using TensorOperations
using Printf
using Test
include("hf_utils.jl")

pyscf = pyimport("pyscf")

mol = pyscf.gto.M(atom = """O          0.00000        0.00000        0.11779
  H          0.00000        0.75545       -0.47116
  H          0.00000       -0.75545       -0.47116""", basis = "631g",)

mf = pyscf.scf.RHF(mol).run()

s = mol.intor_symmetric("int1e_ovlp")
h = mol.intor("int1e_kin") + mol.intor("int1e_nuc")
g = permutedims(mol.intor("int2e", aosym = 1), [4, 3, 2, 1])
e_nuc = mf.energy_nuc()
nelectron = mol.nelectron
n_ao = convert(UInt8, mol.nao_nr())
println(n_ao)

# println(mf.mo_energy)

# Check Fock solve
# pyscf_mo_energy, pyscf_mo_coeff = mf.eig(mf.get_fock(), mf.get_ovlp())
# println(pyscf_mo_coeff[:, 1]' * pyscf_mo_coeff[:, 1])
# println(dot(pyscf_mo_coeff[:, 1], pyscf_mo_coeff[:, 1]))
# println(pyscf_mo_coeff[1, :]' * pyscf_mo_coeff[10, :])
# println(dot(pyscf_mo_coeff[:, 1], pyscf_mo_coeff[:, 10]))

# exit(0)

#
# SCF Check
#

sym_rand = rand(n_ao, n_ao)
sym_rand = dm' + dm
rand_unitary = qr(sym_rand).Q

# Check rand_unitary
# println(diag(rand_unitary' * rand_unitary))

fake_mo_occ = zeros(n_ao, n_ao)
for i = 1:convert(UInt8, nelectron / 2)
  fake_mo_occ[i, i] = 2
end
starting_dm = rand_unitary' * fake_mo_occ * rand_unitary

function mf_scf(h, g, s, dm0, e_nuc)
  dm = copy(dm0)

  δE = 1
  δdm = 1
  last_energy = calc_energy_elec(h, g, dm0)[1]
  mo_occ = zeros(n_ao)
  for i = 1:convert(UInt8, nelectron / 2)
    mo_occ[i] = 2
  end

  @printf "Initial SCF Energy %16.9f\n" last_energy + e_nuc

  # SCF Cycle
  for i = 1:50

    # Solve Roothan-Hall Equations (Helgaker 10.6.17)
    mo_energy, mo_coeff = solve_fock(h, g, dm, s)
    @tensor mo_ovlp[i,j] := mo_coeff[n,i] * mo_coeff[m,j] * s[n,m]
    res = svd(mo_ovlp)
    println(res.S)
    println("MO Energy Error ", norm(mo_energy - mf.mo_energy))
    println("MO COEFF ERROR ", norm(mo_coeff - mf.mo_coeff))
    exit(0)

    # Convergence Check
    new_dm = make_rdm1(mo_occ, mo_coeff)
    energy_elec, e_1b, e_2b = calc_energy_elec(h, g, new_dm)

    # δs
    δE = energy_elec - last_energy
    δdm = norm(new_dm - dm)

    @printf "SCF Energy %16.9f \t δE = %10.3e \t δdm = %10.3e\n" energy_elec + e_nuc δE δdm

    if abs(δE) < 1e-6
      println("CONVERGED")
      break
    end

    # Update DM
    dm = copy(new_dm)
    last_energy = copy(energy_elec)

    exit(0)
  end
end

# mf_scf(h, g, s, starting_dm, e_nuc)
mf_scf(h, g, s, pyscf_rdm1, e_nuc)
println(mf.e_tot)
# function make_fock_ao(hμν, gμνρσ, Dρσ)

#   fμν = h
# end