text:
include("utilities.jl")
import ..Tokens
import ..Tokens: Token, Kind, TokenError
import ..Tokens:  NAME

mutable struct Lexing
    io::IO
    io_startpos::Int

    token_start_row::Int
    token_start_col::Int

    prevpos::Int
    token_startpos::Int

    current_row::Int
    current_col::Int
    current_pos::Int

    last_token::Tokens.Kind
end

Lexing(io::IO) = Lexing(io, position(io), 1, 1, -1, 0, 1, 1, position(io), Tokens.ERROR)
Lexing(str::AbstractString) = Lexing(IOBuffer(str))

"""
    Tokensgraphql(x::String)

Given an origin string, this returns the set of tokens that make up that source.
"""
Tokensgraphql(x::String) = collect(Lexing(x))

caracter_ignored(c::Char) = Base.isspace(c) #'\n   \t'
# Iterator interface
Base.IteratorSize(::Type{Lexing}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{Lexing}) = Base.HasEltype()
Base.eltype(::Type{Lexing}) = Token

function Base.iterate(l::Lexing)
    seek(l.io, l.io_startpos)
    l.token_startpos = position(l)
    l.token_start_row = 1
    l.token_start_col = 1

    l.current_row = 1
    l.current_col = 1
    l.current_pos = l.io_startpos
    t = next_token(l)
    return t, t.kind == Tokens.ENDMARKER
end

function Base.iterate(l::Lexing, isdone::Any)
    isdone && return nothing
    t = next_token(l)
    return t, t.kind == Tokens.ENDMARKER
end

Base.show(io::IO, l::Lexing) = print(io, typeof(l), " at position: ", position(l))

"""
    startpos(l::Lexing)

Return the latest `Token`'s starting position.
"""
startpos(l::Lexing) = l.token_startpos

"""
    prevpos(l::Lexing)

Return the lexer's previous position.
"""
prevpos(l::Lexing) = l.prevpos

"""
    position(l::Lexing)

Returns the current position.
"""
Base.position(l::Lexing) = Base.position(l.io)

"""
    eof(l::Lexing)

Determine whether the end of the lexer's underlying buffer has been reached.
"""
eof(l::Lexing) = eof(l.io)

Base.seek(l::Lexing, pos) = seek(l.io, pos)

"""
    readchar(l::Lexing)

Returns the next character and increments the current position.
"""
function readchar(l::Lexing)
    l.prevpos = position(l)
    c = readchar(l.io)
    return c
end

"""
    emit(l::Lexing, kind::Kind, str::String,
                       err::TokenError=Tokens.NO_ERR)

Returns a `Token` of kind `kind` with contents `str` and starts a new `Token`.
"""
function emit(l::Lexing, kind::Kind, str::Union{String,Int64,Float64}, err::TokenError = Tokens.NO_ERR)
   tok= Token(kind, (l.token_start_row, l.token_start_col),
                (l.current_row, l.current_col),
                startpos(l), position(l) - 1,
                str, err)
    l.last_token = kind
    l.current_col += 1
    return tok
end

"""
    emit_error(l::Lexing, str::String, err::TokenError=Tokens.UNKNOWN)

Throws Diana.GraphQLError if a lexical error is found.
"""
#Returns an `ERROR` token with error `err` and starts a new `Token`.
function emit_error(l::Lexing, str::String, err::TokenError = Tokens.UNKNOWN)
    return throw(GraphQLError("{\"errors\":[{\"locations\": [{\"column\": $(l.current_col),\"line\": $(l.current_row)}],\"message\": \"Syntax Error GraphQL request ($(l.current_row):$(l.current_col)) Unexpected character $(str) \"}]}"))
end

"""
    emit_error(l::Lexing, str::Char, err::TokenError=Tokens.UNKNOWN)

Throws Diana.GraphQLError if a lexical error is found.
"""
function emit_error(l::Lexing, str::Char, err::TokenError = Tokens.UNKNOWN)
    if caracter_ignored(str)
       return throw(GraphQLError("{\"errors\":[{\"locations\": [{\"column\": $(l.current_col),\"line\": $(l.current_row)}],\"message\": \"Syntax Error GraphQL request ($(l.current_row):$(l.current_col)) Unexpected character $(escape_string(string(str))) \"}]}"))
    end

    return throw(GraphQLError("{\"errors\":[{\"locations\": [{\"column\": $(l.current_col),\"line\": $(l.current_row)}],\"message\": \"Syntax Error GraphQL request ($(l.current_row):$(l.current_col)) Unexpected character $(str) \"}]}"))
end

"""
    next_token(l::Lexing)

Returns the next `Token`.
"""
function next_token(l::Lexing)
   c= readchar(l)

   if c == ','
      l.current_col += 1
      c= readchar(l)
   end

    while c == '#'
      while c!='\n'
        l.current_col += 1
        c= readchar(l)
      end
      l.current_row += 1
      l.current_col = 1
     c= readchar(l)
   end

    while caracter_ignored(c)
        if c == '\n'
            l.current_row += 1
            l.current_col = 1
        elseif c== '\t'
           l.current_col += 5
        else
           l.current_col += 1
         end
        c= readchar(l)

        while c == '#'
          while c!='\n'
            l.current_col += 1
            c= readchar(l)
          end
          l.current_row += 1
          l.current_col = 1
         c= readchar(l)
       end
    end

    l.token_startpos = position(l)
    l.token_start_row = l.current_row
    if l.current_col == 1
        l.token_start_col = 1
    else
        l.token_start_col = l.current_col
    end

    if eof(c); return emit(l,Tokens.ENDMARKER,"eof")
    elseif c == '['; return emit(l, Tokens.LSQUARE,"[")
    elseif c == ']'; return emit(l, Tokens.RSQUARE,"]")
    elseif c == '{'; return emit(l, Tokens.LBRACE,"{")
    elseif c == '}'; return emit(l, Tokens.RBRACE,"}")
    elseif c == '('; return emit(l, Tokens.LPAREN,"(")
    elseif c == ')'; return emit(l, Tokens.RPAREN,")")
    elseif c == '|'; return emit(l, Tokens.PIPE,"|")
    elseif c == '@'; return emit(l, Tokens.AT,"@")
    elseif c == '$'; return emit(l, Tokens.DOLLAR,"\$")
    elseif c == '='; return emit(l, Tokens.EQUALS,"=")
    elseif c == '!'; return emit(l, Tokens.BANG,"!")
    elseif c == ':'; return emit(l, Tokens.COLON,":")
    elseif c == '&'; return emit(l, Tokens.AMP,"&")
    elseif c == '"'; return lex_quote(l,c);
    elseif is_identifier_start_char(c); return lex_identifier(l, c)
    elseif c == '.'; return lex_dot(l);
    elseif c == '-'; return lex_digit(l,c);
    elseif isdigit(c); return lex_digit(l,c)
    else return emit_error(l,string(c))
    end

end

function lex_quote(l::Lexing, c::Char)
    multiline= false
    s= ""
    #s*=c
    c= readchar(l)
    l.current_col += 1

    if c=='"'
        c= readchar(l)
        if c=='"'
         c= readchar(l)
         multiline= true
        else
         return emit_error(l,string(c))
        end
    end

    while c!='"'
        if multiline==true

            if c == '\n'
                l.current_row += 1
                l.current_col = 1
            elseif c== '\t'
               l.current_col += 4
            end

            s*=c
            l.current_col += 1
            c= readchar(l)
        else
            if c=='\n'
                return emit_error(l,c)
            else
                s*=c
                l.current_col += 1
                c= readchar(l)
            end
        end
    end
    #s*="\""
    if (multiline==true)
        c= readchar(l)
        if c=='"'
            l.current_col += 1
            c= readchar(l)
            if c=='"'
                l.current_col += 1
            else
                return emit_error(l,string(c))
            end
        else
            return emit_error(l,string(c))
        end
    end

    return emit(l, Tokens.STRING, s)
end

function lex_identifier(l::Lexing, c::Char)
    s= ""
    while is_identifier_char(c)
        s*=c
        c= readchar(l)
        l.current_col += 1
    end
    l.current_col -= 1
    seek(l, prevpos(l))
    return emit(l, NAME,s)
end

function lex_dot(l::Lexing)
    if readchar(l) == '.'
        if readchar(l) == '.'
            return emit(l, Tokens.SPREAD, "...")
        end
    end
    return emit_error(l," ")
end

# A digit has been consumed
function lex_digit(l::Lexing, c::Char)
    s=""
    isfloat=false
    s*=c
    c= readchar(l)
    l.current_col += 1

    while isdigit(c)
        s*=c
        c= readchar(l)
        l.current_col += 1
    end

    if c=='.'
        s*=c
        c= readchar(l)
        l.current_col += 1
        isfloat=true
        if isdigit(c)
            s*=c
            c= readchar(l)
            l.current_col += 1
            while isdigit(c)
                s*=c
                c= readchar(l)
                l.current_col += 1
            end
        else
            return emit_error(l,c)
        end
    end

    if ((c=='e') | (c== 'E'))
        s*=c
        c= readchar(l)
        l.current_col += 1
        isfloat=true
        if ((c=='+') | (c=='-'))
            s*=c
            c= readchar(l)
            l.current_col += 1
        end
        if isdigit(c)
            s*=c
            c= readchar(l)
            l.current_col += 1
            while isdigit(c)
                s*=c
                c= readchar(l)
                l.current_col += 1
            end
        else
            return emit_error(l,c)
        end
    end

    l.current_col -= 1
    seek(l, prevpos(l))

    if isfloat==true
        return emit(l, Tokens.FLOAT,parse(Float64, s))
    end
    return emit(l, Tokens.INT, parse(Int64, s))
end


text:
function kendall(;
  α :: Float64 = 2.0,
  μ :: Float64 = 1.0,
  ν :: Float64 = 0.5)
  m = Network("Kendall's Process")

  m <= Species("X", 5)

  m <= Reaction("birth", α, "X --> X + X")
  m <= Reaction("death", μ, "X --> 0")
  m <= Reaction("immigration", ν, "0 --> X")

  return m
end

model = kendall()

ext: jl


text:
module NonconvexSemidefinite

export SDPBarrierOptions, SDPBarrierAlg, decompress_symmetric

using ChainRulesCore, Reexport, Parameters
@reexport using NonconvexCore
using NonconvexCore: @params, VecModel, AbstractOptimizer, AbstractModel, set_objective
import NonconvexCore: Workspace, optimize!, _optimize_precheck

# Semidefinite programming

# Decompress
function lowertriangind(mat::Matrix)
    indices = [i for i in CartesianIndices(mat) if i[1]>i[2]]
    return LinearIndices(mat)[indices]
end

function rearrange_x(x_L::AbstractVector, x_D::AbstractVector)
    mat_dim = length(x_D)
    L = zeros(mat_dim, mat_dim)
    L[lowertriangind(L)] .= x_L
    D = diagm(x_D)
    return (L, D)
end

function decompress_symmetric(L::Matrix, D::Matrix)
    L + D + L'
end

"""
    decompress_symmetric
    For example:
    a 3*3 positive semidefinite matrix:
    `[  a       b       d;
        b       c       e;
        d       e       f       ]`
    represents by:
    `[  x_D[1]  x_L[3]  x_L[2];
        x_L[1]  x_D[2]  x_L[1];
        x_L[2]  x_L[3]  x_D[3]  ]`
- `x_L::AbstractArray`: representing lower triangular part of a `n*n` matrix, length should be `(n^2-n)÷2`
- `x_D::AbstractArray`: representing diagonal part of a `n*n` matrix, length should be `n`
"""
function decompress_symmetric(x_L::AbstractArray, x_D::AbstractArray)
    L, D = rearrange_x(x_L, x_D)
    return decompress_symmetric(L, D)
end

function ChainRulesCore.rrule(::typeof(rearrange_x), x_L::AbstractVector, x_D::AbstractVector)
    function pullback((ΔL, ΔD))
        Δx_L = ΔL[lowertriangind(ΔL)]
        Δx_D = diag(ΔD)
        NoTangent(), Δx_L, Δx_D
    end
    return rearrange_x(x_L, x_D), pullback
end

"""
    SDPBarrierAlg: a meta-algorithm used for Semidefinite barrier programming
- `sub_alg`: algorithm to solve underlying optimization objective, which should be in consistent with `sub_options` of `SDPBarrierOptions`
"""
struct SDPBarrierAlg{Alg <: AbstractOptimizer} <: AbstractOptimizer
    sub_alg::Alg
end
function SDPBarrierAlg(;sub_alg)
    return SDPBarrierAlg(sub_alg)
end

function _optimize_precheck(model::AbstractModel, optimizer::SDPBarrierAlg, args...; kwargs...)
    nothing
end

