text:
__precompile__(true)

module VertexModels

export RandGenerator,
       Orientation,
       maxorientation,
       twentyvertex,
       picture,
       height,
       issink,
       issource,
       pushdown!,
       pushup!

mutable struct RandGenerator
    D::Dict{Int64,Tuple{Int64,Int64,Bool}}
    m::Int64
    n::Int64
end

function twentyvertex(m,n,maxiter=2^20)
    U = RandGenerator(Dict(),m,n)
    T = 1
    upper0, lower0 = maxorientation(m,n), minorientation(m,n)
    upper,lower = upper0, lower0
    while upper != lower
        upper,lower = copy(upper0), copy(lower0)
        for t = -T:-1
            pushdown!(upper,sample(U,t)...)
            pushdown!(lower,sample(U,t)...)
        end
        T = 2T
        if T >= maxiter
            error("Did not converge")
        end
    end
    return upper
end

function sample(R::RandGenerator,t::Integer)
    if !(t in keys(R.D))
        u = (rand(2:R.m-1),rand(3:R.n-2),rand(Bool))
        R.D[t] = u
        return u
    else
        return R.D[t]
    end
end

mutable struct Orientation
    D::Dict{Tuple{Tuple{Int64,Int64},Tuple{Int64,Int64}},Int64}
    m::Int64
    n::Int64
end

import Base: size, keys, getindex, setindex!, copy, ==, !=, hash

hash(o::Orientation) = hash((o.D,o.m,o.n))
==(o::Orientation,p::Orientation) = o.D == p.D
!=(o::Orientation,p::Orientation) = !(o == p)
size(o::Orientation) = (o.m,o.n)
keys(o::Orientation) = keys(o.D)
getindex(o::Orientation,I::Tuple{Tuple{Int64,Int64},Tuple{Int64,Int64}}) = o.D[I]
setindex!(o::Orientation, i::Int64, k::Tuple{Tuple{Int64,Int64},Tuple{Int64,Int64}}) = setindex!(o.D,i,k)
copy(o::Orientation) = Orientation(copy(o.D),o.m,o.n)

function sinks(h::Orientation)
    m,n = size(h)
    return [(i,j) for i=1:m,j=1:n][[issink(h,i,j) for i=1:m,j=1:n]]
end

function sources(h::Orientation)
    m,n = size(h)
    return [(i,j) for i=1:m,j=1:n][[issource(h,i,j) for i=1:m,j=1:n]]
end

function position(i::Integer,j::Integer,m::Integer,n::Integer)
    h = sqrt(3)/3
    if iseven(m-i) && isodd(j)
        return [j,(m-i)*3h]
    elseif iseven(m-i) && iseven(j)
        return [j,(m-i)*3h-h]
    elseif isodd(m-i) && isodd(j)
        return [j,2h+(m-i-1)*3h]
    else
        return [j,2h+(m-i-1)*3h+h]
    end
end

function neighbors(i::Integer,j::Integer,m::Integer,n::Integer)
    nbs = Array{Int64,1}[]
    t = iseven(j) && isodd(m-i) || isodd(j) && iseven(m-i) ? -1 : 1
    for (a,b) in [(0,1),(0,-1),(t,0)]
        if 1 ≤ i+a ≤ m && 1 ≤ j+b ≤ n
            push!(nbs,[i+a,j+b])
        end
    end
    return nbs
end

"""
function picture(filename::AbstractString,h::Orientation;kwargs...)
    m,n = size(h)
    return Graphics2D.showgraphics(filename,
    [h[e]*Graphics2D.Arrow([position(e[1]...,m,n)'; position(e[2]...,m,n)'];
     arrowsize=0.2,arrowloc=0.6) for e in keys(h)])
end

function picture(h::Orientation;show=true,kwargs...)
    m,n = size(h)
    if show
        return Graphics2D.showgraphics(
	[h[e]*Graphics2D.Arrow([position(e[1]...,m,n)';
	position(e[2]...,m,n)'];
 	arrowsize=0.2,arrowloc=0.6,kwargs...) for e in keys(h)])
    else
        return [h[e]*Graphics2D.Arrow(
	[position(e[1]...,m,n)'; position(e[2]...,m,n)'];
	arrowsize=0.2,arrowloc=0.6,kwargs...) for e in keys(h)]
    end
end
"""

function issource(O::Orientation,i::Integer,j::Integer)
    m,n = size(O)
    for (k,l) in neighbors(i,j,m,n)
        if ((i,j),(k,l)) in keys(O)
            if O[((i,j),(k,l))] == -1
                return false
            end
        elseif ((k,l),(i,j)) in keys(O)
            if O[((k,l),(i,j))] == 1
                return false
            end
        end
    end
    return true
end

function issink(O::Orientation,i::Integer,j::Integer)
    m,n = size(O)
    for (k,l) in neighbors(i,j,m,n)
        if ((i,j),(k,l)) in keys(O)
            if O[((i,j),(k,l))] == 1
                return false
            end
        elseif ((k,l),(i,j)) in keys(O)
            if O[((k,l),(i,j))] == -1
                return false
            end
        end
    end
    return true
end

function height(O::Orientation)
    m,n = size(O)
    D = Dict{Tuple{Int64,Int64},Int64}()
    D[(1,1)] = 0
    while length(D) < m*n
        for t in keys(O.D)
            if t[1] in keys(D)
                if t[2] in keys(D)
                    if D[t[2]] - D[t[1]] != O.D[t]
                        error("height function inconsistent")
                    end
                end
                D[t[2]] = D[t[1]] + O.D[t]
            elseif t[2] in keys(D)
                if t[1] in keys(D)
                    if D[t[1]] - D[t[2]] != -O.D[t]
                        error("height function inconsistent")
                    end
                end
                D[t[1]] = D[t[2]] - O.D[t]
            end
        end
    end
    return [D[i,j] for i=1:m,j=1:n]
end

"""
function picture(D::Dict{Tuple{Int64,Int64},Int64})
    m = maximum([m for (m,n) in keys(D)])
    n = maximum([n for (m,n) in keys(D)])
    grlist = Graphics2D.GraphicElement[]
    for k in keys(D)
        push!(grlist,Graphics2D.GraphicText(position(k...,m,n),string(D[k]);textsize=1.0))
    end
    return grlist
end
"""

function pushup!(h::Orientation,i::Integer,j::Integer)::Bool
    pushdown!(h,i,j,true)
end

function pushdown!(h::Orientation,i::Integer,j::Integer,up=false)::Bool
    m,n = size(h)
    pushed = false
    if up ? issource(h,i,j) : issink(h,i,j)
        pushed = true
        for v in neighbors(i,j,m,n)
            v = tuple(v...)
            if ((i,j),v) in keys(h)
                h[((i,j),v)] *= -1
            else
                h[(v,(i,j))] *= -1
            end
        end
    end
    return pushed
end

function minorientation(m::Integer,n::Integer)::Orientation
    return maxorientation(m,n,down=true)
end

function maxorientation(m::Integer,n::Integer;down=false)::Orientation
    topedges = [((1,j),(1,j+1)) for j=1:n-1]
    topedges = [isodd(i) ? (e[2],e[1]) : e for (i,e) in enumerate(topedges)]
    bottomedges = [((m,j),(m,j+1)) for j=1:n-1][end:-1:1]
    bottomedges = [isodd(i) ? (e[2],e[1]) : e for (i,e) in enumerate(bottomedges)]
    leftedges = vcat([[((i,1),(i+1,1)), ((i+1,1),(i+1,2)),
        ((i+1,2),(i+2,2)), ((i+2,2),(i+2,1))] for i=1:2:m-1]...)[2m-3:-1:1]
    leftedges = [isodd(i) ? (e[2],e[1]) : e for (i,e) in enumerate(leftedges)]
    rightedges = vcat([[((i,n),(i+1,n)), ((i+1,n),(i+1,n-1)),
        ((i+1,n-1),(i+2,n-1)), ((i+2,n-1),(i+2,n))] for i=1:2:m-1]...)[1:2m-3]
    rightedges = [isodd(i) ? (e[2],e[1]) : e for (i,e) in enumerate(rightedges)]
    boundaryedges = [topedges; rightedges; bottomedges; leftedges]
    interioredges = vcat([isodd(i+j) ? ((i,j),(i,j+1)) : ((i,j+1),(i,j)) for j=2:n-2,i=2:m-1]...)
    interioredges = [interioredges;
        vcat([((i+1,j),(i,j)) for i=1:2:m-1,j=3:2:n-2]...);
        vcat([((i+1,j),(i,j)) for i=2:2:m-1,j=4:2:n-2]...);
    ];

    alledges = [boundaryedges; interioredges];

    O = Orientation(Dict(e=>1 for e in alledges),m,n)

    for i=1:10*(m+n)
        for j=2:m-1
            for k=3:n-2
                if down
                    pushdown!(O,j,k)
                else
                    pushup!(O,j,k)
                end
            end
        end
    end
    return O
end

end # module


text:
function AmbiguitySolver(arg0::List, arg1::IntegerLeastSquareSolver, arg2::AmbiguityAcceptance)
    return AmbiguitySolver((List, IntegerLeastSquareSolver, AmbiguityAcceptance), arg0, arg1, arg2)
end

function fix_integer_ambiguities(obj::AmbiguitySolver, arg0::jint, arg1::List, arg2::RealMatrix)
    return jcall(obj, "fixIntegerAmbiguities", List, (jint, List, RealMatrix), arg0, arg1, arg2)
end

function get_all_ambiguity_drivers(obj::AmbiguitySolver)
    return jcall(obj, "getAllAmbiguityDrivers", List, ())
end

function un_fix_ambiguity(obj::AmbiguitySolver, arg0::ParameterDriver)
    return jcall(obj, "unFixAmbiguity", void, (ParameterDriver,), arg0)
end

ext: jl


text:
using DelimitedFiles

using LinearAlgebra
using DSP: conv

