text:
export Codecov

"""
Coverage.Codecov Module

This module provides functionality to push coverage information to the CodeCov.io
web service. It exports the `submit` and `submit_token` methods.
"""
module Codecov
    using HTTP
    using Coverage
    using JSON
    using Compat
    using Compat.LibGit2

    export submit, submit_token, submit_local, submit_generic

    #=
    JSON structure for Codecov.io
    https://codecov.io/api#post-report

    {
      "coverage": {
        "path/to/file.py": [null, 1, 0, null, true, 0, 0, 1, 1],
        "path/to/other.py": [null, 0, 1, 1, "1/3", null]
      },
      "messages": {
        "path/to/other.py": {
          "1": "custom message for line 1"
        }
      }
    }
    =#

    # Turn vector of filename : coverage pairs into a dictionary
    function to_json(fcs::Vector{FileCoverage})
        cov = Dict()
        for fc in fcs
            cov[fc.filename] = vcat(nothing, fc.coverage)
        end
        return Dict("coverage" => cov)
    end

    if VERSION >= v"0.7"
        """
        kwargs provides default values to insert into args_array, only if they are
        not already specified in args_array.
        """
        function set_defaults(args_array; kwargs...)
            defined_names = keys(pairs(args_array))
            is_args_array = Pair{Symbol, Any}[]
            if args_array isa Base.Iterators.Pairs
                is_args_array = vcat(is_args_array, collect(Pair(k, v) for (k,v) in args_array))
            else
                is_args_array = vcat(is_args_array, vec(args_array))
            end
            for kwarg in kwargs
                if !(kwarg[1] in defined_names)
                    push!(is_args_array, Pair(kwarg[1], kwarg[2]))
                end
            end
            return is_args_array
        end
    else
        """
        kwargs provides default values to insert into args_array, only if they are
        not already specified in args_array.
        """
        function set_defaults(args_array; kwargs...)
            defined_names = [k for (k, v) in args_array]
            for kwarg in kwargs
                if !(kwarg[1] in defined_names)
                    push!(args_array, kwarg)
                end
            end
            return args_array
        end
    end

    """
        submit(fcs::Vector{FileCoverage})

    Takes a vector of file coverage results (produced by `process_folder`),
    and submits them to Codecov.io. Assumes that this code is being run
    on TravisCI or AppVeyor. If running locally, use `submit_local`.
    """
    function submit(fcs::Vector{FileCoverage}; kwargs...)
        if lowercase(get(ENV, "APPVEYOR", "false")) == "true"
            appveyor_pr = get(ENV, "APPVEYOR_PULL_REQUEST_NUMBER", "")
            appveyor_job = join(
                [
                    ENV["APPVEYOR_ACCOUNT_NAME"],
                    ENV["APPVEYOR_PROJECT_SLUG"],
                    ENV["APPVEYOR_BUILD_VERSION"],
                ],
                "%2F",
            )
            kwargs = set_defaults(kwargs,
                service      = "appveyor",
                branch       = ENV["APPVEYOR_REPO_BRANCH"],
                commit       = ENV["APPVEYOR_REPO_COMMIT"],
                pull_request = appveyor_pr,
                job          = appveyor_job,
                slug         = ENV["APPVEYOR_REPO_NAME"],
                build        = ENV["APPVEYOR_JOB_ID"],
            )
        elseif lowercase(get(ENV, "TRAVIS", "false")) == "true"
            kwargs = set_defaults(kwargs,
                service      = "travis-org",
                branch       = ENV["TRAVIS_BRANCH"],
                commit       = ENV["TRAVIS_COMMIT"],
                pull_request = ENV["TRAVIS_PULL_REQUEST"],
                job          = ENV["TRAVIS_JOB_ID"],
                slug         = ENV["TRAVIS_REPO_SLUG"],
                build        = ENV["TRAVIS_JOB_NUMBER"],
            )
        elseif lowercase(get(ENV, "CIRCLECI", "false")) == "true"
            circle_slug = join(
                [
                    ENV["CIRCLE_PROJECT_USERNAME"],
                    ENV["CIRCLE_PROJECT_REPONAME"],
                ],
                "%2F",
            )
            kwargs = set_defaults(kwargs,
                service      = "circleci",
                branch       = ENV["CIRCLE_BRANCH"],
                commit       = ENV["CIRCLE_SHA1"],
                pull_request = get(ENV, "CIRCLE_PR_NUMBER", "false"),  # like Travis
                build_url    = ENV["CIRCLE_BUILD_URL"],
                slug         = circle_slug,
                build        = ENV["CIRCLE_BUILD_NUM"],
            )
        elseif lowercase(get(ENV, "JENKINS", "false")) == "true"
            kwargs = set_defaults(kwargs,
                service      = "jenkins",
                branch       = ENV["GIT_BRANCH"],
                commit       = ENV["GIT_COMMIT"],
                job          = ENV["JOB_NAME"],
                build        = ENV["BUILD_ID"],
                build_url    = ENV["BUILD_URL"],
                jenkins_url  = ENV["JENKINS_URL"],
            )
        else
            error("No compatible CI platform detected")
        end

        submit_generic(fcs; kwargs...)
    end

    """
        submit_local(fcs::Vector{FileCoverage}, dir::AbstractString=pwd())

    Take a `Vector` of file coverage results (produced by `process_folder`),
    and submit them to Codecov.io. Assumes the submission is being made from
    a local git installation, rooted at `dir`. A repository token should be specified by a
    `token` keyword argument or the `CODECOV_TOKEN` environment variable.
    """
    function submit_local(fcs::Vector{FileCoverage}, dir::AbstractString=joinpath(pwd(),".."); kwargs...)
        LibGit2.with(LibGit2.GitRepo(dir)) do repo
            LibGit2.with(LibGit2.head(repo)) do headref
                branch_name = LibGit2.shortname(headref) # this function returns a String
                commit_oid  = LibGit2.GitHash(LibGit2.peel(headref))
                kwargs = set_defaults(kwargs,
                    commit = string(commit_oid),
                    branch = branch_name
                    )
            end
        end

        if haskey(ENV, "REPO_TOKEN")
            println("the environment variable REPO_TOKEN is deprecated, use CODECOV_TOKEN instead")
            kwargs = set_defaults(kwargs, token = ENV["REPO_TOKEN"])
        end

        submit_generic(fcs; kwargs...)
    end

    @deprecate submit_token submit_local

    """
        submit_generic(fcs::Vector{FileCoverage})

    Takes a vector of file coverage results (produced by `process_folder`),
    and submits them to a Codecov.io instance. Keyword arguments are converted
    into a generic Codecov.io API uri.  It is essential that the keywords and
    values match the Codecov upload/v2 API specification.
    The `codecov_url` keyword argument or the CODECOV_URL environment variable
    can be used to specify the base path of the uri.
    The `dry_run` keyword can be used to prevent the http request from
    being generated
    """
    function submit_generic(fcs::Vector{FileCoverage}; kwargs...)
        @assert length(kwargs) > 0

        if haskey(ENV, "CODECOV_URL")
            kwargs = set_defaults(kwargs, codecov_url = ENV["CODECOV_URL"])
        end

        if haskey(ENV, "CODECOV_TOKEN")
            kwargs = set_defaults(kwargs, token = ENV["CODECOV_TOKEN"])
        end

        codecov_url = "https://codecov.io"
        dry_run = false
        verbose = true
        for (k,v) in kwargs
            if k == :codecov_url
                codecov_url = v
            end
            if k == :dry_run
                dry_run = true
            end
            if k == :verbose
                verbose = v
            end
        end
        @assert codecov_url[end] != "/" "the codecov_url should not end with a /, given url $(codecov_url)"

        uri_str = "$(codecov_url)/upload/v2?"
        for (k,v) in kwargs
            if k != :codecov_url && k != :dry_run && k != :verbose
                uri_str = "$(uri_str)&$(k)=$(v)"
            end
        end

        if verbose
            println("Codecov.io API URL:")
            println(uri_str)
        end

        if !dry_run
            heads   = Dict("Content-Type" => "application/json")
            data    = to_json(fcs)
            req     = HTTP.post(uri_str; body = JSON.json(data), headers = heads)
            println("Result of submission:")
            println(String(req))
        end
    end

end  # module Codecov


text:
import Ipopt
import PowerModels
import SymbolicAD

function power_model(case::String)
    pm = PowerModels.instantiate_model(
        joinpath(dirname(@__DIR__), "data", case),
        PowerModels.ACPPowerModel,
        PowerModels.build_opf,
    )
    return pm.model
end

model = power_model("pglib_opf_case2853_sdet.m")
SymbolicAD.run_unit_benchmark(model; direct_model = true, atol = 1e-6)
SymbolicAD.run_solution_benchmark(model, Ipopt.Optimizer; atol = 1.0)

# model = power_model("pglib_opf_case5_pjm.m")
# SymbolicAD.run_unit_benchmark(model; direct_model = true, atol = 1e-11)

ext: jl


text:
# Perhaps this could even be a typealias rather than a subtype...
abstract type AbstractArr{T, N} <: Assoc{CartesianIndex{N}, T}; end
const AbstractVec{T} = AbstractArr{T, 1}
const AbstractMat{T} = AbstractArr{T, 2}

#=
@propagate_inbounds function getindex(a::AbstractArr, i...)
    _getindex(IndexStyle(a), a, cartesian(a, i1, i2, is...)]
end

cartesian(a, i::Integer...) = CartesianIndex(i)
cartesian(a, i...) = Base.to_indices(a, i)

getindex(a::AbstractArr, s::Base.Slice) = _getindex(a::)
=#


text:
# ---------------------------------------------------------------------------- #
#
#   mesh3D.jl
#
#   Type for 3D meshes
#   Inherits from the abstract "mesh" type
#
#   λυτέος
#   Fall 2017
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

"""
    Mesh3D

Mesh3D type:
Type for 3D meshes consisting of tetrahedrons. It holds the nodes and
connectivity information.
"""
struct Mesh3D <: Mesh

  dim::Int64              # Dimension of the problem

  porder::Int64           # Polynomial order

  n::Int64                # Number of nodes

  p::Matrix{Float64}      # Nodal locations
  ploc::Matrix{Float64}   # Local nodal locations
  tloc::Matrix{Int64}     # Local tets
  t::Matrix{Int64}        # Tet - node connectivity
  t2f::Matrix{Int64}      # Tet - face connectivity
  f::Matrix{Int64}        # Face - node/tet connectivity
  fb::Matrix{Int64}       # Boundary face info
  nodes::Array{Float64,3} # Nodes on which solution is evaluated
  trorder::Matrix{Int64}  # Different ordering indices for faces

  boundtags::Vector{String} # Boundary names

end

"""
    Mesh3D( name::String, porder_::Int64; N = 5::Int64, M = N, Q = N )

Constructor for one of the default meshes. Currently "cube" is implemented.
"""
function Mesh3D( name::String, porder_::Int64; N = 5::Int64, M = N, Q = N, scl = 1. )

  if name == "cube"
    (p_, t_, bel_, tags_) = makecube( N, M, Q )
  elseif name == "skewcube"
    (p_, t_, bel_, tags_) = makeskewedcube( N, M, Q, 45.0 )
  elseif name == "single"
    p_   = [ 0.0 0.0 0.0;
             1.0 0.0 0.0;
             0.0 1.0 0.0;
             0.0 0.0 1.0 ]
    t_   = [ 1 2 3 4 ]
    bel_ = [ 2 3 4 1;
             1 4 3 2;
             1 2 4 3;
             1 3 2 4 ]
    tags_ = ["diag", "yz", "xz", "xy"]
  elseif name == "four"
    (p_, t_, bel_, tags_) = makefour( N )
  elseif name[end-3:end] == ".su2"
    (p_, t_, bel_, tags_) = readSU2_3D( name )
  elseif name[end-4:end] == ".mesh" # FEFLOA
    (p_, t_, bel_, tags_ ) = readFEFLOA_3D( name )
  else
    error("Unknown mesh type")
  end

  p_ .*= scl

  if porder_ > 3
    error("P>3 not implemented for 3D")
  end

  (f_, t2f_, nodes_, ploc_, tloc_, trorder_, fb_) = genmesh3D( porder_, p_, t_, bel_ )

  n_ = size( p_, 1 )

  Mesh3D( 3, porder_, n_, p_, ploc_, tloc_, t_, t2f_, f_, fb_, nodes_, trorder_, tags_ )

end

"""
    Mesh3D( name::String, porder_::Int64; N = 5::Int64, M = N, Q = N )

Constructor for one of the default meshes. Currently "cube" is implemented.
"""
function Mesh3D( p_::Matrix{Float64}, t_::Matrix{Int64}, bel_::Matrix{Int64}, porder_::Int64 )

  if porder_ > 3
    error("P>3 not implemented for 3D")
  end

  f_, t2f_, nodes_, ploc_, tloc_, trorder_, fb_ = genmesh3D( porder_, p_, t_, bel_ )

  n_ = size( p_, 1 )

  nb = maximum( bel_[:,4] )
  tags = Vector{String}( nb )
  for bb in 1:nb
      tags[bb] = string("bound", bb)
  end

  Mesh3D( 3, porder_, n_, p_, ploc_, tloc_, t_, t2f_, f_, fb_, nodes_, trorder_, tags )

end

"""
    genmesh3D( porder::Int64, p::Matrix{Float64}, t::Matrix{Int64}, bel::Matrix{Int64} )

Generates all necessary connectivity information given nodes locations (`p`),
tetrahedron connectivity (`t`), and boundary element information (`bel`).
"""
function genmesh3D( porder::Int64, p::Matrix{Float64}, t::Matrix{Int64}, bel::Matrix{Int64} )

  f_, t2f_, trorder_, fb_ = genFaces3D( t, bel )
  nodes_, ploc_, tloc_ = genNodes3D( porder, p, t )

  return f_, t2f_, nodes_, ploc_, tloc_, trorder_, fb_

end