"""
    SDPBarrierOptions: options for sdp programming
- `c_init`: initial multiplier `c` for barrier, could be a real number or vector. Vector input is for the case that objective contains multiple matrix and each element of vector for corresponding one of them
- `c_decr`: decreasing rate that multiplies `c` in every iteration. Could be either a real number of a vector, same as `c_init`
- `n_iter`: how many iterations in total for barrier method
- `sub_options`: options to solve underlying optimization target
- `keep_all`: if set to `True`, `SDPBarrierResult` keeps all results in each iteration but not only optima
- `sub_alg`: algorithm to solve underlying optimization objective, which should be in consistent with `sub_options` in `SDPBarrierOptions`
"""
@params mutable struct SDPBarrierOptions
    # Dimension of objective matrix
    # Hyperparameters
    # Initial value of `c` in barrier method:
    # `Real` for using same value for all `sd_constraints`, `AbstractArray` for assign them respectively
    c_init::Union{Real, AbstractArray}
    # Decrease rate of `c` for every epoch, same as above
    c_decr::Union{Real, AbstractArray}
    n_iter::Int
    # sub_option to solve (in)equality constraints
    sub_options
    # Keep all results or not
    keep_all::Bool
end
function SDPBarrierOptions(c_init, c_decr, n_iter; sub_options, keep_all=true)
    @assert all(0 .< c_decr .< 1) "c_decr should be between 0 and 1. "
    @assert all(c_init .> 0) "c_init shoule be larger than 0. "
    SDPBarrierOptions(c_init, c_decr, n_iter, sub_options, keep_all)
end
function SDPBarrierOptions(;sub_options, c_init=1.0, c_decr=0.1, n_iter=10, keep_all=true)
    SDPBarrierOptions(c_init, c_decr, n_iter, sub_options=sub_options, keep_all=keep_all)
end

# Result
@params struct SDPBarrierResult
    minimum
    minimizer
    results
    optimal_ind
end

@params struct SDPBarrierWorkspace <: Workspace
    model::VecModel
    x0::AbstractVector
    options::SDPBarrierOptions
    sub_alg::AbstractOptimizer
end

function Workspace(model::VecModel, optimizer::SDPBarrierAlg, x0, args...; options, kwargs...,)
    @unpack c_init, c_decr = options
    for c in model.sd_constraints.fs
        @assert isposdef(c(x0)) "Initial matrix should be positive definite. "
    end
    if c_init isa AbstractArray
        @assert length(model.sd_constraints.fs) == length(c_init) "c_init should be same length with number of `sd_constraints` when using array. "
    end
    if c_decr isa AbstractArray
        @assert length(model.sd_constraints.fs) == length(c_decr) "c_decr should be same length with number of `sd_constraints` when using array. "
    end
    if c_init isa AbstractArray && c_decr isa AbstractArray
        @assert length(c_init) == length(c_decr) "c_decr should be same length with c_init. "
    end
    return SDPBarrierWorkspace(model, copy(x0), options, optimizer.sub_alg)
end

function sd_objective(objective0, sd_constraints, c::AbstractArray)
    function _objective(args)
        target = objective0(args)
        barrier = sum(c .* -logdet.(map(f -> f(args), sd_constraints.fs)))
        return target + barrier
    end
    return _objective
end

function to_barrier(model, c::AbstractArray)
    sd_constraints, objective0 = model.sd_constraints, model.objective
    _model = set_objective(model, sd_objective(objective0, sd_constraints, c))
    return _model
end

function optimize!(workspace::SDPBarrierWorkspace)
    @unpack model, x0, options, sub_alg = workspace
    @unpack c_init, c_decr, n_iter, sub_options, keep_all = options
    objective0 = model.objective
    x = copy(x0)
    c = c_init isa Real ? ([c_init for _ in 1:length(model.sd_constraints.fs)]) : c_init
    results = []
    for _ in 1:n_iter
        model_i = to_barrier(model, c)
        result_i = optimize(model_i, sub_alg, x, options = sub_options)
        minimizer_i = result_i.minimizer
        push!(results, (objective0(minimizer_i), minimizer_i))
        c = c .* c_decr
        x = copy(minimizer_i)
    end
    optimal_ind = argmin(first.(results))
    minimum, minimizer = results[optimal_ind]
    if keep_all
        return SDPBarrierResult(minimum, minimizer, results, optimal_ind)
    else
        return SDPBarrierResult(minimum, minimizer, nothing, nothing)
    end
end

end


text:
# define main types used in this module

"""
  [`AbstractParamTpe`](@ref) for Navier-Stokes physics.
"""
mutable struct ParamType{Tdim, Tsol, Tres, Tmsh} <: AbstractParamType{Tdim}
  # need this in case we ever want to call node level functions in the Euler
  # module
  euler_params::EulerEquationMod.ParamType{Tdim, :conservative, Tsol, Tres, Tmsh}
  # some of these fields alias those in euler_params so the outside world
  # see a consistent view
  f::BufferedIO
  t::Float64  # current time value (for unsteady simulations)
  order::Int  # accuracy of elements (p=1,2,3...)
  time::Timings
  isViscous::Bool  # can this ever be false?
  penalty_relaxation::Float64
  const_tii::Float64
  Re::Float64  # free stream Reynolds number

  function ParamType{Tdim, Tsol, Tres, Tmsh}(euler_eqn::EulerEquationMod.EulerData, mesh, sbp, opts, order::Integer) where {Tdim, Tsol, Tres, Tmsh}

    # incomplete initialization
    params = new{Tdim, Tsol, Tres, Tmsh}()

    params.euler_params = euler_eqn.params
    params.f = params.euler_params.f  # alias
    params.t = 0.0
    params.order = order
    params.time = params.euler_params.time  # alias

    params.penalty_relaxation = opts["Cip"]
    params.isViscous = opts["isViscous"]

    if params.isViscous
      params.const_tii = calcTraceInverseInequalityConst(sbp, mesh.sbpface)
    end

    params.Re = opts["Re"]
    return params
  end

end  # end ParamType definition

"""
  Useful alias for 2D ParamType
"""
const ParamType2 = ParamType{2}

"""
  Useful alias for 3D ParamType
"""
const ParamType3 = ParamType{3}

"""
  Datatype for holding all the data required to solve the Naver-Stokes equations

  In particular, it contains the field `euler_eqn`, which is an
  [`EulerData`](@ref).  This object is used to evaluate the inviscid terms
  and generally allows the `NavierStokesMod` to build on top of
  `EulerEquationMod`.

  Note that many of the corresponding fields between `eqn_euler` and
  the `NSData` are aliased, including:

    * `q`
    * `q_vec`
    * `res`
    * `res_vec`
    * `res_edge`
    * `Minv`
    * `M`
"""
mutable struct NSData_{Tsol, Tres, Tdim, Tmsh} <: NSData{Tsol, Tres, Tdim}

  params::ParamType{Tdim, Tsol, Tres, Tmsh}
  comm::MPI.Comm
  commsize::Int
  myrank::Int

  # the nested EulerData object
  euler_eqn::EulerData_{Tsol, Tres, Tdim, Tmsh, :conservative}

  # required arrays
  q::Array{Tsol,3}        # holds conservative variables for all nodes
  q_vec::Array{Tres,1}    # initial condition in vector form
  res::Array{Tres, 3}     # element-based residual
  res_vec::Array{Tres, 1} # vector form of residual (aliases res in DG case)
  res_edge::Array{Tres, 4}       # edge based residual (not used)

  Minv3D::Array{Float64, 3}       # inverse mass matrix for application to res, not res_vec
  Minv::Array{Float64, 1}         # inverse mass matrix
  M::Array{Float64, 1}            # mass matrix

  # required functions
  multiplyA0inv::Function         # multiply an array by inv(A0), where A0
                                  # is the coefficient matrix of the time derivative
  majorIterationCallback::Function # called before every major (Newton/RK) itr

  # physics-specific stuff
  src_func::SRCType  # functor for the source term
  viscous_flux_func::FluxType  # functor for the viscous flux numerical flux function

  q_face::Array{Tsol, 4}     # q interpolated to face
  q_bndry::Array{Tsol, 3}    # q interpolated to boundary
  bndryflux::Array{Tres, 3}  # inviscid-type boundary flux
  flux_face::Array{Tres, 3}  # inviscid-type face flux
  area_sum::Array{Tmsh, 1}		  # the wet area of each element
  # vecflux_face::Array{Tres, 4}    # stores (u+ - u-)nx*, (numDofs, numNodes, numFaces)
  vecflux_faceL::Array{Tres, 4}     # stores (u+ - u-)nx*, (numDofs, numNodes, numFaces)
  vecflux_faceR::Array{Tres, 4}     # stores (u+ - u-)nx*, (numDofs, numNodes, numFaces)
  vecflux_bndry::Array{Tres, 4}     # stores (u+ - u-)nx*, (numDofs, numNodes, numFaces)

  file_dict::Dict{String, IO}  # dictionary of all files used for logging

  function NSData_{Tsol, Tres, Tdim, Tmsh}(mesh::AbstractMesh, sbp::AbstractOperator, opts) where {Tsol, Tres, Tdim, Tmsh}

    eqn = new{Tsol, Tres, Tdim, Tmsh}()  # incomplete initialization

    EulerEquationMod.checkOptions(opts)  # get any default values Euler needs
    eqn.euler_eqn = EulerEquationMod.EulerData_{Tsol, Tres, Tdim, Tmsh,
                                                :conservative}(mesh, sbp, opts)

    eqn.params = ParamType{Tdim, Tsol, Tres, Tmsh}(eqn.euler_eqn, mesh, sbp, opts, mesh.order)
    eqn.comm = eqn.euler_eqn.comm
    eqn.commsize = eqn.euler_eqn.commsize
    eqn.myrank = eqn.euler_eqn.myrank

    # alias the arrays to avoid having to update them manually
    eqn.q = eqn.euler_eqn.q
    eqn.q_vec = eqn.euler_eqn.q_vec
    eqn.res = eqn.euler_eqn.res
    eqn.res_vec = eqn.euler_eqn.res_vec
    eqn.res_edge = eqn.euler_eqn.res_edge
    eqn.q_face = eqn.euler_eqn.q_face
    eqn.q_bndry = eqn.euler_eqn.q_bndry
    eqn.bndryflux = eqn.euler_eqn.bndryflux
    eqn.flux_face = eqn.euler_eqn.flux_face
    eqn.Minv3D = eqn.euler_eqn.Minv3D
    eqn.Minv = eqn.euler_eqn.Minv
    eqn.M = eqn.euler_eqn.M
    eqn.file_dict = eqn.euler_eqn.file_dict

    # because NSData is solving in the same variables as Euler, we can reuse
    # their A0inv function
    eqn.multiplyA0inv = eqn.euler_eqn.multiplyA0inv

    eqn.majorIterationCallback = majorIterationCallback

   if eqn.params.isViscous
     numfacenodes = mesh.numNodesPerFace
     numfaces = mesh.numInterfaces
     numBndFaces = mesh.numBoundaryFaces
     numvars  = mesh.numDofPerNode
     # eqn.vecflux_face = zeros(Tsol, Tdim, numvars, numfacenodes, numfaces)
     eqn.vecflux_faceL = zeros(Tsol, Tdim, numvars, numfacenodes, numfaces)
     eqn.vecflux_faceR = zeros(Tsol, Tdim, numvars, numfacenodes, numfaces)
     eqn.vecflux_bndry = zeros(Tsol, Tdim, numvars, numfacenodes, numBndFaces)
     eqn.area_sum = zeros(Tmsh, mesh.numEl)
     calcElemSurfaceArea(mesh, sbp, eqn)
   else
     # eqn.vecflux_face  = Array{Tsol}(0, 0, 0, 0)
     eqn.vecflux_faceL = Array{Tsol}(0, 0, 0, 0)
     eqn.vecflux_faceR = Array{Tsol}(0, 0, 0, 0)
     eqn.vecflux_bndry = Array{Tsol}(0, 0, 0, 0)
     eqn.area_sum = Array{Tsol}(0)
   end

   return eqn
 end  # end constructor

end  # end NSData_ type definition

import PDESolver: copyParameters

function copyParameters(eqn_old::NSData, eqn_new::NSData)

  EulerEquationMod.copyParameters(eqn_old.euler_params, eqn_new.euler_params)
  eqn_new.params.order = eqn_old.params.order
  eqn_new.params.isViscous = eqn_old.params.isViscous
  eqn_new.params.penalty_relaxation = eqn_old.params.penalty_relaxation
  eqn_new.params.const_tii = eqn_old.params.const_tii
  eqn_new.params.Re = eqn_old.params.Re

  return nothing
end


text:
# module Ulti
include("ai.jl")

#TODO: separate UI/UX from engine - getLicit() -> licit(bemondas), etc.
#TODO: parallel console and WebSocket implementation (interactive = true/false)
#TODO: websocketd? Julia websockets?

#TODO UI: auto-40-20 - if clicked, no need to click it every time.
#TODO: separate state types by *, such as S*..., L*..., G*... sent to parseState