function derivative(signal, order, halflen, dt)
    @assert 2*halflen ≥ order ≥ 0
    @assert iseven(halflen)
    S = (-halflen:halflen)' .^ (0:order);
    F = qr(S')
    filter = Matrix(F.Q)/Matrix(F.R')
    der = conv(signal, 1/-dt*filter[:, 2])
    return der[halflen+1:end-halflen]
end

function bayeslinreg(dΩ, Ω, U, σ, μβ, invΣβ)
    p = Ω[:,1]; q = Ω[:,2]; r = Ω[:,3]
    ψ = hcat(p.^2, q.^2, r.^2, p.*q, p.*r, q.*r, U)'
    invA = inv(1/σ^2*ψ*ψ' + invΣβ)
    μ(y) = 1/σ^2*invA*(ψ*y + invΣβ*μβ)
    μ.(eachcol(dΩ)), invA
end

σ = 1.0
μβ = zeros(9)
invΣβ = inv(diagm(0=>ones(9)))
Ω = readdlm("/tmp/omega.txt")
U = readdlm("/tmp/controls.txt")
dΩ = hcat(derivative.(eachcol(Ω), 4, 30, 0.01)...)

using MAT
file = matopen("/home/arun/Documents/shared-scripts/learning.mat", "w")
write(file, "cflie", hcat(dΩ, Ω, U))
close(file)

# slice some inbetween data to train
#  Ω = Ω[251:500,:]
#  U = U[251:500,:]
#  dΩ = dΩ[251:500,:]

postμ, postσ = bayeslinreg(dΩ, Ω, U, σ, μβ, invΣβ)

## for testing only
function predictivedist(Ωx, Ux, postμ, postσ, σ)
    p, q, r = Ωx
    ψx = [p^2, q^2, r^2, p*q, p*r, q*r, Ux...]
    μ = Ref(ψx').*postμ
    σ = ψx'*postσ*ψx .+ σ^2
    μ, σ
end
pred = predictivedist.(eachrow(Ω), eachrow(U), Ref(postμ), Ref(postσ), Ref(σ))
predx = [p[2][2] for p in pred]
plot(dΩ[:,2]); plot!(predx)


text:
include("tokens.jl")
include("data_structures.jl")
grammar_map = Dict{Symbol,Int}()
grammar = Production[]

function addProduction(id::Symbol, body_::Production)
  body = body_
  if !haskey(grammar_map,id)
    #if !isdefined(id)
      body.enum =length(grammar)+1
      body.lexem = string(id)
      push!(grammar, body)
      grammar_map[id] = length(grammar)
    #else
    # error(string("The element ",id," is already defined"))
    #end
  else
    error("Production already exists")
  end
end

function addProduction(id::Symbol,body_)
    addProduction(id,Production(body_))
end
function addProd(id::Symbol,body_)
    addProduction(id,Production(body_))
end

function getProd_idx(s::Union{Symbol,Int})
   try
    return grammar[grammar_map[s]].enum
   catch

    return false
   end
end

function getProd(s::Union{Symbol,Int})
   try
    return grammar[grammar_map[s]]
   catch
    if isinteger(s)
      return grammar[s]
    else
      error("Wrong type input for getProd")
     end
   end
end

#=
   normal:
    S -> TYPE 'id' PARAM '{' ALL_ITER '}' S
      -> EPSILON
=#
addProduction(:S, [[:IDT_TYPE_OR_VOID,ID, :PARAM, O_C_BRCKT, :ALL_ITER, C_C_BRCKT, :S],
                   [EOF],
                   [EPSILON]])

#=
   normal:
     TYPE  -> CONST 'int'
           -> CONST 'char'
           -> CONST 'float'
           -> CONST 'string'

     CONST -> 'const'
           -> EPSILON

   fectoring:
     TYPE  -> CONST IDT_TYPE

     CONST -> 'const'
           -> EPSILON

     IDT_TYPE    -> 'int'
           -> 'char'
           -> 'float'
           -> 'string'
=#
addProduction(:TYPE, [[:CONST_R, :IDT_TYPE]])
addProduction(:CONST_R, [[CONST],
                         [EPSILON]])
addProduction(:IDT_TYPE, [[IDT_INT],
                    [IDT_CHAR],
                    [IDT_FLOAT],
                    [IDT_STRING],
                    [IDT_BOOL]])

addProduction(:IDT_TYPE_OR_VOID, [[:TYPE], [IDT_VOID]])

#=
   normal:
     PARAM  -> '(' ')'
            -> '(' PR ')'

     PR     -> TYPE ID_OR_VEC
            -> TYPE ID_OR_VEC ',' PR
   factoring to remove ambiguity:
     PARAM  -> '(' PRH

     PRH    -> ')'
            -> PR ')'

     PR     -> TYPE ID_OR_VEC PRL

     PRL    -> ',' PR
            -> EPSILON
=#
addProduction(:PARAM, [[O_BRCKT, :PRH]])
addProduction(:PRH, [[C_BRCKT],
                     [:PR, C_BRCKT]])
addProduction(:PR, [[:TYPE, :ID_OR_VEC, :PRL]])
addProduction(:PRL, [[COMMA, :PR],
                     [EPSILON]])

#=
   normal:
     ATTR  -> 'IDT_INT' ID_OR_VEC OPR_ATR EXPR_NUM
           -> 'IDT_STRING' OPR_ATR EXPR_STRING
           -> 'IDT_CHAR' ID_OR_VEC OPR_ATR EXPR_NUM
           -> 'IDT_FLOAT' ID_OR_VEC OPR_ATR EXPR_NUM

     ID_OR_VEC -> ID
           -> ID, VEC_IN, :EXPR_NUM

   factoring to remove ambiguity:

     ATTR  -> 'IDT_INT' ID_OR_VEC OPR_ATR EXPR_NUM
           -> 'IDT_STRING' OPR_ATR EXPR_STRING
           -> 'IDT_CHAR' ID_OR_VEC OPR_ATR EXPR_NUM
           -> 'IDT_FLOAT' ID_OR_VEC OPR_ATR EXPR_NUM

     ID_OR_VEC -> 'id' IDT
     IDT   -> '::' EXPR_NUM | EPSILON
=#

#=
addProduction(:ATTR, [[IDT_INT, :ID_OR_VEC, OPR_ATR, :EXPR_NUM],
                      [IDT_STRING, :ID_OR_VEC, OPR_ATR, :EXPR_STRING],
                      [IDT_FLOAT, :ID_OR_VEC, OPR_ATR, :EXPR_NUM],
                      [IDT_BOOL, :ID_OR_VEC, OPR_ATR, :EXPR_BOOL],
                      [IDT_CHAR, :ID_OR_VEC, OPR_ATR, :EXPR_NUM]])
=#

addProduction(:CT_BOOL, [[CT_TRUE],
                         [CT_FALSE]])

addProduction(:CT, [[CT_STRING],
                    [CT_FLOAT],
                    [CT_INT],
                    [:CT_BOOL]])

addProduction(:FCALL_OR_ATRIB, [[:ID_OR_VEC, :FCALL_OR_ATRIB_F]])

addProduction(:ATRIB, [[OPR_ATR, :EXPR_NUM]] )

addProduction(:FCALL_OR_ATRIB_F, [[:ATRIB],
                                   [:FN_CH]])

addProduction(:VAR_DCLR, [[:IDT_TYPE, :ID_OR_VEC, OPR_ATR, :DATA]])

addProduction(:DATA, [[:CT],
                      [:ID_OR_VEC]])

addProduction(:ID_OR_VEC, [[ID, :VEC_INDX]])

addProduction(:INT_OR_ID, [[ID],
                           [CT_INT]])

addProduction(:VEC_INDX, [[VEC_IN, :VEC_INDX_F],
                     [EPSILON]])

addProduction(:VEC_INDX_F, [[O_BRCKT, :EXPR_NUM,C_BRCKT],
                            [:INT_OR_ID]])

#=
normal:
   EXPR_STRING -> EXPR_STRING '+' EXPR_STRING
   EXPR_STRING -> (EXPR_STRING)
   EXPR_STRING -> 'ct_string'
   EXPR_STRING -> 'id'

factoring:
  EXPR_STRING -> ESH '+' EXPR_STRING

  ESH      -> 'ct_string'
           -> 'id'
           -> (EXPR_STRING)

addProduction(:EXPR_STRING, [[:ESH, OPR_PM, :EXPR_STRING]])
addProduction(:ESH, [[CT_STRING],
                     [:FN_H_STR],
                     [O_BRCKT, :EXPR_STRING, C_BRCKT]])

addProduction(:FN_H_STR, [[ID, :FN_H_V_STR]])
addProduction(:FN_H_STR_V, [[VEC_IN, :EXPR_NUM],
                            [:FN_CH],
                            [EPSILON]])

   normal:
     EXPR_NUM  -> EXPR_NUM + EXPR_NUM
               -> EXPR_NUM * EXPR_NUM
               -> '-' EXPR_NUM
               -> '(' EXPR_NUM ')'
               -> 'ct_int'
               -> 'ct_float'
   precedence:
     EXPR_NUM -> EXPR_NUM_K
              -> EXPR_NUM_K + EXPR_NUM

     EXPR_NUM_K        -> EXPR_NUM_G
              -> EXPR_NUM_K * EXPR_NUM_G

     EXPR_NUM_G        -> 'ct_int'
              -> 'ct_float'
              -> '(' EXPR_NUM ')'
   removed left recursion:
     EXPR_NUM -> EXPR_NUM_K
              -> EXPR_NUM_K '+' EXPR_NUM

     EXPR_NUM_K        -> EXPR_NUM_G EXPR_NUM_KH

     EXPR_NUM_KH       -> '*' EXPR_NUM_G EXPR_NUM_KH
              -> EPSILON

     EXPR_NUM_G        -> 'ct_int'
              -> 'ct_float'
              -> '(' EXPR_NUM ')'
   ambiguity:
     EXPR_NUM -> EXPR_NUM_K EXPR_NUM_KR
     EXPR_NUM_KR       -> '+' EXPR_NUM
              -> EPSILON

     EXPR_NUM_K        -> EXPR_NUM_G EXPR_NUM_KH

     EXPR_NUM_KH       -> '*' EXPR_NUM_G EXPR_NUM_KH
              -> EPSILON

     EXPR_NUM_G        -> 'ct_int'
              -> 'ct_float'
              -> '(' EXPR_NUM ')'
=#

#=
addProduction(:NUM_OPRS, [[OPR_PM],
                          [OPR_DM]]
                           )

addProduction(:EXPR_NUM, [[O_C_BRCKT,:EXPR_NUM,C_C_BRCKT],
                          [:NUM_OPRS, :EXPR_NUM],
                          [:DATA, :EXPR_NUM],
                          [EPSILON]])

addProduction(:EXPR_NUM_R, [[:DATA, :NUM_OPRS, :EXPR_NUM_R]
                            [EPSILON]])
=#

addProduction(:EXPR_NUM, [[:EXPR_NUM_K, :EXPR_NUM_KR]])
addProduction(:EXPR_NUM_KR, [[OPR_PM ,:EXPR_NUM],
                    [EPSILON]])
addProduction(:EXPR_NUM_K, [[:EXPR_NUM_G,:EXPR_NUM_KH]])
addProduction(:EXPR_NUM_KH, [[OPR_DM,:EXPR_NUM_G, :EXPR_NUM_KH],
                    [EPSILON]])
addProduction(:EXPR_NUM_G, [[CT_FLOAT],
                   [CT_INT],
                   [:ID_OR_VEC],
                   [O_BRCKT, :EXPR_NUM, C_BRCKT]])

#addProduction(:FN_H_NUM, [[ID, :FN_H_V_NUM]])
#addProduction(:FN_H_NUM_V, [[VEC_IN, :EXPR_NUM],
#                            [:FN_CH],
#                            [EPSILON]])

#=
   normal:
     EXPR_BOOL -> EXPR_BOOL 'or' EXPR_BOOL
               -> EXPR_BOOL 'and' EXPR_BOOL
               -> 'not' EXPR_BOOL
               -> '(' EXPR_BOOL ')'
               -> EXPR_NUM 'rel' EXPR_NUM
               -> 'true'
               -> 'false'
   precedence:
     EXPR_BOOL -> EXPR_BOOL_T
               -> EXPR_BOOL 'or' EXPR_BOOL_T
               -> 'not' EXPR_BOOL_T

     EXPR_BOOL_T         -> EXPR_BOOL_F
               -> EXPR_BOOL_T 'and' EXPR_BOOL_F

     EXPR_BOOL_F         -> '(' EXPR_BOOL ')'
               -> EXPR_NUM 'rel' EXPR_NUM
               -> 'true'
               -> 'false'

   removed left recursion:
     EXPR_BOOL  -> EXPR_BOOL_T EXPR_BOOL_H
                -> 'not' EXPR_BOOL_T EXPR_BOOL_H

     EXPR_BOOL_H -> 'or' EXPR_BOOL_T EXPR_BOOL_H
                -> EPSILON

     EXPR_BOOL_T          -> EXPR_BOOL_F EXPR_BOOL_TH

     EXPR_BOOL_TH         -> 'and'EXPR_BOOL_F EXPR_BOOL_TH
                -> EPSILON

     EXPR_BOOL_F          -> '(' EXPR_BOOL ')'
                -> EXPR_NUM OPRLR_REL EXPR_NUM
                -> 'true'
                -> 'false'
=#
addProduction(:EXPR_BOOL, [[:EXPR_BOOL_T, :EXPR_BOOL_H],
                           [OPRL_NOT, :EXPR_BOOL_T, :EXPR_BOOL_H]])
addProduction(:EXPR_BOOL_H, [[OPRLR_OR, :EXPR_BOOL_T, :EXPR_BOOL_H],
                            [EPSILON]])

addProduction(:EXPR_BOOL_T, [[:EXPR_BOOL_F, :EXPR_BOOL_TH]])
addProduction(:EXPR_BOOL_TH, [[OPRLR_AND, :EXPR_BOOL_F, :EXPR_BOOL_TH],
                    [EPSILON]])
addProduction(:EXPR_BOOL_F, [[O_BRCKT, :EXPR_BOOL, C_BRCKT],
                   [CT_INT, :OPRLR_EQQ_REL , :EXPR_NUM] ,
                   [CT_FLOAT, :OPRLR_EQQ_REL , :EXPR_NUM] ,
                   [:FN_H_BOOL, :OPRLR_EQQ_REL , :EXPR_NUM] ,
                   [CT_FALSE],
                   [CT_TRUE]])

addProduction(:OPRLR_EQQ_REL , [[OPR_ATR, OPR_ATR],[OPRLR_REL]])

addProduction(:FN_H_BOOL, [[ID, :FN_H_BOOL_V]])
addProduction(:FN_H_BOOL_V, [[VEC_IN, :EXPR_NUM],
                            [:FN_CH],
                            [EPSILON]])
#=
   normal:
     ALL_ITER -> ITER_IF ALL_ITER
               -> ATTR ALL_ITER
               -> ITER_WHILE ALL_ITER
               -> ITER_FOR ALL_ITER
               -> ITER_CTRL ALL_ITER
=#
addProduction(:ALL_ITER, [[:ITER_IF, :ALL_ITER],
                           [:VAR_DCLR , :ALL_ITER],
                           [:ITER_WHILE, :ALL_ITER],
                           [:ITER_FOR, :ALL_ITER],
                           [:ITER_CTRL, :ALL_ITER],
                           [:FCALL_OR_ATRIB, :ALL_ITER],
                           [EPSILON]])

#=
   normal:
     ITER_IF  -> 'if' '(' EXPR_BOOL ')' '{' ALL_ITER '}' ITER_IF_R
     ITER_IF_R -> 'else' ITER_IF
          -> 'else' '{' ALL_ITER '}'
          -> EPSILON

   factoring:
     ITER_IF  -> 'if' '(' EXPR_BOOL ')' '{' ALL_ITER '}' ITER_IF_R

     ITER_IF_R -> 'else' ITER_IF_RR
          -> EPSILON

     ITER_IF_RR -> ITER_IF
          -> '{' ALL_ITER '}'
=#
addProduction(:ITER_IF,[[IF, O_BRCKT, :EXPR_BOOL, C_BRCKT,
                     O_C_BRCKT, :ALL_ITER,C_C_BRCKT, :ITER_IF_R]])
addProduction(:ITER_IF_R, [[ELSE, :ITER_IF_RR],
                      [EPSILON]])
addProduction(:ITER_IF_RR, [[:ITER_IF],
                      [O_C_BRCKT, :ALL_ITER, C_C_BRCKT]])
#=
   normal:
      ITER_WHILE -> 'while' '(' EXPR_BOOL ')' '{' ALL_ITER '}'
=#
addProduction(:ITER_WHILE, [[WHILE, O_BRCKT, :EXPR_BOOL, O_BRCKT, O_C_BRCKT,
                         :ALL_ITER, C_C_BRCKT]])

#=
  normal:
   ITER_FOR   -> 'for' '(' ATTR_I ',' EXPR_NUM ',' EXPR_NUM ') '{' '{' ALL_ITER '}'
   ATTR_I -> 'int' 'id' '=' :EXPR_NUM
          -> 'id' ATTR_IH

  ATTR_IH -> '=' EXPR_NUM
          -> EPSILON
=#
addProduction(:FOR_EXP_LIMITER,[[:ID_OR_VEC, OPRLR_REL,:EXPR_NUM]])

addProduction(:ITER_FOR, [[FOR, O_BRCKT, :VAR_DCLR , COMMA,:FOR_EXP_LIMITER, COMMA,
                       :ID_OR_VEC, OPR_ATR, :EXPR_NUM,C_BRCKT, O_C_BRCKT, :ALL_ITER, C_C_BRCKT]])

#=
   normal:
    ITER_CTRL -> 'continue'
          -> 'break'
          -> 'return' RH

   RH     -> EXPR_STRING
          -> EXPR_NUM
          -> EXPR_BOOL
          -> EPSILON
=#
addProduction(:ITER_CTRL, [[CONTINUE], [BREAK], [RETURN,:DATA]])
#addProduction(:RH, [#[:EXPR_STRING],
#                    [ID],
#                    [EPSILON]])
#=
addProduction(:FN_H_RET, [[ID, :FN_H_V_RET]])
addProduction(:FN_H_RET_V, [[VEC_IN, :EXPR_NUM],
                            [:FN_CH],
                            [EPSILON]])
=#

#=
  normal:
    FN_CALL -> 'id' FN_CH

    FN_CH   -> '(' FN_CHP

    FN_CHP  ->  FN_PR ')'
    FN_CHP  ->  ')'

    FN_PR   -> 'id' FN_PRE
            -> 'CT' FN_PRE

    FN_PRE  ->  ',' FN_PR
            ->  EPSILON
    CT -> 'ct_int'
          -> 'ct_float'
          -> 'ct_string'
          -> 'ct_bool'
=#

addProduction(:FN_CH, [[O_BRCKT, :FN_CHP]])
addProduction(:FN_CHP, [[:FN_PR, C_BRCKT],
                        [C_BRCKT]])
addProduction(:FN_PR, [[ID, :FN_PRE],
                       [:CT, :FN_PRE]])
addProduction(:FN_PRE, [[ID, COMMA, :FN_PR],
                        [EPSILON]])

#=
normal:
EXPR    -> EPXR  '+'  EXPR
        -> EXPR  '*'  EXPR
        -> EXPR 'or'  EXPR
        -> EXPR 'and' EXPR
        -> EXPR 'rel' EXPR
        -> '(' EXPR ')'

precedence:
EXPR    -> EPXR_T  '+'  EXPR
        -> EXPR_T

EXPR_T  -> EXPR_T  '*'  EXPR_H
        -> EXPR_B

EXPR_B  -> EXPR_I 'or'  EXPR
        -> EXPR_I 'and' EXPR
        -> EXPR_I

EXPR_I  -> EXPR_I 'rel' EXPR_J
        -> EXPR_J

EXPR_J  -> CT
        -> FN_CALL
        -> ID_OR_VEC
        -> '(' EXPR ')'

left recursion:

EXPR    -> EPXR_T  '+'  EXPR
        -> EXPR_T

EXPR_T  -> EXPR_H EXPR_TL

EXPR_TL -> '*' EXPR_B EXPR_TL
        -> 'EPSILON'

EXPR_B  -> EXPR_I 'or'  EXPR
        -> EXPR_I 'and' EXPR
        -> EXPR_I

EXPR_I  -> EXPR_J EXPR_IL
EXPR_IL -> 'rel' EXPR EXPR_IL
        -> 'EPSILON'

EXPR_J  -> CT
        -> FN_CALL
        -> ID_OR_VEC
        -> '(' EXPR ')'

ambiguity:

EXPR    -> EPXR_T EXPRL

EXPRL   -> '+' EXPR
        -> 'EPSILON'

EXPR_T  -> EXPR_H EXPR_TL

EXPR_TL -> '*' EXPR_B EXPR_TL
        -> 'EPSILON'

EXPR_B  -> EXPR_I EXPR_BL

EXPR_BL -> 'or'  EXPR
        -> 'and' EXPR
        -> 'EPSILON'

EXPR_I  -> EXPR_J EXPR_IL
EXPR_IL -> 'rel' EXPR EXPR_IL
        -> 'EPSILON'

EXPR_J  -> CT
        -> FN_CALL
        -> ID_OR_VEC
        -> '(' EXPR ')'

addProduction(:EXPR, [[:EXPR_T, :EXPRL]])
addProduction(:EXPRL, [[OPR_PM, :EXPR],
                       [EPSILON]])
addProduction(:EXPR_T, [[:EXPR_H, :EXPR_TL]])
addProduction(:EXPR_TL, [[OPR_DM, :EXPR_B, :EXPR_TL],
                         [EPSILON]])
addProduction(:EXPR_B, [[:EXPR_I, :EXPR_BL]])
addProduction(:EXPR_BL, [[OPRLR_OR, :EXPR],
                         [OPRLR_AND, :EXPR],
                         [EPSILON]])
addProduction(:EXPR_I, [[:EXPR_J, :EXPR_IL]])

addProduction(:EXPR_IL, [[OPRLR_REL, :EXPR_J, :EXPR_IL],
                         [EPSILON]])
addProduction(:EXPR_J, [[:CT],
                        [:FN_H_EXPR],
                        [O_BRCKT, :EXPR, C_BRCKT]])

addProduction(:FN_H_EXPR, [[ID, :FN_H_EXPR_RET]])
addProduction(:FN_H_EXPR_V, [[VEC_IN, :EXPR_NUM],
                             [:FN_CH],
                             [EPSILON]])

=#

max_stars_repo_name: NaelsonDouglas/namelesscompiler

max_issues_repo_name: NaelsonDouglas/namelesscompiler

max_forks_repo_name: NaelsonDouglas/namelesscompiler


text:
using POMDPs
using FiniteHorizonPOMDPs
using POMDPModelTools

#####################
# MDP and State types
#####################
struct PyramidState
    position::Int64
    epoch::Int64
end

# initial state constructor
PyramidState(position::Int64)::PyramidState = PyramidState(position, 0)

struct PyramidMDP <: MDP{PyramidState, Symbol} # Note that our MDP is parametarized by the state and the action
    _horizon::Int64
    actions::Vector{Symbol}
    action_cost::Float64
    actions_impact::Base.ImmutableDict{Symbol, Int64}
    reward_states::Vector{Vector{Int64}}
    reward::Float64
    discount_factor::Float64 # discount factor
    noise::Float64
end

function PyramidMDP(;_horizon=10, action_cost=1., reward_states=[], reward=-10., discount_factor=1., noise=.3)
    actions = [:l, :r]
    actions_impact = Base.ImmutableDict(:l => 0, :r => 1)
    return PyramidMDP(_horizon, actions, action_cost, actions_impact, reward_states, reward, discount_factor, noise)
end

FiniteHorizonPOMDPs.HorizonLength(::Type{<:PyramidMDP}) = FiniteHorizon()
FiniteHorizonPOMDPs.horizon(mdp::PyramidMDP) = mdp._horizon

###################################
# changed elements of POMDPs interface
###################################

# Creates (horizon(mdp) + 1) * mdp.no_states states to be evaluated and mdp.no_states sink states
function POMDPs.states(mdp::PyramidMDP)
    return [PyramidState(i, e) for e in 1:horizon(mdp) + 1 for i in 1:e]
end

POMDPs.stateindex(mdp::PyramidMDP, ss::PyramidState)::Int64 = sum(1:ss.epoch - 1) + ss.position
POMDPs.isterminal(mdp::PyramidMDP, ss::PyramidState) = FiniteHorizonPOMDPs.stage(mdp, ss) > horizon(mdp) || (ss.epoch, ss.position) in mdp.reward_states

# returns transition distributions - works only for 1D Gridworld with possible moves to left and to right
function POMDPs.transition(mdp::PyramidMDP, ss::PyramidState, a::Symbol)::SparseCat
    sp = (  PyramidState(ss.position + mdp.actions_impact[a], ss.epoch + 1),
            PyramidState(ss.position + mdp.actions_impact[a == :l ? :r : :l], ss.epoch + 1))
    prob = (1. - mdp.noise, mdp.noise)

    return SparseCat(sp, prob)
end

POMDPs.actions(mdp::PyramidMDP)::Vector{Symbol} = mdp.actions
POMDPs.actions(mdp::PyramidMDP, ss::PyramidState) = mdp.actions
POMDPs.actionindex(mdp::PyramidMDP, a::Symbol)::Int64 = findfirst(x->x==a, POMDPs.actions(mdp))

###############################
# FiniteHorizonPOMDPs interface
###############################
FiniteHorizonPOMDPs.stage(mdp::PyramidMDP, ss::PyramidState) = ss.epoch

function FiniteHorizonPOMDPs.stage_states(mdp::PyramidMDP, stage::Int)
    return (PyramidState(i, stage) for i in 1:stage)
end

FiniteHorizonPOMDPs.stage_stateindex(mdp::PyramidMDP, ss::PyramidState) = ss.position

###############################
# Forwarded parts of POMDPs interface
###############################
function isreward(mdp::PyramidMDP, ss::PyramidState)::Bool
    return (ss.epoch, ss.position) in mdp.reward_states
end

function POMDPs.reward(mdp::PyramidMDP, ss::PyramidState, a::Symbol, sp::PyramidState)::Float64
    isreward(mdp, sp) ? mdp.reward : mdp.action_cost
end

POMDPs.discount(mdp::PyramidMDP)::Number = mdp.discount_factor


text:
"""
Basic types and functions to support `Ballot`s, a fundamental unit of data exchange
for successful Paxos execution.
"""
module Ballots

export Ballot,
  BallotNumber, Command, InstanceID, Operation, Request, RequestID, SequenceNumber, Promise

using UUIDs

using ..Nodes
using ..Configurations

"""
An `InstanceID` uniquely names a specific instance of the Paxos algorithm.
Ballots occur within a specific instance, with the ultimate goal of
choosing a ballot by consensus.

`InstanceID`s are ordered, with "earlier" instances coming before later ones.
The ordering is intended to support ordering instances (and their chosen ballot)
into an ordered ledger of chosen commands.
"""
struct InstanceID
  value::UInt128
end

"""
Return `true` if the `left` id refers to an instance that occurs earlier
than the `right`; otherwise, return `false`
"""
function before(left::InstanceID, right::InstanceID)
  left.value < right.value
end

"""
  Return `true` if the left id refers to an instance that occurs later than
  the `right`; otherwise, return true
"""
function after(left, right)
  before(right, left)
end

struct SequenceNumber
  value::UInt128
end

struct Operation
  action::Any
end

Command = Union{Operation,Configuration}

struct RequestID
  id::UUID
  clientID::NodeID
  clientSequenceNumber::SequenceNumber
end

struct Request
  id::RequestID
  command::Command
end

struct BallotNumber
  instanceID::InstanceID
  sequenceNumber::SequenceNumber
  leaderID::NodeID
end
BallotNumber(original::BallotNumber) =
  BallotNumber(original.instanceID, original.sequenceNumber + 1, original.leaderID)

struct Ballot
  number::BallotNumber
  request::Request
end

"""
A `Promise` is a message from a member committing to a specific
`Ballot`, and agreeing not to make a similar promise for any ballot
in the same instance but with a lower sequence number.
"""
struct Promise
  memberID::NodeID
  ballotNumber::BallotNumber
end

Ballot(original::Ballot) =
  Ballot(BallotNumber(original.number), original.request)

end


text:
@with_kw  struct _color_types
    sequential::Dict{Symbol,ColorScheme} = filter(kv -> occursin("sequential", kv[2].notes), pairs(colorschemes))
    diverging::Dict{Symbol,ColorScheme} = filter(kv -> occursin("diver", kv[2].notes), pairs(colorschemes))
    cyclical::Dict{Symbol,ColorScheme} = filter(kv -> occursin("cycl", kv[2].notes), pairs(colorschemes))
    discrete::Dict{Symbol,ColorScheme} = filter(kv -> occursin("quali", kv[2].notes), pairs(colorschemes))
    all::Dict{Symbol,ColorScheme} = colorschemes
end

Base.getindex(ct::_color_types, k::Symbol) = ct.all[k]
Base.getproperty(ct::_color_types, k::Symbol) = hasfield(_color_types, k) ?  getfield(ct, :all) : ct[k]
Base.propertynames(ct::_color_types) = vcat(collect(fieldnames(typeof(ct))), collect(keys(ct.all)))

const colors = _color_types()


text:
using MindReader
using Documenter

DocMeta.setdocmeta!(MindReader, :DocTestSetup, :(using MindReader); recursive=true)

makedocs(;
    modules=[MindReader],
    authors="DanielRivasMD <danielrivasmd@gmail.com> and contributors",
    repo="https://github.com/DanielRivasMD/MindReader.jl/blob/{commit}{path}#{line}",
    sitename="MindReader.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://DanielRivasMD.github.io/MindReader.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/DanielRivasMD/MindReader.jl",
)

max_stars_repo_name: DanielRivasMD/MindReader

max_issues_repo_name: DanielRivasMD/MindReader

max_forks_repo_name: DanielRivasMD/MindReader


text:
module ModiaLogging

logTranslationDefault = false # if logging of translation should be made (in a file in the folder testresults or in the terminal window).
global logOnFile # if logging should made on file or in terminal window
global logName = "Modia_log"

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

export openLogModia, logModia, loglnModia, closeLogModia, logFileModia, setOptions, setDefaultLogName, resetTestStatus, setTestStatus, increaseLogCategory, printTestStatus

global logTranslation = logTranslationDefault

@static if VERSION < v"0.7.0-DEV.2005"
    global defaultOutput = STDOUT
else
    global defaultOutput = stdout
end

global logFileModia = defaultOutput
global defaultLogName
setDefaultLogName(name) = global defaultLogName = name

function setOptions(options)
    global logTranslation = logTranslationDefault
    if haskey(options, :logTranslation)
        global logTranslation = options[:logTranslation]
        @show logTranslation
    end

    global logOnFile = true
    if haskey(options, :logOnFile)
        global logOnFile = options[:logOnFile]
        @show logOnFile
    end

    global logName = defaultLogName
    if haskey(options, :logName)
        global logName = options[:logName]
        @show logName
    end
end

function openLogModia()
    if logOnFile && logTranslation
        logDirectory = homedir() * "/" * "ModiaResults"

        if !isdir(logDirectory)
            mkdir(homedir() * "/" * "ModiaResults")
        end

        logFileName = homedir() * "/" * "ModiaResults/" * logName * ".txt"
        global logFileModia = open(logFileName, "w")
        println("Log file: ", homedir(), "/", logFileName)
    else
        global logFileModia = defaultOutput
    end
    # write(logFileModia, "\r\n") # Trying to enable viewing log in Notepad
end

logModia(args...) = if logTranslation; print(logFileModia, args...) else end
loglnModia(args...) = if logTranslation; println(logFileModia, args...) else end
closeLogModia() = if logOnFile && logTranslation; close(logFileModia) else end

global nOK = 0
global nNOTOK = 0
global logCategories = Dict()

function resetTestStatus()
    global nOK = 0
    global nNOTOK = 0
    global logCategories = Dict()
    return
end

function setTestStatus(OK)
    if OK
        global nOK += 1
    else
        global nNOTOK += 1
    end
end

function increaseLogCategory(category)
    if haskey(logCategories, category)
        global logCategories[category] += 1
    else
        global logCategories[category] = 1
    end
end

@static if VERSION < v"0.7.0-DEV.2005"
    printstyled(s, i; bold=false, color=:black) = print_with_color(color, s, i, bold=bold)
end

function printTestStatus()
    global nOK
    global nNOTOK
    global logCategories
    println()
    println("\n----------------------\n")
    println()
    printstyled("Number of simulations OK    : ", nOK, bold=true, color=:green); println()
    printstyled("Number of simulations NOT OK: ", nNOTOK, bold=true, color=:red); println()
    println()
    println("Log category statistics:")
    for (cat, count) in logCategories
        println(cat, ": ", count)
    end
    println("\n----------------------\n")
    println()
    resetTestStatus()
end

end


text:
# To run:
# include("DI_design.jl")

# DESCRIPTION
# Trims and linearizes the rotorcraft model at incremental flight speeds, and
# generates the gains of the dynamic inversion controller. Always run after
# modifying the aircraft parameters in the AIRCRAFTNAME_constants.jl file.
#
# ------------------------------------------------------------------------------

# include packages
#using DifferentialEquations
using LinearAlgebra
using SparseArrays
using ControlSystems
using MAT
using Interpolations
using ApproxFun
using Printf
using Plots
using JLD

# include functions
include("H60_mixing.jl")
include("table_lookup.jl")
include("H60_StabSched.jl")
include("atan2.jl")
include("eqnmot.jl")
include("interpp1.jl")
include("linearize.jl")
include("interpp2.jl")
include("wrapper.jl")
include("H60_trimmer.jl")
include("ModRed.jl")
include("ModRed8.jl")
include("H60.jl")
#include("H60_DI.jl")
include("rk4.jl")
include("finp.jl")
include("simulate.jl")
# load H60 constants
include("H60_constants.jl")
# vector of trim speeds [kts]
VXTRIM_vec=[0.01,20.0,40.0,60.0,80.0,100.0,120.0,140.0,160.0]
# initialize DI matrices
CATAB = zeros(3,3,length(VXTRIM_vec))
CBinvTAB = zeros(3,3,length(VXTRIM_vec))
x0_mat = zeros(NSTATES,length(VXTRIM_vec))
u0_mat = zeros(NCTRLS,length(VXTRIM_vec))
# trim at different speeds
for iv=1:length(VXTRIM_vec)
    # load H60 constants
    include("H60_constants.jl")
    # forward, lateral, and velocities in trim [ft/s]
    VXTRIM=VXTRIM_vec[iv]*1.688
    VYTRIM=0.0
    VZTRIM=0.0
    # turn rate in trim [rad/s]
    PSIDTRIM=0.0
    # initial heading [rad]
    PSITRIM = 0.0*D2R
    # initial position [ft]
    ALTTRIM = 0.0
    XNTRIM = 0.0
    YETRIM = 0.0
    # initialize state vectors and initial guess
    # xf = [u v w p q r phi theta psi x y z]'
    xf0=zeros(NFSTATES)
    # xr = [MBC_flap, dMBC_flap, MBC_lag, dMBC_lag, MBC_inflow, azimuth, dyn twist]'
    xr0=zeros(NRSTATES)
    # xtr = TR_inflow
    xtr0=zeros(NTRSTATES)
    # initial guess for controls [%]
    u0=[50.0; 50.0; 50.0; 50.0]
    # set non-zero inflow states
    xr0[17]=0.05
    # initialize dynamic twist blade loads value [lb]
    xr0[21]=WEIGHT/NB
    # tail rotor inflow
    xtr0=0.05
    # initial guess for trim solution
    xf0[1]=VXTRIM
    xf0[2]=VYTRIM
    xf0[3]=VZTRIM
    xf0[9]=PSITRIM
    xf0[10]=XNTRIM
    xf0[11]=YETRIM
    xf0[12]=-ALTTRIM
    x0=[xf0; xr0; xtr0]
    # set trim targets
    xfdot0=zeros(NFSTATES)
    xfdot0[9]=PSIDTRIM
    xfdot0[10:12]=[VXTRIM*cos(PSITRIM)-VYTRIM*cos(PSITRIM); VXTRIM*sin(PSITRIM)+
                   VYTRIM*cos(PSITRIM); VZTRIM]
    xrdot0=zeros(NRSTATES)
    xtrdot0=zeros(NTRSTATES)
    xdot0=[xfdot0; xrdot0; xtrdot0]
    xdot_targ=xdot0[TRIMTARG]
    # Option to trim in zero sideselip or zero bank angle
    if (VXTRIM>=60*1.688)
        if PSIDTRIM==0
            # Trim heading for zero bank angle above 60 kts
            x0[7]=0.0;
            TRIMVARS[7:8]=[8; 9]
        else
            x0[9]=atan2(VYTRIM,VXTRIM)
            TRIMVARS[7:8]=[7; 8]
            # initial guess for phi (using turn coord)
            x0[7]=atan(PSIDTRIM*sqrt(VXTRIM^2+VYTRIM^2+VZTRIM^2)/G)
            # initial guess for r
            x0[6]=G/sqrt(VXTRIM^2+VYTRIM^2+VZTRIM^2)*sin(x0[7])
            # initial guess for q
            x0[5]=x0[6]*tan(x0[7])
        end
    else
        # otherwise set heading to flight path
        x0[9]=atan2(VYTRIM,VXTRIM)
        TRIMVARS[7:8]=[7; 8]
    end
    # time [s]
    t=0.0
    # trim aircraft
    x0, u0 = H60_trimmer(H60!,x0,u0,xdot_targ,t)
    # store trim state and controls
    x0_mat[:,iv] = x0
    u0_mat[:,iv] = u0
    # linearized aircraft dynamics
    A, B = linearize(H60!,x0,u0,xdot0,t)
    # reduce to 8-state model
    sysRed = ModRed8(A,B)
    A8=sysRed[1]
    B8=sysRed[2]
    eigs8=eigvals(A8)
    # dynamic inverse matrices
    C_DI=[1 0 0
          0 1 0
          0 0 1]
    CATAB[:,:,iv]=C_DI*A8[4:6,4:6]
    CBinvTAB[:,:,iv]=inv(C_DI*B8[4:6,[1,2,4]])
end
# command filters time constants
tau_p=1/3.5
tau_q=1/4.5
tau_r=0.5
# washout filter for controls
tau_u=5.0;
# disturbance rejection frequency and damping and integrator poles
wnroll_d=3.5
dmproll_d=1.0
wnpitch_d=4.5
dmppitch_d=1.0
wnyaw_d=2.;
dmpyaw_d=1.;
# feedback gains
kp_roll=2*wnroll_d*dmproll_d
ki_roll=wnroll_d^2
kp_pitch=2*wnpitch_d*dmppitch_d
ki_pitch=wnpitch_d^2
kp_yaw=2*dmpyaw_d*wnyaw_d
ki_yaw=wnyaw_d^2;
# pilot gains
k_lat=pi/100
k_lon=pi/100
k_ped=0.644
# save DI matrices and gains into a file
save("DI.jld","CATAB",CATAB,"CBinvTAB",CBinvTAB,"kp_roll",kp_roll,"ki_roll",
    ki_roll,"kp_pitch",kp_pitch,"ki_pitch",ki_pitch,"kp_yaw",kp_yaw,"ki_yaw",
    ki_yaw,"k_lat",k_lat,"k_lon",k_lon,"k_ped",k_ped,"tau_p",tau_p,
    "tau_q",tau_q,"tau_r",tau_r,"tau_u",tau_u)

max_stars_repo_name: umbertosaetti/J-GenHel

max_issues_repo_name: umbertosaetti/J-GenHel

max_forks_repo_name: umbertosaetti/J-GenHel


text:
using JuMP
include("utils.jl")

function variable_branch_flow_slacks(pm::AbstractPowerModel; nw::Int=pm.cnw)
    ### "line flow discrepancies (d), with separated positive and negative components"
    pd_br = var(pm, nw)[:pd_br] = JuMP.@variable(pm.model,
            [(l,i,j) in ref(pm,nw,:arcs),k=0:1], base_name="$(nw)_pd_br",
            lower_bound = 0,
            start = 0
    )
    qd_br = var(pm, nw)[:qd_br] = JuMP.@variable(pm.model,
            [(l,i,j) in ref(pm,nw,:arcs),k=0:1], base_name="$(nw)_qd_br",
            lower_bound = 0,
            start = 0
    )

    ### "line flow targets (t), with separated positive and negative components"
    pt_br = var(pm, nw)[:pt_br] = JuMP.@variable(pm.model,
            [(l,i,j) in ref(pm,nw,:arcs),k=0:1], base_name="$(nw)_pt_br",
            lower_bound = 0,
            start = 0
    )
    qt_br = var(pm, nw)[:qt_br] = JuMP.@variable(pm.model,
            [(l,i,j) in ref(pm,nw,:arcs),k=0:1], base_name="$(nw)_qt_br",
            lower_bound = 0,
            start = 0
    )
    for l in ids(pm,nw,:branch)
        branch = ref(pm, nw, :branch, l)
        f_bus,t_bus = branch["f_bus"],branch["t_bus"]
        f_idx,t_idx = (l, f_bus, t_bus),(l, t_bus, f_bus)
        if haskey(var(pm,nw),:p)
            p_fr = var(pm, nw,:p, f_idx)
            p_to = var(pm, nw,:p, t_idx)
            if typeof(p_fr)==VariableRef
                if JuMP.has_upper_bound(p_fr)
                    @assert JuMP.upper_bound(p_fr) >= 0
                    JuMP.set_upper_bound(pt_br[f_idx,1], JuMP.upper_bound(p_fr))
                    JuMP.delete_upper_bound(p_fr)
                end
                if JuMP.has_lower_bound(p_fr)
                    @assert JuMP.lower_bound(p_fr) <= 0
                    JuMP.set_upper_bound(pt_br[f_idx,0], abs(JuMP.lower_bound(p_fr)))
                    JuMP.delete_lower_bound(p_fr)
                end
            else
                #TODO
            end
            if typeof(p_to)==VariableRef
                if JuMP.has_upper_bound(p_to)
                    @assert JuMP.upper_bound(p_to) >= 0
                    JuMP.set_upper_bound(pt_br[t_idx,1], JuMP.upper_bound(p_to))
                    JuMP.delete_upper_bound(p_to)
                end
                if JuMP.has_lower_bound(p_to)
                    @assert JuMP.lower_bound(p_to) <= 0
                    JuMP.set_upper_bound(pt_br[t_idx,0], abs(JuMP.lower_bound(p_to)))
                    JuMP.delete_lower_bound(p_to)
                end
            else
                #TODO
            end
        end
        if haskey(var(pm,nw),:q)
            q_fr = var(pm, nw, :q, f_idx)
            q_to = var(pm, nw, :q, t_idx)
            if typeof(q_fr)==VariableRef
                if JuMP.has_upper_bound(q_fr)
                    @assert JuMP.upper_bound(q_fr) >= 0
                    JuMP.set_upper_bound(qt_br[f_idx,1], JuMP.upper_bound(q_fr))
                    JuMP.delete_upper_bound(q_fr)
                end
                if JuMP.has_lower_bound(q_fr)
                    @assert JuMP.lower_bound(q_fr) <= 0
                    JuMP.set_upper_bound(qt_br[f_idx,0], abs(JuMP.lower_bound(q_fr)))
                    JuMP.delete_lower_bound(q_fr)
                end
            else
                #TODO
            end
            if typeof(q_to)==VariableRef
                if JuMP.has_upper_bound(q_to)
                    @assert JuMP.upper_bound(q_to) >= 0
                    JuMP.set_upper_bound(qt_br[t_idx,1], JuMP.upper_bound(q_to))
                    JuMP.delete_upper_bound(q_to)
                end
                if JuMP.has_lower_bound(q_to)
                    @assert JuMP.lower_bound(q_to) <= 0
                    JuMP.set_upper_bound(qt_br[t_idx,0], abs(JuMP.lower_bound(q_to)))
                    JuMP.delete_lower_bound(q_to)
                end
            else
                #TODO
            end
        end
    end
    replaceThermalLineLimits(pm)
end

function variable_ordering_auxiliary(pm::AbstractPowerModel; nw::Int=pm.cnw)
    undecided_branches = filter(l->!(l in pm.data["protected_branches"] || l in pm.data["inactive_branches"]), ids(pm,nw,:branch))
    u_ord_aux = var(pm, nw)[:u_ord_aux] = JuMP.@variable(pm.model,
        [b in undecided_branches], base_name="$(nw)_u_ord_aux",
        lower_bound = 0,
        start = 0
    )
    u_K = var(pm,nw)[:u_K] = JuMP.@variable(pm.model, base_name="$(nw)_u_K",lower_bound=0,start=0)
end


text:
export polyinterp!

function polyinterp!(ydes::VecO{Ty}, xdes::VecO{Tx}, ysrc::VecI, xsrc::VecI; order::Int=3) where {Tx<:Real,Ty<:Real}
    lx = 1                     # left index
    rx = sz = order + 1        # right index
    ia = sz >> 1 + 1           # anchor index
    mx = size(xsrc, 1)         # upper bound of window
    xb = Vector{Tx}(undef, sz) # xsrc buffer
    yb = Vector{Ty}(undef, sz) # ysrc buffer
    tb = Vector{Ty}(undef, sz) # interpolation buffer
    @inbounds for i in 1:sz
        xb[i] = xsrc[i]
        yb[i] = ysrc[i]
    end
    @inbounds xanchor = xb[ia]
    @inbounds for i in eachindex(xdes)
        xi = xdes[i]
        #### moving interpolation window
        @inbounds while xi > xanchor && rx < mx
            lx += 1
            rx += 1
            ix  = 1
            jx  = 2
            while ix < sz
                @inbounds xb[ix] = xb[jx]
                @inbounds yb[ix] = yb[jx]
                ix += 1
                jx += 1
            end
            xb[ix]  = xsrc[rx]
            yb[ix]  = ysrc[rx]
            xanchor = xb[ia]
        end
        ydes[i] = polyinterp!(xi, xb, yb, tb, sz)
    end
    return nothing
end

function polyinterp!(x::Real, xv::VecI, yv::VecI, bv::VecB, n::Int)
    one2n = eachindex(1:n)
    @inbounds for i in one2n
        x == xv[i] && return yv[i]
    end
    Δx = Inf
    yp = 0.0
    @inbounds for i in one2n
        δx = abs(x - xv[i])
        δx < Δx && (Δx = δx; yp = yv[i])
    end
    @simd for i in one2n
        @inbounds bv[i] = yv[i] - yp
    end
    @inbounds for k in 1:n-1, i in 1:n-k
        bv[i] += (bv[i] - bv[i+1]) * (x - xv[i]) / (xv[i] - xv[i+k])
    end
    @inbounds return yp + bv[1]
end


text:
# RTG model(s)
using ModelingToolkit

"""
```julia
hil(x, k, n=1)
```

Hill function

```math
hil(x, k, n) := \frac{x^n}{x^n + k^n}
```
"""
hil(x, k) = x / (x + k)
hil(x, k, n) = hil(x^n, k^n)

"""Constant damage signal"""
CONSTANT_SIGNAL(s) = t -> s

"""Damage signal = 0"""
ZERO_SIGNAL = CONSTANT_SIGNAL(0)

"""Damage signal = 1"""
ONE_SIGNAL = CONSTANT_SIGNAL(1)

# Stressed protein levels
const STRESSED = Dict(
    ΣRtg1 => 17.2068,
    ΣRtg2 => 221.9025,
    ΣRtg3 => 25.7976,
    ΣMks => 47.8059,
    ΣBmh => 1611.5946
)

# Unstressed protein levels
const UNSTRESSED = Dict(
    ΣRtg1 => 8.0878,
    ΣRtg2 => 80.7769,
    ΣRtg3 => 50.7607,
    ΣMks => 33.158,
    ΣBmh => 1239.5666
)

"""
```julia
rtgMTK(S = DEFAULT_SIGNAL; name, proteinlevels=STRESSED)
```

Make a RTG signalling ModelingToolkit (mtk) system with a time-dependent damage signal function `S` and `proteinlevels` for Rtg1, Rtg2, Rtg3, Bmh, and Mks proteins.
"""
function RtgMTK(S=ZERO_SIGNAL; name,
    proteinlevels=STRESSED, simplify=true)

    D = Differential(t)
    @variables v[1:13](t)

    # MTK does not support elimination in overdetermined systems
    # ODEs are eliminated manually
    eqs = [
        s ~ S(t),
        # Rtg2I_c <=> Rtg2A_c
        v[1] ~ ksV * hil(s * mul_S, ksD, n_S) * Rtg2I_c - k2I * Rtg2A_c,
        # Rtg2A_c + Mks <=> Rtg2Mks_c
        v[2] ~ k2M * Rtg2A_c * Mks - kn2M * Rtg2Mks_c,
        # Bmh + Mks <=> BmhMks
        v[3] ~ kBM * Bmh * Mks - knBM * BmhMks,
        # Rtg13A_c => Rtg13I_c
        v[4] ~ (k13I + k13IV * hil(BmhMks, k13ID)) * Rtg13A_c,
        # Rtg3I_c <=> Rtg3A_c
        v[5] ~ k3A_c * Rtg3I_c - k3I_c * Rtg3A_c,
        # Rtg3A_n => Rtg3I_n
        v[6] ~ k3I_n * Rtg3A_n,
        # Rtg1_c + Rtg3A_c <=> Rtg13A_c
        v[7] ~ k13_c * Rtg1_c * Rtg3A_c - kn13_c * Rtg13A_c,
        # Rtg1_c + Rtg3I_c <=> Rtg13I_c
        v[8] ~ k13_c * Rtg1_c * Rtg3I_c - kn13_c * Rtg13I_c,
        # Rtg1_n + Rtg3A_n <=> Rtg13A_n
        v[9] ~ k13_n * Rtg1_n * Rtg3A_n - kn13_n * Rtg13A_n,
        # Rtg1_n + Rtg3I_n <=> Rtg13I_n
        v[10] ~ k13_n * Rtg1_n * Rtg3I_n - kn13_n * Rtg13I_n,
        # Rtg1_c <=> Rtg1_n
        v[11] ~ k1in * Rtg1_c - k1out * Rtg1_n,
        # Rtg3A_c <=> Rtg3A_n
        v[12] ~ k3inA * Rtg3A_c - k3outA * Rtg3A_n,
        # Rtg3I_c <=> Rtg3I_n
        v[13] ~ k3inI * Rtg3I_c - k3outI * Rtg3I_n,

        # D(Rtg2I_c) ~ -v[1],
        D(Rtg2A_c) ~ v[1] - v[2],
        # D(Mks) ~ -v[2] - v[3],
        D(Rtg2Mks_c) ~ v[2],
        # D(Bmh) ~ -v[3],
        D(BmhMks) ~ v[3],
        D(Rtg13A_c) ~ -v[4] + v[7],
        D(Rtg13I_c) ~ v[4] + v[8],
        # D(Rtg3I_c) ~ -v[5] - v[8] - v[13],
        D(Rtg3A_c) ~ v[5] - v[7] - v[12],
        D(Rtg3A_n) ~ -v[6] - v[9] + v[12],
        D(Rtg3I_n) ~ v[6] - v[10] + v[13],
        # D(Rtg1_c) ~ -v[7] - v[8] - v[11],
        D(Rtg1_n) ~ -v[9] - v[10] + v[11],
        D(Rtg13A_n) ~ v[9],
        D(Rtg13I_n) ~ v[10],

        # Conservation relationships
        ΣRtg1 ~ Rtg1_c + Rtg1_n + Rtg13A_c + Rtg13I_c + Rtg13A_n + Rtg13I_n,
        ΣRtg2 ~ Rtg2I_c + Rtg2A_c + Rtg2Mks_c,
        ΣRtg3 ~ Rtg13A_c + Rtg13I_c + Rtg13A_n + Rtg13I_n + Rtg3I_c + Rtg3A_c + Rtg3A_n + Rtg3I_n,
        ΣMks ~ Mks + Rtg2Mks_c + BmhMks,
        ΣBmh ~ Bmh + BmhMks,
    ]

    sys = ODESystem(eqs, t; name, defaults=proteinlevels)

    if simplify
        sys = structural_simplify(sys)
    end
    return sys
end

"""
Make resting initial conditions fot the system. Default to all zeros.
"""
resting_u0(sys, val=0) = Dict(states(sys) .=> val)


text:
"""
    judge(pkg::String,
          [target]::Union{String, BenchmarkConfig},
          baseline::Union{String, BenchmarkConfig};
          kwargs...)

**Arguments**:

- `pkg` - The package to benchmark.
- `target` - What do judge, given as a git id or a [`BenchmarkConfig`](@ref). If skipped, use the current state of the package repo.
- `baseline` - The commit / [`BenchmarkConfig`](@ref) to compare `target` against.

**Keyword arguments**:

- `f` - Estimator function to use in the [judging](https://github.com/JuliaCI/BenchmarkTools.jl/blob/master/doc/manual.md#trialratio-and-trialjudgement).
- `judgekwargs::Dict{Symbol, Any}` - keyword arguments to pass to the `judge` function in BenchmarkTools

The remaining keyword arguments are passed to [`benchmarkpkg`](@ref)

**Return value**:

Returns a [`BenchmarkJudgement`](@ref)
"""
function BenchmarkTools.judge(pkg::String, target::Union{BenchmarkConfig,String}, baseline::Union{BenchmarkConfig,String};
                              f=minimum, judgekwargs=Dict(), kwargs...)

    target, baseline = BenchmarkConfig(target), BenchmarkConfig(baseline)

    group_target = benchmarkpkg(pkg, target; kwargs...)
    group_baseline = benchmarkpkg(pkg, baseline; kwargs...)

    return judge(group_target, group_baseline, f; judgekwargs=judgekwargs)
end

function BenchmarkTools.judge(pkg::String, baseline::Union{BenchmarkConfig,String}; kwargs...)
    judge(pkg, BenchmarkConfig(), baseline; kwargs...)
end

"""
    judge(target::BenchmarkResults, baseline::BenchmarkResults, f;
          judgekwargs = Dict())

Judges the two [`BenchmarkResults`](@ref) in `target` and `baseline` using the function `f`.

**Return value**

Returns a [`BenchmarkJudgement`](@ref)
"""
function BenchmarkTools.judge(target::BenchmarkResults, baseline::BenchmarkResults, f = minimum; judgekwargs = Dict())
    judged = judge(f(benchmarkgroup(target)), f(benchmarkgroup(baseline)); judgekwargs...)
    return BenchmarkJudgement(target, baseline, judged)
end


text:
using Roots

export ng_zfactor
"""
    ng_zfactor(p,T;CH4=0.965,N2=0.003,CO2=0.006,C2H6=0.018,C3H8=0.045,H2O=0.,H2S=0.,
    H2=0.,CO=0.,O2=0.,iC4H10=0.001,nC4H10=0.001,iC5H12=0.0005,nC5H12=0.0003,nC6H14=0.0007,
    nC7H16=0.,nC8H18=0.,nNonane=0.,nDecane=0.,He=0.,Ar=0.)
Calculate the compress factor of specific pressure and temperature nature gas

# Auguments
...
- "p::Float64" : thre pressure,unit MPa
- "T::Float64" : the temperature,unit K
...
# Example
```
julia>ng_z_factor(6,300;CO2=0.006,N2=0.003,CH4=0.965,C2H6=0.018,C3H8=0.0045,
iC4H10=0.001,nC4H10=0.001,iC5H12=0.0005,nC5H12=0.0003,nC6H14=0.0007))
0.8953514530758864
```
"""
function ng_zfactor(p,T;
    CH4=0.965,
    N2=0.003,
    CO2=0.006,
    C2H6=0.018,
    C3H8=0.0045,
    H2O=0.,
    H2S=0.,
    H2=0.,
    CO=0.,
    O2=0.,
    iC4H10=0.001,
    nC4H10=0.001,
    iC5H12=0.0005,
    nC5H12=0.0003,
    nC6H14=0.0007,
    nC7H16=0.,
    nC8H18=0.,
    nNonane=0.,
    nDecane=0.,
    HE=0.,
    AR=0.)
    #参数检查
    #=
    ul=[12,338,1,0.2,0.2,0.1,0.035,0.015,0.005,0.001,0.0005,0.0005,0.1,0.03,0.005,0.00015]
    ll=[0;263;0.7;zeros(13)]
    pars=["pressure"=>p,"Temperature"=>T,"CH4"=>CH4,"N2"=>N2,"CO2"=>CO2,"C2H6"=>C2H6,"C3H8"=>C3H8,
          "C4H10"=>iC4H10+nC4H10,"C5H12"=>iC5H12+nC5H12,"C6"=>nC6H14,"C7"=>nC7H16,"C8+"=>nC8H18+nNonane+nDecane,
          "H2"=>H2,"CO"=>CO,"HE"=>HE,"H2O"=>H2O]
    for i=1:16
        !(ll[i]<pars[i].second<ul[i]) && error("$(pars[i].first) should be in range of $(ll[i]) - $(ul[i])")
    end
    =#

    x=[CH4,N2,CO2,C2H6,C3H8,H2O,H2S,H2,CO,O2,iC4H10,nC4H10,iC5H12,nC5H12,nC6H14,nC7H16,nC8H18,nNonane,nDecane,HE,AR]
    x=x./sum(x)
    #sum(x)!=100. && error("all componnet sum should be 100")
    N=21
    R=8.31451e-3
    #状态参数值
    a=[  0.153832600, 1.341953000,-2.998583000, -0.048312280, 0.375796500, -1.589575000,-0.053588470, 0.886594630,
        -0.710237040,-1.471722000, 1.321850350, -0.786659250, 0.2291290e-08,0.157672400,-0.436386400,-0.044081590,
        -0.003433888, 0.032059050, 0.024873550,  0.073322790,-0.001600573,  0.642470600,-0.416260100,-0.066899570,
         0.279179500,-0.696605100,-0.002860589, -0.008098836, 3.150547000,  0.007224479,-0.705752900, 0.534979200,
        -0.079314910,-1.418465000,-0.5999050e-16,0.105840200, 0.034317290, -0.007022847, 0.024955870, 0.042968180,
         0.746545300,-0.291961300, 7.294616000, -9.936757000,-0.005399808, -0.243256700, 0.049870160, 0.003733797,
         1.874951000, 0.002168144,-0.658716400,  0.000205518, 0.009776195, -0.020487080, 0.015573220, 0.006862415,
        -0.001226752, 0.002850908]
    b=Float64[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,5,5,5,5,5,6,6,7,7,8,8,8,9,9]
    k=Float64[0,0,0,0,0,0,0,0,0,0,0,0,3,2,2,2,4,4,0,0,2,2,2,4,4,4,4,0,1,1,2,2,3,3,4,4,4,0,0,2,2,2,4,4,0,2,2,4,4,0,2,0,2,1,2,2,2,2]
    c=Float64[0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,0,1,1,1,1,0,1,0,1,1,1,1,1,1]
    g=Float64[0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,1,0,0]
    f=Float64[0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
    q=Float64[0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,1,0,0,1,0,0,0,0,0,1]
    s=Float64[0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
    w=Float64[0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
    u=Float64[0,0.5 ,1,3.5 ,-0.5 ,4.5 ,0.5 ,7.5 ,9.5 ,6,12,12.5 ,-6,2,3,2,2,11,-0.5 ,0.5 ,0,4,6,21,23,22,-1,-0.5 ,7,-1,6,
              4,1,9,-13,21,8,-0.5 ,0,2,7,9,22,23,1,9,3,8,23,1.5 ,5,-0.5 ,4,7,3,0,1,0]
    #特征参数值
    M=[16.0430 , 28.0135 , 44.0100 , 30.0700 , 44.0970 , 18.0153 , 34.0820 ,  2.0159 , 28.0100 , 31.9988 ,
       58.1230 , 58.1230 , 72.1500 , 72.1500 , 86.1770 ,100.2040 ,114.2310 ,128.2580 ,142.2850 ,  4.0026 ,
       39.9480 ]
    E=[151.318300 ,  99.737780 , 241.960600 , 244.166700 , 298.118300 , 514.015600 , 296.355000 ,  26.957940 ,
       105.534800 , 122.766700 , 324.068900 , 337.638900 , 365.599900 , 370.682300 , 402.636293 , 427.722630 ,
       450.325022 , 470.840891 , 489.558373 ,   2.610111 , 119.629900 ]
    K=[ 0.4619255 ,  0.4479153 , 0.4557489 , 0.5279209 , 0.5837490 , 0.3825868 , 0.4618263 , 0.3514916 ,
        0.4533894 ,  0.4186954 , 0.6406937 , 0.6341423 , 0.6738577 , 0.6798307 , 0.7175118 , 0.7525189 ,
        0.7849550 ,  0.8152731 , 0.8437826 , 0.3589888 , 0.4216551 ]
    G=[0, 0.027815 , 0.189065 , 0.079300 , 0.141239 ,  0.332500 , 0.088500 , 0.034369 , 0.038953 , 0.021000 ,
          0.256692 , 0.281835 , 0.332267 , 0.366911 , 0.289731 , 0.337542 , 0.383381 , 0.427354 , 0.469659 , 0, 0]
    Q=[zeros(2);0.69; zeros(2); 1.06775; 0.633276; zeros(14)]
    F=[zeros(7);1.0;zeros(13,1)]
    S=[zeros(5); 1.5822; 0.390; zeros(14)]
    W=[zeros(5); 1; zeros(15)]

    #交互作用参数值
    Ex=[1.0 0.971640 0.960644 1.0 0.994635 0.708218 0.931484 1.170520 0.990126 1.0 1.019530 0.989844  1.002350   0.999268   1.107274   0.880880  0.880973   0.881067   0.881161   1.0  1.0;
    	ones(1,2) 1.022740   0.970120   0.945939   0.746954   0.902271  1.086320   1.005710   1.021000   0.946914   0.973384  0.959340   0.945520  ones(1,7);
        ones(1,3) 0.925053   0.960237   0.849408   0.955052   1.281790  1.5   1  0.906849   0.897362   0.726255  0.859764   0.855134   0.831229   0.808310   0.786323  0.765171  ones(1,2);
        ones(1,4) 1.022560   0.693168   0.946871  1.164460   ones(1,3)  1.013060   1  1.00532  ones(1,7);
        ones(1,7) 1.034787   ones(1,3)  1.0049  ones(1,9);
        ones(1,21);
        ones(1,14) 1.008692   1.010126   1.011501  1.012821   1.014089  ones(1,2);
        ones(1,8) 1.1   1  1.3   1.3  ones(1,9);ones(13,21)]

    Gx=ones(21,21);Gx[1,3]=0.807653;Gx[1,8]=1.957310;Gx[2,3]=0.982746; Gx[3,4]=0.370296; Gx[3,6]=1.673090

    Ux=[1.0 0.886106   0.963827   1  0.990877   1.0  0.736833  1.156390   ones(1,3)  0.992291   1  1.003670   1.302576  1.191904   1.205769   1.219634   1.233498  ones(1,2);
    ones(1,2) 0.835058   0.816431   0.915502   1.0 0.993476   0.408838   ones(1,3)  0.993556  ones(1,9);
    ones(1,3) 0.969870   ones(1,2)  1.045290   1.0  0.9  ones(1,5)  1.066638   1.077634   1.088178   1.098291  1.108021  ones(1,2);
    ones(1,4) 1.065173   1.0  0.971926   1.616660   ones(1,2)  ones(1,4)*1.25  ones(1,7);
    ones(2,21);
    ones(1,14)  1.028973   1.033754  1.038338   1.042735   1.046966  ones(1,2);
    ones(14,21)]

    Kx=[ones(1,1) 1.003630   0.995933   1  1.007619   1  1.000080  1.023260   ones(1,3)  0.997596   1  1.002529   0.982962  0.983565   0.982707   0.981849   0.980991  ones(1,2);
    ones(1,2) 0.982361   1.007960   1  1  0.9425960  1.032270  ones(1,13);
    ones(1,3) 1.008510   ones(1,2)  1.00779  ones(1,7) 0.910183   0.895362   0.881152   0.867520   0.854406  ones(1,2);
    ones(1,4) 0.986893   1  0.999969   1.020340  ones(1,13);
    ones(2,21);
    ones(1,16)  0.968130   0.962870   0.957828   0.952441   0.948338;
    ones(14,21)]

    Eij=[Ex[i,j]*sqrt(E[i]*E[j]) for i=1:N,j=1:N]
    Gij=[Gx[i,j]*(G[i]+G[j])/2.0 for i=1:N,j=1:N]
    Bnij=[(Gij[i,j]+1.0-g[n])^g[n]*
          (Q[i]*Q[j]+1.0-q[n])^q[n]*
          (sqrt(F[i]*F[j])+1.0-f[n])^f[n]*
          (S[i]*S[j]+1.0-s[n])^s[n]*
          (W[i]*W[j]+1.0-w[n])^w[n] for n=1:18,i=1:N,j=1:N]
    B=sum(a[n]*T^(-u[n])*sum(x[i]*x[j]*Bnij[n,i,j]*Eij[i,j]^u[n]*(K[i]*K[j])^1.5 for i=1:N,j=1:N) for n=1:18)

    U0=(sum(x.*(E.^2.5))^2.0+sum(x[i]*x[j]*(Ux[i,j]^5.0-1.0)*(E[i]*E[j])^2.5 for i=1:N-1 for j=i+1:N))^0.2
    G0=sum(x.*G)+sum(x[i]*x[j]*(Gx[i,j]-1.0)*(G[i]+G[j]) for i=1:N-1 for j=i+1:N)
    Q0=sum(x.*Q)
    F0=sum(x.^2.0.*F)
    K0=(sum(x.*K.^2.5)^2.0+2*sum(x[i]*x[j]*(Kx[i,j]^5.0-1.0)*(K[i]*K[j])^2.5  for i=1:N-1 for j=i+1:N))^0.2

    Cn=[a[n]*(G0+1-g[n])^g[n]*(Q0^2+1.0-q[n])^q[n]*(F0+1.0-f[n])^f[n]*U0^u[n]*T^(-u[n]) for n=13:58]
    Cn=[zeros(12);Cn]

    function fun(ρm)
        ρr=K0^3.0*ρm
        z=1.0+B*ρm-ρr*sum(Cn[n] for n=13:18)+sum(Cn[n]*(b[n]-c[n]*k[n]*ρr^k[n])*ρr^b[n]*exp(-c[n]*ρr^k[n]) for n=13:58)
        ρm*R*T*z-p
    end

    ρm=find_zero(fun,40,rtol=1e-6)
    p/(ρm*R*T)
end

max_stars_repo_name: hzgzh/NG

max_issues_repo_name: hzgzh/NG

max_forks_repo_name: hzgzh/NG


text:
mutable struct WeightedBagNode{T <: Union{Missing, AbstractNode}, B <: AbstractBags, W, C} <: AbstractBagNode
    data::T
    bags::B
    weights::Vector{W}
    metadata::C

    function WeightedBagNode(d::T, b::B, w::Vector{W}, m::C=nothing) where {T <: Union{Missing, AbstractNode}, B <: AbstractBags, W, C}
        ismissing(d) && any(_len.(b.bags) .!= 0) && error("WeightedBagNode with nothing in data cannot have a non-empty bag")
        new{T, B, W, C}(d, b, w, m)
    end
end

WeightedBagNode(x::T, b::Vector, weights::Vector{W}, metadata::C=nothing) where {T, W, C} =
WeightedBagNode(x, bags(b), weights, metadata)

mapdata(f, x::WeightedBagNode) = WeightedBagNode(mapdata(f, x.data), x.bags, x.weights, x.metadata)
function Base.getindex(x::WeightedBagNode{T, B, W}, i::VecOrRange) where {T, B, W}
    nb, ii = remapbag(x.bags, i)
    isempty(ii) && return(WeightedBagNode(missing, nb, W[], nothing))
    WeightedBagNode(subset(x.data,ii), nb, subset(x.weights, ii), subset(x.metadata, ii))
end

function reduce(::typeof(catobs), as::Vector{T}) where {T <: WeightedBagNode}
    data = filter(!ismissing, [x.data for x in as])
    metadata = filter(!isnothing, [x.metadata for x in as])
    bags = vcat((d.bags for d in as)...)
    WeightedBagNode(isempty(data) ? missing : reduce(catobs, data),
                    bags, vcat((a.weights for a in as)...),
            isempty(metadata) ? nothing : reduce(catobs, metadata))
end

removeinstances(a::WeightedBagNode, mask) = WeightedBagNode(subset(a.data, findall(mask)), adjustbags(a.bags, mask), subset(a.weights, findall(mask)), a.metadata)

function dsprint(io::IO, n::WeightedBagNode{T}; pad=[], s="", tr=false) where T
    c = COLORS[(length(pad)%length(COLORS))+1]
    m = T <: Missing ? " missing " : " "
    paddedprint(io, "WeightedNode with $(length(n.bags))$(m)bag(s) and weights Σw = $(sum(n.weights))$(tr_repr(s, tr))\n", color=c)
    paddedprint(io, "  └── ", color=c, pad=pad)
    dsprint(io, n.data, pad = [pad; (c, "      ")], s=s * encode(1, 1), tr=tr)
end


text:
const ESC = "\u001B" # escape code
print("$ESC[?25l")                       # hide the cursor
print("Enter anything, press RETURN: ")  # prompt shown
input = readline()                       # but no cursor
print("$ESC[0H$ESC[0J$ESC[?25h")         # reset, visible again
sleep(3)
println()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
import Distributions: logpdf, Poisson, Normal, MvNormal

################################
# poisson convenience wrappers #
################################

function elementwise_poisson(lam)
    C, H, W = size(lam)
    samp    = Array(Float64, size(lam))
    for c in 1:C, h in 1:H, w in 1:W
        samp[c, h, w] = float(rand(Distributions.Poisson(lam[c, h, w])))
    end
    return samp
end

# output a single poisson log likelihood
function poisson_lnpdf(data, lam)
    #lldist = float(logpdf(Distributions.Poisson(lam), data))
    return data*log(lam) - lam - lgamma(data+1.)
end

# constrain checking
function inrange(val, a, b)
    if (val <= a) || (val >= b)
        return false
    end
    return true
end

#######################################################
# save images + plotting stuff for numpy/matplotlib   #
#######################################################

function save_image_array(img_stack, fname)
    npzwrite(fname, img_stack)
end

function save_blob_images(images, fname)
    imgs = permutedims(cat(3, [i.pixels for i in mod_images]...), (3, 1, 2))
    npzwrite(fname, imgs)
end

function save_plot(domain, range, fname; true_val=nothing)
    domain_name = @sprintf "%s-%s" "fin" fname
    range_name  = @sprintf "%s-%s" "fval" fname
    npzwrite(domain_name, domain)
    npzwrite(range_name, range)

    if true_val!=nothing
      tname = @sprintf "%s-%s" "ftrue" fname
      npzwrite(tname, [true_val])
    end
end

"""
Given a function (lnpdf), that takes in a D dimensional vector (th),
this function fixes D-1 of the dimensions of th, and returns a function
handle that only varies with respect to dimension d.
"""
function curry_dim(lnpdf::Function, th0::Vector{Float64}, dim::Int)
    thd0 = th0[dim]

    function clnpdf(thd::Float64)
        th0[dim] = thd
        return lnpdf(th0)
    end

    return clnpdf
end

"""
Function that mimics the AccuracyBenchmark module function
variational_parameters_to_data_frame_row, for direct comparison of ground
truth parameters to VB fit parameters.
"""
function catalog_to_data_frame_row(catalog_entry; objid="truth")
    ce = catalog_entry
    df = DataFrame()
    fluxes = ce.is_star ? ce.star_fluxes : ce.gal_fluxes
    df[:objid]                          = [objid]
    df[:ra]            = [ce.pos[1]]
    df[:dec]                = [ce.pos[2]]
    df[:is_star]                        = [ce.is_star]
    df[:gal_frac_dev]  = [ce.gal_frac_dev]
    df[:gal_axis_ratio]         = [ce.gal_axis_ratio]
    df[:gal_radius_px]           = [ce.gal_radius_px]
    df[:gal_angle_deg]                      = [ce.gal_angle]
    df[:flux_r_nmgy]       = [fluxes[3]]
    df[:log_flux_r_stderr] = [NaN]
    df[:color_ug]             = [log(fluxes[2]) - log(fluxes[1])]
    df[:color_gr]             = [log(fluxes[3]) - log(fluxes[2])]
    df[:color_ri]             = [log(fluxes[4]) - log(fluxes[3])]
    df[:color_iz]             = [log(fluxes[5]) - log(fluxes[4])]
    df[:color_ug_stderr]      = [NaN]
    df[:color_gr_stderr]      = [NaN]
    df[:color_ri_stderr]      = [NaN]
    df[:color_iz_stderr]      = [NaN]
    return df
end

function samples_to_dataframe(chain; is_star=true)
    """ Turns MCMC samples (star or galaxy chain) into a dataframe of samples
    with parameters we can compare to Photo and VB inferences

    Conversions performed:
      - lnfluxes  : lnfluxes => reference band lnflux and colors
      - gal_angle : radians [0, pi] => degrees [0, 180]
      - half_light_radius_px : the Celeste parameterization is uncoupled from
         the minor_major_axis_ratio --- in the coadd+primary catalogs the
         :half_light_radius_px is scaled by the sqrt of the
         :minor_major_axis_ratio, e.g. px => px / sqrt(gal_ab)
    """
    df = DataFrame()
    # reference band flux (+ log flux) and colors
    df[:log_flux_r]  = chain[:,3]
    df[:flux_r_nmgy] = exp.(chain[:,3])
    df[:color_ug]    = chain[:,2] .- chain[:,1]
    df[:color_gr]    = chain[:,3] .- chain[:,2]
    df[:color_ri]    = chain[:,4] .- chain[:,3]
    df[:color_iz]    = chain[:,5] .- chain[:,4]
    df[:ra]          = chain[:, 6]
    df[:dec]         = chain[:, 7]
    if !is_star
      df[:gal_frac_dev]   = chain[:, 8]
      df[:gal_axis_ratio] = chain[:, 9]
      df[:gal_angle_deg]  = chain[:, 10] * 360 / (2*pi) # rad => deg
      df[:gal_radius_px]  = chain[:, 11] .* sqrt.(df[:gal_axis_ratio])
    end
    return df
end

function samples_to_dataframe_row(sampdf; is_star=true)
    """ summarize a set of samples into a single dataframe row """
    df = DataFrame()
    df[:ra]                = [mean(sampdf[:ra])]
    df[:dec]               = [mean(sampdf[:dec])]
    df[:is_star]           = [true]
    df[:gal_frac_dev]      = [NaN]
    df[:gal_axis_ratio]    = [NaN]
    df[:gal_radius_px]     = [NaN]
    df[:gal_angle_deg]     = [NaN]
    df[:flux_r_nmgy]       = [mean(sampdf[:flux_r_nmgy])]
    df[:log_flux_r]        = [mean(sampdf[:log_flux_r])]
    df[:log_flux_r_stderr] = [std(sampdf[:log_flux_r])]
    df[:color_ug]          = [mean(sampdf[:color_ug])]
    df[:color_gr]          = [mean(sampdf[:color_gr])]
    df[:color_ri]          = [mean(sampdf[:color_ri])]
    df[:color_iz]          = [mean(sampdf[:color_iz])]
    df[:color_ug_stderr]   = [std(sampdf[:color_ug])]
    df[:color_gr_stderr]   = [std(sampdf[:color_gr])]
    df[:color_ri_stderr]   = [std(sampdf[:color_ri])]
    df[:color_iz_stderr]   = [std(sampdf[:color_iz])]
    if !is_star
        df[:is_star]                        = [false]
        df[:gal_frac_dev]   = [mean(sampdf[:gal_frac_dev])]
        df[:gal_axis_ratio] = [mean(sampdf[:gal_axis_ratio])]
        df[:gal_radius_px]  = [mean(sampdf[:gal_radius_px])]
        df[:gal_angle_deg]  = [mean(sampdf[:gal_angle_deg])]
    end
    return df
end

# run multiple chains
function run_multi_chain_mcmc(lnpdf::Function, th0s)
  chains, logprobs = [], []
  for c in 1:num_chains
    # initialize a star state (around the existing catalog entry)
    println("---- chain ", c, " ----")
    th0 = init_star_params(star_state)
    samples, lls =
      run_mh_sampler_with_warmup(star_logpdf, th0,
                                 num_samples,
                                 num_warmup;
                                 print_skip=print_skip,
                                 keep_warmup=true,
                                 warmup_prop_scale=warmup_prop_scale,
                                 chain_prop_scale=chain_prop_scale)

    push!(chains, samples[(num_warmup+1):end, :])
    push!(logprobs, lls[(num_warmup+1):end])

    # report PSRF after c = 3 chains
    if c > 2
      psrfs = MCMC.potential_scale_reduction_factor(chains)
      println(" potential scale red factor", psrfs)
    end
  end
end

######################################################################
# unit flux image utils --- writes unit flux images to pixel array   #
# Very similar to Synthetic.jl functions                             #
######################################################################

"""
Add a star to a matrix of pixels.  Defaults to unit flux.
"""
function write_star_unit_flux(pos::Array{Float64, 1},
                              psf::Array{Model.PsfComponent,1},
                              wcs::WCS.WCSTransform,
                              iota::Float64,
                              pixels::Matrix{Float64};
                              offset::Array{Float64, 1}=[0., 0.],
                              flux::Float64=1.)
    # write the unit-flux scaled, correctly offset psf
    for k in 1:length(psf)
        # mean in pixels space
        the_mean = SVector{2}(WCS.world_to_pix(wcs, pos) - offset) + psf[k].xiBar
        the_cov  = psf[k].tauBar
        intensity = flux * iota * psf[k].alphaBar
        write_gaussian(the_mean, the_cov, intensity, pixels)
    end
end

"""
Add a galaxy model image to a matrix of pixels.  Defaults to unit flux.
"""
function write_galaxy_unit_flux(pos::Array{Float64, 1},
                                psf::Array{Model.PsfComponent,1},
                                wcs::WCS.WCSTransform,
                                iota::Float64,
                                gal_frac_dev::Float64,
                                gal_axis_ratio::Float64,
                                gal_angle::Float64,
                                gal_scale::Float64,
                                pixels::Matrix{Float64};
                                offset::Array{Float64, 1}=[0., 0.],
                                flux::Float64=1.)
    # write the unit-flux scaled, correctly offset galaxy shape +
    # psf convolution
    e_devs = [gal_frac_dev, 1 - gal_frac_dev]
    XiXi = Model.get_bvn_cov(gal_axis_ratio, gal_angle, gal_scale)
    for i in 1:2
        for gproto in galaxy_prototypes[i]
            for k in 1:length(psf)
                the_mean = SVector{2}(WCS.world_to_pix(wcs, pos) - offset) +
                           psf[k].xiBar
                the_cov = psf[k].tauBar + gproto.nuBar * XiXi
                intensity = flux * iota * psf[k].alphaBar * e_devs[i] * gproto.etaBar
                write_gaussian(the_mean, the_cov, intensity, pixels)
            end
        end
    end
end

function write_galaxy_unit_flux_pixel(px_pos::Array{Float64, 1},
                                      psf::Array{Model.PsfComponent,1},
                                      iota::Float64,
                                      gal_frac_dev::Float64,
                                      gal_axis_ratio::Float64,
                                      gal_angle::Float64,
                                      gal_scale::Float64,
                                      pixels::Matrix{Float64};
                                      flux::Float64=1.)
    # write the unit-flux scaled, correctly offset galaxy shape +
    # psf convolution
    e_devs = [gal_frac_dev, 1 - gal_frac_dev]
    XiXi = Model.get_bvn_cov(gal_axis_ratio, gal_angle, gal_scale)
    for i in 1:2
        for gproto in galaxy_prototypes[i]
            for k in 1:length(psf)
                the_mean = px_pos + psf[k].xiBar
                the_cov = psf[k].tauBar + gproto.nuBar * XiXi
                intensity = flux * iota * psf[k].alphaBar * e_devs[i] * gproto.etaBar
                write_gaussian(the_mean, the_cov, intensity, pixels)
            end
        end
    end
end

"""
Generate a model image on a patch, according to that image/patch psf
"""
function render_patch_nmgy(img::Image, patch::ImagePatch, n_bodies::Vector{CatalogEntry})
    # create sky noise background image in nmgy
    patch_pixels = ones(Float32, size(img.sky)) .* img.sky

    # write star/gal model images onto patch_pixels
    for body in n_bodies
        if body.is_star
            Model.write_star_nmgy!(body.pos, body.star_fluxes[img.b], patch,
                patch_pixels; write_to_patch=true)
        else
            Model.write_galaxy_nmgy!(body.pos, body.gal_fluxes[img.b],
                body.gal_frac_dev, body.gal_axis_ratio, body.gal_angle,
                body.gal_radius_px, img.psf, [patch][:,:], patch_pixels;
                write_to_patch=true)
        end
    end
    return Array{Float64}(patch_pixels)
end

"""
Write a gaussian bump on a pixel array
"""
function write_gaussian(the_mean, the_cov, intensity, pixels)
    the_precision = inv(the_cov)
    c = sqrt(det(the_precision)) / 2pi

    H, W = size(pixels)
    w_range, h_range = get_patch(the_mean, H, W)

    for w in w_range, h in h_range
        y = @SVector [the_mean[1] - h, the_mean[2] - w] # Maybe not hard code Float64
        ypy = dot(y,  the_precision * y)
        pdf_hw = c * exp(-0.5 * ypy)
        pixel_rate = intensity * pdf_hw
        pixels[h, w] += pixel_rate
    end
    pixels
end

function get_patch(the_mean::SVector{2,Float64}, H::Int, W::Int)
    radius = 50
    hm = round(Int, the_mean[1])
    wm = round(Int, the_mean[2])
    w11 = max(1, wm - radius):min(W, wm + radius)
    h11 = max(1, hm - radius):min(H, hm + radius)
    return (w11, h11)
end

"""
Potential Scale Reduction Factor --- Followed the formula from the
following website:
http://blog.stata.com/2016/05/26/gelman-rubin-convergence-diagnostic-using-multiple-chains/
"""
function potential_scale_reduction_factor(chains)
    # each chain has to be size N x D, we have M chains
    N, D  = size(chains[1])
    M     = length(chains)

    # mean and variance of each chain
    means = vcat([mean(s, 1) for s in chains]...)
    vars  = vcat([var(s, 1)  for s in chains]...)

    # grand mean
    gmu   = mean(means, 1)

    # between chain variance:w
    B = float(N)/(float(M)-1)*sum( broadcast(-, means, gmu).^2, 1)

    # average within chain variance
    W = mean(vars, 1)

    # compute PRSF ratio
    Vhat = (float(N)-1.)/float(N) * W + (float(M)+1)/float(N*M) * B
    psrf = Vhat ./ W
    return psrf
end

"""
Convert `chains` (list of sample arrays), `logprobs` (list of
log-likelihood traces) and `colnames` (list of parameter names for chains)
into a DataFrame object --- for saving and plotting
"""
function chains_to_dataframe(chains, logprobs, colnames)
    # sample, lls, chain identifier all in arrays
    samples    = vcat(chains...)
    sample_lls = vcat(logprobs...)
    chain_id   = vcat([s*ones(Integer, size(chains[s], 1))
                       for s in 1:length(chains)]...)

    # convert to data frame, rename variables
    sdf = convert(DataFrame, samples)
    rename!(sdf, names(sdf), [Symbol(s) for s in colnames])

    # add lls, chain index
    sdf[:lls]   = sample_lls
    sdf[:chain] = chain_id
    return sdf
end

"""
creates a new image from an existing larger (field) and a patch object
"""
function patch_to_image(patch::ImagePatch, img::Image; round_pixels_to_int=true)
    # subselect patch pixels from image
    Hr = patch.bitmap_offset[1]:(patch.bitmap_offset[1] +
                                size(patch.active_pixel_bitmap)[1] - 1)
    Wr = patch.bitmap_offset[2]:(patch.bitmap_offset[2] +
                                size(patch.active_pixel_bitmap)[2] - 1)
    patch_pixels = img.pixels[Hr, Wr]
    if round_pixels_to_int
        patch_pixels = round.(patch_pixels)
    end
    H, W = length(Hr), length(Wr)

    # create sub sky intensity object
    sky_small   = fill(median(img.sky.sky_small), H, W)
    calibration = img.sky.calibration[Hr]
    sky_x       = img.sky.sky_x[Hr]
    sky_y       = img.sky.sky_y[Wr]
    sky         = SDSSBackground(sky_small, sky_x, sky_y, calibration)
    nelec_per_nmgy = img.nelec_per_nmgy[Hr]

    # TODO create an adjusted WCS object
    #wcs = deepcopy(img.wcs)
    #wcs[:crpix] = wcs[:crpix] - patch.bitmap_offset

    # instantiate a smaller patch image
    patch_image = Image(patch_pixels, img.b, img.wcs,
                        patch.psf,
                        sky, nelec_per_nmgy, img.psfmap)
    return patch_image
end


text:
module RIdb

using Interpolations, HDF5

export aluminum, air, bk7, chrome, dummy, glass, gold, silicon, silicontemperature, silver, sno2f, h2o, etoh

# Workig directory
file = "RefractiveIndexDB.h5"

function air(λ::AbstractArray{T,M}) where {T<:Number, M}
    N::Array{ComplexF64} = 1.00029 * ones.(lastindex(λ)) + im .* zeros.(lastindex(λ))
end

function dummy(λ::AbstractArray{T,N}, x::S, y::S) where {T<:Number, N, S<:Number}
    Nc::Array{ComplexF64} = ones.(length(λ)) .* (x + im*y)
end # EOF dummy(...)

#http://www.cvilaser.com/PublicPages/Pages/TechnicalTips.aspx
function glass(λ::AbstractArray{T,M}) where {T<:Number, M}
    N::Array{ComplexF64} = ones.(lastindex(λ)) .* (1.52 + im * 0)
end # EOF glass(...)

#source: http://refractiveindex.info/
function etoh(λ)
    λ = λ * 1e6
    C1 = 1.35265
    C2 = 0.00306
    C3 = 0.00002
    N::Array{ComplexF64} = C1 + C2./(λ.^2) + C3./(λ.^4) + im * zero(λ)
end

#source: http://www-swiss.ai.mit.edu/~jaffer/FreeSnell/nk.html
function aluminum(λ::AbstractArray{T,M}) where {T<:Number, M}
    readf = h5open(file, "r") do file
        read(file, "aluminum")
    end

    knots = (sort(vec(readf["lambda"])).*1e9,)
    spl_n = interpolate(knots, vec(readf["n"]), Gridded(Linear()))
    spl_k = interpolate(knots, vec(readf["k"]), Gridded(Linear()))
    N::Array{ComplexF64} = spl_n(vec(λ)) + im*spl_k(vec(λ))
    return N
end # EOF aluminum(...)

function bk7(λ::AbstractArray{T,M}) where {T<:Number, M}
    readf = h5open(file, "r") do file
        read(file, "bk7")
    end
    knots = (sort(vec(readf["lambda"]).*1e9),)
    spl_n = interpolate(knots, vec(readf["n"]), Gridded(Linear()))
    spl_k = interpolate(knots, vec(readf["k"]), Gridded(Linear()))
    N::Array{ComplexF64} = spl_n(vec(λ)) + im*spl_k(λ)
    return N
end # EOF bk7(...)

#source: http://www-swiss.ai.mit.edu/~jaffer/FreeSnell/nk.html
function chrome(λ::AbstractArray{T,M}) where {T<:Number, M}
    readf = h5open(file, "r") do file
        read(file, "chrome")
    end
    knots = (sort(vec(readf["lambda"])),)
    spl_n = interpolate(knots, vec(readf["n"]), Gridded(Linear()))
    spl_k = interpolate(knots, vec(readf["k"]), Gridded(Linear()))
    N::Array{ComplexF64} = spl_n(vec(λ)) + im*spl_k(vec(λ))
    return N
end # EOF chrome(...)

#source: people.csail.mit.edu/jaffer/FreeSnell/nk.html#nk
function gold(λ::AbstractArray{T,M}) where {T<:Number, M}
    readf = h5open(file, "r") do file
        read(file, "gold")
    end
    knots = (sort(vec(readf["lambda"])),)
    spl_n = interpolate(knots, vec(readf["n"]), Gridded(Linear()))
    spl_k = interpolate(knots, vec(readf["k"]), Gridded(Linear()))
    N::Array{ComplexF64} = spl_n(vec(λ)) + im*spl_k(vec(λ))
    return N
end # EOF gold(...)

#source:
function silicon(λ::AbstractArray{T,M}) where {T<:Number, M}
    readf = h5open(file, "r") do file
        read(file, "silicon")
    end
    knots = (sort(vec(readf["lambda"]))*1e9,)
    spl_n = interpolate(knots, vec(readf["n"]), Gridded(Linear()))
    spl_k = interpolate(knots, vec(readf["k"]), Gridded(Linear()))
    N::Array{ComplexF64} = spl_n(vec(λ)) + im*spl_k(vec(λ))
    return N
end # EOF silicon(...)

#source: http://refractiveindex.info/, temperature in C (valid range: 20-450)
function silicontemperature(λ::AbstractArray{T,M}, t::S) where {T<:Number, M, S<:Number}

    if ( t > 450) | (t < 20 )
        error("Temperature range is invalid (20-450 C).")
    else
        readf = h5open(file, "r") do file
            read(file, "silicontemperature")
        end
        # λ = vec(λ)
        lambda = sort(vec(readf["lambda"]))
        # https://discourse.julialang.org/t/a-question-on-extrapolation-with-interpolations-jl/3669/4
        lambda = [prevfloat(lambda[1]);lambda[2:end-1];nextfloat(lambda[end])]
        knots = (lambda,)
        spl_n20 = interpolate(knots, vec(readf["n20"]), Gridded(Linear()))
        aux_n20 = spl_n20(λ)
        spl_k20 = interpolate(knots, vec(readf["k20"]), Gridded(Linear()))
        aux_k20 = spl_k20(λ)
        spl_n450 = interpolate(knots, vec(readf["n450"]), Gridded(Linear()))
        aux_n450 = spl_n450(λ)
        spl_k450 = interpolate(knots, vec(readf["k450"]), Gridded(Linear()))
        aux_k450 = spl_k450(λ)
        dndt = (aux_n450-aux_n20)/(450-20)
        dkdt = (aux_k450-aux_k20)/(450-20)
        if t < 215
            spl_n = interpolate(knots, vec(readf["n20"]), Gridded(Linear()))
            aux = spl_n(λ).*(1 .+ dndt*(t-20.))
            spl_k = interpolate(knots, vec(readf["k20"]), Gridded(Linear()))
            aux2 = spl_k(λ).*(1 .+ dkdt*(t-20.))
        else
            spl_n = interpolate(knots, vec(readf["n450"]), Gridded(Linear()))
            aux = spl_n(λ).*(1 .+ dndt*(t-20.))
            spl_k = interpolate(knots, vec(readf["k450"]), Gridded(Linear()))
            aux2 = spl_k(λ).*(1 .+ dkdt*(t-20.))
        end
        N::Array{ComplexF64} = aux + im*aux2
        return N
    end
end # EOF silicontemperature(...)

#source:
function silver(λ::AbstractArray{T,M}) where {T<:Number, M}
    readf = h5open(file, "r") do file
        read(file, "silver")
    end
    knots = (sort(vec(readf["lambda"])),)
    spl_n = interpolate(knots, vec(readf["n"]), Gridded(Linear()))
    spl_k = interpolate(knots, vec(readf["k"]), Gridded(Linear()))
    N::Array{ComplexF64} = spl_n(vec(λ)) + im*spl_k(vec(λ))
    return N
end # EOF silver(...)

#source: Mater. Res. Soc. Symp. Proc., 426, (1996) 449
function sno2f(λ::AbstractArray{T,M}) where {T<:Number, M}
    readf = h5open(file, "r") do file
        read(file, "sno2f")
    end
    knots = (sort(vec(readf["lambda"])*1e9),)
    spl_n = interpolate(knots, vec(readf["n"]), Gridded(Linear()))
    spl_k = interpolate(knots, vec(readf["k"]), Gridded(Linear()))
    N::Array{ComplexF64} = spl_n(vec(λ)) + im*spl_k(vec(λ))
    return N
end # EOF sno2f(...)

function h2o(λ::AbstractArray{T,M}) where {T<:Number, M}
    readf = h5open(file, "r") do file
        read(file, "h2o")
    end
    knots = (sort(vec(readf["lambda"])*1000),)
    spl_n = interpolate(knots, vec(readf["n"]), Gridded(Linear()))
    spl_k = interpolate(knots, vec(readf["k"]), Gridded(Linear()))
    N::Array{ComplexF64} = spl_n(vec(λ)) + im*spl_k(vec(λ))
    return N
end # EOF h2o(...)

end # module RIdb


text:
using MassInstallAction
using Test

@testset "MassInstallAction.jl" begin
    @testset "public.jl" begin
    end

    @testset "repo_url.jl" begin
        @test MassInstallAction._get_repo_url_with_auth("MYORG", "MYPKG"; token = "MYTOKEN") == "https://x-access-token:MYTOKEN@github.com/MYORG/MYPKG.jl"
        @test MassInstallAction._get_repo_url_with_auth("MYORG", "MYPKG"; token = nothing) == "https://github.com/MYORG/MYPKG.jl"
    end
end


text:
@testset "24.swap-nodes-in-pairs.jl" begin
    @test swap_pairs(ListNode{Int}([1, 2, 3, 4])) == ListNode{Int}([2, 1, 4, 3])
    @test swap_pairs(ListNode{Int}()) == ListNode{Int}()
    @test swap_pairs(ListNode{Int}([1])) == ListNode{Int}([1])
end


text:
#############################################################################
#############################################################################
#
# This file implements polynomial division
#
#############################################################################
#############################################################################

"""  Modular algorithm.
f divide by g

f = q*g + r

p is a prime
"""
function divide(num::Polynomial, den::Polynomial)
    function division_function(p::Int)
        f, g = mod(num,p), mod(den,p)
        degree(f) < degree(num) && return nothing
        iszero(g) && throw(DivideError())
        q = Polynomial()
        prev_degree = degree(f)
        while degree(f) ≥ degree(g)
            h = Polynomial( (leading(f) ÷ leading(g))(p) )  #syzergy
            f = mod((f - h*g), p)
            q = mod((q + h), p)
            prev_degree == degree(f) && break
            prev_degree = degree(f)
        end
        @assert iszero( mod((num  - (q*g + f)),p))
        return q, f
    end
    return division_function
end

"""
The quotient from polynomial division. Returns a function of an integer.
"""
÷(num::Polynomial, den::Polynomial)  = (p::Int) -> first(divide(num,den)(p))

"""
The remainder from polynomial division. Returns a function of an integer.
"""
rem(num::Polynomial, den::Polynomial)  = (p::Int) -> last(divide(num,den)(p))


text:
"""
Subtypes of `DistributionStrategy` define how work should be distributed among
all available processing devices.
"""
abstract type DistributionStrategy end

const TILESIZE = 32

include("strategies/serial.jl")
include("strategies/threads.jl")
include("strategies/tasks_local.jl")
include("strategies/tasks_distributed.jl")
include("strategies/procs_local.jl")
include("strategies/procs_distributed.jl")


text:
# Copyright 2021 The CLVR Authors
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""Read a libsvm binary dataset into a sparse data matrix transposed (i.e. y *. X^T)."""
function read_libsvm_into_yXT_sparse(filepath::String, dim_dataset::Int, num_dataset::Int)

    train_indices = Array{Int}([])
    feature_indices = Array{Int}([])
    values = Array{Float64}([])

    open(filepath) do f

        line = 1
        while !eof(f)
            s = readline(f)
            split_line = split(s, " ", keepempty=false)

            label = nothing
            for v in split_line
                if isnothing(label)
                    label = parse(Int, v)
                    continue
                end

                _index, _value = split(v, ":")
                index = parse(Int, _index)
                value = label * parse(Float64, _value)  # value = b_i * x_{i j}

                push!(train_indices, line)
                push!(feature_indices, index)
                push!(values, value)
            end
            line += 1
        end
    end

    sparse(feature_indices, train_indices, values, dim_dataset, num_dataset)
end

max_stars_repo_name: ericlincc/Efficient-GLP

max_issues_repo_name: ericlincc/Efficient-GLP

max_forks_repo_name: ericlincc/Efficient-GLP


text:
using LinearAlgebra
using TimerOutputs

const to = TimerOutput()

include("activations.jl")

"""
    initialize_parameters(layer_dims::Array{Int})::Dict{String, Array{Float64}}

Create and initialize a parameters dict to store parameters of a L-layer feedforward neural network.
The first argument to this function is an array representing number of units in each layer.
The second argument is a placeholder to infer the type of Arrays to initialize - either a normal array or a GPU based CuArray.

Parameters Wi (weights) and bi (biases) are initialized for every layer other than the input layer.
Weights are intialized randomly using Xavier initialization method (to normalize them) and
biases are initialized to zeros.

# Examples
```jldoctest
julia> parameters = initialize_parameters([5 10 1], [])
Dict{String,Array{Float32,N} where N} with 4 entries:
  "W2" => Float32[0.374419 -0.0859404 … -0.203984 -0.303371]
  "W1" => Float32[0.217445 -0.593327 … -0.199775 0.211567; -0.324463 0.611785 … -0.595275 0.668538; … ; -0.105528 0.425086 … -0.706531 0.53281; 0.225039 -0.08…
  "b2" => Float32[0.0]
  "b1" => Float32[0.0; 0.0; … ; 0.0; 0.0]

julia> for (key, value) in parameters
           println(key, " ", size(value))
       end
W2 (1, 10)
W1 (10, 5)
b2 (1, 1)
b1 (10, 1)
```
"""
function initialize_parameters(layer_dims::Array{Int}, Y::Array)::Dict{String, Array{Float32}}
    parameters = Dict{String, Array{Float32}}()
    for i = 2:length(layer_dims)
        # randomly initialize weights and reduce their magnitude (leads to small weights, prevents gradient explosion)
        @timeit to "Init weights" parameters[string("W", i-1)] = randn(Float32, layer_dims[i], layer_dims[i-1]) / sqrt(layer_dims[i-1]) # Xavier initialization
        # initialize biases to zero
        @timeit to "Init biases" parameters[string("b", i-1)] = zeros(Float32, layer_dims[i], 1)
    end
    return parameters
end

function forward_prop(X::Array{Float32}, parameters::Dict{String, Array{Float32}}, activations::Tuple)::Tuple{Array{Float32}, Dict{String, Array{Float32}}}
    num_layers = length(parameters) ÷ 2  # number of layers in the network

    caches = Dict{String, Array{Float32}}()  # init dict to store Z and A values for backprop

    # A0 is X (input)
    caches[string("A", 0)] = X

    Ai = Nothing
    for i = 1:num_layers
        # Z = W * A_prev + b
        @timeit to "calculate Zi" begin
            Wi = parameters[string("W", i)]
            Ai_prev = caches[string("A", i-1)]
            Zi = zeros(Float32, size(Wi)[1], size(Ai_prev)[2])
            mul!(Zi, Wi, Ai_prev)
            Zi .+= parameters[string("b", i)]
        end
        # A = activation(Z)
        @timeit to "calculate Ai" Ai = activations[i].(Zi)
        # store A and Z for use in backprop
        caches[string("Z", i)] = Zi
        caches[string("A", i)] = Ai
    end

    return Ai, caches
end

"""
Simple log binary loss
"""
function cost_binary(Y::Array{Float32}, Ŷ::Array{Float32})::Float32
    @assert length(Y) == length(Ŷ)
    m = length(Y)

    cost = - sum(Y .* log.(Ŷ) .+ (1 .- Y) .* log.(1 .- Ŷ)) / m
    return cost
end

function backward_prop(Y::Array{Float32},
                       Ŷ::Array{Float32},
                       parameters::Dict{String, Array{Float32}},
                       caches::Dict{String, Array{Float32}},
                       layer_dims::Array{Int},
                       activations::Tuple)::Dict{String, Array{Float32}}

    num_layers = length(layer_dims)
    @assert length(Y) == length(Ŷ) # verify that predictions and targets have same length
    m = size(Y)[2]

    # formula to get initial gradient
    @timeit to "Initial dA" dA = (.- Y ./ Ŷ .+ (1 .- Y) ./ (1 .- Ŷ))
    if all(isnan.(dA))
        println("dA was NaN!")
        dA = randn(Float32) # TODO: Remove this hack
    end
    # println("dA: ", dA)

    # initialize dict to store gradients
    grads = Dict{String, Array{Float32}}()

    for l in num_layers-1:-1:1
        @timeit to "calculate dZ" dZ = dA .* activations[l].(caches[string("Z", l)])  # dZ = dA * activation(Z)
        @timeit to "calculate dw" grads[string("dw", l)] = 1/m .* (dZ * transpose(caches[string("A", l-1)]))  # dW = 1/m (dZ * A(l-1)')
        @timeit to "calculate db" grads[string("db", l)] = 1/m .* sum(dZ, dims=2)
        @timeit to "calculate dA" dA = transpose(parameters[string("W", l)]) * dZ  # dA = Wl' * dZ
    end

    return grads
end

function update_parameters(parameters::Dict{String, Array{Float32}}, grads::Dict{String, Array{Float32}}, layer_dims::Array{Int}, learning_rate::Number)::Dict{String, Array{Float32}}
    num_layers = length(layer_dims)
    for l = 1:num_layers-1
        # Simple gradient descent
        @timeit to "update dW" parameters[string("W", l)] -= learning_rate .* grads[string("dw", l)]
        @timeit to "update db" parameters[string("b", l)] -= learning_rate .* grads[string("db", l)]
    end
    return parameters
end

"""
Helper functions to get backward activations from activations
"""
function get_back_activations(activations)
    activations_back = []
    for activation in activations
        push!(activations_back, @eval ($(Symbol("$activation", "_back"))))
    end
    activations_back = Tuple(activations_back)
    return activations_back
end

# Prevent Rank 0 array and use Float32 for better consistency
function reshape_Y(Y::Array)
    Y = convert(Array{Float32, ndims(Y)}, Y)
    Y = reshape(Y, 1, :)
    return Y
end

"""
    neural_network_dense(X, Y, layer_dims::Array{Int}, num_iterations::Int, learning_rate::Number; activations=Nothing, print_stats=false, parameters=nothing, resume=false, checkpoint_steps=100)

Build and train a dense neural network for binary classification. Performs batch gradient descent for optimization and uses binary logistic loss for cost function.

Also supports resuming training from previously trained parameters.

Parameters:

- `X`: training inputs (the first value in size of this should be the same as the first value in layer_dims)

- `Y`: training outputs (for binary classification, first value of size should be 1)

- `layer_dims`: Number of neurons in each layer. First layer size should be equal to the number of features in input. Output layer should be 1 in case of binary classification

- `num_iterations`: Number of iterations to train for.

- `learning_rate`: for gradient descent optimizer

- `activations` (optional): A tuple of activation functions for every layer other than the input layer (default: (relu, relu...., sigmoid))

- `print_stats` (optional): Whether to print mean and variance of parameters for every `checkpoint_steps` steps (Statistics package must be imported to use this) (default: false)

- `checkpoint_steps` (optional): the loss (and stats, if print_stats is true) is printed every `checkpoint_steps` steps (default: 100)

- `resume` (optional): whether to resume training by using the given parameters (if parameters are not given, they are initialized again) (default: false)

- `parameters` (optional): parameters for resuming training from checkpoint. Used only if resume is true (default: nothing)

Returns:
- `parameters`: trained parameters
- `activations`: activations used in training (for passing them to `predict` function)
"""
function neural_network_dense(X, Y,
                              layer_dims::Array{Int},
                              num_iterations::Int,
                              learning_rate::Number;
                              activations=Nothing,
                              print_stats=false,
                              parameters=nothing,
                              resume=false,
                              checkpoint_steps=100)
    @timeit to "calculate num layers" num_layers = length(layer_dims) # calculate number of layers

    @timeit to "reshape Y" Y = reshape_Y(Y)
    @assert ndims(Y) == 2

    @timeit to "get back activations" begin
        # if activations are not given, assume that all hidden layers have relu and output layer has sigmoid
        if activations === Nothing
            activations = Array{Function}(undef, num_layers-1)
            for i = 1:num_layers-2
                activations[i] = relu
            end
            activations[num_layers-1] = sigmoid
        end
        activations = Tuple(activations)
        activations_back = get_back_activations(activations)
    end

    @timeit to "Initialize params" begin
        # Check if training has to resume or start over form beginning
        init_params = false
        if !resume
            init_params=true
        elseif (resume && parameters==nothing)
            println("Cannot resume without parameters, pass parameters=parameters to resume training. Reinitializing parameters")
            init_params=true
        end

        # initialize params if it has to
        if init_params
            parameters = initialize_parameters(layer_dims, Y)
        end
    end

    if print_stats
        for i in eachindex(parameters)
            println("\tInitial Mean of parameter ", i, " is ", mean(parameters[i]))
            println("\tInitial Variance of parameter ", i, " is ", var(parameters[i]))
        end
    end

    # pass through the whole dataset num_iterations times
    @timeit to "Actual training" begin
        for iteration = 1:num_iterations
            @timeit to "forward prop" Ŷ, caches = forward_prop(X, parameters, activations)
            @timeit to "backprop" grads = backward_prop(Y, Ŷ, parameters, caches, layer_dims, activations_back)
            @timeit to "update params" parameters = update_parameters(parameters, grads, layer_dims, learning_rate)

            @timeit to "print stats" begin
                # print stats every few steps
                if iteration % checkpoint_steps == 0
                    cost = cost_binary(Y, Ŷ)
                    println("Cost at iteration $iteration is $cost")
                    if print_stats
                        for i in eachindex(parameters)
                            println("\tMean of parameter ", i, " is ", mean(parameters[i]))
                            println("\tVariance of parameter ", i, " is ", var(parameters[i]))
                        end
                    end
                end
            end
        end
    end

    return parameters, activations
end

"""
    predict(X, Y, parameters, activations::Tuple)

Predict using the trained parameters and calculate the accuracy.

Parameters:

- `X`: testing inputs
- `Y`: outputs to check with
- `parameters`: trained parameters which is taken from the output of `neural_network_dense`
- `activations`: also given by the outputs of `neural_network_dense`

Returns:

- `predicts`: predictions from the NN
- `accuracy`: accuracy of predictions
"""
function predict(X, Y, parameters, activations::Tuple)
    m = size(X)[2]
    n = length(parameters)
    predicts = zeros((1, m))

    # Copy Y to CPU
    if Y != nothing
        @timeit to "copy to CPU" Y = Array(Y)
        Y = reshape_Y(Y)
    end

    @timeit to "forward prop" probas, caches = forward_prop(X, parameters, activations)
    @timeit to "calculate probas" begin
        probas = Array(probas)

        for i = 1 : m
            if probas[1, i] > 0.5f0
                predicts[1, i] = 1
            else
                predicts[1, i] = 0
            end
        end
    end

    accuracy = nothing
    if Y != nothing
        @timeit to "calculate accuracy" accuracy = sum(predicts .== Y) / m
        println("Accuracy is ", accuracy*100, "%")
    end

    return predicts, accuracy
end