"""
    genFaces3D( t::Matrix{Int64}, bel::Matrix{Int64} )

Generates face connectivity given tetrahedron connectivity (`t`), and
boundary element information (`bel`).
"""
function genFaces3D( t::Matrix{Int64}, bel::Matrix{Int64} )

  nt = size(t,1)

  # 2 3 4 - 1 4 3  - 1 2 4 - 1 3 2 : First ind is the vertices of face which does not included current vertices
  faces = vcat( t[:,[2,3,4]], t[:,[1,4,3]], t[:,[1,2,4]], t[:,[1,3,2]] ) # This holds all faces (but multiple copies)
  tets  = 1:nt
  tets  = vcat( tets, tets, tets, tets )

  nf = size(faces,1)

  orderind = [1 2 3;
              1 3 2;
              2 1 3;
              2 3 1;
              3 1 2;
              3 2 1]
  orderflip = [1, -1, -1, 1, 1, -1]

  tflip     = fill( 1::Int64, nf )
  typeorder = fill( 1::Int64, nf )

  # sort in ascending order
  for ii in 1:nf
    indorder      = sortperm(  faces[ii,:] )
    ind           = find( all( indorder' .== orderind, 2 ) )
    tflip[ii]     = orderflip[ ind[1] ]
    typeorder[ii] = ind[1]
    faces[ii,:]   = faces[ii, orderind[ind[1],:] ]
  end

  # This index links to the first unique index in the array faces
  ix  = groupslices( faces, 1 )

  t2f     = fill(0::Int64, nt, 8)
  tforder = fill(0::Int64, nt, 3)

  # Find unique vector index
  # Need to find an index that links from faces to the unique faces (without gaps)
  jx = Array{Int64}( size(ix) )
  jx[1] = 0
  lind  = 0 # maximum unique index
  mind  = 0 # maximum index in faces that had unique index
  nf    = 0  # number of unique elements

  indUni = fill( false, size(ix,1) ) # index of unique elements

  for ii = 1:length(ix)
    temp    = ix[ii] - lind - 1
    temp2   = ix[ii] - mind - 1
    jx[ii]  = temp
    mind    = max(mind,ix[ii])

    if temp2 >= 0 # Found a unique index
      lind = lind + 1
      nf  += 1
      indUni[ii] = true
    end

    if temp2 < 0 # Need to link to the unique index
      jx[ii] = jx[ix[ii]]
    end

  end

  f = fill( 0::Int64, nf, 5 )

  f[:,1:3] = faces[indUni,:]

  # make t2f and complete f
  for jj = 1:length(tflip)

    ind  = tets[jj]
    indt = Int64(ceil(jj / nt))

    indf = ix[jj] - jx[jj]

    if tflip[jj] == 1
      f[indf,4] = ind
    else
      f[indf,5] = ind
    end

    #t2f[ ind, indt ] = tflip[jj] * indf # NOTE: JUST A MINUS SIGN IS NOT SUFFICIENT ANYMORE, NEED TO KNOW THE ORDERING EXACTLY (I.E. ONE OF SIX DIFFERENT ORDERINGS)
    t2f[ ind, indt ]     = indf
    t2f[ ind, indt + 4 ] = tflip[jj] * typeorder[jj] # Sign indicates whether outward normal or not (outward: pos)

  end

  ### Boundary Elements
  # include information from 'bel' into 'f': the index of the boundaries

  bel[:,1:3] = sort(bel[:,1:3],2) # sort such that we compare against f
  fb = fill( 0, size(bel,1), 2 )
  nfb = 1

  # Look in fourth column of array
  indb = find( f[:,4] .== 0 )
  for jj = 1:length(indb)
    indf = indb[jj]

    # Find corresponding element in boundary element vector
    indel = find( all( bel[:,1:3] .== f[indf,1:3]', 2) )

    # Ensure boundary element is oriented counter-clockwise
    #   flip orientation of face in t2f
    indt     = f[indf,5]
    indfint  = find( t2f[indt,1:4] .== indf )
    orderorg = orderind[ abs.(t2f[indt,4+indfint]), : ]
    # flip nodes to get outward normal
    temp        = orderorg[3]
    orderorg[3] = orderorg[2]
    orderorg[2] = temp
    indorder = find( all( orderorg .== orderind, 2 ) )
    t2f[indt,4+indfint] = indorder
    #   flip nodes in f
    temp = f[indf,2]
    f[indf,2] = f[indf,3]
    f[indf,3] = temp
    f[indf,4] = f[indf,5]

    # Mark correct boundary
    f[indf,5] = - bel[indel[1],4]
    fb[nfb,1] = indf
    fb[nfb,2] = bel[indel[1],4]
    nfb += 1
  end

  # Look in fourth array
  indb = find( f[:,5] .== 0 )
  for jj = 1:length(indb)
    indf = indb[jj]

    # Find corresponding element in boundary element vector
    indel = find( all( bel[:,1:3] .== f[indf,1:3]', 2) )

    # Mark correct boundary
    f[indf,5] = - bel[indel[1],4]
    fb[nfb,1] =   indf
    fb[nfb,2] =   bel[indel[1],4]
    nfb += 1
  end

  return f, t2f, orderind, fb

end

"""
    genNodes3D( porder::Int64, p::Matrix{Float64}, t::Matrix{Int64} )

Generates the higher order nodes inside each element.
"""
function genNodes3D( porder::Int64, p::Matrix{Float64}, t::Matrix{Int64} )

  nt     = size(t,1)
  # (plocal, tlocal) = genlocal( porder )
  (plocal,tlocal) = genlocal3D( porder )

  npl    = size(plocal,1)
  nodes  = fill( 0.0, npl, 3, nt )

  # We use barycentric coordinates for this
  # for interpolation with barycentric coordinates at the point (xi,yi,zi) with
  # barycentric coordinates (li1, li2, li3, li3), the following holds:
  #   f(xi,yi) = li1 * f(x1,y1) + li2 * f(x2,y2) + li3 * f(x3,y3) + li4 * f(x4,y4)

  for ii in 1:nt, jj in 1:npl
    for kk in 1:4
      nodes[jj,1,ii] += plocal[jj,kk] * p[ t[ii,kk], 1 ]
      nodes[jj,2,ii] += plocal[jj,kk] * p[ t[ii,kk], 2 ]
      nodes[jj,3,ii] += plocal[jj,kk] * p[ t[ii,kk], 3 ]
    end
  end

  return nodes, plocal, tlocal

end

"""
    genlocal3D( porder::Int64 )

Generates the barycentric coordinates for the higher order nodes inside
the master element (tetrahedron).
"""
function genlocal3D( porder::Int64 )

  # The barycentric coordinates are implemented as (1 - x - y - z, x, y, z)

  n   = Int64( round( 1/6 * (porder+1) * (porder+2) * (porder+3) ) )
  n2d = Int64( round( 0.5 * (porder+1) * (porder+2) ) )

  plocal = Array{Float64}( n, 3 )
  qq = 1
  for kk = 1:porder+1
    for jj = 1:porder+2-kk
      for ii = 1:(porder + 3 - jj - kk)
        plocal[qq,:] = [ (ii-1)/porder, (jj-1)/porder, (kk-1)/porder ]
        qq += 1
      end
    end
  end

  temp = 1 - sum(plocal,2)

  plocal = hcat( temp, plocal )

  # Switch to have corners first, then edges
  plocal2 = fill( 0.0, size(plocal) )
  #COORD(3, s) = {0, 1, 0, 0, 0 , 0 , O3, T3, T3, O3, 0 , 0 , 0 , 0 , O3, T3, O3, 0 , O3, O3};
  #COORD(3, t) = {0, 0, 1, 0, T3, O3, 0 , 0 , O3, T3, T3, O3, 0 , 0 , 0 , 0 , O3, O3, 0 , O3};
  #COORD(3, u) = {0, 0, 0, 1, O3, T3, T3, O3, 0 , 0 , 0 , 0 , O3, T3, 0 , 0 , O3, O3, O3, 0 };
  indmap = [1,2,3,4]
  if porder == 2
    indmap = [1,porder+1,n2d,n,9,8,5,4,7,2]
  elseif porder == 3
    indmap = [1,porder+1,n2d,n,16,19,18,13,7,9,8,5,11,17,2,3,15,14,12,6]
  end
  plocal2 = plocal[ indmap, :]

  # generate local tetrahedra -- used for plotting
  if porder == 1
    tlocal      = fill( 0, 1, 4)
    tlocal[1,:] = [1 2 3 4]
  elseif porder == 2
    tlocal = [10     6     2     7
               4     9     5     6
               8     5     9    10
               3     5     8     7
               6     9     5    10
               8     5    10     7
               6    10     5     7
               8     9     1    10]
  elseif porder == 3
    tlocal = [2    16     8     9
             18     5    17    20
             20     5    17    10
             18    15    12    20
             18    19    15    20
             20    17     9    10
              3     5    11    10
             17     7    19     8
             18    17    19    20
             11     5    20    10
             16    17     8     9
              6     4    14     7
             11     5    18    20
             11    18    12    20
             20    17    16     9
             18     5     6    17
             15    20    19    16
             16    17    19     8
             20    17    19    16
             18    14    13    19
             19     6    14     7
             19    17     6     7
             18    17     6    19
             18     6    14    19
             18    13    15    19
             18    13    12    15
              1    12    13    15]
  end

  return (plocal2, tlocal)

end

"""
    makecube( n::Int64, m::Int64, q::Int64 )

Generates the nodes locations, tetrahedron connectivity and boundary element
information for a cube [0,1] x [0,1] x [0,1].
"""
function makecube( n::Int64, m::Int64, q::Int64 )

  # boundary 1 is z=0
  # boundary 2 is y=1
  # boundary 3 is z=1
  # boundary 4 is y=0
  # boundary 5 is x=0
  # boundary 6 is x=1

  p   = Array{Float64}( n*m*q, 3 )
  hex = Array{Int64}( (n-1)*(m-1)*(q-1), 8 ) # Map to ordering per hex cell
  t   = Array{Int64}( 6*(n-1)*(m-1)*(q-1), 4 )
  bel = Array{Int64}( 4*(n-1)*(q-1) + 4*(m-1)*(q-1) + 4*(m-1)*(n-1), 4 )

  # nodes
  for ii = 1:n, jj = 1:m, kk = 1:q
    p[ii + (jj-1) * n + (kk-1) * n * m, :] = [ 1/(n-1) * (ii - 1),
      1/(m-1) * (jj - 1), 1/(q-1) * (kk - 1) ]
  end
  # hex nodes
  kk     = 1
  hexind = 1
  for qq = 1:(q-1)
    for jj = 1:(m-1)
      for ii = 1:(n-1)

        hex[ kk, : ] = [hexind,     hexind+1,     hexind+n,     hexind+n+1,
                        hexind+n*m, hexind+n*m+1, hexind+n+n*m, hexind+n+1+n*m]

        hexind += 1
        kk += 1

      end
      hexind += 1
    end
    hexind += n
  end

  #         y
  #  3----------4
  #  |\     ^   |\
  #  | \    |   | \
  #  |  \   |   |  \
  #  |   7------+---8
  #  |   |  +-- |-- | -> x
  #  1---+---\--2   |
  #   \  |    \  \  |
  #    \ |     \  \ |
  #     \|      z  \|
  #      5----------6

  # tet index within hex
  tethexind = [1 2 3 5;
               2 5 4 3;
               7 3 4 5;
               2 5 6 4;
               5 7 6 4;
               8 7 4 6]

  # tetrahedra
  ntet_level = (n-1) * (m-1)
  for ii = 1:(n-1), jj = 1:(m-1), kk = 1:(q-1)

    hexind   = (ii-1) + (jj-1) * (n-1) + (kk-1) * (n-1) * (m-1)

    for pp in 1:6
      t[ 6*hexind + pp, :] = hex[ hexind+1, tethexind[pp,:] ]
    end

  end

  # triangle index for boundary
  indyzpos  = [2 4 6;
               4 8 6]
  indyzneg  = [1 5 3;
               5 7 3]
  indxzneg  = [6 5 2;
               5 1 2]
  indxzpos  = [8 4 7;
               4 3 7]
  indxypos  = [5 6 7;
               6 8 7]
  indxyneg  = [4 2 3;
               2 1 3]

  # boundary elements
  kk = 1
  for ii in 1:n-1, jj in 1:m-1 # xy neg
    indhex = (ii-1) + (jj-1)*(n-1)
    bel[kk,:] = [ hex[ indhex+1, indxyneg[1,:] ]'  1 ]
    kk = kk + 1
    bel[kk,:] = [ hex[ indhex+1, indxyneg[2,:] ]'  1 ]
    kk = kk + 1
  end
  for jj in 1:m-1, qq in 1:q-1 # yz pos
    indhex = (n-2) + (jj-1)*(n-1) + (qq-1)*(n-1)*(m-1)
    bel[kk,:] = [ hex[ indhex+1, indyzpos[1,:] ]'  2 ]
    kk = kk + 1
    bel[kk,:] = [ hex[ indhex+1, indyzpos[2,:] ]'  2 ]
    kk = kk + 1
  end
  for ii in 1:n-1, jj in 1:m-1 # xy pos
    indhex = (ii-1) + (jj-1)*(n-1) + (q-2)*(n-1)*(m-1)
    bel[kk,:] = [ hex[ indhex+1, indxypos[1,:] ]'  3 ]
    kk = kk + 1
    bel[kk,:] = [ hex[ indhex+1, indxypos[2,:] ]'  3 ]
    kk = kk + 1
  end
  for jj in 1:m-1, qq in 1:q-1 # yz neg
    indhex = (jj-1)*(n-1) + (qq-1)*(n-1)*(m-1)
    bel[kk,:] = [ hex[ indhex+1, indyzneg[1,:] ]'  4 ]
    kk = kk + 1
    bel[kk,:] = [ hex[ indhex+1, indyzneg[2,:] ]'  4 ]
    kk = kk + 1
  end
  for ii in 1:n-1, qq in 1:q-1 # xz neg
    indhex = (ii-1) + (qq-1)*(n-1)*(m-1)
    bel[kk,:] = [ hex[ indhex+1, indxzneg[1,:] ]'  5 ]
    kk = kk + 1
    bel[kk,:] = [ hex[ indhex+1, indxzneg[2,:] ]'  5 ]
    kk = kk + 1
  end
  for ii in 1:n-1, qq in 1:q-1 # xz pos
    indhex = (ii-1) + (m-2)*(n-1) + (qq-1)*(n-1)*(m-1)
    bel[kk,:] = [ hex[ indhex+1, indxzpos[1,:] ]'  6 ]
    kk = kk + 1
    bel[kk,:] = [ hex[ indhex+1, indxzpos[2,:] ]'  6 ]
    kk = kk + 1
  end

  # Add boundary tags
  tags = Vector{String}( 6 )
  tags[1] = "xyneg"
  tags[2] = "yzpos"
  tags[3] = "xypos"
  tags[4] = "yzneg"
  tags[5] = "xzneg"
  tags[6] = "xzpos"

  return p, t, bel, tags

end

"""
    makesymcube( n::Int64, m::Int64, q::Int64 )

Generates the nodes locations, tetrahedron connectivity and boundary element
information for a cube [0,1] x [0,1] x [0,1].
"""
function makesymcube( n::Int64, m::Int64, q::Int64 )

  # boundary 1 is z=0
  # boundary 2 is y=1
  # boundary 3 is z=1
  # boundary 4 is y=0
  # boundary 5 is x=0
  # boundary 6 is x=1

  p   = Array{Float64}( n*m*q, 3 )
  hex = Array{Int64}( (n-1)*(m-1)*(q-1), 8 ) # Map to ordering per hex cell
  t   = Array{Int64}( 6*(n-1)*(m-1)*(q-1), 4 )
  bel = Array{Int64}( 4*(n-1)*(q-1) + 4*(m-1)*(q-1) + 4*(m-1)*(n-1), 4 )

  # nodes
  for ii = 1:n, jj = 1:m, kk = 1:q
    p[ii + (jj-1) * n + (kk-1) * n * m, :] = [ 1/(n-1) * (ii - 1),
      1/(m-1) * (jj - 1), 1/(q-1) * (kk - 1) ]
  end
  # hex nodes
  kk     = 1
  hexind = 1
  for qq = 1:(q-1)
    for jj = 1:(m-1)
      for ii = 1:(n-1)

        hex[ kk, : ] = [hexind,     hexind+1,     hexind+n,     hexind+n+1,
                        hexind+n*m, hexind+n*m+1, hexind+n+n*m, hexind+n+1+n*m]

        hexind += 1
        kk += 1

      end
      hexind += 1
    end
    hexind += n
  end

  #         y
  #  3----------4
  #  |\     ^   |\
  #  | \    |   | \
  #  |  \   |   |  \
  #  |   7------+---8
  #  |   |  +-- |-- | -> x
  #  1---+---\--2   |
  #   \  |    \  \  |
  #    \ |     \  \ |
  #     \|      z  \|
  #      5----------6

  # tet index within hex
  tethexindorg = [1 2 3 5;
                  2 5 4 3;
                  7 3 4 5;
                  2 5 6 4;
                  5 7 6 4;
                  8 7 4 6]

  negxmap = [ [5,5],
              [6,6],
              [1,1],
              [3,3],
              [8,15],
              [7,16],
              [2,13],
              [4,14] ]

  negzmap = [ [1,1],
              [2,2],
              [3,3],
              [4,4],
              [8,10],
              [7,12],
              [6,11],
              [5,9] ]

  negyzmap = [ [1,1],
               [3,3],
               [2,13],
               [4,14],
               [5,9],
               [6,11],
               [8,17],
               [7,18] ]

  tethexind = fill( 0, 4*6 )
  tethexind[1:6,:] = tethexindorg

  texhextemp = copy(tethexindorg)
  for jj in 1:8
    tethextemp[ tethexindorg .== negxmap[jj][1] ] = negxmap[jj][2]
  end
  tethexind[7:12,:] = tethexindorg

  texhextemp = copy(tethexindorg)
  for jj in 1:8
    tethextemp[ tethexindorg .== negzmap[jj][1] ] = negxmap[jj][2]
  end
  tethexind[13:18,:] = tethexindorg

  texhextemp = copy(tethexindorg)
  for jj in 1:8
    tethextemp[ tethexindorg .== negyzmap[jj][1] ] = negxmap[jj][2]
  end
  tethexind[19:24,:] = tethexindorg

  # tetrahedra
  ntet_level = (n-1) * (m-1)
  for ii = 1:(n-1), jj = 1:(m-1), kk = 1:(q-1)

    hexind   = (ii-1) + (jj-1) * (n-1) + (kk-1) * (n-1) * (m-1)

    for pp in 1:6
      t[ 6*hexind + pp, :] = hex[ hexind+1, tethexind[pp,:] ]
    end

  end

  # triangle index for boundary
  indyzpos  = [2 4 6;
               4 8 6]
  indyzneg  = [1 5 3;
               5 7 3]
  indxzneg  = [6 5 2;
               5 1 2]
  indxzpos  = [8 4 7;
               4 3 7]
  indxypos  = [5 6 7;
               6 8 7]
  indxyneg  = [4 2 3;
               2 1 3]

  # boundary elements
  kk = 1
  for ii in 1:n-1, jj in 1:m-1 # xy neg
    indhex = (ii-1) + (jj-1)*(n-1)
    bel[kk,:] = [ hex[ indhex+1, indxyneg[1,:] ]'  1 ]
    kk = kk + 1
    bel[kk,:] = [ hex[ indhex+1, indxyneg[2,:] ]'  1 ]
    kk = kk + 1
  end
  for jj in 1:m-1, qq in 1:q-1 # yz pos
    indhex = (n-2) + (jj-1)*(n-1) + (qq-1)*(n-1)*(m-1)
    bel[kk,:] = [ hex[ indhex+1, indyzpos[1,:] ]'  2 ]
    kk = kk + 1
    bel[kk,:] = [ hex[ indhex+1, indyzpos[2,:] ]'  2 ]
    kk = kk + 1
  end
  for ii in 1:n-1, jj in 1:m-1 # xy pos
    indhex = (ii-1) + (jj-1)*(n-1) + (q-2)*(n-1)*(m-1)
    bel[kk,:] = [ hex[ indhex+1, indxypos[1,:] ]'  3 ]
    kk = kk + 1
    bel[kk,:] = [ hex[ indhex+1, indxypos[2,:] ]'  3 ]
    kk = kk + 1
  end
  for jj in 1:m-1, qq in 1:q-1 # yz neg
    indhex = (jj-1)*(n-1) + (qq-1)*(n-1)*(m-1)
    bel[kk,:] = [ hex[ indhex+1, indyzneg[1,:] ]'  4 ]
    kk = kk + 1
    bel[kk,:] = [ hex[ indhex+1, indyzneg[2,:] ]'  4 ]
    kk = kk + 1
  end
  for ii in 1:n-1, qq in 1:q-1 # xz neg
    indhex = (ii-1) + (qq-1)*(n-1)*(m-1)
    bel[kk,:] = [ hex[ indhex+1, indxzneg[1,:] ]'  5 ]
    kk = kk + 1
    bel[kk,:] = [ hex[ indhex+1, indxzneg[2,:] ]'  5 ]
    kk = kk + 1
  end
  for ii in 1:n-1, qq in 1:q-1 # xz pos
    indhex = (ii-1) + (m-2)*(n-1) + (qq-1)*(n-1)*(m-1)
    bel[kk,:] = [ hex[ indhex+1, indxzpos[1,:] ]'  6 ]
    kk = kk + 1
    bel[kk,:] = [ hex[ indhex+1, indxzpos[2,:] ]'  6 ]
    kk = kk + 1
  end

  # Add boundary tags
  tags = Vector{String}( 6 )
  tags[1] = "xyneg"
  tags[2] = "yzpos"
  tags[3] = "xypos"
  tags[4] = "yzneg"
  tags[5] = "xzneg"
  tags[6] = "xzpos"

  return p, t, bel, tags

end

"""
    makeskewedcube( n::Int64, m::Int64, q::Int64, th::Float64 )

Generates the nodes locations, tetrahedron connectivity and boundary element
information for a skewed cube [0,1] x [0,1] x [0,1].
"""
function makeskewedcube( n::Int64, m::Int64, q::Int64, th::Float64 )

  # boundary 1 is z=0
  # boundary 2 is x=1
  # boundary 3 is z=1
  # boundary 4 is x=0
  # boundary 5 is y=0
  # boundary 6 is y=1

  p   = Array{Float64}( n*m*q, 3 )
  hex = Array{Int64}( (n-1)*(m-1)*(q-1), 8 ) # Map to ordering per hex cell
  t   = Array{Int64}( 6*(n-1)*(m-1)*(q-1), 4 )
  bel = Array{Int64}( 4*(n-1)*(q-1) + 4*(m-1)*(q-1) + 4*(m-1)*(n-1), 4 )

  # nodes
  for ii = 1:n, jj = 1:m, kk = 1:q
    p[ii + (jj-1) * n + (kk-1) * n * m, :] = [ 1/(n-1) * (ii - 1),
      1/(m-1) * (jj - 1), 1/(q-1) * (kk - 1) ]
  end

  # skew cube
  p[:,1] = p[:,1] + tand(th) * p[:,2]
  p[:,2] = p[:,2] + tand(th) * p[:,3]

  # hex nodes
  kk     = 1
  hexind = 1
  for qq = 1:(q-1)
    for jj = 1:(m-1)
      for ii = 1:(n-1)

        hex[ kk, : ] = [hexind,     hexind+1,     hexind+n,     hexind+n+1,
                        hexind+n*m, hexind+n*m+1, hexind+n+n*m, hexind+n+1+n*m]

        hexind += 1
        kk += 1

      end
      hexind += 1
    end
    hexind += n
  end

  #         y
  #  3----------4
  #  |\     ^   |\
  #  | \    |   | \
  #  |  \   |   |  \
  #  |   7------+---8
  #  |   |  +-- |-- | -> x
  #  1---+---\--2   |
  #   \  |    \  \  |
  #    \ |     \  \ |
  #     \|      z  \|
  #      5----------6

  # tet index within hex
  tethexind = [1 2 3 5;
               2 5 4 3;
               7 3 4 5;
               2 5 6 4;
               5 7 6 4;
               8 7 4 6]

  # tetrahedra
  ntet_level = (n-1) * (m-1)
  for ii = 1:(n-1), jj = 1:(m-1), kk = 1:(q-1)

    hexind   = (ii-1) + (jj-1) * (n-1) + (kk-1) * (n-1) * (m-1)

    for pp in 1:6
      t[ 6*hexind + pp, :] = hex[ hexind+1, tethexind[pp,:] ]
    end

  end

  # triangle index for boundary
  indyzpos  = [2 4 6;
               4 8 6]
  indyzneg  = [1 5 3;
               5 7 3]
  indxzneg  = [6 5 2;
               5 1 2]
  indxzpos  = [8 4 7;
               4 3 7]
  indxypos  = [5 6 7;
               6 8 7]
  indxyneg  = [4 2 3;
               2 1 3]

  # boundary elements
  kk = 1
  for ii in 1:n-1, jj in 1:m-1 # xy neg
    indhex = (ii-1) + (jj-1)*(n-1)
    bel[kk,:] = [ hex[ indhex+1, indxyneg[1,:] ]'  1 ]
    kk = kk + 1
    bel[kk,:] = [ hex[ indhex+1, indxyneg[2,:] ]'  1 ]
    kk = kk + 1
  end
  for jj in 1:m-1, qq in 1:q-1 # yz pos
    indhex = (n-2) + (jj-1)*(n-1) + (qq-1)*(n-1)*(m-1)
    bel[kk,:] = [ hex[ indhex+1, indyzpos[1,:] ]'  2 ]
    kk = kk + 1
    bel[kk,:] = [ hex[ indhex+1, indyzpos[2,:] ]'  2 ]
    kk = kk + 1
  end
  for ii in 1:n-1, jj in 1:m-1 # xy pos
    indhex = (ii-1) + (jj-1)*(n-1) + (q-2)*(n-1)*(m-1)
    bel[kk,:] = [ hex[ indhex+1, indxypos[1,:] ]'  3 ]
    kk = kk + 1
    bel[kk,:] = [ hex[ indhex+1, indxypos[2,:] ]'  3 ]
    kk = kk + 1
  end
  for jj in 1:m-1, qq in 1:q-1 # yz neg
    indhex = (jj-1)*(n-1) + (qq-1)*(n-1)*(m-1)
    bel[kk,:] = [ hex[ indhex+1, indyzneg[1,:] ]'  4 ]
    kk = kk + 1
    bel[kk,:] = [ hex[ indhex+1, indyzneg[2,:] ]'  4 ]
    kk = kk + 1
  end
  for ii in 1:n-1, qq in 1:q-1 # xz neg
    indhex = (ii-1) + (qq-1)*(n-1)*(m-1)
    bel[kk,:] = [ hex[ indhex+1, indxzneg[1,:] ]'  5 ]
    kk = kk + 1
    bel[kk,:] = [ hex[ indhex+1, indxzneg[2,:] ]'  5 ]
    kk = kk + 1
  end
  for ii in 1:n-1, qq in 1:q-1 # xz pos
    indhex = (ii-1) + (m-2)*(n-1) + (qq-1)*(n-1)*(m-1)
    bel[kk,:] = [ hex[ indhex+1, indxzpos[1,:] ]'  6 ]
    kk = kk + 1
    bel[kk,:] = [ hex[ indhex+1, indxzpos[2,:] ]'  6 ]
    kk = kk + 1
  end

  # Add boundary tags
  tags = Vector{String}( 6 )
  tags[1] = "xyneg"
  tags[2] = "yzpos"
  tags[3] = "xypos"
  tags[4] = "yzneg"
  tags[5] = "xzneg"
  tags[6] = "xzpos"

  return p, t, bel, tags

end

function makefour( n::Int64 )

    if n == 2
        p = [ 0.0 0.0         0.0;
              1.0 0.0         0.0;
              0.5 0.5*sqrt(3) 0.0;
              0.5 0.5/sqrt(3) 0.5*sqrt(3);
              0.5 0.5/sqrt(3) 0.25 ]
        t = [ 5     1     4     3
              1     5     2     3
              2     5     4     3
              4     1     5     2 ]
        bel = [1 3 2 1;
               1 2 4 2;
               3 4 2 3;
               1 4 3 4]
    else
        error("n is not 2")
    end

    # Add boundary tags
    tags = Vector{String}( 4 )
    tags[1] = "bound1"
    tags[2] = "bound2"
    tags[3] = "bound3"
    tags[4] = "bound4"

    return p, t, bel, tags

end


text:
module ImageUtils

using AxisArrays
using ColorSchemes
using ColorVectorSpace
using FFTW
using ImageMagick
using Interpolations
using LinearAlgebra
using Random
using Reexport
using Roots
using Statistics
using TestImages
using Unitful

@reexport using Colors
@reexport using FixedPointNumbers
@reexport using Images

const axes = Base.axes

export shepp_logan

include("Arrays.jl")
include("Coloring.jl")
include("Slicing.jl")
include("Transformation.jl")
include("PerfusionMaps.jl")
include("FileExport.jl")
include("Permutation.jl")
include("AnalyzeNifti.jl")

shepp_logan(args...) = Float64.(TestImages.shepp_logan(args...))

end # module


text:
@auto_hash_equals struct Point{I,C,F}
    ec::C   # the curve containing this point
    x::F
    y::F
end

const ConcretePoint{C,F} = Point{false,C,F}
const IdealPoint{C,F}    = Point{true, C,F}

point(ec::Curve{F}, x, y) where {F} = point(ec, convert(F, x), convert(F, y))
point(ec::Curve{F}, x::F, y::F, valid = false) where {F} = begin
    valid || contains(ec, x, y) || error("($(x), $(y)) is not a valid point on $(ec)")
    ConcretePoint{typeof(ec),F}(ec, x, y)
end

point(ec::Curve{F}) where {F} = (z = zero(ring(ec)); IdealPoint{typeof(ec),F}(ec, z, z))
ideal(ec::Curve) = point(ec)

# http://sagenb.org/src/schemes/elliptic_curves/ell_generic.py: is_x_coord(self, x)
# http://sagenb.org/src/schemes/elliptic_curves/ell_generic.py: lift_x(self, x, all=False)
point(ec::WNFCurve{T}, x::T) where {T} = begin
    y = x
    try
        r = sqrt(x^3 + ec.a * x + ec.b)
        y = convert(T, r)
    catch
    end

    !contains(ec, x, y) && error("no point with x-coordinate $(x) on $(ec)")
    point(ec, x, y, true)
end

rand(ec::Curve{T}) where {T<:FinFieldElem} = begin
    R = field(ec)
    q = order(R)

    # the following allows the ideal point to be picked
    rand() <= 1 / convert(BigInt, q + 1) && return ideal(ec)

    for _ in 1:100
        x = rand(R)
        try
            p = point(ec, x)
            return p
        catch e
            !isa(e, ErrorException) && rethrow()
        end
    end

    error("no random point on $(ec) could be found")
end

show(io::IO, p::ConcretePoint) = print(io, "($(p.x), $(p.y)) on elliptic curve $(curve(p))")
show(io::IO, p::IdealPoint) = print(io, "𝒪 on elliptic curve $(curve(p))")

curve(p::Point) = p.ec
in(p::Point, ec::Curve) = (ec == curve(p)) && contains(ec, p.x, p.y)

samecurve(p::Point{I,C,F}, q::Point{J,C,F}) where {I,J,C,F} = curve(p) == curve(q)
samecurve(p::Point, q::Point) = false

isideal(::IdealPoint) = true
isideal(::ConcretePoint) = false

-(p::IdealPoint) = p
-(p::ConcretePoint{C}) where {C<:WNFCurve} = point(curve(p), p.x, -p.y, true)
-(p::ConcretePoint{C}) where {C<:GWNFCurve} = point(curve(p), p.x, -p.y - curve(p).a₁ * p.x - curve(p).a₃, true)
-(p::Point, q::Point) = p + (-q)

+(p::Point, q::IdealPoint) = (samecurve(p, q) || error("$(p) and $(q) have different curves"); p)
+(p::IdealPoint, q::ConcretePoint) = (samecurve(p, q) || error("$(p) and $(q) have different curves"); q)
+(p::ConcretePoint{C}, q::ConcretePoint{C}) where {C<:WNFCurve} = begin
    samecurve(p, q) || error("$(p) and $(q) have different curves")

    ec = curve(p)
    if p == q
        p.y == 0 && return ideal(ec)        # vertical line
        m = divexact(3*p.x^2 + ec.a, 2*p.y) # slope of the tangent line
    else
        p.x == q.x && return ideal(ec)      # vertical line
        m = divexact(q.y - p.y, q.x - p.x)  # slope of the secant line
    end

    x = m^2 - q.x - p.x                     # using Vieta's formula for the sum of the roots
    y = m * (x - p.x) + p.y
    point(ec, x, -y, true)                  # do the reflection to get the sum of the two points
end
+(p::ConcretePoint{C}, q::ConcretePoint{C}) where {C<:GWNFCurve} = begin
    samecurve(p, q) || error("$(p) and $(q) have different curves")

    ec = curve(p)
    if p.x == q.x
        p.y + q.y + ec.a₁*p.x + ec.a₃ == 0 && return ideal(ec)

        t = 2*p.y + ec.a₁*p.x + ec.a₃
        c = divexact(3*p.x^2  + 2*ec.a₂*p.x +   ec.a₄ - ec.a₁*p.y, t)
        d = divexact(-(p.x^3) +   ec.a₄*p.x + 2*ec.a₆ - ec.a₃*p.y, t)
     else
        c = divexact(q.y - p.y,         q.x - p.x)
        d = divexact(p.y*q.x - q.y*p.x, q.x - p.x)
    end

    Σx = c^2 + ec.a₁*c - ec.a₂ - p.x - q.x
    Σy = -(c + ec.a₁) * Σx - d - ec.a₃
    point(ec, Σx, Σy, true)
end

*(n::Integer, p::IdealPoint) = p
*(p::IdealPoint, n::Integer) = p
*(n::Integer, p::ConcretePoint) = p * n
*(p::ConcretePoint, n::Integer) = begin
    n == 0 && return ideal(curve(p))
    n < 0 && return -p * -n

    q = p
    r = (n & 1 == 1) ? p : ideal(curve(p))

    i = oftype(n, 2)
    while i <= n
      q = q + q
      if n & i == i
        r = q + r
      end
      i = i << 1
    end

    r
end


text:
#using HDF5
#using LegendHDF5IO: readdata, writedata
using JLD2
using SolidStateDetectors
using Unitful

using Plots;
pyplot(fmt = :png)

T = Float32

function SaveSimulation(filename::String, simulation::Simulation{T})::Nothing
    @info "Saving simulation"
    if !isfile(filename) mkpath(dirname(filename)) end
    jldopen(filename, "w") do f
        write(f, "Simulation", NamedTuple(simulation))
    end
end

function ReadSimulation(filename)::Union{Nothing, Simulation{T}}
    if isfile(filename)
        @info "Reading simulation"
        simulation = jldopen(filename, "r") do f
            ConvertNamedTuple(read(f, "Simulation"))
        end
        # SetChargeDriftModel!(simulation)

        return simulation
    end
    return nothing
end

function ConvertNamedTuple(nt::NamedTuple)::Union{Nothing, Simulation{T}}
    if !ismissing( nt.detector_json_string )
        det = SolidStateDetector{T}( Dict(nt.detector_json_string) )
        sim = Simulation(det)

        if hasproperty(nt, :weighting_potentials)
            for contact in sim.detector.contacts
                if !ismissing( values( nt.weighting_potentials[contact.id] )[1] )
                    sim.weighting_potentials[contact.id] = WeightingPotential( nt.weighting_potentials[contact.id] )
                end
            end
        end

        d = Dict(:ρ => EffectiveChargeDensity,
                 :ρ_fix => EffectiveChargeDensity,
                 :ϵ => DielectricDistribution,
                 :point_types => PointTypes,
                 :electric_potential => ElectricPotential,
                 :electric_field => ElectricField,
                 :electron_drift_field => ElectricField,
                 :hole_drift_field => ElectricField)

        for (prop, func) in d
            if !ismissing( nt[prop][1] ) setproperty!( sim, prop, func( nt[prop] ) ) end
        end
        return sim
    end
    return nothing
end

function CalculateDefaultDetectorFields!(sim::Simulation{T},
    driftConfigFile::Union{Missing, String} = missing,
    max_ref = 4)

    CalculateElectricPotential!(sim, max_ref)
    CalculateElectricField!(sim)
    CalculateWeightingPotential!(sim, max_ref)
    SetChargeDriftModel!(sim, driftConfigFile)
end

function SetupDefaultSimulation(geomConfigFile::String, CCDName::String="")::Simulation{T}
    simulation = Simulation{T}(geomConfigFile)

    if CCDName != ""
        # include(CCDName) #put this line in the main script you are running
        simulation.detector.semiconductors[1].charge_density_model = ccdm
    end

    apply_initial_state!(simulation, ElectricPotential)
    for c in simulation.detector.contacts
        apply_initial_state!(simulation, WeightingPotential, c.id)
    end
    return simulation
end

function SetBiasVoltage!(sim, contact, voltage)
    sim.detector.contacts[contact].potential = voltage
end

function CalculateElectricPotential!(simulation::Simulation{T}, max_ref = 4, verbose = false)
    calculate_electric_potential!(simulation, depletion_handling = true, max_refinements = max_ref)
    if verbose
        println("Active volume: ", get_active_volume(simulation.point_types))
    end
end

function CalculateWeightingPotential!(simulation::Simulation{T}, max_ref = 2)
    # TODO Add options
    for contact in simulation.detector.contacts
        calculate_weighting_potential!(simulation, contact.id, max_refinements = max_ref, n_points_in_φ = 2, verbose = false)
    end
end

function CalculateElectricField!(simulation::Simulation{T}, points_in_φ = 72, verbose = true)
    calculate_electric_field!(simulation, n_points_in_φ = points_in_φ)
    if verbose
        println("Capacitance: ", calculate_capacitance(simulation))
    end
end

function SetChargeDriftModel!(simulation::Simulation{T}, driftConfigFile::Union{Missing, String})
    charge_drift_model = ADLChargeDriftModel(driftConfigFile, T=T)
    set_charge_drift_model!(simulation, charge_drift_model)
    calculate_drift_fields!(simulation)
end

#-------------------------------------------------------------------------------
# Plotting
#

function PlotDetector(simulation::Simulation{T})
    p = plot(simulation.detector, size = (700, 700), title = "Detector geometry")
    return p
end

function PlotElectricPotential(simulation::Simulation{T})
    p = plot(
    plot(simulation.electric_potential), # initial electric potential (boundary conditions)
    plot(simulation.point_types), # map of different point types: fixed point / inside or outside detector volume / depleted/undepleted
    plot(simulation.ρ), # charge density distribution
    plot(simulation.ϵ), # dielectric distribution
    layout = (1, 4), size = (1400, 700)
    )
    return p
end

function PlotWeightingPotential(simulation::Simulation{T}, contact::Int64)
    p = plot(
    plot(simulation.weighting_potentials[contact]),
    plot(simulation.point_types),
    plot(simulation.electric_potential),
    layout = (1, 3), size = (1000, 700)
    )
    return p
end

max_stars_repo_name: leenderthayen/NDP

max_issues_repo_name: leenderthayen/NDP

max_forks_repo_name: leenderthayen/NDP


text:
using KhepriBase
using Test

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


text:
#=
language_detection.jl; this script shows the application of hypervectors to a language_detection problem.
=#
using Pkg
Pkg.activate(".")
using CSV, DataFrames, HyperdimensionalComputing, MLJ

"""
This function makes all letters in a text lowercase, removes accents and all-round 'normalizes' language-specific variations on letters to a 26-letter alphabet.
"""
function language_preprocessing(text::String)
        #make text lowercase
        text = lowercase(text)

        #remove odd accents and peculiarities
        text = replace(text, "é"=>"e")
        text = replace(text, "ê"=>"e")
        text = replace(text, "è"=>"e")
        text = replace(text, "ë"=>"e")
        text = replace(text, "ę"=>"e")
        text = replace(text, "á"=>"a")
        text = replace(text, "à"=>"a")
        text = replace(text, "â"=>"a")
        text = replace(text, "ä"=>"a")
        text = replace(text, "ą"=>"a")
        text = replace(text, "ç"=>"c")
        text = replace(text, "ć"=>"c")
        text = replace(text, "î"=>"i")
        text = replace(text, "ï"=>"i")
        text = replace(text, "í"=>"i")
        text = replace(text, "ł"=>"l")
        text = replace(text, "ñ"=>"n")
        text = replace(text, "ń"=>"n")
        text = replace(text, "ô"=>"o")
        text = replace(text, "ö"=>"o")
        text = replace(text, "ó"=>"o")
        text = replace(text, "ś"=>"s")
        text = replace(text, "ß"=>"ss")
        text = replace(text, "û"=>"u")
        text = replace(text, "ü"=>"u")
        text = replace(text, "ú"=>"u")
        text = replace(text, "ź"=>"z")
        text = replace(text, "ż"=>"z")

        #remove punctuation and double spaces
        text = string([t for t in text if t ∈ vcat('a':'z', ' ')]...)
        text = replace(text, "  "=>" ")

        return text
end

# load in and preprocess training paragraphs from project Gutenberg
labelled_languages_df = DataFrame(CSV.File("data/labelled_languages_train.csv"))
insertcols!(labelled_languages_df, 2, :"PREPROC"=>Vector{String}(undef, size(labelled_languages_df, 1)))
labelled_languages_df[!, "PREPROC"] = [language_preprocessing(labelled_languages_df[i, "TEXT"]) for i in 1:size(labelled_languages_df, 1)]

# load in and preprocess testing sentences from tatoeba.org
examples_df = DataFrame(CSV.File("data/labelled_languages_test.csv"))
insertcols!(examples_df, 2, :"PREPROC"=>Vector{String}(undef, size(examples_df, 1)))
examples_df[!, "PREPROC"] = [language_preprocessing(examples_df[i, "TEXT"]) for i in 1:size(examples_df, 1)]

for vec in [BipolarHDV, BinaryHDV, GradedHDV, RealHDV]
    println(vec)

    # encode both train and test data
    alphabet = Dict(c=>vec() for c in vcat('a':'z', ' '))
    encoded_references = [sequence_embedding(ref, alphabet, 3) for ref in labelled_languages_df[!, "PREPROC"]]
    encoded_examples = [sequence_embedding(ex, alphabet, 3) for ex in examples_df[!, "PREPROC"]]

    # predict
    centers = train(labelled_languages_df[!, "LANGUAGE"], encoded_references)
    test_predictions = [HyperdimensionalComputing.predict(hv, centers) for hv in encoded_examples]

    # evaluate performance
    println(string("weighted F1-score: ", multiclass_f1score(test_predictions, examples_df[!, "LANGUAGE"])))
end


text:
const loadedModules = Dict{String, Module}()
const loadedModulesTimes = Dict{String, Float64}()

function loadModule(fn::String, force::Bool=false)
    # Explicit check for .jl extensions
    # Modules from zip will have .jl.from_zip to prevent loading via this method

    if !hasExt(fn, ".jl")
        return false
    end

    try
        fileMtime = stat(fn).mtime

        if fileMtime > get(loadedModulesTimes, fn, 0) || force
            loadedModules[fn] = include(fn)
            loadedModulesTimes[fn] = fileMtime

            return true

        else
            return false
        end

    catch e
        println(Base.stderr, "! Failed to load \"$fn\"")
        println(Base.stderr, e)

        return false
    end
end

hasModuleField(modul::Module, funcname::String) = isdefined(modul, Symbol(funcname))
hasModuleField(modul::String, funcname::String) = haskey(loadedModules, modul) && isdefined(loadedModules[modul], Symbol(funcname))

getModuleField(modul::Module, funcname::String) = getfield(modul, Symbol(funcname))
getModuleField(modul::String, funcname::String) = getfield(loadedModules[modul], Symbol(funcname))

function eventToModule(modul::Module, funcname::String, args...)
    if hasModuleField(modul, funcname)
        # Get function and call with arguments
        func = getModuleField(modul, funcname)

        if isa(func, Function) && applicable(func, args...)
            try
                res = func(args...)

                return res

            catch e
                # Error running the function
                println(Base.stderr, "Exception running function $funcname for $modul")

                for (exc, bt) in Base.catch_stack()
                    showerror(Base.stderr, exc, bt)
                    println(Base.stderr, "")
                end
            end

        else
            # Can't call the function with the arguments we have
            # This is normal, event is not consumed

            return false
        end

    else
        # Module doesn't have the method
        # Return false as this cannot consume the event

        return false
    end
end

eventToModule(v::Nothing, funcname::String, args...) = false

function eventToModule(fn::String, funcname::String, args...)
    # Make sure the module is properly loaded
    # If it isn't we can just return nothing
    if haskey(loadedModules, fn)
        return eventToModule(loadedModules[fn], funcname, args...)
    end
end

function eventToModules(funcname::String, args...)
    for (filename, modul) in loadedModules
        moduleRes = eventToModule(modul, funcname, args...)

        # False and nothing as only return value means we are not handing it and it should be taken to another function
        if moduleRes == false || moduleRes === nothing
            continue
        end

        return moduleRes
    end

    return false
end

function eventToModules(moduleNames::Array{String, 1}, funcname::String, args...)
    for moduleName in moduleNames
        moduleRes = eventToModule(moduleName, funcname, args...)

        # False and nothing as only return value means we are not handing it and it should be taken to another function
        if moduleRes == false || moduleRes === nothing
            continue
        end

        return moduleRes
    end

    return false
end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
include("../src/geom_module.jl")

module gtest
using BenchmarkTools
using ..geom
using PyPlot
PyPlot.pygui(true)

function test_for_concrete_points()

    x_ = [ 0.5715937461408043
    0.8916206497257309
    0.8203136478441682
    0.25164289137427875
    0.4403285873605831
    0.43681263870560505
    0.2217704001238916
    0.20004462640858622
    0.5542972315770416
    0.9659375830665178
    0.1
    1]

    y_ = [ 0.712953799846491
    0.8172755930467877
    0.8094867906788563
    0.5920712401620924
    0.7383511283992972
    0.4897454560982253
    0.34283329277824315
    0.39353333699051496
    0.7313215009997691
    0.19999430363560622
    0.3
    0.3]

    point_array = map(geom.Point, x_, y_)
    (p_x, p_y) = geom.discretize_for_plot(point_array)

    ch_graham   = geom.graham_scan!(point_array)
    ch_jarvis   = geom.jarvis_march!(point_array)
    ch_jarvis2  = geom.jarvis2!(point_array)
    ch_jarvis3  = geom.jarvis3!(point_array)
    ch_incr     = geom.incremental_hull!(point_array)

    (ch_graham_x, ch_graham_y)      = geom.discretize_for_plot(ch_graham, true)
    (ch_jarvis_x, ch_jarvis_y)      = geom.discretize_for_plot(ch_jarvis, true)
    (ch_jarvis2_x, ch_jarvis2_y)    = geom.discretize_for_plot(ch_jarvis2, true)
    (ch_jarvis3_x, ch_jarvis3_y)    = geom.discretize_for_plot(ch_jarvis3, true)
    (ch_incr_x, ch_incr_y)          = geom.discretize_for_plot(ch_incr, true)

    PyPlot.figure()
    PyPlot.grid()

    PyPlot.plot(p_x, p_y,
        color = "k",
        alpha = 0.7,
        linestyle = "none",
        marker = ".",
        markersize = 10)

    PyPlot.plot(ch_graham_x, ch_graham_y,
        linewidth = 2,
        marker = "o")

    PyPlot.figure()
    PyPlot.grid()

    PyPlot.plot(p_x, p_y,
        color = "k",
        alpha = 0.7,
        linestyle = "none",
        marker = ".",
        markersize = 10)

    PyPlot.plot(ch_jarvis_x, ch_jarvis_y,
        linewidth = 2,
        marker = "o")

    PyPlot.figure()
    PyPlot.grid()

    PyPlot.plot(p_x, p_y,
        color = "k",
        alpha = 0.7,
        linestyle = "none",
        marker = ".",
        markersize = 10)

    PyPlot.plot(ch_incr_x, ch_incr_y,
        linewidth = 2,
        marker = "o")
end

function test_random_test()

    for kk = 1:3# 100
        N_point = 10
        point_array = geom.random_point(N_point)

        point_array = map(p -> geom.Point(p.x + p.y, p.x - p.y), point_array)

        ch_graham = geom.graham_scan!(point_array)
        # ch = geom.jarvis_march!(point_array)
        # ch = geom.jarvis2!(point_array)
        # ch = geom.jarvis3!(point_array)
        ch = geom.incremental_hull!(point_array)
        reduced_point_array = geom.reduce_array_for_convex_hull!(point_array)

        # println("The graham length $(length(ch_graham))")
        # @show length(ch)
        # @show all( ch_graham .== ch )

        L_incr = length(ch)
        L_graham = length(ch_graham)
        if L_incr == L_graham
            b = all(ch_graham .== ch)
            if b == false
                # they are not identical because of the permutation of the points
                println("stop here")

                (p_x, p_y) = geom.discretize_for_plot(point_array)
                (ch_graham_x, ch_graham_y) = geom.discretize_for_plot(ch_graham, true)
                (ch_x, ch_y) = geom.discretize_for_plot(ch, true)

                PyPlot.figure()
                PyPlot.grid()
                PyPlot.plot(p_x, p_y,
                    color = "k",
                    alpha = 0.7,
                    linestyle = "none",
                    marker = ".",
                    markersize = 10)
                PyPlot.plot(ch_x, ch_y,
                    linewidth = 2,
                    marker = "o")
                PyPlot.plot(ch_graham_x, ch_graham_y,
                    linewidth = 2,
                    marker = "*")

            end
            println( b )

        else
            if L_incr < L_graham
                println("graham is longer  Lincr $(L_incr), Lgr $(L_graham)")
            else
                println("incr is longer Lincr $(L_incr), Lgr $(L_graham)")
            end
        end
    end

end

test_random_test()

function test_graham_scan()

    N_point = 2000
    point_array = geom.random_point(N_point)

    ch = geom.naive_graham_scan!(point_array)

    (ch_x, ch_y) = geom.discretize_for_plot(ch, true)
    (xs, ys) = geom.discretize_for_plot(point_array)

    PyPlot.figure()
    PyPlot.grid()
    PyPlot.plot(xs, ys, linestyle = :none, marker=:., markersize = 10)

    PyPlot.plot(ch_x, ch_y, linewidth=2)

    # n_new = geom.reduce_array_for_convex_hull(point_array)

end

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Dipole model for the Earth geomagnetic field.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] http://helios.fmi.fi/~juusolal/geomagnetism/Lectures/Chapter3_dipole.pdf
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export geomag_dipole

"""
    function geomag_dipole(r_e::AbstractVector, pole_lat::Number, pole_lon::Number, m::Number)

Compute the geomagnetic field [nT] using the simplified dipole model at position
`r_e` (ECEF reference frame). This function considers that the latitude of the
South magnetic pole (which lies in the North hemisphere) is `pole_lat` [rad] and
the longitude is `pole_lon` [rad]. Furthermore, the dipole moment is considered
to be `m` [A.m²].

    function geomag_dipole(r_e::AbstractVector, year::Number = 2019)

Compute the geomagnetic field [nT] using the simplified dipole model at position
`r_e` (ECEF reference frame). This function uses the year `year` to obtain the
position of the South magnetic pole (which lies in the North hemisphere) and the
dipole moment. If `year` is omitted, then it will be considered as 2019.

# Remarks

In both functions, the output vector will be represented in the ECEF reference
frame.

"""
function geomag_dipole(r_e::AbstractVector, year::Number = 2019)
    d2r = π/180

    # Obtain the pole location and the dipole moment.
    pole_lat = _itp_dipole_lat(year)*d2r
    pole_lon = _itp_dipole_lon(year)*d2r
    m        = _itp_dipole_mag(year)*1e22

    return geomag_dipole(r_e, pole_lat, pole_lon, m)
end

function geomag_dipole(r_e::AbstractVector, pole_lat::Number, pole_lon::Number,
                       m::Number)

    # DCM that convertes the ECEF into the geomagnetic coordinates.
    Dge = angle_to_dcm(pole_lon, π/2-pole_lat, 0, :ZYX)

    # Compute the dipole momentum represented in the ECEF reference frame.
    k₀_e = 1e-7m*(Dge'*SVector{3}(0,0,-1))

    # Compute the distance from the Earth center of the desired point.
    r = norm(r_e)

    # Compute the unitary vector that points to the desired direction.
    eᵣ_e = SVector{3}(r_e/r)

    # Compute the geomagnetic field vector [nT].
    B_e = 1/r^3*( 3eᵣ_e*eᵣ_e' - I )*k₀_e*1e9

    return B_e
end


text:
"""
    Stateful{T, S}
    Stateful(schedule::T; advance = state -> true)

Create a stateful iterator around `schedule`.
Pass in a predicate, `advance(state)`, to conditionally control iteration.
See also [`ParameterSchedulers.next!`](#) and [`ParameterSchedulers.reset!`](#).
"""
mutable struct Stateful{T, S<:Integer, R}
    schedule::T
    state::S
    advance::R
end
Stateful(schedule; advance = state -> true) = Stateful(schedule, 1, advance)

"""
    next!(iter::Stateful)

Advance `iter` by one iteration
(if `iter.advance(state) == true`) and return the next value.
See also [`ParameterSchedulers.Stateful`](#).
"""
function next!(iter::Stateful)
    val = iter.schedule(iter.state)
    if iter.advance(iter.state)
        iter.state += 1
    end

    return val
end

"""
    reset!(iter::Stateful)

Reset `iter` to its initial state.
See also [`ParameterSchedulers.Stateful`](#).
"""
function reset!(iter::Stateful)
    iter.state = 1

    return iter
end

at(x::Number, t) = x
at(f, t) = f(t)

"""
    Sequence{T, S}
    Sequence(schedules, step_sizes)
    Sequence(schedule1 => step1, schedule2 => step2, ...)

A sequence of schedules.
The output of this schedule is the concatenation of `schedules` where each
schedule is evaluated for each step size in `step_sizes`.

Note that `schedules` can also be a vector of numbers (not just schedules).

# Arguments
- `schedules`: a vector of schedules or numbers
- `step_sizes`: a vector of iteration lengths for each schedule
"""
struct Sequence{T, S}
    schedules::T
    step_sizes::S
end
Sequence(stages::Pair...) = Sequence(first.(stages), last.(stages))

function (schedule::Sequence)(t)
    accum_steps = cumsum(schedule.step_sizes)
    i = findlast(x -> t > x, accum_steps)
    i = isnothing(i) ? 1 :
            (i >= length(schedule.schedules)) ? length(schedule.schedules) : i + 1
    toffset = (i > 1) ? t - accum_steps[i - 1] : t

    return at(schedule.schedules[i], toffset)
end

Base.IteratorEltype(::Type{<:Sequence}) = Base.EltypeUnknown()
Base.IteratorSize(::Type{<:Sequence}) = Base.SizeUnknown()

function Base.iterate(schedule::Sequence, state = (1, 1, 1))
    t, i, t0 = state
    if (i < length(schedule.step_sizes)) && (t >= t0 + schedule.step_sizes[i])
        # move onto next step range
        i += 1
        t0 = t
    end

    return at(schedule.schedules[i], t - t0 + 1), (t + 1, i, t0)
end

"""
    Loop{T, S<:Integer}
    Loop(f, period)

Create a schedule that loops `f` every `period` iterations.
`f` must be callabe (a function or schedule).

# Arguments
- `f`: the schedule to loop
- `period::Integer`: how often to loop
"""
struct Loop{T, S<:Integer}
    f::T
    period::S
end
Loop(f, period) = Loop(f, period)

(schedule::Loop)(t) = schedule.f(mod1(t, schedule.period))

Base.IteratorEltype(::Type{<:Loop{T}}) where T = Base.IteratorEltype(T)
Base.eltype(::Type{<:Loop{T}}) where T = eltype(T)
Base.IteratorSize(::Type{<:Loop}) = Base.IsInfinite()

Base.iterate(schedule::Loop, t = 1) = schedule(t), t + 1

Base.axes(::Loop) = (OneToInf(),)

"""
    reverse(f, period)

Return a reverse function such that `reverse(f, period)(t) == f(period - t)`.
"""
reverse(f, period) = t -> f(period - t)
"""
    symmetric(f, period)

Return a symmetric function such that for `t ∈ [1, period / 2)`,
the symmetric function evaluates to `f(t)`, and when `t ∈ [period / 2, period)`,
the symmetric functions evaluates to `f(period - t)`.
"""
symmetric(f, period) = t -> (t < period / 2) ? f(t) : f(period - t)


text:
# Keeps compatibility with bariess code moved to Base/stdlib on older releases

using LinearAlgebra
using SparseArrays
using SparseArrays: AbstractSparseMatrixCSC, getcolptr

macro swap(a, b)
    esc(:(($a, $b) = ($b, $a)))
end

# https://github.com/JuliaLang/julia/pull/42678
@static if VERSION > v"1.8.0-DEV.762"
    import Base: swaprows!
else
    function swaprows!(a::AbstractMatrix, i, j)
        i == j && return
        rows = axes(a,1)
        @boundscheck i in rows || throw(BoundsError(a, (:,i)))
        @boundscheck j in rows || throw(BoundsError(a, (:,j)))
        for k in axes(a,2)
            @inbounds a[i,k],a[j,k] = a[j,k],a[i,k]
        end
    end
    function Base.circshift!(a::AbstractVector, shift::Integer)
        n = length(a)
        n == 0 && return
        shift = mod(shift, n)
        shift == 0 && return
        reverse!(a, 1, shift)
        reverse!(a, shift+1, length(a))
        reverse!(a)
        return a
    end
    function Base.swapcols!(A::AbstractSparseMatrixCSC, i, j)
        i == j && return

        # For simplicitly, let i denote the smaller of the two columns
        j < i && @swap(i, j)

        colptr = getcolptr(A)
        irow = colptr[i]:(colptr[i+1]-1)
        jrow = colptr[j]:(colptr[j+1]-1)

        function rangeexchange!(arr, irow, jrow)
            if length(irow) == length(jrow)
                for (a, b) in zip(irow, jrow)
                    @inbounds @swap(arr[i], arr[j])
                end
                return
            end
            # This is similar to the triple-reverse tricks for
            # circshift!, except that we have three ranges here,
            # so it ends up being 4 reverse calls (but still
            # 2 overall reversals for the memory range). Like
            # circshift!, there's also a cycle chasing algorithm
            # with optimal memory complexity, but the performance
            # tradeoffs against this implementation are non-trivial,
            # so let's just do this simple thing for now.
            # See https://github.com/JuliaLang/julia/pull/42676 for
            # discussion of circshift!-like algorithms.
            reverse!(@view arr[irow])
            reverse!(@view arr[jrow])
            reverse!(@view arr[(last(irow)+1):(first(jrow)-1)])
            reverse!(@view arr[first(irow):last(jrow)])
        end
        rangeexchange!(rowvals(A), irow, jrow)
        rangeexchange!(nonzeros(A), irow, jrow)

        if length(irow) != length(jrow)
            @inbounds colptr[i+1:j] .+= length(jrow) - length(irow)
        end
        return nothing
    end
    function swaprows!(A::AbstractSparseMatrixCSC, i, j)
        # For simplicitly, let i denote the smaller of the two rows
        j < i && @swap(i, j)

        rows = rowvals(A)
        vals = nonzeros(A)
        for col = 1:size(A, 2)
            rr = nzrange(A, col)
            iidx = searchsortedfirst(@view(rows[rr]), i)
            has_i = iidx <= length(rr) && rows[rr[iidx]] == i

            jrange = has_i ? (iidx:last(rr)) : rr
            jidx = searchsortedlast(@view(rows[jrange]), j)
            has_j = jidx != 0 && rows[jrange[jidx]] == j

            if !has_j && !has_i
                # Has neither row - nothing to do
                continue
            elseif has_i && has_j
                # This column had both i and j rows - swap them
                @swap(vals[rr[iidx]], vals[jrange[jidx]])
            elseif has_i
                # Update the rowval and then rotate both nonzeros
                # and the remaining rowvals into the correct place
                rows[rr[iidx]] = j
                jidx == 0 && continue
                rotate_range = rr[iidx]:jrange[jidx]
                circshift!(@view(vals[rotate_range]), -1)
                circshift!(@view(rows[rotate_range]), -1)
            else
                # Same as i, but in the opposite direction
                @assert has_j
                rows[jrange[jidx]] = i
                iidx > length(rr) && continue
                rotate_range = rr[iidx]:jrange[jidx]
                circshift!(@view(vals[rotate_range]), 1)
                circshift!(@view(rows[rotate_range]), 1)
            end
        end
        return nothing
    end
end

function bareiss_update!(zero!, M::StridedMatrix, k, swapto, pivot, prev_pivot::Base.BitInteger)
    flag = zero(prev_pivot)
    prev_pivot = Base.MultiplicativeInverses.SignedMultiplicativeInverse(prev_pivot)
    @inbounds for i in k+1:size(M, 2)
        Mki = M[k,i]
        @simd ivdep for j in k+1:size(M, 1)
            M[j,i], r = divrem(M[j,i]*pivot - M[j,k]*Mki, prev_pivot)
            flag = flag | r
        end
    end
    iszero(flag) || error("Overflow occurred")
    zero!(M, k+1:size(M, 1), k)
end

function bareiss_update!(zero!, M::StridedMatrix, k, swapto, pivot, prev_pivot)
    @inbounds for i in k+1:size(M, 2), j in k+1:size(M, 1)
        M[j,i] = exactdiv(M[j,i]*pivot - M[j,k]*M[k,i], prev_pivot)
    end
    zero!(M, k+1:size(M, 1), k)
end

@views function bareiss_update!(zero!, M::AbstractMatrix, k, swapto, pivot, prev_pivot)
    if prev_pivot isa Base.BitInteger
        prev_pivot = Base.MultiplicativeInverses.SignedMultiplicativeInverse(prev_pivot)
    end
    V = M[k+1:end, k+1:end]
    V .= exactdiv.(V .* pivot .- M[k+1:end, k] * M[k, k+1:end]', prev_pivot)
    zero!(M, k+1:size(M, 1), k)
end

function bareiss_update_virtual_colswap!(zero!, M::AbstractMatrix, k, swapto, pivot, prev_pivot)
    if prev_pivot isa Base.BitInteger
        prev_pivot = Base.MultiplicativeInverses.SignedMultiplicativeInverse(prev_pivot)
    end
    V = @view M[k+1:end, :]
    V .= @views exactdiv.(V .* pivot .- M[k+1:end, swapto[2]] * M[k, :]', prev_pivot)
    zero!(M, k+1:size(M, 1), swapto[2])
end

bareiss_zero!(M, i, j) = M[i,j] .= zero(eltype(M))

function find_pivot_col(M, i)
    p = findfirst(!iszero, @view M[i,i:end])
    p === nothing && return nothing
    idx = CartesianIndex(i, p + i - 1)
    (idx, M[idx])
end

function find_pivot_any(M, i)
    p = findfirst(!iszero, @view M[i:end,i:end])
    p === nothing && return nothing
    idx = p + CartesianIndex(i - 1, i - 1)
    (idx, M[idx])
end

const bareiss_colswap = (Base.swapcols!, swaprows!, bareiss_update!, bareiss_zero!)
const bareiss_virtcolswap = ((M,i,j)->nothing, swaprows!, bareiss_update_virtual_colswap!, bareiss_zero!)

"""
    bareiss!(M, [swap_strategy])

Perform Bareiss's fraction-free row-reduction algorithm on the matrix `M`.
Optionally, a specific pivoting method may be specified.

swap_strategy is an optional argument that determines how the swapping of rows and coulmns is performed.
bareiss_colswap (the default) swaps the columns and rows normally.
bareiss_virtcolswap pretends to swap the columns which can be faster for sparse matrices.
"""
function bareiss!(M::AbstractMatrix{T}, swap_strategy=bareiss_colswap;
                  find_pivot=find_pivot_any, column_pivots=nothing) where T
    swapcols!, swaprows!, update!, zero! = swap_strategy
    prev = one(eltype(M))
    n = size(M, 1)
    pivot = one(T)
    column_permuted = false
    for k in 1:n
        r = find_pivot(M, k)
        r === nothing && return (k - 1, pivot, column_permuted)
        (swapto, pivot) = r
        if column_pivots !== nothing && k != swapto[2]
            column_pivots[k] = swapto[2]
            column_permuted |= true
        end
        if CartesianIndex(k, k) != swapto
            swapcols!(M, k, swapto[2])
            swaprows!(M, k, swapto[1])
        end
        update!(zero!, M, k, swapto, pivot, prev)
        prev = pivot
    end
    return (n, pivot, column_permuted)
end

function nullspace(A; col_order=nothing)
    column_pivots = collect(1:size(A, 2))
    B = copy(A)
    (rank, d, column_permuted) = bareiss!(B; column_pivots)
    reduce_echelon!(B, rank, d)

    # The first rank entries in col_order are columns that give a basis
    # for the column space. The remainder give the free variables.
    if col_order !== nothing
        resize!(col_order, size(A,2))
        col_order .= 1:size(A,2)
        for (i,cp) in enumerate(column_pivots)
            @swap(col_order[i],col_order[cp])
        end
    end

    N = ModelingToolkit.reduced_echelon_nullspace(rank, B)
    apply_inv_pivot_rows!(N, column_pivots)
end

function apply_inv_pivot_rows!(M, ipiv)
    for i in size(M, 1):-1:1
        swaprows!(M, i, ipiv[i])
    end
    M
end

###
### Modified from AbstractAlgebra.jl
###
### https://github.com/Nemocas/AbstractAlgebra.jl/blob/4803548c7a945f3f7bd8c63f8bb7c79fac92b11a/LICENSE.md
function reduce_echelon!(A::AbstractMatrix{T}, rank, d) where T
    m, n = size(A)
    for i = rank + 1:m
        for j = 1:n
            A[i, j] = zero(T)
        end
    end
    if rank > 1
        t = zero(T)
        q = zero(T)
        d = -d
        pivots = zeros(Int, n)
        np = rank
        j = k = 1
        for i = 1:rank
            while iszero(A[i, j])
                pivots[np + k] = j
                j += 1
                k += 1
            end
            pivots[i] = j
            j += 1
        end
        while k <= n - rank
            pivots[np + k] = j
            j += 1
            k += 1
        end
        for k = 1:n - rank
            for i = rank - 1:-1:1
                t = A[i, pivots[np + k]] * d
                for j = i + 1:rank
                    t += A[i, pivots[j]] * A[j, pivots[np + k]] + q
                end
                A[i, pivots[np + k]] = exactdiv(-t, A[i, pivots[i]])
            end
        end
        d = -d
        for i = 1:rank
            for j = 1:rank
                if i == j
                    A[j, pivots[i]] = d
                else
                    A[j, pivots[i]] = zero(T)
                end
            end
        end
    end
    return A
end

function reduced_echelon_nullspace(rank, A::AbstractMatrix{T}) where T
    n = size(A, 2)
    nullity = n - rank
    U = zeros(T, n, nullity)
    if rank == 0
        for i = 1:nullity
            U[i, i] = one(T)
        end
    elseif nullity != 0
        pivots = zeros(Int, rank)
        nonpivots = zeros(Int, nullity)
        j = k = 1
        for i = 1:rank
            while iszero(A[i, j])
                nonpivots[k] = j
                j += 1
                k += 1
            end
            pivots[i] = j
            j += 1
        end
        while k <= nullity
            nonpivots[k] = j
            j += 1
            k += 1
        end
        d = -A[1, pivots[1]]
        for i = 1:nullity
            for j = 1:rank
                U[pivots[j], i] = A[j, nonpivots[i]]
            end
            U[nonpivots[i], i] = d
        end
    end
    return U
end


text:
module XMLconvert

using LightXML
using DataStructures

export xml2json, xml2dict, xmlASCII2file, parse_file, parse_string,
       root, show_key_structure, flatten, attributes_to_elements!

include("xml_to_json.jl")
include("xml_to_dict.jl")
include("xml_to_file.jl")
include("utils.jl")
include("flatten_multidict.jl")

end # module


text:
mutable struct FPGroupIter{S, T, GEl}
    seen::S
    seen_iter_state::T
    current::GEl
    gen_idx::Int
    u_tmp::GEl
    v_tmp::GEl
end

function Base.iterate(G::AbstractFPGroup)
    seen = OrderedSet([one(G)])
    current, seen_state = iterate(seen)
    gr_iter = FPGroupIter(seen, seen_state, current, 1, one(G), one(G))
    return one(G), gr_iter
end

function _next_elt!(itr::FPGroupIter)
    res = iterate(itr.seen, itr.seen_iter_state)
    res === nothing && return true

    itr.current = first(res)
    itr.seen_iter_state = last(res)
    itr.gen_idx = 1

    return false
end

function Base.iterate(G::AbstractFPGroup, state)
    if state.gen_idx > length(alphabet(G))
        finished = _next_elt!(state)
        finished && return nothing
    end

    elt = let u = state.u_tmp, v = state.v_tmp, current = state.current
        copyto!(v, current)
        push!(word(v), state.gen_idx)
        _setnormalform!(v, false)
        _setvalidhash!(v, false)
        @assert !isnormalform(v)

        resize!(word(u), 0)
        normalform!(u, v)
    end

    state.gen_idx += 1

    if elt in state.seen
        return iterate(G, state)
    else
        @assert isnormalform(elt)
        push!(state.seen, deepcopy(elt))
        return elt, state
    end
end

# Groups.Core default:
# Base.IteratorSize(::Type{<:AbstractFPGroup}) = Base.SizeUnknown()


text:
"""
Collections of utility functions for data simulation and evaluations

Author: Wayne Wang
Last modified: 10/06/2021
"""

# Precision estimation methods type
struct SG_PALM end
struct TERALASSO end
struct KGLASSO end
struct GLASSO end
# Covariance estimation methods type
struct KPCA end
# Customized precision matrix structure types
struct AR1 end
struct SB end
struct UNIF end
# Customized Kalman Filter process types
struct POISSON end
struct POISSON_AR end
struct CONVECTION_DIFFUSION end
# Customized Kalman Filter observation types
struct IDENTITY end
struct LINEAR_PERM end
struct LINEAR_PERM_MISS end

function method_str_to_type(method::AbstractString)
    if method == "sg_palm"
        method = SG_PALM()
    elseif method == "teralasso"
        method = TERALASSO()
    elseif method == "kglasso"
        method = KGLASSO()
    elseif method == "glasso"
        method = GLASSO()
    elseif method == "kpca"
        method = KPCA()
    else
        print("Precision/covariance matrix estimation method unsupported!")
    end
    return method
end

function method_type_to_str(method::Union{SG_PALM, TERALASSO, KGLASSO, GLASSO, KPCA})
    if method == SG_PALM()
        method = "sg_palm"
    elseif method == TERALASSO()
        method = "teralasso"
    elseif method == KGLASSO()
        method = "kglasso"
    elseif method == GLASSO()
        method = "glasso"
    elseif method == KPCA()
        method = "kpca"
    else
        print("Precision/covariance matrix estimation method unsupported!")
    end
    return method
end

function inv_cov_model(model_type::AR1, p::Int; ρ::Real = 0.5)
    """
    Create sparse precision matrix according to one of Autoregressive1 (AR),
        Uniform distributed weights (UNIF), and Start-Block (SB) models
    In:
        - p::Int: dimension of the cov matrix
        - ρ::Real: weight param between 0 and 1
    Out:
        - inv_cov::SparseMatrixCSC
    """
    inv_cov = sparse(Diagonal(repeat([1 + ρ^2], p)))
    for i = 1:(p-1)
        inv_cov[i, i+1] = -ρ
        inv_cov[i+1, i] = -ρ
    end

    return inv_cov
end

function inv_cov_model(model_type::UNIF, p::Int; dens::Real = 0.1, min_eigen::Real = 1e-3)
    """
    Create sparse precision matrix according to one of Autoregressive (AR),
        Uniform distributed weights (UNIF), and Start-Block (SB) models
    In:
        - p: the desired dimension of the (inv)covariance matrix
        - dens: percentage of non-zeros of the final inverse covariance
        - min_eigen: the minimum bound for the eigenvalues of the desired (inv)covariance matrix
    Out:
        - inv_cov::SparseMatrixCSC
    """
    U = sprand(p, p, dens)
    for i = 1:p
        for j = 1:p
            if U[i, j] != 0
                U[i, j] = rand([-1.0, 1.0], 1)[1]
            end
        end
    end
    U .= copy(U') * U
    Udiag = Diagonal(U) .+ Diagonal(ones(p))
    Uoff = max.(min.(U .- Diagonal(U), 1.0), -1.0)
    inv_cov = Uoff .+ Udiag
    eigs, _ = eigen(collect(inv_cov))
    inv_cov .+= Diagonal(repeat([max(-1.2 * minimum(eigs), min_eigen)], p))

    return inv_cov
end

function inv_cov_model(model_type::SB, p::Int; ρ::Real = 0.5, num_subgraph::Int = 4,
    tol::Real = 1e-5)
    """
    Create sparse precision matrix according to one of Autoregressive (AR),
        Uniform distributed weights (UNIF), and Start-Block (SB) models
    In:
        - p::Int: dimension of the cov matrix
        - ρ::Real: weight param between 0 and 1
        - num_subgraph::Int: number of star structures
        - tol::Real: thresholding value for zeros
    Out:
        - inv_cov::SparseMatrixCSC
    """
    p_subgraph = Int(floor(p / num_subgraph))
    uneq_graph = (p != p_subgraph * num_subgraph)

    list_inv_subgraph = []

    for i = 1:num_subgraph
        if i == num_subgraph && uneq_graph
            p_subgraph = p - p_subgraph * (i - 1)
        end
        central_node = Int(rand(Tuple(collect(1:p_subgraph))))
        subgraph = repeat([ρ^2], p_subgraph, p_subgraph)
        subgraph[:, central_node] .= ρ
        subgraph[central_node, :] .= ρ
        subgraph .= subgraph .- Diagonal(subgraph) .+ Diagonal(ones(p_subgraph))
        inv_subgraph = inv(subgraph)
        push!(list_inv_subgraph, sparse(inv_subgraph))
    end

    inv_cov = SparseArrays.blockdiag(list_inv_subgraph...)
    inv_cov[abs.(inv_cov) .< tol] .= 0

    return inv_cov
end

function metric_edge(true_edge, est_edge)
    """
    Computes FNR, FPR, MCC, etc
    In
    - true_edge::AbstractArrary: true sparse precision matrix
    - est_edge::AbstractArray: estimated sparse precision matrix
    Out
    - metric_res::Array: FP, FN, TP, TN, FPR, FNR, Precision, Recall, MCC
    """
    p = size(true_edge, 1)
    nz_pos = (true_edge .!= 0)
    z_pos = (true_edge .== 0)

    FP = (sum((est_edge[z_pos] .!= 0)))/2
    FN = (sum((est_edge[nz_pos] .== 0)))/2
    TP = (sum((est_edge[nz_pos] .!= 0)) - p)/2
    TN = (sum((est_edge[z_pos] .== 0)))/2

    FPR = FP / (FP + TN)
    FNR = FN / (TP + FN)
    Precision = TP / (TP + FP)
    Recall = TP / (TP + FN)
    MCC = (TP * TN - FP * FN) / sqrt((TP + FP) * (TP + FN) * (TN + FP) * (TN + FN))

    metric_res = [FP, FN, TP, TN, FPR, FNR, Precision, Recall, MCC]

    return metric_res
end

function gen_kronecker_data(model::AbstractString, inv_cov_model_type::AbstractString, K::Int,
    N::Int, d_list::AbstractVector; Ψ_list::AbstractArray = Array{Array{Float64, 2}, 1}(),
    tensorize_out::Bool = true)
    """
    Generate multivariate Gaussian data with precision matrix satisfying either the
        Kronecker product, Kronecker sum, or Sylvester structre.
    """
    if inv_cov_model_type == "ar1"
        inv_cov_model_type = AR1()
    elseif inv_cov_model_type == "sb"
        inv_cov_model_type = SB()
    elseif inv_cov_model_type == "unif"
        inv_cov_model_type = UNIF()
    else
        print("Unsupported inverse covariance model type")
    end

    # create each factor matrix if not provided
    if length(Ψ_list) == 0
        for k = 1:K
            push!(Ψ_list, inv_cov_model(inv_cov_model_type, d_list[k]))
        end
        d = prod(d_list)
    else
        K = length(Ψ_list)
        d_list = size.(Ψ_list, 1)
        d = prod(d_list)
    end

    # eigen decomposition for each Ψ_k
    U_list = []; Λ_list = []
    for k = 1:K
        Λk, Uk = eigen(collect(Ψ_list[k]))
        push!(Λ_list, Diagonal(Λk))
        push!(U_list, Uk)
    end
    reverse!(U_list)
    eigvecsSigmaSqrt = kron(U_list...)

    if model == "kp"
        eigsSigmaSqrt = inv(Diagonal(kron(Λ_list...)).^(1/2))
    elseif model == "ks"
        eigsSigmaSqrt = inv(Diagonal(kroneckersum_list(Λ_list)).^(1/2))
    elseif model == "sylvester"
        eigsSigmaSqrt = inv(Diagonal(kroneckersum_list(Λ_list)))
    else
        print("Unsupported model type")
    end

    # generate fake data
    X = zeros(d, N)
    z = zeros(d)
    xtilde = zeros(d)
    x = zeros(d)
    for i = 1:N
        randn!(z)
        mul!(xtilde, eigsSigmaSqrt, z)
        mul!(x, eigvecsSigmaSqrt, xtilde)
        X[:, i] .= x
    end

    if tensorize_out
        return matten(copy(X'), K + 1, push!(d_list, N))
    end

    return X
end

function gen_kalmanfilter_data(dynamic_type::AbstractString, obs_type::AbstractString,
    T::Int, px::Tuple, py::Tuple,
    obs_noise::Real, process_noise::Real,
    add_process_noise::Bool)
    """
    Generate Kalman Filtering obervations
    Input:
        - T: number of time steps
        - px: 2-tuple for process dimension
        - py: 2-tuple for observation dimension
        - obs_noise: obervations will be subject to Gaussian noise N(0,obs_noise*I)
        - process_noise: processes will be subject to Gaussian noise N(0,process_noise*I)
        - add_process_noise: whether to introduce process noise
    Output:
        - (X, Y): T true processes and observations
    """
    # initialization
    X0 = rand(MatrixNormal(2 * randn(px), ScalMat(px[1], 2), ScalMat(px[2], 2)))
    X = zeros((prod(px), T + 1))
    Y = (obs_type == "linear_perm_miss") ? zeros((Int(0.5 * prod(py)), T)) : zeros((prod(py), T))
    X[:, 1] .= X0[:]

    # swtich observation types
    if obs_type == "identity"
        obs_type = IDENTITY()
    elseif obs_type == "linear_perm"
        obs_type = LINEAR_PERM()
    elseif obs_type == "linear_perm_miss"
        obs_type = LINEAR_PERM_MISS()
    else
        print("Observation type type unsupported!")
    end
    H = measure_operator(obs_type, py, px)

    # evolution
    for t = 1:T
        # dynamics update
        X[:, t + 1] .= kalmanfilter_dynamic_update(dynamic_type, X[:, t], px, add_process_noise, process_noise)
        # observations update
        Y[:, t] .= kalmanfilter_observation_update(H, X[:, t + 1], py, px, obs_noise)
    end

    return X, Y, H
end

function kalmanfilter_dynamic_update(dynamic_type::AbstractString, X_curr::AbstractArray,
    px::Tuple,
    add_process_noise::Bool,
    process_noise::Real;
    α::Real = 0.05,
    ϵ::Real = 0.01,
    Δx::Real = 0.005,
    Δt::Real = 0.0005)
    # step sizes for 2D convection-diffusion PDE discretization
    # convergence/stability criteria: Δt <= (Δx^2)/(4*α)
    # diffusion constant α = 0.5
    # convection constant ϵ = 0.5
    # assume same spatial steps in x&y Δx = 1/200
    # time steps Δt = 0.0005

    # swtich observation types
    if dynamic_type == "poisson"
        dynamic_type = POISSON()
    elseif dynamic_type == "poisson_ar"
        dynamic_type = POISSON_AR()
    elseif dynamic_type == "convection_diffusion"
        dynamic_type = CONVECTION_DIFFUSION()
    else
        print("Dynamic type unsupported!")
    end

    # dynamics update, i.e., solve the Sylvester equation
    X_new = similar(X_curr)
    sylv_eqn_solver!(X_new, X_curr, px, dynamic_type; α = α, ϵ = ϵ, Δx = Δx, Δt = Δt)

    if add_process_noise
        w = zeros(prod(px))
        rand!(MvNormal(prod(px), process_noise), w)
        X_new .+= w
    end

    return X_new
end

function sylv_eqn_solver!(X_new::AbstractArray, X_curr::AbstractArray, px::Tuple,
    type::CONVECTION_DIFFUSION;
    α::Real = 0.05,
    ϵ::Real = 0.01,
    Δx::Real = 0.005,
    Δt::Real = 0.0005)
    # step sizes for 2D convection-diffusion PDE discretization
    # convergence/stability criteria: Δt <= (Δx^2)/(4*α)
    # diffusion constant α = 0.5
    # convection constant ϵ = 0.5
    # assume same spatial steps in x&y Δx = 1/200
    # time steps Δt = 0.0005
    γ = (α * Δt) / Δx^2
    # 2D oncvection-diffusion equation operator
    A = spdiagm(-1 => [-(1 + ϵ / (2 * α) * Δx) * γ for _ = 1:(px[1] - 1)],
                0 => [(4 + 1 / γ) * γ / 2 for _ = 1:px[1]],
                1 => [-(1 - ϵ / (2 * α) * Δx) * γ for _ = 1:(px[1] - 1)])
    # X = sylvc(A,B,C)
    # Solve the continuous Sylvester matrix equation AX + XB = C
    X_new .= sylvc(A, A, reshape(X_curr, px))[:]
    return nothing
end

function sylv_eqn_solver!(X_new::AbstractArray, X_curr::AbstractArray, px::Tuple,
    type::POISSON_AR;
    α::Real = 0.05,
    ϵ::Real = 0.01,
    Δx::Real = 1/200,
    Δt::Real = 0.0005)
    # 2D oncvection-diffusion equation operator
    A = spdiagm(-1 => [-1.0 for _ = 1:(px[1] - 1)],
                0 => [2.0 for _ = 1:px[1]],
                1 => [-1.0 for _ = 1:(px[1]-1)])
    B = spdiagm(-1 => [-1.0 for _ = 1:(px[1] - 1)],
                0 => [2.0 for _ = 1:px[1]],
                1 => [-1.0 for _ = 1:(px[1] - 1)])
    # X = gsylv(A,B,E)
    # Solve the generalized Sylvester matrix equation AXB  = E.
    X_new .= gsylv(A, B, randn(px))[:]
    return nothing
end

function sylv_eqn_solver!(X_new::AbstractArray, X_curr::AbstractArray, px::Tuple,
    type::POISSON;
    α::Real = 0.05,
    ϵ::Real = 0.01,
    Δx::Real = 1/200,
    Δt::Real = 0.0005)
    # 2D oncvection-diffusion equation operator
    A = spdiagm(-1 => [-1.0 for _ = 1:(px[1] - 1)],
                0 => [2.0 for _ = 1:px[1]],
                1 => [-1.0 for _ = 1:(px[1] - 1)])
    # X = sylvc(A,B,C)
    # Solve the continuous Sylvester matrix equation AX + XB = C
    # X_new .= sylvc(A, A, (Δx^2) * randn(px))[:]
    X_new .= sylvc(A, A, randn(px))[:]
    return nothing
end

function kalmanfilter_observation_update(H::AbstractArray, X::AbstractArray,
    py::Tuple,
    px::Tuple,
    obs_noise::Real)
    # measurement operation plus noise
    v = zeros(size(H, 1))
    rand!(MvNormal(length(v), obs_noise), v)
    Y = H * X .+ v
    return Y
end

function measure_operator(type::IDENTITY,
    py::Tuple,
    px::Tuple)
    # measurement operator
    H = 1.0 * I(prod(px))
    return H
end

function measure_operator(type::LINEAR_PERM,
    py::Tuple,
    px::Tuple)
    # measurement operator
    H = zeros(prod(py), prod(px))
    perm = randperm(prod(px))
    for i = 1:prod(px)
        H[i, perm[i]] = 1.0
    end
    return H
end

function measure_operator(type::LINEAR_PERM_MISS,
    py::Tuple,
    px::Tuple)
    # measurement operator
    H = zeros(prod(py), prod(px))
    perm = randperm(prod(px))
    for i = 1:prod(px)
        H[i, perm[i]] = 1.0
    end
    H = H[randperm(prod(px))[1:Int(prod(px)/2)], :] # half of the rows missing
    return H
end

# # gif plot all time steps
# function plot_state_gif(X::AbstractArray, px::Tuple;
#     out_dir::AbstractString = "",
#     clim::Tuple = ())
#     """
#     Plot all time steps of the state images
#     """
#     anim = @animate for i = 2:(T + 1)
#         plot_state_img(X, px, i, clim = clim)
#     end

#     if length(out_dir) != 0
#         gif(anim, out_dir, fps = 5)
#     else
#         gif(anim, fps = 5)
#     end

#     return nothing
# end

# function plot_state_img(X::AbstractArray, px::Tuple, time_step::Int;
#     clim::Tuple = ())
#     """
#     Plot a single state img
#     """
#     if length(clim) != 0
#         display(plot(jim(reshape(X[:, time_step], px), clim = clim),
#                     title = string("Time step: ", time_step)))
#     end

#     display(plot(jim(reshape(X[:, time_step], px)),
#                 title = string("Time step: ", time_step)))

#     return nothing
# end

function plot_state_img(X::AbstractArray, px::Tuple, time_step::Int;
    vmin = nothing, vmax = nothing, cmap::AbstractString = "inferno")
    """
    Plot a single state img *using PyPlot*
    """
    fig = figure()

    if !isnothing(vmin) && !isnothing(vmax)
        PyPlot.imshow(reshape(X[:, time_step], px), cmap = cmap, vmin = vmin, vmax = vmax)
    else
        PyPlot.imshow(reshape(X[:, time_step], px), cmap = cmap)
    end

    PyPlot.axis("off")
    PyPlot.title(string("Time step: ", time_step))
    PyPlot.colorbar()
    PyPlot.display(fig)

    return nothing
end

# plot nrmse
function plot_nrmse!(NRMSEs::AbstractArray, method::AbstractString)
    T, N = size(NRMSEs)
    ## NRMSEs over time steps
    ## showing the mean NRMSEs and the 90% CI
    Plots.plot!(1:T, mean(NRMSEs, dims=2),
            yaxis = :log,
            ribbon = sort(NRMSEs, dims = 2),
            fillalpha = .2,
            markerstrokecolor = :auto,
            label = string("Mean RMSE ", method)
        )

    return nothing
end

# compute nrmse
function compute_nrmse(X::AbstractArray, Xhat::AbstractArray)
    ## compute NRMSEs for all ensemble members
    T = size(X, 2)
    N = size(Xhat, 3)
    NRMSEs = zeros((T - 1, N))

    for t = 2:T
        for i = 1:N
            NRMSEs[t - 1, i] = sqrt(mean((Xhat[t, :, i] .- X[:, t]).^2)) #/(maximum(X[:])-minimum(X[:]))
        end
    end

    return NRMSEs
end


text:
function test_tntp()
  data_dir = download_tntp()

  # Test
  for d in readdir(data_dir)
    if isdir(joinpath(data_dir, d))
      try
        read_ta_network(d)
        load_ta_network(d)
        @info "Network '$d' is OK."
      catch e
        @show e
        @warn "Network '$d' is not usable."
      end
    end
  end
end

function summarize_ta_data()
  data_dir = download_tntp()

  # Test
  dic = Dict()
  for d in readdir(data_dir)
    if isdir(joinpath(data_dir, d))
      try
        network_data_file, trip_table_file = read_ta_network(d)
        number_of_zones, number_of_links, number_of_nodes = read_ta_summary(network_data_file)
        dic[d] = (number_of_zones, number_of_links, number_of_nodes)
      catch

      end

    end
  end

  max_len = [0, 0, 0, 0]
  for net in sort(dic)
    max_len[1] = max( length(net[1]), max_len[1] )
    max_len[2] = max( length(digits(net[2][1])), max_len[2] )
    max_len[3] = max( length(digits(net[2][2])), max_len[3] )
    max_len[4] = max( length(digits(net[2][3])), max_len[4] )
  end

  function format(mlen, val)
    len = 0
    str = ""
    if isa(val, String)
      len = length(val)
      str = val * " "^(mlen - len + 1 )
    elseif isa(val, Number)
      len = length(digits(val))
      str = " "^(mlen - len + 1) * string(val)
    end
    return str
  end

  println("-"^(17+sum(max_len)))
  println("| $(format(max_len[1],"Network")) | $(format(max_len[2],"Zones")) | $(format(max_len[3],"Links")) | $(format(max_len[4],"Nodes")) |")
  println("| $(format(max_len[1],":---")) | $(format(max_len[2],"---:")) | $(format(max_len[3],"---:")) | $(format(max_len[4],"---:")) |")
  for net in sort(dic)
    println("| $(format(max_len[1],net[1])) | $(format(max_len[2],net[2][1])) | $(format(max_len[3],net[2][2])) | $(format(max_len[4],net[2][3])) |")
  end
  println("-"^(17+sum(max_len)))

  return dic
end

function read_ta_summary(network_data_file)
  @assert ispath(network_data_file)

  number_of_zones = 0
  number_of_links = 0
  number_of_nodes = 0
  first_thru_node = 0

  n = open(network_data_file, "r")

  search_sc(s,c) = something(findfirst(isequal(c), s), 0)

  while (line=readline(n)) != ""
      if occursin("<NUMBER OF ZONES>", line)
          number_of_zones = parse(Int, line[ search_sc(line, '>')+1 : end ] )
      elseif occursin("<NUMBER OF NODES>", line)
          number_of_nodes = parse(Int, line[ search_sc(line, '>')+1 : end ] )
      elseif occursin("<FIRST THRU NODE>", line)
          first_thru_node = parse(Int, line[ search_sc(line, '>')+1 : end ] )
      elseif occursin("<NUMBER OF LINKS>", line)
          number_of_links = parse(Int, line[ search_sc(line, '>')+1 : end ] )
      elseif occursin("<END OF METADATA>", line)
          break
      end
  end

  return number_of_zones, number_of_links, number_of_nodes
end


text:
module PyQHA

__precompile__() # this module is safe to precompile

using PyCall: PyNULL, PyError, PyObject, PyAny, pyimport, pycall
# Extending methods
import PyCall

const qha = PyNULL()
const qha_calculator = PyNULL()
const qha_basicio_out = PyNULL()
const qha_settings = PyNULL()
const qha_plotting = PyNULL()
const qha_input_maker = PyNULL()

@eval import Conda
Conda.pip_interop(true)
Conda.pip("install", "git+https://github.com/MineralsCloud/qha.git@v1.0.20")

function __init__()
    copy!(qha, pyimport("qha"))
    copy!(qha_calculator, pyimport("qha.calculator"))
    copy!(qha_basicio_out, pyimport("qha.basic_io.out"))
    copy!(qha_settings, pyimport("qha.settings"))
    copy!(qha_plotting, pyimport("qha.plotting"))
    copy!(qha_input_maker, pyimport("qha.basic_io.input_maker"))
    # Code from https://github.com/JuliaPy/PyPlot.jl/blob/caf7f89/src/init.jl#L168-L173
    vers = qha.__version__
    global version = try
        VersionNumber(vers)
    catch
        v"0.0.0" # fallback
    end
end

macro pyinterface(T)
    T = esc(T)
    return quote
        # Code from https://github.com/JuliaPy/PyPlot.jl/blob/6b38c75/src/PyPlot.jl#L54-L62
        PyCall.PyObject(f::$T) = getfield(f, :o)
        Base.convert(::Type{$T}, o::PyObject) = $T(o)
        Base.:(==)(f::$T, g::$T) = PyObject(f) == PyObject(g)
        Base.:(==)(f::$T, g::PyObject) = PyObject(f) == g
        Base.:(==)(f::PyObject, g::$T) = f == PyObject(g)
        Base.hash(f::$T) = hash(PyObject(f))
        PyCall.pycall(f::$T, args...; kws...) = pycall(PyObject(f), args...; kws...)
        (f::$T)(args...; kws...) = pycall(PyObject(f), PyAny, args...; kws...)
        Base.Docs.doc(f::$T) = Base.Docs.doc(PyObject(f))
        # Code from https://github.com/JuliaPy/PyPlot.jl/blob/6b38c75/src/PyPlot.jl#L65-L71
        Base.getproperty(f::$T, s::Symbol) = getproperty(PyObject(f), s)
        Base.getproperty(f::$T, s::AbstractString) = getproperty(PyObject(f), s)
        Base.setproperty!(f::$T, s::Symbol, x) = setproperty!(PyObject(f), s, x)
        Base.setproperty!(f::$T, s::AbstractString, x) = setproperty!(PyObject(f), s, x)
        PyCall.hasproperty(f::$T, s::Symbol) = hasproperty(PyObject(f), s)
        Base.propertynames(f::$T) = propertynames(PyObject(f))
        Base.haskey(f::$T, x) = haskey(PyObject(f), x)
    end
end

include("Settings.jl")
include("Calculators.jl")
include("io.jl")
include("run.jl")
include("plot.jl")
include("convert.jl")

end


text:
function Attitude(arg0::AbsoluteDate, arg1::Frame, arg2::AngularCoordinates)
    return Attitude((AbsoluteDate, Frame, AngularCoordinates), arg0, arg1, arg2)
end

function Attitude(arg0::AbsoluteDate, arg1::Frame, arg2::Rotation, arg3::Vector3D, arg4::Vector3D)
    return Attitude((AbsoluteDate, Frame, Rotation, Vector3D, Vector3D), arg0, arg1, arg2, arg3, arg4)
end

function Attitude(arg0::Frame, arg1::TimeStampedAngularCoordinates)
    return Attitude((Frame, TimeStampedAngularCoordinates), arg0, arg1)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_date(obj::Attitude)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function get_orientation(obj::Attitude)
    return jcall(obj, "getOrientation", TimeStampedAngularCoordinates, ())
end

function get_reference_frame(obj::Attitude)
    return jcall(obj, "getReferenceFrame", Frame, ())
end

function get_rotation(obj::Attitude)
    return jcall(obj, "getRotation", Rotation, ())
end

function get_rotation_acceleration(obj::Attitude)
    return jcall(obj, "getRotationAcceleration", Vector3D, ())
end

function get_spin(obj::Attitude)
    return jcall(obj, "getSpin", Vector3D, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function interpolate(obj::Attitude, arg0::AbsoluteDate, arg1::Stream)
    return jcall(obj, "interpolate", Attitude, (AbsoluteDate, Stream), arg0, arg1)
end

function interpolate(obj::TimeInterpolable, arg0::AbsoluteDate, arg1::Collection)
    return jcall(obj, "interpolate", TimeInterpolable, (AbsoluteDate, Collection), arg0, arg1)
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

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

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

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end

function with_reference_frame(obj::Attitude, arg0::Frame)
    return jcall(obj, "withReferenceFrame", Attitude, (Frame,), arg0)
end


text:
export AdaGrad

"""
    AdaGrad

AdaGrad Optimizer.

# References
* Duchi t al., ["Adaptive Subgradient Methods for Online Learning and Stochastic Optimization"](http://www.jmlr.org/papers/volume12/duchi11a/duchi11a.pdf), JMLR 2011.
"""
mutable struct AdaGrad
    alpha::Float64
    states::IdDict
end

AdaGrad(alpha::Float64) = AdaGrad(alpha, IdDict())

function (opt::AdaGrad)(value::Array{T}, grad::Array{T}) where T
    state = get!(opt.states, value, nothing)
    if state == nothing
        sqgrad = zeros(T, length(value))
        opt.states[value] = sqgrad
    else
        sqgrad = state::Array{T}
    end
    @inbounds for i = 1:length(grad)
        sqgrad[i] += grad[i] * grad[i]
        if abs(sqgrad[i]) > T(1e-8)
            value[i] -= T(opt.alpha) * grad[i] / sqrt(sqgrad[i])
        end
    end
    fill!(grad, T(0.0))
end
(opt::AdaGrad)(x::Var) = opt(x.data, x.grad)