##############
#Engine runtime loop
##############

#RNG
SEED = 0
rng = MersenneTwister(SEED)
# rng = RandomDevice() #truly random

function getLicit(player::Player, currentLicit::Contract)
    while true
        try
            lct = chomp(readline())
            contract = parseContract(lct, player, currentLicit)

            if contract != nocontract && contract.totalvalue <= currentLicit.totalvalue
                println("Contract value must be higher than the current one $(currentLicit.totalvalue) > $(contract.totalvalue)")
            # elseif #TODO: kontrazni csak az eddigi bemondast lehet

            else
                return contract
            end
        catch e
            if isa(e, InterruptException) || isa(e, SystemError)
                throw(e)
            end
            println("unable to parse, please try again")
        end

    end
end

function mainLoop()
    println("Enter the number of players (default: 3)")
    numberOfPlayers = 3
    try
      numberOfPlayers = parse(Int64, readline())
    catch ArgumentError
    end

    println("Enter player names separated by comma. Default: Alfa, Béla, Gamma, ...")
    try
      nameStrings = readline() |> x -> chomp(x) |> x -> split(x, ',')
      assert(length(nameStrings) == numberOfPlayers)

      for (i, name) in enumerate(nameStrings)
        playerNames[Player(i)] = name
      end
    catch
        #NOP
    end

    #create sessionState
    session = sessionState(Vector{Tuple{Player, Int}}())
    firstplayer = 0x01 #elso kezd eloszor

    while true
        #Create licitState
        licitState = LicitState(numberOfPlayers, firstplayer)
        dealInitialCards!(licitState)

        #Licit phase
        while true
            print(stdout, licitState, false)
            lct = getLicit(licitState.currentPlayer, licitState.contract)
            licitState = licit!(licitState, lct)
            #check if licit is over
            #either no contract for a round (korpassz)
            #or no contract since current player made one
            if (length(licitState.contractHistory) == 0 &&
                    licitState.currentPlayer == licitState.firstPlayer)
                #korpassz, uj osztas
                break
            elseif  (length(licitState.contractHistory) > 0 &&
                    licitState.contractHistory[end].felvevo == licitState.currentPlayer)
                #create gameState
                gameState = startGamePhase(licitState)
                print(stdout, gameState, false)

                #Lejatszas
                result = lejatszas(GameState)
                # while true

                # end
                break
            end
        end

        #Update scores and history

        #Rotate firstplayer
        firstplayer = nextPlayer(firstplayer, numberOfPlayers)
    end
end

function lejatszas(g::GameState)
  while ~isempty(ps(g, g.currentPlayer).hand)
    print(stdout, g, false)
    # value, bestMoves = alfabeta(g, -1, typemin(Int), typemax(Int))
    value, bestMoves = alfabeta(g, -1, -g.contract.totalvalue, g.contract.totalvalue)
    chosenMove = dealCards!(1, [x for x in bestMoves])
    println(stdout, "\nChosen move: $chosenMove with value $value. Best moves: $bestMoves. Enter card to play or nothing to accept")
    _l, manualCard = readline() |> chomp |> parseCards
    if length(manualCard) > 0
      chosenMove = manualCard
    end
    g = newState(g, chosenMove)
  end
  gameScore = score(g)
  print(stdout, "$g\nGame over. Score: $gameScore", false)
end

##############
#Tests
##############

function testGame()
    g = GameState(
        #TODO BUG kontra parti nem jelenik meg
      Contract(0x01, p, [ContractElement(elolrol, ulti, Kontrak(), 96), ContractElement(elolrol, parti, Kontrak([KontraElement(pX, [KH])]), 48)], 144), #elolrol piros ulti (+ passz)
      (
        PlayerState(
          Player(1),
          CardSet([p9,pT,pK,pU,p7, mA,mK,m8, tA, tU, t9, z7]), #felvevo lapja talonnal
          # CardSet([pA,pT,pK,pU,p7, mA,mK,m8, tA, z7]), #felvevo lapja
          CardSet(), UInt8(0),UInt8(0), #nincs meg utese, 20 vagy 40
        ),
        PlayerState(
          Player(2),
          CardSet([pA,pF,p8, mT,mF,m7, zK, zF, zA, tF]), #2. jatekos lapja
          CardSet(), UInt8(0),UInt8(0), #nincs meg utese, 20 vagy 40
        ),
        PlayerState(
          Player(3),
          CardSet([mU,m9, z9,zT,zU,z8, tT,tK,t8,t7 ]), #3. jatekos lapja
          CardSet(), UInt8(0),UInt8(0), #nincs meg utese, 20 vagy 40
        )
      ),
      CardSet(), #pakli ures (kiosztva)
      CardSet(), #asztal ures
      CardSet(), #talon
      # CardSet([t9, tU]), #talon
      p1 #current player
    )

    # g=newState(g, rand(rng, validMoves(g)))
    # alfabeta(g, -1, typemin(Int), typemax(Int))
    alfabeta(g, -1, -g.contract.totalvalue, g.contract.totalvalue)
    # end #module
end

max_stars_repo_name: albapa/ulti

max_issues_repo_name: albapa/ulti

max_forks_repo_name: albapa/ulti


text:
Base.show(io::IO, A::AbstractStarAlgebra) =
    print(io, "*-algebra of $(object(A))")
Base.show(io::IO, ::Type{<:StarAlgebra{O,T}}) where {O,T} =
    print(io, "StarAlgebra{$O, $T, …}")

__prints_with_minus(::Any) = false
__prints_with_minus(x::Real) = x < 0
__needs_parens(::Any) = false
__needs_parens(a::AlgebraElement) = true

function _coeff_elt_print(io, c, elt)
    print(io, c, "·")
    __needs_parens(elt) && print(io, "(")
    print(io, elt)
    __needs_parens(elt) && print(io, ")")
    return
end

function Base.show(io::IO, a::AlgebraElement)
    A = parent(a)
    if iszero(a)
        T = eltype(a)
        if hasbasis(A)
            _coeff_elt_print(io, zero(T), first(basis(A)))
        else
            print(io, zero(T))
        end
    elseif hasbasis(A)
        elts = String[]
        nzeros = findall(!iszero, coeffs(a))
        for (counter, idx) in enumerate(nzeros)
            c, elt = coeffs(a)[idx], basis(A)[idx]
            if counter == 1
                _coeff_elt_print(io, c, elt)
            else
                print(io, ' ')
                __prints_with_minus(c) || print(io, '+')
                _coeff_elt_print(io, c, elt)
                counter == length(nzeros) || print(io, ' ')
            end
        end
    else
        println(io, "algebra element without defined basis")
        show(io, MIME("text/plain"), a.coeffs)
    end
end

function Base.show(io::IO, ::MIME"text/plain", mstr::TrivialMStructure)
    Tw = _istwisted(mstr)
    l = length(basis(mstr))
    Tw && print(io, "twisted ")
    print(io, "TrivialMStructure over basis with $l elements")
    return
end


text:
export centerofmass, findmid

# find center of mass in 2D image
function centerofmass(img)
    x_pos = []
    y_pos = []
    for i = 1:size(img, 3)
        tmp = img[:, :, i]
        x_data = mean(tmp, dims = 2)
        push!(x_pos, findmid(x_data))
        y_data = mean(tmp, dims = 1)
        push!(y_pos, findmid(y_data))
    end
    return x_pos, y_pos
end

# find weighted center in axis
function findmid(data)
    weighted_sum = 0
    int_sum = 0
    for i = 1:length(data)
        weighted_sum += data[i] * i
        int_sum += data[i]
    end
    return weighted_sum / int_sum
end


text:
module KnuthBendix

using ProgressMeter

export Word, Alphabet, RewritingSystem, LenLex, knuthbendix

include("abstract_words.jl")
include("words.jl")
include("bufferwords.jl")
include("alphabets.jl")
include("orderings.jl")
include("rules.jl")

include("rewriting.jl")
include("automata.jl")
include("helper_structures.jl")
include("derive_rule.jl")
include("force_confluence.jl")
include("kbs.jl")

include("parsing.jl")
end


text:
###################################################################
#  Scale related definitions & functions
###################################################################

config_scale_spec = Dict(
:type     => AbstractString,
:domain   => Vector,
:range    => Vector,
:round    => Bool,
:bandSize => Real,
:padding  => Real,
:clamp    => Bool,
:nice     => Union{Bool, AbstractString},
:exponent => Real,
:zero     => Bool
)

scale_help = """

- `type` (`AbstractString`) :  ScaleType, one of "linear", "log", "pow", "sqrt", "quantile", "quantize", "ordinal", "time", "utc".
- `domain` (`AbstractString[]`) : The domain of the scale, representing the set of data values. For quantitative data, this can take the form of a two-element array with minimum and maximum values. For ordinal/categorical data, this may be an array of valid input values. The domain may also be specified by a reference to a data source.
- `range` (`AbstractString[]`) : The range of the scale, representing the set of visual values. For numeric values, the range can take the form of a two-element array with minimum and maximum values. For ordinal or quantized data, the range may by an array of desired output values, which are mapped to elements in the specified domain. For ordinal scales only, the range can be defined using a DataRef: the range values are then drawn dynamically from a backing data set.
- `round` (`Bool`) :  If true, rounds numeric output values to integers. This can be helpful for snapping to the pixel grid.

Ordinal

- `bandSize` (`Real`) : minimum 0
- `padding` (`Real`) : Applies spacing among ordinal elements in the scale range. The actual effect depends on how the scale is configured. If the __points__ parameter is `true`, the padding value is interpreted as a multiple of the spacing between points. A reasonable value is 1.0, such that the first and last point will be offset from the minimum and maximum value by half the distance between points. Otherwise, padding is typically in the range [0, 1] and corresponds to the fraction of space in the range interval to allocate to padding. A value of 0.5 means that the range band width will be equal to the padding width. For more, see the [D3 ordinal scale documentation](https://github.com/mbostock/d3/wiki/Ordinal-Scales).

Typical

- `clamp` (`Bool`) : If true, values that exceed the data domain are clamped to either the minimum or maximum range value
- `nice` (`Union{Bool, AbstractString}`) : If specified, modifies the scale domain to use a more human-friendly value range. If specified as a true boolean, modifies the scale domain to use a more human-friendly number range (e.g., 7 instead of 6.96). If specified as a string, modifies the scale domain to use a more human-friendly value range. For time and utc scale types only, the nice value should be a string indicating the desired time interval ("second", "minute", "hour", "day", "week", "month", "year").
- `exponent` (`Real`) : Sets the exponent of the scale transformation. For pow scale types only, otherwise ignored.
- `zero` (`Bool`) : If true, ensures that a zero baseline value is included in the scale domain. This option is ignored for non-quantitative scales.
"""

########### helper function for encoding()  ########
"""
Sets scale properties for the scale argument within `encoding()` :

```julia
encoding_x(..., scale=scale(round=true, bandSize=12), ...)
```

$scale_help
"""
scale(;properties...) = _mkdict(config_scale_spec, properties)

########### config function ########
"""
Sets scale properties

`config_scale(;properties...)`

$scale_help
"""
config_scale(;properties...) = _mkvis((:config, :scale), config_scale_spec, properties)


text:
#**************************************************************************************
# Input_Batch.jl
# =============== part of the GeoEfficiency.jl package.
#
# all the input from the csv files or for batch work to the package is handled by some function here.
#
#**************************************************************************************

#------------------ consts - globals - imports ----------------------------

using .MathConstants, DelimitedFiles

isconst(@__MODULE__, :dataFolder) 	||	const dataFolder = string(@__MODULE__)
isconst(@__MODULE__, :dataDir)		||	const dataDir    = joinpath(homedir(), dataFolder)
mkpath(dataDir)

const Detectors  = "Detectors.csv";
const srcHeights = "srcHeights.csv";
const srcRhos    = "srcRhos.csv";
const srcRadii   = "srcRadii.csv";
const srcLengths = "srcLengths.csv";

@enum SrcType srcUnknown = -1 srcPoint = 0 srcLine = 1 srcDisk = 2 srcVolume = 3 srcNotPoint = 4
global srcType = srcUnknown

#------------------ typeofSrc --------------------------------------

"""

    typeofSrc()::SrcType

return the current value of the global `GeoEfficiency.srcType`.

"""
typeofSrc()::SrcType = srcType  # srcType != SrcType, but

"""

    typeofSrc(x::Int)::SrcType

set and return the value of the global `GeoEfficiency.srcType` corresponding to `x`.

  *  srcUnknown = -1 also any negative integer treated as so,
  *  srcPoint   = 0,
  *  srcLine    = 1,
  *  srcDisk    = 2,
  *  srcVolume  = 3,
  *  srcNotPoint = 4 also any greater than 4 integer treated as so.

"""
function typeofSrc(x::Int)::SrcType
   	global srcType = if x < 0
   					SrcType(-1)
				elseif x > 4
   					SrcType(4)
				else
   					SrcType(x)
				end
end #function

#------------------------ setSrcToPoint ---------------------------

"""

    setSrcToPoint()::Bool

return whether the source type is a point or not.
"""
setSrcToPoint()::Bool = srcType === srcPoint

"""

    setSrcToPoint(yes::Bool)::Bool

return whether the source type is a point or not after setting `srcType` to `srcPoint` if
`yes` = `true` else if `yes` = `false` setting it to `srcNotPoint` if it was not already
set to other non-point type (`srcDisk`, `srcLine`, `srcVolume`).

!!! note
    *  The user can use this function to change the source type any time.
    *  The source type is set the fist time asked for source.

**see also:** [`typeofSrc(::Int)`](@ref).

"""
function setSrcToPoint(yes::Bool)::Bool
   	global srcType = if yes
  						srcPoint
				elseif srcType in [srcUnknown, srcPoint]
  						srcNotPoint
				else
  						srcType
				end
   	return srcType === srcPoint
end

"""

	setSrcToPoint(prompt::AbstractString)::Bool

return whether the source type is a point or not. only prompt the user to set the source
type if it were not already set before.

**see also:** [`typeofSrc(::Int)`](@ref), [`setSrcToPoint(::Bool)`](@ref).

"""
setSrcToPoint(prompt::AbstractString)::Bool = srcType != srcUnknown ?	setSrcToPoint() :
											setSrcToPoint(input(prompt) |> lowercase != "n")

#---------------------------- detector_info_from_csvFile ------------------------------

"""# UnExported

	 detector_info_from_csvFile(detectors::AbstractString = Detectors,
                                      datadir::AbstractString = dataDir)
return a vector{Detector} based on information in the file of name `detectors` found in the
directory `datadir`.

!!! note
    *  if no path is given the second argument `datadir` is default to `$(dataDir)` as set by
       the constant `dataDir`.
    *  if no file name is specified the name of the predefined file `$Detectors` as set by
       the constant `Detectors`.
    *  the no argument method is the most useful; other methods are mainly for `test propose`.

"""
function detector_info_from_csvFile(detectors::AbstractString = Detectors,
                                      datadir::AbstractString = dataDir)
    detector_info_array::Matrix{Float64} = Matrix{Float64}(undef, 0, 0)
    @info("opening '$(detectors)'......", _file=nothing)
    try
        detector_info_array = readdlm(joinpath(datadir, detectors), ',', header = true)[1];
        return getDetectors(detector_info_array)

    catch err
        if isa(err, SystemError)
      		    @error("detector_info_from_csvFile: Some thing went wrong, may be the file '$(joinpath(datadir, detectors))' can't be found")
      		end
        rethrow()

    end #try

end #function

#---------------------------- read_from_csvFile --------------------------------

"""# UnExported

	read_from_csvFile(csv_data::AbstractString,
                       datadir::AbstractString = dataDir)::Vector{Float64}

return Vector{Float64} based on data in csv file named `csv_data`. directory `datadir` point to
where the file is located default to ``$(dataDir)`` as set by the constant `dataDir`.

"""
function read_from_csvFile(csv_data::AbstractString, datadir::AbstractString = dataDir)::Vector{Float64}
   	@info("Opening `$(csv_data)`......", _file=nothing)
   	try
      		indata = readdlm(joinpath(datadir, csv_data), ',',  header = true)[1][:,1]
      		return float(indata) |> sort;

   	catch err
   	    if isa(err, SystemError)
      		    @error("Some thing went wrong, may be `$(csv_data)` can't be found in `$(datadir)`", _file=nothing)

      		else
      		    @error("Some thing went wrong, may be `$(csv_data)` in `$(datadir)` format is bad or empty", _file=nothing)

      		end
      		return Float64[0.0]

   	end #try
end #function

#--------------------------- read_batch_info ------------------------------------

"""# UnExported

	read_batch_info()

read `detectors` and `sources` parameters from the predefined csv files.

Return a tuple
	   (detectors_array,
		srcHeights_array,
		srcRhos_array,
		srcRadii_array,
		srcLengths_array,
		GeoEfficiency_isPoint)

"""
read_batch_info() = read_batch_info(dataDir,
                                  Detectors,
								 srcHeights,
								    srcRhos,
								   srcRadii,
								 srcLengths)

"""# UnExported

	read_batch_info(datadir::AbstractString,
                  detectors::AbstractString,
                 srcHeights::AbstractString,
                    srcRhos::AbstractString,
                   srcRadii::AbstractString,
                 srcLengths::AbstractString)

read `detectors` and `sources` parameters from the location given in the argument list.

Return a tuple

	   (detectors_array,
		srcHeights_array,
		srcRhos_array,
		srcRadii_array,
		srcLengths_array,
		isPoint)

"""
function read_batch_info(datadir::AbstractString,
                       detectors::AbstractString,
					  srcHeights::AbstractString,
					     srcRhos::AbstractString,
					    srcRadii::AbstractString,
					  srcLengths::AbstractString)

   	@info("Starting the Batch Mode ....", _file=nothing)
   	isPoint = setSrcToPoint("\n Is it a point source {Y|n} ?")

   	@info("Reading data from `CSV files` at $datadir .....", _file=nothing)
   	detectors_array::Vector{Detector} = try
								detector_info_from_csvFile(detectors, datadir)
				catch err
								getDetectors()
			 end
   	srcHeights_array::Vector{Float64} = read_from_csvFile(srcHeights, datadir)
   	srcRhos_array::Vector{Float64} = [0.0]
   	srcRadii_array::Vector{Float64} = [0.0]
   	srcLengths_array::Vector{Float64} = [0.0]

   	function batchfailure(err::AbstractString)
      		@warn(err, ", transfer to direct data input via the `console`......", _file=nothing)
      		sleep(3); src = source()
      		srcHeights_array, srcRhos_array, srcRadii_array, srcLengths_array   =
		[src[1].Height], [src[1].Rho], [src[2]], [src[3]]
      		return nothing
   	end #function

   	if srcHeights_array == [0.0]
      		batchfailure("`$(srcHeights)` is not found in `$(datadir)`)")

   	elseif isPoint
      		srcRhos_array =	read_from_csvFile(srcRhos, datadir)

   	else
      		srcRadii_array = read_from_csvFile(srcRadii, datadir)
      		if srcRadii_array == [0.0]
         			batchfailure("`$(srcRadii)` is not found in `$(datadir)`)")

      		else
         			srcLengths_array = read_from_csvFile(srcLengths, datadir)

      		end #if
   	end #if
	#println("\n Results log\n=============")
   	return (detectors_array,
		srcHeights_array,
		srcRhos_array,
		srcRadii_array,
		srcLengths_array,
		isPoint,)
end #function

#------------------------- getDetectors -------------------------------------

"""

    getDetectors(detectors_array::Vector{<:Detector} = Detector[])::Vector{Detector}

return the `detectors_array` as Vector{Detector} extended by the entered detectors and sorted according to the
detector volume.
prompt the user to input detector parameters from the `console`.

!!! note
    If no array received in the input an empty array will be created to receive the converted detectors.

"""
function getDetectors(detectors_array::Vector{<:Detector} = Detector[])::Vector{Detector}
   	Vector{Detector}(detectors_array); @info("Please, input the detector information via the console", _file=nothing)
   	while(true)
      		try
         			push!(detectors_array, Detector());

      		catch err
         			println(err); @warn("Please: Enter a New Detector", _file=nothing)
         			continue
      		end #try

      		lowercase(input("""\n
    	                - To add a new detector press return\n
    	                - To quit press 'q'|'Q' then return\n
			\n\t your Choice: """, :blue))  == "q" && break
   	end #while
   	return detectors_array |> sort
end #function

"""

	getDetectors(detector_info_array::Matrix{<:Real},
					 detectors_array::Vector{<:Detector} = Detector[];
					 						   console_FB=true)::Vector{Detector}

return `detectors_array` as Vector{Detector}, after extending it with the successfully converted detectors. while,
attempt to convert detectors from the information in `detector_info_array`.

!!! note
    if `console_FB` argument is set to true , the function will call `getDetectors()` to take input
    from the `console` if the `detector_info_array` is empty or contain no numerical element.

"""
function getDetectors(detector_info_array::Matrix{<:Real},
						  detectors_array::Vector{<:Detector} = Detector[];
						  							console_FB = true)::Vector{Detector}

   	if isempty(detector_info_array)
      		if console_FB
         			@info("The new detectors information may entered via the console", _file=nothing)
         			return getDetectors(detectors_array)
      		else
        		 	error("getDetectors: Empty `detector_info_array`", _file=nothing)
      		end

   	else
      		Vector{Detector}(detectors_array)
      		for i_th_line = 1:size(detector_info_array)[1]
         			try
            				push!(detectors_array, Detector((detector_info_array[i_th_line,:])...))

         			catch err
         			end #try
      		end #for

      		return detectors_array |> sort
   	end   #if
end #function


text:
@testset "linear algebra" begin
    @test issquare(reshape(rand(Double32,3*3), 3, 3))

    # issue #77
    M = Double64.([1 0; 1 1])
    @test_throws ErrorException exp(M)
    M = Complex{Double64}.([1 0; 1 1])
    @test_throws ErrorException exp(M)

    t=[Complex{Double64}(1.0,0.0) Complex{Double64}(0.0,1.0);
       Complex{Double64}(0.0,-1.0) Complex{Double64}(1.0,0.0)]
    @test ishermitian(t)

    @test norm(Complex{Double64}[2, 2im]) ≈ norm(ComplexF64[2, 2im])
    @test norm(Complex{Double64}[2, 2im]) isa Double64

    @test norm(Complex{Double64}[2, 2im], 3.0) ≈ norm(ComplexF64[2, 2im], 3.0)
    @test norm(Complex{Double64}[2, 2im], 3.0) isa Double64

    @test norm(Double64[2, 1]) ≈ norm(Float64[2, 1])
    @test norm(Double64[2, 1], 3.0) ≈ norm(Float64[2, 1], 3.0)

    # issue #105
    for T in [Double16, Double32, Double64]
        for p in [0, 1.0, 2.0, Inf, -Inf]
            @test norm(T[], p) == zero(T)
            @test norm(Complex{T}[], p) == zero(T)

            @test normalize(T[]) == T[]
            @test normalize(Complex{T}[]) == Complex{T}[]
        end
    end
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 3f2a7e7c-7e03-11eb-3e84-316469e7d063
begin
	using HDF5
	using Plots
	using PlutoUI
	using LaTeXStrings
	using Printf
	import Pda
end

# ╔═╡ 4713e678-7e03-11eb-3b02-bb688d4e9390
file = "./U5.000tauR1.000.h5"

# ╔═╡ 682d1eba-7e03-11eb-304f-e9733f0a44ce
tauR =1

# ╔═╡ 53f76e3c-7e03-11eb-187c-55c93b0ead58
begin
	plot(Pda.get_h5data(file,"x/t")/tauR,
		Pda.get_h5data(file, "x/m"),
		framestyle=:box,
		label=nothing)
	xlabel!(L"t/\tau_R")
	ylabel!(L"<x>")
end

# ╔═╡ 8714fc80-7e03-11eb-03f2-cd28a1f1af3a
begin
	plot(Pda.get_h5data(file,"x/t")/tauR,
		Pda.get_h5data(file, "x/v"),
		framestyle=:box,
		label=nothing,
		xscale=:log10,
		yscale=:log10)
	xlabel!(L"t/\tau_R")
	ylabel!("variance")
end

# ╔═╡ 0b0af80a-7e04-11eb-2351-090a0530f61b
 @bind frame PlutoUI.Slider(10:10:100, show_value = true)

# ╔═╡ 099e3036-7e04-11eb-36bc-fb3772b3578d
begin
	t = Pda.get_h5data(file, "config/$(frame)/t")
	y = Pda.get_h5data(file, "config/$(frame)/x")
	loc, n = Pda.density1d(y, bins=100, normalize=true)
end

# ╔═╡ 49247ae4-7e04-11eb-13f5-cbb82ebbd241
begin
	scatter(loc[2:end-1],n[2:end-1], framestyle=:box, label=nothing)#,yscale=:log10)
	xlabel!(L"y/H")
	ylabel!("number density")
	time_label = @sprintf "t/\\tau_R=%.2f" t/tauR
	title!(latexstring(time_label))
end

# ╔═╡ Cell order:
# ╠═3f2a7e7c-7e03-11eb-3e84-316469e7d063
# ╠═4713e678-7e03-11eb-3b02-bb688d4e9390
# ╠═682d1eba-7e03-11eb-304f-e9733f0a44ce
# ╠═53f76e3c-7e03-11eb-187c-55c93b0ead58
# ╠═8714fc80-7e03-11eb-03f2-cd28a1f1af3a
# ╠═0b0af80a-7e04-11eb-2351-090a0530f61b
# ╠═099e3036-7e04-11eb-36bc-fb3772b3578d
# ╠═49247ae4-7e04-11eb-13f5-cbb82ebbd241


text:
module colour

using ..lib.nc
using ..util
using ..util.string

using Color, Color.names
using Compat

colour_config = Dict(
:FOREGROUND => :whitesmoke,
:BACKGROUND => :black,
:HIGHLIGHT => :yellow,
:BOLD => :white,
:LINKS => :blue,
:SELECTION => :blue,
:SELECTED_TEXT => :brightblue,
:BADGE => :brightblue,
:BORDER => :grey,
:CURSOR => :grey,
:CURSOR_TEXT => :darkgrey,
:DEBUG => :darkgray,
:INFO => :white,
:WARN => :yellow,
:EXCEPTION => :yellow,
:ERROR => :red,
:PROMPT => :white,
:ANSI_BLACK => :black,
:ANSI_BRIGHT_BLACK => :black,
:ANSI_RED => :red,
:ANSI_BRIGHT_RED => :red,
:ANSI_GREEN => :green,
:ANSI_BRIGHT_GREEN => :green,
:ANSI_YELLOW => :yellow,
:ANSI_BRIGHT_YELLOW  => :yellow,
:ANSI_BLUE => :blue,
:ANSI_BRIGHT_BLUE => :blue,
:ANSI_MAGENTA => :magenta,
:ANSI_BRIGHT_MAGENTA => :magenta,
:ANSI_CYAN => :cyan,
:ANSI_BRIGHT_CYAN  => :cyan,
:ANSI_WHITE=> :white,
:ANSI_BRIGHT_WHITE => :white
)

type TColor
    name::Symbol
    color_name::ASCIIString
    red::Int
    green::Int
    blue::Int
    idx::Int
    TColor(name::ASCIIString, index::Int, r::Int, g::Int, b::Int) =
    new(string2symbol(name), name, r, g, b, index)
end

const colorNames = Dict{Symbol, TColor}()
const colorIndexes = Dict{Int, TColor}()

function rgb2nc(r::Int, g::Int, b::Int)
    local red, green, blue
    local rat = 1000/255
    red = r*rat
    green = g*rat
    blue = b*rat
    int(red), int(green), int(blue)
end

function load_color_symbols()

    idx = 1
    for color_name in keys(Color.color_names)
        r, g, b = Color.color_names["$color_name"]
        r, g, b = rgb2nc(r,g,b)
        col = TColor("$color_name", idx, r, g, b)
        setindex!(colorNames, col, col.name)
        setindex!(colorIndexes, col, idx)
        idx=idx+1
    end

    for color in keys(colour_config)
        # println("$colorNames[colour_config[color]]")
        # println(colorNames[colour_config[color]])
    end

end

function color(color_symbol::Symbol)
    colorNames[color_symbol].idx
end

function color(idx::Int)
    colorIndexes[idx].name
end

    # function getcolor(name::ASCIIString)
        # return color_names[name]
    # end

    function init_nc_colors()
        for idx in keys(colorIndexes)
            init_color(idx, colorIndexes[idx].red, colorIndexes[idx].green, colorIndexes[idx].blue)
        end
    end

    function __init__()
        load_color_symbols()
    end

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************

function add_load_balance_constraints(m, p; _n="")

	##Constraint (8a): Electrical Load Balancing with Grid
    if isempty(p.etariff.export_bins)
        conrefs = @constraint(m, [ts in p.time_steps_with_grid],
            sum(p.production_factor[t, ts] * p.levelization_factor[t] * m[Symbol("dvRatedProduction"*_n)][t,ts] for t in p.elec_techs) +
            sum( m[Symbol("dvDischargeFromStorage"*_n)][b,ts] for b in p.storage.types ) +
            m[Symbol("dvGridPurchase"*_n)][ts] ==
            sum( sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for b in p.storage.types)
                + m[Symbol("dvCurtail"*_n)][t, ts] for t in p.elec_techs)
            + sum(m[Symbol("dvGridToStorage"*_n)][b, ts] for b in p.storage.types)
            + p.elec_load.loads_kw[ts]
        )
    else
        conrefs = @constraint(m, [ts in p.time_steps_with_grid],
            sum(p.production_factor[t, ts] * p.levelization_factor[t] * m[Symbol("dvRatedProduction"*_n)][t,ts] for t in p.elec_techs) +
            sum( m[Symbol("dvDischargeFromStorage"*_n)][b,ts] for b in p.storage.types ) +
            m[Symbol("dvGridPurchase"*_n)][ts] ==
            sum(  sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for b in p.storage.types)
                + sum(m[Symbol("dvProductionToGrid"*_n)][t, u, ts] for u in p.export_bins_by_tech[t])
                + m[Symbol("dvCurtail"*_n)][t, ts]
            for t in p.elec_techs)
            + sum(m[Symbol("dvGridToStorage"*_n)][b, ts] for b in p.storage.types)
            + p.elec_load.loads_kw[ts]
        )
    end

	for (i, cr) in enumerate(conrefs)
		JuMP.set_name(cr, "con_load_balance"*_n*string("_t", i))
	end

	##Constraint (8b): Electrical Load Balancing without Grid
	@constraint(m, [ts in p.time_steps_without_grid],
		sum(p.production_factor[t,ts] * p.levelization_factor[t] * m[Symbol("dvRatedProduction"*_n)][t,ts] for t in p.elec_techs) +
		sum( m[Symbol("dvDischargeFromStorage"*_n)][b,ts] for b in p.storage.types )  ==
        sum( sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for b in p.storage.types) +
        m[Symbol("dvCurtail"*_n)][t, ts] for t in p.elec_techs) +
        p.elec_load.critical_loads_kw[ts]
	)
end

function add_production_constraints(m, p; _n="")
	# Constraint (4d): Electrical production sent to storage or export must be less than technology's rated production
    if isempty(p.etariff.export_bins)
        @constraint(m, [t in p.elec_techs, ts in p.time_steps_with_grid],
            sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for b in p.storage.types)
          + m[Symbol("dvCurtail"*_n)][t, ts]
         <= p.production_factor[t, ts] * p.levelization_factor[t] * m[Symbol("dvRatedProduction"*_n)][t, ts]
        )
    else
        @constraint(m, [t in p.elec_techs, ts in p.time_steps_with_grid],
            sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for b in p.storage.types)
          + m[Symbol("dvCurtail"*_n)][t, ts]
          + sum(m[Symbol("dvProductionToGrid"*_n)][t, u, ts] for u in p.export_bins_by_tech[t])
         <= p.production_factor[t, ts] * p.levelization_factor[t] * m[Symbol("dvRatedProduction"*_n)][t, ts]
        )
    end

	# Constraint (4e): Electrical production sent to storage or grid must be less than technology's rated production - no grid
	@constraint(m, [t in p.elec_techs, ts in p.time_steps_without_grid],
		sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for b in p.storage.types)  <=
		p.production_factor[t, ts] * p.levelization_factor[t] * m[Symbol("dvRatedProduction"*_n)][t, ts]
	)

end

max_stars_repo_name: NLaws/REoptLite

max_issues_repo_name: NLaws/REoptLite

max_forks_repo_name: NLaws/REoptLite


text:
module LightningBreakerBox

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Lightning Breaker Box" => Ahorn.EntityPlacement(
        Maple.LightningBreakerBox
    )
)

sprite = "objects/breakerBox/Idle00"

function Ahorn.selection(entity::Maple.LightningBreakerBox)
    x, y = Ahorn.position(entity)
    scaleX = get(entity, "flipX", false) ? -1 : 1

    return Ahorn.getSpriteRectangle(sprite, x, y, sx=scaleX, jx=0.25, jy=0.25)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.LightningBreakerBox, room::Maple.Room)
    scaleX = get(entity, "flipX", false) ? -1 : 1

    Ahorn.drawSprite(ctx, sprite, 0, 0, sx=scaleX, jx=0.25, jy=0.25)
end

end


text:
module CoreExtended

export debug

global _debug = false
debug(msg::String) = if _debug info(msg) end

# -------------------------------------------------------------

export EMPTY_FUNCTION
const EMPTY_FUNCTION = () -> nothing

export iscallable

iscallable(f) = !isempty(methods(f))

# -------------------------------------------------------------

export invoke
export stabilize
export execute

execute(m::Module, s::Symbol, args...) = isdefined(m,s) ? execute(:($m.$s),args...) : nothing #(if isdefined(m,s) return execute(eval(m,s),args...); end; nothing)
#execute(o::Any, s::Symbol, args...) = isdefined(o,s) ? execute(:($o.$s),args...) : nothing
execute(e::Expr, args...) = execute(catchException(()->eval(e)),args...)
execute(f::Function, args...) = iscallable(f) ? invoke(f,args...) : nothing
#execute(f::Function, args...) = (debug("execute function("*string(args...)*")"); if iscallable(f) return invoke(f, args...); end; nothing)
execute(t::Tuple{Bool,Any}, args...) = execute(t[1]?t[2]:nothing,args...)
execute(r::Any, args...) = (debug("result "*string(typeof(r))); r)
execute(r::Void, args...) = (warn("Cannot execute nothing"); r)

function invoke(f::Function, args...)
	result = nothing
	catchException(()	-> result = f(args...)) #@eval
	result
end

stabilize(f::Function) = (args...) -> invoke(f, args...)

# -------------------------------------------------------------

export AbstractObjectReference
export EMPTY_OBJECT

abstract AbstractObjectReference

type EmptyObject <: AbstractObjectReference
	EmptyObject() = new()
end

const EMPTY_OBJECT = EmptyObject()

# -------------------------------------------------------------

export createSortedDict

import DataStructures
using DataStructures.SortedDict
export SortedDict

SortedDict(K, D, Ord=Base.Order.ForwardOrdering) = SortedDict{K,D,Ord}
SortedDict(D) = SortedDict{Symbol,D,Base.Order.ForwardOrdering}

createSortedDict(K, D, Ord=Base.Order.ForwardOrdering) = SortedDict(K,D,Ord)()
createSortedDict(D) = SortedDict(D)()
# -------------------------------------------------------------

export catchException
export linkToException

#Base.show_backtrace(STDOUT,backtrace())

OnException = (x)->nothing

function linkToException(f::Function)
	global OnException = f
end

function backTraceException(ex::Exception)
	println("--- [ BACKTRACE ] ---")
	Base.showerror(STDERR, ex, catch_backtrace())
	println("\n---------------------")
end

function catchException(f::Function, exf=OnException)
	try return f()
	catch ex exf(ex)
	end
	nothing
end

linkToException(backTraceException)

# -------------------------------------------------------------

export hasVal
export replace
export update

function hasVal(a::AbstractArray, getindex::Function)
	i=0; for v in a
		i+=1
		if getindex(v) return (i,v) end
	end
	#found=find(f,a)
	(0,nothing)
end

function replace(a::AbstractArray, f::Function, v::Any)
	found=hasVal(f,a)
	if found[1] == 0 push!(a, v)
	else a[found[1]]=v
	end
end

#found=find(x -> x.name == p.name,program.properties)
#!haskey() push!(program.properties, p)

function update(a::AbstractArray, getindex::Function, f::Function)
	found=hasVal(getindex,a)
	if found[1] == 0 push!(a, f((false,nothing)))
	else a[found[1]]=f((true,found[2]))
	end
end

function update(dict::Dict, index::Any, f::Function)
	v=nothing
	try
		v=(true,dict[index])
	catch error
		if isa(error, KeyError) v=(false,nothing) end
	end
	if v != nothing dict[index]=f(v) end
end

# -------------------------------------------------------------

export presetManager

presetManager(mod::Module,T::DataType) = presetManager(mod,T,T)

function presetManager(mod::Module,T::DataType, S::DataType)
	#mname=Base.replace(string(S),r"\..*","")
	#mod=Module(Symbol(mname))
	#typ="EMPTY_"*uppercase(Base.replace(string(S),r"[^\.]+\.",""))
	reset = isdefined(mod,:resets) ? mod.resets : ()->nothing
	link = isdefined(mod,:link) ? mod.link : (x)->nothing
	unlink = isdefined(mod,:unlinks) ? mod.unlinks : ()->nothing

	eval(mod,
		Expr(:toplevel,
			:(
				const Typ = Union{Void,$T};
				list = createSortedDict($S);

				function create(k::Symbol);
					if !haskey(list,k);
						e=$S(k);
						list[k]=e;
						setCurrent(e);
						if isdefined($mod,:init) init(e) end;
					else
						e=list[k];
						setCurrent(e);
					end;
					e;
				end;

				current = nothing;
				getCurrent() = current;
				setCurrent(obj::Typ) = global current = obj;
				#setCurrent(obj::Symbol) = setCurrent(create(obj));

				function reset();
					unlink();
					$reset();
				end;

				isInvalid(obj::Typ) = obj == nothing || !isa(obj,$S);
				isLinked(obj::Typ) =	!isInvalid(obj) && obj == getCurrent();

				function unlink() setCurrent(nothing); $unlink(); end;

				function link();
					obj=getCurrent();
					if !isInvalid(obj); $link(obj) else unlink() end;
				end;

				function linkTo(obj::Typ);
					other=getCurrent();
					if other == obj return end;
					setCurrent(obj);
					link();
				end;
			)
		)
	)
end

end # CoreExtended

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
# Function to merge eventlist and write merged event file
# Author: Gulab Dewangan
# Date : 29 October 2018
#

using DataFrames, FITSIO, FITSIO.Libcfitsio

"""
    sxt_l2evtlist_merge(l2evtfilelist, merged_evtfile)

Merge the level2 orbit-wise event files in FITS format from the AstroSat/SXT instrument.

...
# Arguments
- `l2evtfilelist::String`: An ascii file listing the level2 event files for different orbits of a given SXT observation.
- `merged_evtfile::String`: Name of the merged event file to be created.
...
"""
function sxt_l2evtlist_merge(l2evtfilelist::String, merged_evtfile::String)
   evtfiles = readlines(l2evtfilelist)
   # Read fits event files
   f=Array{FITS}(undef, length(evtfiles))
   for i=1:length(evtfiles)
     f[i]=FITS(evtfiles[i], "r+")
     write_key(f[i][2],"TFORM10","1I")
     write_key(f[i][4],"TFORM5","1I")
     close(f[i])
   end
   f=FITS.(evtfiles, "r+")
   evtlist=Array{TableHDU}(undef, length(evtfiles))
   badpixlist=Array{TableHDU}(undef,length(evtfiles))
   gtilist=Array{TableHDU}(undef,length(evtfiles))
   evt_df=Array{DataFrame}(undef, length(evtfiles))
   badpix_df=Array{DataFrame}(undef, length(evtfiles))
   gti_df=Array{DataFrame}(undef, length(evtfiles))
   for i=1:length(f)
     evtlist[i]=f[i][2]
     gtilist[i]=f[i][3]
     badpixlist[i]=f[i][4]
     evt_df[i]=sxtevtlist2df(evtlist[i])
     badpix_df[i]=sxtbadpixlist2df(badpixlist[i])
     gti_df[i]=gti2df(gtilist[i])
     write_key(f[i][2],"TFORM10","16X")
     write_key(f[i][4],"TFORM5","16X")
    # close(f)
   end
   close.(f)
   # Merge event lists and badpixel lists
   evtlist_merged=vcat(evt_df...)
   gtilist_merged = vcat(gti_df...)
   gtilist_unique = gti_unique_or(gtilist_merged)
   badpixlist_merged=vcat(badpix_df...)
   # Retain only unique events or badpixels
   sort!(evtlist_merged, [:TIME, :RAWX, :RAWY])
   unique!(evtlist_merged, [:CCDFrame,:RAWX, :RAWY])
   sort!(badpixlist_merged, [:RAWX, :RAWY])
   unique!(badpixlist_merged)

   # Convert dataframes to vectors to be able to write to fits files
   (evtlist_col_names, evtlist_col_vals) = sxtevt_df2vec(evtlist_merged)
   (badpix_col_names,badpix_col_vals) = sxtbadpix_df2vec(badpixlist_merged)
   (gti_start, gti_stop) = (gtilist_unique[:,:START], gtilist_unique[:,:STOP])
   exposure_sec = sum(gti_stop .- gti_start)

   # Create GTIs
   # evttimes = evtlist_merged[:TIME]
   # (gti_start,gti_stop) = buildgti(evttimes)
   # exptime_sec=sum(gti_stop .- gti_start)

   gti_col_names=["START","STOP"]
   gti_col_vals=[gti_start,gti_stop]

   # Read headers from the first evtfile
   # ff=FITS(evtfiles[1])
   # primary_hdr=read_header(ff[1])
   # evt_hdr=read_header(ff[2])
   # gti_hdr=read_header(ff[3])
   # badpix_header=read_header(ff[4])
   # close(ff)
   #
   # # Update header keywords for merged eventfile
   # evt_hdr["EXPOSURE"]=exptime_sec
   # evt_hdr["ONTIME"]=exptime_sec
   # evt_hdr["LIVETIME"]=exptime_sec

   # Write the merged event file

  ff=fits_create_file(merged_evtfile)
  evtlist_coldef=[("TIME", "1D", "seconds"),("CCDFrame", "1J", " "),("X", "1I", "pixel"),("Y", "1I", "pixel"),("RAWX", "1I", "pixel"),("RAWY", "1I", "pixel"),("DETX", "1I", "pixel"),("DETY", "1I", "pixel"),("PHAS", "9I", "channel"),("STATUS", "16X", " "), ("PHA", "1J", " "), ("GRADE", "1I", " "),  ("PixsAbove", "1I", " "), ("PI", "1J", "channel"), ("RA", "D", " "),("Dec", "D", " "), ("PHASO", "9I", " ")]
  fits_create_binary_tbl(ff,0,evtlist_coldef,"EVENTS")
  for i=1:length(evtlist_col_names)
    fits_write_col(ff,i,1,1,evtlist_col_vals[i])
  end

  gti_coldef=[("START", "1D", "s"),("STOP", "1D", "s")]
  fits_create_binary_tbl(ff,0,gti_coldef,"STDGTI")
  for i=1:length(gti_col_names)
    fits_write_col(ff,i,1,1,gti_col_vals[i])
  end

  badpix_coldef=[("RAWX", "1I", " "),("RAWY", "1I"," "),("TYPE", "1I", " "),("YEXTENT", "1I"," "), ("BADFLAG", "16I", " ")]
  fits_create_binary_tbl(ff,0,badpix_coldef,"BADPIX ")
  for i=1:length(badpix_col_names)
    fits_write_col(ff,i,1,1,badpix_col_vals[i])
  end
fits_close_file(ff)

# Determine start-stop keywords from the list of event files
  ef = Array{FITS}(undef, length(evtfiles))
  ef_evtheader = Array{FITSHeader}(undef, length(evtfiles))
  ef = FITS.(evtfiles, "r")
  for n=1:length(ef)
     ef_evtheader[n] = read_header(ef[n][2])
     close(ef[n])
   end

   # Determine keyword values
   tstart = Array{Float64}(undef,length(evtfiles))
   tstop = Array{Float64}(undef,length(evtfiles))
   date_obs = Array{String}(undef,length(evtfiles))
   date_end = Array{String}(undef,length(evtfiles))
   time_obs = Array{String}(undef,length(evtfiles))
   time_end = Array{String}(undef,length(evtfiles))
   mjd_obs = Array{Float64}(undef,length(evtfiles))
  for i=1:length(evtfiles)
    tstart[i] = ef_evtheader[i]["TSTART"]
    tstop[i] = ef_evtheader[i]["TSTOP"]
    date_obs[i] = ef_evtheader[i]["DATE-OBS"]
    date_end[i] = ef_evtheader[i]["DATE-END"]
    time_obs[i] = ef_evtheader[i]["TIME-OBS"]
    time_end[i] = ef_evtheader[i]["TIME-END"]
    mjd_obs[i] = ef_evtheader[i]["MJD-OBS"]
  end
  tstart_merged = minimum(tstart)
  tstop_merged = maximum(tstop)
  mjd_obs_merged = minimum(mjd_obs)
  telapse_merged = tstop_merged - tstart_merged
  date_obs_merged = date_obs[findall(x ->x==tstart_merged, tstart)][1]
  date_end_merged = date_end[findall(x ->x==tstop_merged, tstop)][1]
  time_obs_merged = time_obs[findall(x ->x==tstart_merged, tstart)][1]
  time_end_merged = time_end[findall(x ->x==tstop_merged, tstop)][1]

  # Write header keywords from the first eventfile
  ef1=FITS(readlines(l2evtfilelist)[1])
  mef=FITS(merged_evtfile,"r+")
  # ef1h1=read_header(ef1[1])
  # ef1h2=read_header(ef1[2])
  # ef1h3=read_header(ef1[3])
  # ef1h4=read_header(ef1[4])
  ef1_header=[read_header(ef1[i]) for i=1:4]
  for i=1:4
    for j=1:length(ef1_header[i].keys)
      if ef1_header[i].keys[j] in ["SIMPLE", "BITPIX", "NAXIS", "EXTEND", "NAXIS1", "NAXIS2", "PCOUNT", "GCOUNT", "TFIELDS", "XTENSION", "BITPIX", "TFORM1", "TFORM2", "TFORM3", "TFORM4", "TFORM5", "TFORM6", "TFORM7", "TFORM8",  "TFORM9", "TFORM10", "TFORM11", "TFORM12", "TFORM13", "TFORM14","TFORM15", "TFORM16", "TFORM17", "TOTAL2K", "ERROR2K", "FILIN001", "EXPOSURE"]
#        println(ef1_header[i].keys[j] * "Not updated")
      else
        write_key(mef[i],ef1_header[i].keys[j],ef1_header[i].values[j],ef1_header[i].comments[j])
#        println(ef1_header[i].keys[j] * " updated.")
      end
      write_key(mef[i],"EXPOSURE",exposure_sec,"Exposure time")
      write_key(mef[i],"LIVETIME",exposure_sec,"On source time")
      write_key(mef[i],"ONTIME",exposure_sec,"On source time")
      write_key(mef[i], "TSTART",tstart_merged,"time start")
      write_key(mef[i], "TSTOP",tstop_merged ,"time stop")
      write_key(mef[i], "DATE-OBS",date_obs_merged,"start date of observation")
      write_key(mef[i], "DATE-END",date_end_merged,"end date of observation")
      write_key(mef[i], "TIME-OBS",time_obs_merged,"start time of observation")
      write_key(mef[i], "TIME-END",time_end_merged,"end time of observation")
      write_key(mef[i], "TELAPSE", telapse_merged, "Total elapsed time")
      write_key(mef[i], "MJD-OBS", mjd_obs_merged, "MJD of observation start")
  end
end
close(ef1)
close(mef)
println("Merged event file: "  * merged_evtfile * " successfully written!")
return 1
end

max_stars_repo_name: gulabd/SXTMerger

max_issues_repo_name: gulabd/SXTMerger


text:
# Multivariate Normal distribution

export MvNormal

@parameterized MvNormal(μ,σ)

# MvNormal(;kwargs...) = MvNormal(kwargs)

@kwstruct MvNormal(μ)
@kwstruct MvNormal(σ)
@kwstruct MvNormal(ω)
@kwstruct MvNormal(μ,σ)
@kwstruct MvNormal(μ,ω)

supportdim(d::MvNormal) = supportdim(params(d))

proxy(d::MvNormal) = affine(params(d), Normal() ^ supportdim(d))
logdensity(d::MvNormal, x) = logdensity(proxy(d), x)
basemeasure(d::MvNormal) = basemeasure(proxy(d))

rand(rng::AbstractRNG, ::Type{T}, d::MvNormal) where {T} = rand(rng, T, proxy(d))

# function MvNormal(nt::NamedTuple{(:μ,)})
#     dim = size(nt.μ)
#     affine(nt, Normal() ^ dim)
# end

# function MvNormal(nt::NamedTuple{(:σ,)})
#     dim = colsize(nt.σ)
#     affine(nt, Normal() ^ dim)
# end

# function MvNormal(nt::NamedTuple{(:ω,)})
#     dim = rowsize(nt.ω)
#     affine(nt, Normal() ^ dim)
# end

# function MvNormal(nt::NamedTuple{(:μ, :σ,)})
#     dim = colsize(nt.σ)
#     affine(nt, Normal() ^ dim)
# end

# function MvNormal(nt::NamedTuple{(:μ, :ω,)})
#     dim = rowsize(nt.ω)
#     affine(nt, Normal() ^ dim)
# end


text:
function name2index(names::Vector{RString})
    n2i = Dict{RString,Int64}()
    i2n = Dict{Int64,RString}()
    @inbounds for (i,k) in enumerate(names)
        if k != ""
            n2i[k] = i
            i2n[i] = k
        end
    end
    n2i, i2n
end

"""
    Container that mimics R vector behaviour.
    Elements could be accessed either by indices as a normal vector,
    or (optionally) by string keys as a dictionary.
"""
immutable DictoVec{T}
    data::T
    name2index::Dict{RString, Int}
    index2name::Dict{Int, RString}

    @compat function (::Type{DictoVec}){T}(data::T, names::Vector{RString} = Vector{RString}())
        n2i, i2n = name2index(names)
        new{T}(data, n2i, i2n)
    end
end

Base.length(dict::DictoVec) = length(dict.data)
Base.isempty(dict::DictoVec) = isempty(dict.data)

Base.haskey(dict::DictoVec, key) = haskey(dict.name2index, key)

function Base.setindex!(dict::DictoVec, value, key)
    ix = get(dict.name2index, key, 0)
    if ix > 0
      setindex!(dict.data, value, ix)
    else
      dict.name2index[key] = length(dict.data)+1
      dict.index2name[ix] = key
      push!(dict.data, value)
    end
end

function Base.setindex!(dict::DictoVec, value, index::Int64)
    setindex!(dict.data, value, index)
end

Base.getindex(dict::DictoVec, index::Int) = getindex(dict.data, index)

function Base.getindex(dict::DictoVec, key)
    ix = get(dict.name2index, key, 0)
    if ix > 0
        return getindex(dict.data, ix)
    else
        throw(KeyError(key))
    end
end

Base.get(dict::DictoVec, index::Int, default) = get(dict.data, index, default)

function Base.get(dict::DictoVec, key, default)
    ix = get(dict.name2index, key, 0)
    return ix > 0 ? dict.data[ix] : default
end

Base.get(f::Function, dict::DictoVec, index::Int) = get(f, dict.data, index)

function Base.get(f::Function, dict::DictoVec, key)
    ix = get(dict.name2index, key, 0)
    return ix > 0 ? dict.data[ix] : f()
end

function Base.keys(dict::DictoVec)
    return keys(dict.name2index)
end

function Base.values(dict::DictoVec)
    return dict.data
end

function Base.show(io::IO, dict::DictoVec)
    if isempty(dict.name2index)
        # no keys
        show(dict.data)
    else
        first = true
        print(io, "DictoVec(")
        n = 0
        for i in eachindex(dict.data)
            first || print(io, ',')
            first = false
            k = get(dict.index2name, i, "")
            if k != ""
                show(io, k)
                print(io, "=>")
            end
            show(io, dict.data[i])
            n += 1
            # limit && n >= 10 && (print(io, "…"); break)
        end
        print(io, ")")
    end
end

function Base.convert(::Type{Dict{RString,Any}}, dv::DictoVec)
    res = Dict{RString,Any}()
    for (k,v) in dv.name2index
        res[k] = dv.data[v]
    end
    res
end


text:
module TestDistributions
using Test
using Distributions
using EvidentialDeepLearning
using SpecialFunctions: loggamma

@testset "NormalInverseGamma" begin
    @testset "logpdf" begin
        for _ in 1:100
            μ = rand(Normal(0,10))
            ν = rand(Uniform(1,10))
            α = rand(Uniform(1,10))
            β = rand(Uniform(1,10))
            nig = NIG(μ,ν,α,β)
            x = rand(Normal(μ, 5))
            σ² = rand(Uniform(0,10))
            # We have
            # x |σ²,μ,ν ~ Normal(μ, √(σ²/ν)), x)
            # σ²|α,β    ~ InverseGamma(α,β)
            @test logpdf(nig, (;x,σ²)) ≈
                logpdf(Normal(μ, √(σ²/ν)), x) + logpdf(InverseGamma(α,β), σ²)
        end
    end
end

function S26(yᵢ,γ,ν,α,β)
    # taken from [Deep Evidential Regression](https://arxiv.org/abs/1910.02600)
    Ω = 2*β*(1+ν)
    1//2*log(π/ν) - α*log(Ω) + (α + 1//2) * log((yᵢ-γ)^2*ν+Ω) + loggamma(α) - loggamma(α+1//2)
end

@testset "Deep Evidential Regression (S26)" begin
    for _ in 1:100
        yᵢ = rand(Uniform(-10,10))
        γ = rand(Uniform(-5,5))     # Normal mean
        ν = rand(Uniform(0,100))    # nobs mean
        α = rand(Uniform(1, 50))    # nobs variance
        β = rand(Uniform(0,10)) * α # halfsum of square deviations

        nig = NIG(γ,ν,α,β)
        pp = posterior_predictive(nig)
        @test nll(pp, yᵢ) ≈ S26(yᵢ,γ,ν,α,β)
    end
end

end#module


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

# Stage-decomposition optimization #
# ================================ #
function optimize!(structure::StageDecompositionStructure, optimizer::AbstractStructuredOptimizer, x₀::AbstractVector)
    # Sanity checks
    supports_structure(optimizer, structure) || throw(UnsupportedStructure{typeof(optimizer), typeof(structure)}())
    check_loadable(optimizer, structure)
    # Load structure
    load_structure!(optimizer, structure, x₀)
    # Run structure-exploiting optimization procedure
    MOI.optimize!(optimizer)
    return nothing
end

function set_master_optimizer!(structure::StageDecompositionStructure, optimizer)
    set_optimizer(structure.first_stage, optimizer)
    return nothing
end

function set_master_optimizer_attribute!(structure::StageDecompositionStructure, attr::MOI.AbstractOptimizerAttribute, value)
    MOI.set(backend(structure.first_stage), attr, value)
    return nothing
end

function set_subproblem_optimizer!(structure::StageDecompositionStructure, optimizer)
    set_optimizer!(scenarioproblems(structure), optimizer)
    return nothing
end

function set_subproblem_optimizer_attribute!(structure::StageDecompositionStructure, attr::MOI.AbstractOptimizerAttribute, value)
    MOI.set(scenarioproblems(structure), attr, value)
    return nothing
end

function cache_solution!(stochasticprogram::StochasticProgram{2}, structure::StageDecompositionStructure{2}, optimizer::MOI.AbstractOptimizer)
    cache = solutioncache(stochasticprogram)
    # Cache main solution
    variables = decision_variables_at_stage(stochasticprogram, 1)
    constraints = decision_constraints_at_stage(stochasticprogram, 1)
    cache[:solution] = SolutionCache(optimizer, variables, constraints)
    # Cache first-stage solution
    cache[:node_solution_1] = SolutionCache(backend(structure.first_stage), variables, constraints)
    # Cache scenario-dependent solutions (Skip if more than 100 scenarios for performance)
    if num_scenarios(stochasticprogram) <= 1e2
        variables = decision_variables_at_stage(stochasticprogram, 2)
        constraints = decision_constraints_at_stage(stochasticprogram, 2)
        cache_solution!(cache, scenarioproblems(structure), optimizer, 2, variables, constraints)
    end
    return nothing
end


text:
# Helper functions
include("adapt/adapt.jl")
include("support/hmc_core.jl")

# Concrete algorithm implementations.
include("hmcda.jl")
include("nuts.jl")
include("sghmc.jl")
include("sgld.jl")
include("hmc.jl")
@init @require DynamicHMC="bbc10e6e-7c05-544b-b16e-64fede858acb" @eval begin
    include("dynamichmc.jl")
end
include("mh.jl")
include("is.jl")
include("smc.jl")
include("pgibbs.jl")
include("pmmh.jl")
include("ipmcmc.jl")
include("gibbs.jl")

## Fallback functions

# utility funcs for querying sampler information
require_gradient(spl::Sampler) = false
require_particles(spl::Sampler) = false

assume(spl::Sampler, dist::Distribution) =
error("Turing.assume: unmanaged inference algorithm: $(typeof(spl))")

observe(spl::Sampler, weight::Float64) =
error("Turing.observe: unmanaged inference algorithm: $(typeof(spl))")

## Default definitions for assume, observe, when sampler = nothing.
##  Note: `A<:Union{Nothing, SampleFromPrior, HamiltonianRobustInit}` can be
##   simplified into `A<:Union{SampleFromPrior, HamiltonianRobustInit}` when
##   all `spl=nothing` is replaced with `spl=SampleFromPrior`.
function assume(spl::A,
    dist::Distribution,
    vn::VarName,
    vi::VarInfo) where {A<:Union{Nothing, SampleFromPrior, HamiltonianRobustInit}}

    if haskey(vi, vn)
        r = vi[vn]
    else
        r = isa(spl, HamiltonianRobustInit) ? init(dist) : rand(dist)
        push!(vi, vn, r, dist, 0)
    end
    # NOTE: The importance weight is not correctly computed here because
    #       r is genereated from some uniform distribution which is different from the prior
    # acclogp!(vi, logpdf_with_trans(dist, r, istrans(vi, vn)))

    r, logpdf_with_trans(dist, r, istrans(vi, vn))

end

function assume(spl::A,
    dists::Vector{T},
    vn::VarName,
    var::Any,
    vi::VarInfo) where {T<:Distribution, A<:Union{Nothing, SampleFromPrior, HamiltonianRobustInit}}

    @assert length(dists) == 1 "Turing.assume only support vectorizing i.i.d distribution"
    dist = dists[1]
    n = size(var)[end]

    vns = map(i -> copybyindex(vn, "[$i]"), 1:n)

    if haskey(vi, vns[1])
        rs = vi[vns]
    else
        rs = isa(spl, HamiltonianRobustInit) ? init(dist, n) : rand(dist, n)

        if isa(dist, UnivariateDistribution) || isa(dist, MatrixDistribution)
            for i = 1:n
                push!(vi, vns[i], rs[i], dist, 0)
            end
            @assert size(var) == size(rs) "Turing.assume: variable and random number dimension unmatched"
            var = rs
        elseif isa(dist, MultivariateDistribution)
            for i = 1:n
                push!(vi, vns[i], rs[:,i], dist, 0)
            end
            if isa(var, Vector)
                @assert length(var) == size(rs)[2] "Turing.assume: variable and random number dimension unmatched"
                for i = 1:n
                    var[i] = rs[:,i]
                end
            elseif isa(var, Matrix)
                @assert size(var) == size(rs) "Turing.assume: variable and random number dimension unmatched"
                var = rs
            else
                @error("Turing.assume: unsupported variable container"); error()
            end
        end
    end

    # acclogp!(vi, sum(logpdf_with_trans(dist, rs, istrans(vi, vns[1]))))

    var, sum(logpdf_with_trans(dist, rs, istrans(vi, vns[1])))

end

function observe(spl::A,
    dist::Distribution,
    value::Any,
    vi::VarInfo) where {A<:Union{Nothing, SampleFromPrior, HamiltonianRobustInit}}

    vi.num_produce += one(vi.num_produce)
    @debug "dist = $dist"
    @debug "value = $value"

    # acclogp!(vi, logpdf(dist, value))
    logpdf(dist, value)

end

function observe(spl::A,
    dists::Vector{T},
    value::Any,
    vi::VarInfo) where {T<:Distribution, A<:Union{Nothing, SampleFromPrior, HamiltonianRobustInit}}

    @assert length(dists) == 1 "Turing.observe only support vectorizing i.i.d distribution"
    dist = dists[1]
    @assert isa(dist, UnivariateDistribution) || isa(dist, MultivariateDistribution) "Turing.observe: vectorizing matrix distribution is not supported"
    if isa(dist, UnivariateDistribution)  # only univariate distributions support broadcast operation (logpdf.) by Distributions.jl
        # acclogp!(vi, sum(logpdf.(Ref(dist), value)))
        sum(logpdf.(Ref(dist), value))
    else
        # acclogp!(vi, sum(logpdf(dist, value)))
        sum(logpdf(dist, value))
    end

end

##############
# Utilities  #
##############

# VarInfo to Sample
@inline Sample(vi::VarInfo) = begin
  value = Dict{Symbol, Any}() # value is named here because of Sample has a field called value
  for vn in keys(vi)
    value[sym(vn)] = vi[vn]
  end

  # NOTE: do we need to check if lp is 0?
  value[:lp] = getlogp(vi)

  if ~isempty(vi.pred)
    for sym in keys(vi.pred)
      # if ~haskey(sample.value, sym)
        value[sym] = vi.pred[sym]
      # end
    end
    # TODO: check why 1. 2. cause errors
    # TODO: which one is faster?
    # 1. Using empty!
    # empty!(vi.pred)
    # 2. Reassign an enmtpy dict
    # vi.pred = Dict{Symbol,Any}()
    # 3. Do nothing?
  end

  Sample(0.0, value)
end

# VarInfo, combined with spl.info, to Sample
@inline Sample(vi::VarInfo, spl::Sampler) = begin
  s = Sample(vi)

  if haskey(spl.info, :wum)
    s.value[:epsilon] = getss(spl.info[:wum])
  end

  if haskey(spl.info, :lf_num)
    s.value[:lf_num] = spl.info[:lf_num]
  end

  if haskey(spl.info, :eval_num)
    s.value[:eval_num] = spl.info[:eval_num]
  end

  return s
end


text:
module DiracQC

import Automa
import Automa.RegExp: @re_str
const re = Automa.RegExp

using DataFrames
using LinearAlgebra

using Formatting

# https://physics.nist.gov/cgi-bin/cuu/Value?auedm
const ea₀ = 8.4783536255e-30 # e*Bohr
# https://en.wikipedia.org/wiki/Debye
const Debye = 1e-21/299792458 # C*m
const Debye_au = Debye/ea₀

function read_until(fun::Function, io::IO, pred::Function)
    l = ""
    while !eof(io)
        l = readline(io)
        pred(l) && break
        fun(l)
    end
    l
end

read_until(fun::Function, io::IO, s::String) =
    read_until(fun, io, l -> occursin(s, l))

read_until(io::IO, pred) =
    read_until(l -> nothing, io, pred)

# * Automa

dec      = re"[-+]?[0-9]+"
prefloat = re"[-+]?([0-9]+\.[0-9]*|[0-9]*\.[0-9]+)"
float    = prefloat | re.cat(prefloat | re"[-+]?[0-9]+", re"[eE][-+]?[0-9]+")

dec.actions[:enter] = [:mark]
dec.actions[:exit] = [:dec]
float.actions[:enter] = [:mark]
float.actions[:exit] = [:float]

# * Energies

orbital_energies_machine = let
    energy = re.cat(re" *", float, re" +",
                    "(", re" *", dec, ")")
    energies = re.rep(energy)

    Automa.compile(energies)
end

orbital_energies_actions = Dict(
    :mark => :(mark = p),
    :float => :(push!(energies, parse(Float64, data[mark:p-1]))),
    :dec => :(push!(gs, parse(Int, data[mark:p-1])))
)

context = Automa.CodeGenContext()
@eval function read_orbital_energies(data::String)
    energies = Vector{Float64}()
    gs = Vector{Int}()
    mark = 0

    $(Automa.generate_init_code(context, orbital_energies_machine))

    # p_end and p_eof were set to 0 and -1 in the init code,
    # we need to set them to the end of input, i.e. the length of `data`.
    p_end = p_eof = lastindex(data)

    $(Automa.generate_exec_code(context, orbital_energies_machine, orbital_energies_actions))

    # We need to make sure that we reached the accept state, else the
    # input did not parse correctly
    iszero(cs) || error("failed to parse on byte ", p)

    (energies=energies, gs=gs)
end

function read_orbital_block(io::IO, header)
    symmetry = last(split(first(split(header, ":")), " "))
    parity = last(symmetry) == 'g' ? 1 : -1
    Ω = Rational(parse.(Int, split(last(split(header, " ")), "/"))...)

    energies = Vector{Float64}()
    gs = Vector{Int}()
    occupation = Vector{Float64}()
    f = 1.0
    read_until(io, l -> isempty(strip(l))) do l
        if occursin("*", l)
            f = parse(Float64, last(split(l, "=")))
        else
            e,g = read_orbital_energies(l)
            append!(energies, e)
            append!(gs, g)
            append!(occupation, fill(f, length(e)))
        end
    end
    df = DataFrame(E=energies, g=gs, occupation=occupation)
    (symmetry=symmetry, parity=parity, Ω=Ω,
     orbitals=df)
end

function read_energies(io::IO)
    energies = (;)
    read_until(io, "TOTAL ENERGY")
    read_until(io, "Eigenvalues") do l
        if occursin(":", l)
            k,v = strip.(split(l, ":"))
            energies = merge(energies, (Symbol(replace(k, " " => "_")) => parse(Float64, v),))
        end
    end
    for i = 1:3 readline(io) end
    h = readline(io)
    blocks = []
    while occursin("* Block", h)
        push!(blocks, read_orbital_block(io, h))
        h = readline(io)
    end
    dfs = map(blocks) do b
        m = size(b.orbitals,1)
        hcat(b.orbitals, DataFrame(Symmetry=fill(b.symmetry, m),
                                   Parity=fill(b.parity, m),
                                   Ω=fill(b.Ω, m)))
    end

    energies = merge(energies, (:orbitals => blocks,
                                :orbitals_combined => reduce(vcat, dfs)))

    energies
end

# * Interface

function load(io::IO)
    (energies=read_energies(io),)
end

load(filename::String) = open(load, filename)

end


text:
using FSimPlots
using FSimZoo
using FlightSims
const FS = FlightSims
using Plots  # you should use `Plots` as well

gr()
# ENV["GKSwstype"] = "nul"
using LinearAlgebra
using DiffEqCallbacks
using Transducers
using UnPack
using ReferenceFrameRotations

function gen_gif()
    multicopter = LeeHexacopter()
    @unpack m, g = multicopter
    x0 = State(multicopter)()
    anim = Animation()
    function affect!(integrator)
        state = copy(integrator.u)
        fig = plot(multicopter, state;
                   xlim=(-1, 1), ylim=(-1, 1), zlim=(-1, 20),
                   camera=(45, 45),
                  )
        frame(anim)
    end
    Δt = 0.1
    tf = 10.0
    cb = PeriodicCallback(affect!, Δt)
    function Dynamics!(multicopter::Multicopter)
        FSimZoo.__Dynamics!(multicopter)
    end
    prob, df = sim(x0,
                   apply_inputs(Dynamics!(multicopter);
                                f=(state, p, t) -> m*g + (0.5*tf-t),
                                M=zeros(3),
                               );
                   tf=tf,
                   callback=cb,
                  )
    gif(anim, "figures/anim.gif", fps=60)
    nothing
end

function topview()
    multicopter = LeeHexacopter()
    x0 = State(multicopter)()
    fig = plot(multicopter, x0;
               ticks=nothing, border=:none,
               # background_color=:transparent,
               xlabel="", ylabel="", zlabel="",
               camera=(0, 90),
               dpi=300,
              )
    savefig(fig, "figures/topview.png")
end

function model_description()
    multicopter = LeeHexacopter()
    x0 = State(multicopter)()
    length_param = 0.5
    x0.p += [length_param, length_param, -length_param]  # NED
    ϕ, θ, ψ = 0, deg2rad(30), 0
    x0.R = ReferenceFrameRotations.angle_to_dcm(ψ, θ, ϕ, :ZYX)
    @unpack p, R = x0
    p_enu = ned2enu(p)
    # plot
    fig = plot3d(;)
    opacity_axes = 0.5
    opacity_pos_vec = 0.5
    # inertial frame
    origin_length = -length_param
    axis_length = length_param
    origin_pos = origin_length*ones(3)
    Naxis_pos = origin_pos + axis_length*ned2enu([1, 0, 0])
    Eaxis_pos = origin_pos + axis_length*ned2enu([0, 1, 0])
    Daxis_pos = origin_pos + axis_length*ned2enu([0, 0, 1])
    Naxis_1 = LinRange(origin_pos[1], Naxis_pos[1], 2)
    Naxis_2 = LinRange(origin_pos[2], Naxis_pos[2], 2)
    Naxis_3 = LinRange(origin_pos[3], Naxis_pos[3], 2)
    Eaxis_1 = LinRange(origin_pos[1], Eaxis_pos[1], 2)
    Eaxis_2 = LinRange(origin_pos[2], Eaxis_pos[2], 2)
    Eaxis_3 = LinRange(origin_pos[3], Eaxis_pos[3], 2)
    Daxis_1 = LinRange(origin_pos[1], Daxis_pos[1], 2)
    Daxis_2 = LinRange(origin_pos[2], Daxis_pos[2], 2)
    Daxis_3 = LinRange(origin_pos[3], Daxis_pos[3], 2)
    plot!(fig,
          Naxis_1, Naxis_2, Naxis_3;
          color=:red,
          label=nothing,
          opacity=opacity_axes,
         )
    plot!(fig,
          Eaxis_1, Eaxis_2, Eaxis_3;
          color=:green,
          label=nothing,
          opacity=opacity_axes,
         )
    plot!(fig,
          Daxis_1, Daxis_2, Daxis_3;
          color=:blue,
          label=nothing,
          opacity=opacity_axes,
         )
    # body frame
    Xaxis_pos = p_enu + ned2enu(R'*enu2ned(Naxis_pos-origin_pos))
    Yaxis_pos = p_enu + ned2enu(R'*enu2ned(Eaxis_pos-origin_pos))
    Zaxis_pos = p_enu + ned2enu(R'*enu2ned(Daxis_pos-origin_pos))
    Xaxis_1 = LinRange(p_enu[1], Xaxis_pos[1], 2)
    Xaxis_2 = LinRange(p_enu[2], Xaxis_pos[2], 2)
    Xaxis_3 = LinRange(p_enu[3], Xaxis_pos[3], 2)
    Yaxis_1 = LinRange(p_enu[1], Yaxis_pos[1], 2)
    Yaxis_2 = LinRange(p_enu[2], Yaxis_pos[2], 2)
    Yaxis_3 = LinRange(p_enu[3], Yaxis_pos[3], 2)
    Zaxis_1 = LinRange(p_enu[1], Zaxis_pos[1], 2)
    Zaxis_2 = LinRange(p_enu[2], Zaxis_pos[2], 2)
    Zaxis_3 = LinRange(p_enu[3], Zaxis_pos[3], 2)
    plot!(fig,
          Xaxis_1, Xaxis_2, Xaxis_3;
          color=:red,
          label=nothing,
          opacity=opacity_axes,
         )
    plot!(fig,
          Yaxis_1, Yaxis_2, Yaxis_3;
          color=:green,
          label=nothing,
          opacity=opacity_axes,
         )
    plot!(fig,
          Zaxis_1, Zaxis_2, Zaxis_3;
          color=:blue,
          label=nothing,
          opacity=opacity_axes,
         )
    # position vector
    plot!(fig,
          [origin_pos[1], p_enu[1]], [origin_pos[2], p_enu[2]], [origin_pos[3], p_enu[3]];
          color=:black,
          label=nothing,
          opacity=opacity_pos_vec,
         )
    # BE CAREFUL; figures w.r.t. ENU
    plot!(fig, multicopter, x0;
          ticks=nothing, border=:none,
          # background_color=:transparent,
          xlabel="", ylabel="", zlabel="",
          xlim=(-2*length_param, 2*length_param),
          ylim=(-2*length_param, 2*length_param),
          zlim=(-2*length_param, 2*length_param),
          camera=(45, 45),
          # dpi=300,
         )
    savefig(fig, "figures/hexacopter_description.png")
    display(fig)
end

function prob_description()
    multicopter = LeeHexacopter()
    x0 = State(multicopter)()
    Δt = 0.1
    traj_des_func = (t) -> [-0.75, -0.75, 0.75]*(1-t)^3 + (-0.5)*ones(3)*(1-t)^2*t + (-1)*ones(3)*(1-t)*t^2
    η = 0.2
    x0.p = traj_des_func(η)
    x0.R = ReferenceFrameRotations.angle_to_dcm(0, deg2rad(-10), deg2rad(30), :ZYX)
    fig = plot(;
               legend=:bottomright,
               legendfont=font(14),
              )
    plot!(fig,
          multicopter, x0;
          ticks=nothing, border=:none,
          # background_color=:transparent,
          xlabel="", ylabel="", zlabel="",
          xlim=(-1.0, 0.1), ylim=(-1.0, 0.1), zlim=(-1.0, 0.1),
          dpi=300,
         )
    plot!(fig,
          [0], [0], [0];
          st=:scatter, color=:red,
          markersize=6,
          label="set point",
         )
    traj_des_enu_before = 0:Δt:η |> Map(t -> ned2enu(traj_des_func(t))) |> collect
    traj_des_e_before = traj_des_enu_before |> Map(enu -> enu[1]) |> collect
    traj_des_n_before = traj_des_enu_before |> Map(enu -> enu[2]) |> collect
    traj_des_u_before = traj_des_enu_before |> Map(enu -> enu[3]) |> collect
    traj_des_enu_after = η:Δt:1.0 |> Map(t -> ned2enu(traj_des_func(t))) |> collect
    traj_des_e_after = traj_des_enu_after |> Map(enu -> enu[1]) |> collect
    traj_des_n_after = traj_des_enu_after |> Map(enu -> enu[2]) |> collect
    traj_des_u_after = traj_des_enu_after |> Map(enu -> enu[3]) |> collect
    plot!(fig,
          traj_des_e_before, traj_des_n_before, traj_des_u_before;
          color=:black,
          label="trajectory",
         )
    plot!(fig,
          traj_des_e_after, traj_des_n_after, traj_des_u_after;
          color=:black,
          ls=:dash,
          label="",
         )
    savefig(fig, "figures/prob_description.png")
    display(fig)
end

gen_gif()
topview()
model_description()
prob_description()


text:
module DataCreationDriver

    using CSVFiles
    using Base.Filesystem

    using CareData:
        DataSettings,
        load_all,
        create_base_datasets,
        default_data_settings
    using GlobalDecls

    export createdata

    #settings.reaching_18s_source = DFE
    #settings.num_iterations = 200
    # push!(settings.targets, "E10000030" )
    #settings.agglevel = national

    function create_data( settings :: DataSettings )
        run_data_dir = DATADIR*"/populations/$(settings.datayear)/$(settings.dataset)/"
        mkpath( run_data_dir )
        alldata = load_all( settings.datayear )
        for k in 1:settings.num_iterations
            println( "creating iteration $k")
            created = create_base_datasets(
                alldata.ofdata,
                settings );
            CSVFiles.save( run_data_dir*"/yp_data_$k.csv", created.yp_data )
            CSVFiles.save( run_data_dir*"/carer_data_$k.csv", created.carer_data )
        end
    end # create data

end # module