text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

@when("transaction commits") do context
    commit_req = g.commit_req()
    g.execute(context[:transaction], commit_req)
end

@when("session opens transaction of type: read") do context
    transaction = g.transaction(context[:session], trans_read)
    @expect transaction !== nothing
    context[:transaction] = transaction
end

@when("session opens transaction of type: write") do context
    transaction = g.transaction(context[:session], trans_write)
    @expect transaction !== nothing
    context[:transaction] = transaction
end

@then("session transaction is null: false") do context
    trans_isempty = isempty(transactions(context))
    @expect trans_isempty === false
end

@then("session transaction is open: true") do context
    sess_trans = transactions(context)
    is_open = map(item->g.is_open(item),sess_trans)
    all_open = all(is_open)
    @expect all_open === true
end

@then("session transaction has type: read") do context
    is_open = g.is_open(context[:transaction])
    @expect is_open === true
    delete_all_databases(context[:client])
end

@then("session transaction has type: write") do context
    transaction_write = context[:transaction].type
    @expect transaction_write == trans_write
    delete_all_databases(context[:client])
end

@then("session transaction commits") do context
    commit_req = g.TransactionRequestBuilder.commit_req()
    g.execute(context[:transaction], commit_req)
end

@then("session transaction commits; throws exception") do context
    commit_req = g.TransactionRequestBuilder.commit_req()
    try
        g.execute(context[:transaction], commit_req)
    catch ex
        @expect ex !== nothing
    end
    delete_all_databases(context[:client])
end

# Scenario: one database, one session, re-committing transaction throws
@when("for each session, open transaction of type: write") do context
    for session in sessions(context)
        g.transaction(session, trans_write)
    end
end

@then("for each session, transaction commits") do context
    for session in sessions(context)
        for trans in transactions(session)
            commit_trans = g.TransactionRequestBuilder.commit_req()
            g.execute(trans, commit_trans)
        end
    end
end

@then("for each session, transaction commits; throws exception") do context
    for session in sessions(context)
        for trans in transactions(session)
            commit_trans = g.TransactionRequestBuilder.commit_req()
            try
                g.execute(trans, commit_trans)
            catch ex
                @expect ex !== nothing
            end
        end
    end
    delete_all_databases(context[:client])
end

# Scenario: one database, one session, transaction close is idempotent
@then("for each session, transaction closes") do context
    for session in sessions(context)
        close.(collect(values(session.transactions)))
    end
end

@then("for each session, transaction is open: false") do context
    for session in sessions(context)
        @expect isempty(values(session.transactions))
    end
    delete_all_databases(context[:client])
end

@when("for each session, open transactions of type:") do context
    read = g.Proto.Transaction_Type.READ
    types_of_read = [row[1] for row in context.datatable]
    sessions = sessions(context)
    for session in sessions
        for type_of_read in types_of_read
            type_of_read == "read" && g.transaction(session,read)
        end
    end
end

@then("for each session, transactions are null: false") do context
    for session in sessions(context)
        transactions = collect(values(session.transactions))
        @expect isempty(transactions) === false
    end
end

@then("for each session, transactions are open: true") do context
    for session in sessions(context)
        transactions = collect(values(session.transactions))
        for transaction in transactions
            @expect g.is_open(transaction) === true
        end
    end
end

@then("for each session, transactions have type:") do context
    types_of_read = [row[1] for row in context.datatable]
    transactions = transactions(context)
    for nr in 1:length(types_of_read)
        type_transaction = types_of_read[nr] == "read" ? trans_read : trans_write
        @expect transactions[nr].type == type_transaction
    end
    delete_all_databases(context[:client])
end

@when("for each session, open transactions in parallel of type:") do context
    trans_types = [row[1] for row in context.datatable]
    for session in sessions(context)
        for nr in 1:length(trans_types)
            trans = trans_types[nr] == "read" ? g.transaction(session, trans_read) : g.transaction(session, trans_write)
            result = trans !== nothing
            @expect result === true
        end
        expectation = length(session.transactions) == length(trans_types)
        @expect expectation === true
    end
end

@then("for each session, transactions in parallel are null: false") do context
    for session in sessions(context)
        @expect isempty(transactions(session)) === false
    end
end

@then("for each session, transactions in parallel are open: true") do context
    for session in sessions(context)
        for trans in transactions(session)
            @expect g.is_open(trans) === true
        end
    end
end

@then("for each session, transactions in parallel have type:") do context
    trans_types = [row[1] for row in context.datatable]
    changed_trans_types = map(x-> x=="read" ? trans_read : trans_write, trans_types)
    res_trans_types = group_count_items(changed_trans_types)
    for session in sessions(context)
        sess_trans_type = [item.type for item in transactions(session)]
        res_trans_sess = group_count_items(sess_trans_type)
        @expect res_trans_types == res_trans_sess
    end
    delete_all_databases(context[:client])
end

@given("connection open sessions for database:") do context
    dbs = [row[1] for row in context.datatable]
    for db in dbs
        g.CoreSession(context[:client], db , g.Proto.Session_Type.DATA, request_timout=Inf)
    end
    count_result = length(context[:client].sessions) == length(dbs)
    @expect count_result === true
end

@when("for each session, open transaction of type: read") do context
    for session in sessions(context)
        g.transaction(session, trans_read)
        test_trans = length(session.transactions) == 1
        @expect test_trans === true
    end
end

@then("for each session, transaction is null: false") do context
    result_empty = map(x->!isempty(x.transactions) ,sessions(context))
    result_expect = all(result_empty)
    @expect result_expect === true
end

@then("for each session, transaction is open: true") do context
    trans = transactions.(sessions(context))
    sess_trans = reduce(vcat, trans)
    result_expect = all(g.is_open.(sess_trans))
    @expect result_expect === true
end

@then("for each session, transaction has type: read") do context
    trans = transactions.(sessions(context))
    sess_trans = reduce(vcat, trans)
    result_expect = all([trans.type == trans_read for trans in sess_trans])
    @expect result_expect === true
    delete_all_databases(context[:client])
end

@then("for each session, transaction has type: write") do context
    trans = transactions.(sessions(context))
    sess_trans = reduce(vcat, trans)
    result_expect = all([trans.type == trans_write for trans in sess_trans])
    @expect result_expect === true
    delete_all_databases(context[:client])
end

# Scenario: write in a read transaction throws
@then("graql define; throws exception containing \"schema writes when transaction type does not allow\"") do context
    define_string =   "define person sub entity;"
    try
        g.define(context[:transaction], define_string)
    catch ex
        res_comparisson = occursin("schema writes when transaction type does not allow", string(ex.error_message))
        @expect res_comparisson === true
    end
    delete_all_databases(context[:client])
end

@then("transaction commits; throws exception") do context
    try
        g.commit(context[:transaction])
    catch ex
        @expect typeof(ex) == g.TypeDBClientException
    end
    delete_all_databases(context[:client])
end


text:
function next_gaussian(obj::RandomGenerator)
    return jcall(obj, "nextGaussian", jdouble, ())
end

function next_int(obj::RandomGenerator, arg0::jint)
    return jcall(obj, "nextInt", jint, (jint,), arg0)
end

function next_long(obj::RandomGenerator, arg0::jlong)
    return jcall(obj, "nextLong", jlong, (jlong,), arg0)
end

function set_seed(obj::RandomGenerator, arg0::jint)
    return jcall(obj, "setSeed", void, (jint,), arg0)
end

function set_seed(obj::RandomGenerator, arg0::jlong)
    return jcall(obj, "setSeed", void, (jlong,), arg0)
end

ext: jl


text:
# test initilization functions

X = rand(8, 12)

# randinit

W, H = NMF.randinit(X, 5)
@test size(W) == (8, 5)
@test size(H) == (5, 12)
@test all(W .>= 0.0)
@test all(H .>= 0.0)

W, H = NMF.randinit(X, 5; zeroh=true)
@test size(W) == (8, 5)
@test size(H) == (5, 12)
@test all(W .>= 0.0)
@test all(H .== 0.0)

W, H = NMF.randinit(X, 5; normalize=true)
@test size(W) == (8, 5)
@test size(H) == (5, 12)
@test all(W .>= 0.0)
@test all(H .>= 0.0)
@test vec(sum(W, 1)) ≈ ones(5)

## nndsvd

W, H = NMF.nndsvd(X, 5)
@test size(W) == (8, 5)
@test size(H) == (5, 12)
@test all(W .>= 0.0)
@test all(H .>= 0.0)
@test vec(sum(abs2, W, 1)) ≈ ones(5) atol=1.0e-8

W2, H2 = NMF.nndsvd(X, 5; zeroh=true)
@test size(W) == (8, 5)
@test size(H) == (5, 12)
@test all(W2 .== W)
@test all(H2 .== 0.0)

W, H = NMF.nndsvd(X, 5; variant=:ar)
@test all(W .> 0.0)
# NMF.printf_mat(W)


text:
using Distributions

p = 6
λ = 9

function K(t)
    return 2λ*t / (1-2t) - p / 2 * log(1-2t)
end

function D1K(t)
    return ( 4λ*t ) / (1-2t)^2 + ( 2λ + p ) / ( 1 - 2t )
end

function D2K(t)
    return 2(p*(1-2t) + 4λ) / (1-2t)^3
end

function logf(n, t)
    return n * (K(t) - t*D1K(t)) + 0.5log(D2K(t))
end

function logg(n, t)
    return -1n * D2K(0) * t^2/2
end

function mh_saddle(T::Int = 10000; n::Int = 1)
    Z = zeros(T)
    g = Normal(0, 1 / sqrt(n*D2K(0)))
    for t = 1:T-1
        z = rand(g)
        logr = logf(n, z) - logf(n, Z[t]) + logg(n, Z[t]) - logg(n, z)
        if log(rand()) < logr
            Z[t+1] = z
        else
            Z[t+1] = Z[t]
        end
    end
    return Z
end

Z = mh_saddle(n = 1)

function tau(x)
    return ( -1p + 2x - sqrt(p^2 + 8λ * x) ) / (4x)
end

println(sum(Z .> tau(36.225)) / 10000)
println(sum(Z .> tau(40.542)) / 10000)
println(sum(Z .> tau(49.333)) / 10000)

# n = 10
Z = mh_saddle(n = 10)

println(sum(Z .> tau(113.6667)) / 10000)
println(sum(Z .> tau(102.063)) / 10000)
println(sum(Z .> tau(96.19335)) / 10000)

# n = 100
Z = mh_saddle(n = 100)

println(sum(Z .> tau(25.18054)) / 10000)
println(sum(Z .> tau(25.52361)) / 10000)
println(sum(Z .> tau(26.17395)) / 10000)

max_stars_repo_name: szcf-weiya/MCMC

max_issues_repo_name: szcf-weiya/MCMC

max_forks_repo_name: szcf-weiya/MCMC


text:
all_examples = [
    ("literate_src/1_pubmed_search_and_save.jl",    " Running Example: Search and Save"),
    # TODO how can we run UMLS examples?
    # ("literate_src/2_umls_map_and_filter.jl",       " Running Example: MeSH/UMLS Map and Filter"),
    ("literate_src/4_pubmed_export_citations.jl",   " Running Example: Export Citations"),
    ("literate_src/5_load_medline.jl", " Running Example: Load MEDLINE")
    ]

println("Running examples:")

for (example, str) in all_examples
    println("-----------------------------------------")
    println("-----------------------------------------")
    println(str)
    println("-----------------------------------------")
    println("-----------------------------------------")

    include(example)
end


text:
using LinearAlgebra: diagind

using NQCModels: NQCModels
using NQCDynamics: Calculators, RingPolymers
using .Calculators:
    AbstractDiabaticCalculator,
    DiabaticCalculator,
    RingPolymerDiabaticCalculator

abstract type StateType end
"Singleton type for labelling states as diabatic."
struct Diabatic <: StateType end
"Singleton type for labelling states as adiabatic."
struct Adiabatic <: StateType end

"""
    ElectronicDistribution{S<:StateType} <: NonadiabaticDistribution

Abstract type for distributions of electronic degrees of freedom only.
"""
abstract type ElectronicDistribution{S<:StateType} <: NonadiabaticDistribution end

"""
    SingleState{S} <: ElectronicDistribution{S}

Electronic distribution for representing a system confined to a single state.
"""
struct SingleState{S} <: ElectronicDistribution{S}
    state::Int
    statetype::S
end
SingleState(state) = SingleState(state, Diabatic())

"""
    ElectronicPopulation{T,S} <: ElectronicDistribution{S}

Electronic distribution for representing a mixed state with non-zero population in multiple states.
"""
struct ElectronicPopulation{T,S} <: ElectronicDistribution{S}
    populations::Vector{T}
    statetype::S
end
ElectronicPopulation(state) = ElectronicPopulation(state, Diabatic())

function initialise_adiabatic_density_matrix(
    electronics::ElectronicDistribution{Diabatic},
    calculator::AbstractDiabaticCalculator,
    r
)

    diabatic_density = initialise_density_matrix(electronics, calculator)
    return transform_density!(diabatic_density, calculator, r, :to_adiabatic)
end

function initialise_adiabatic_density_matrix(
    electronics::ElectronicDistribution{Adiabatic},
    calculator::AbstractDiabaticCalculator,
    r
)

    adiabatic_density = initialise_density_matrix(electronics, calculator)
    return adiabatic_density
end

function initialise_density_matrix(
    electronics::ElectronicDistribution, calculator::AbstractDiabaticCalculator
)

    n = NQCModels.nstates(calculator.model)
    density = Matrix{eltype(calculator)}(undef, n, n)
    fill_density!(density, electronics)
    return density
end

function fill_density!(density::AbstractMatrix, electronics::SingleState)
    fill!(density, zero(eltype(density)))
    density[electronics.state, electronics.state] = 1
    return density
end

function fill_density!(density::AbstractMatrix, electronics::ElectronicPopulation)
    fill!(density, zero(eltype(density)))
    density[diagind(density)] .= electronics.populations
    return density
end

function transform_density!(
    density::AbstractMatrix, calculator::AbstractDiabaticCalculator, r, direction
)
    U = evaluate_transformation(calculator, r)
    if direction === :to_diabatic
        U = U'
    elseif !(direction === :to_adiabatic)
        throw(ArgumentError("`direction` $direction not recognised."))
    end
    density .= U' * density * U
    return density
end

function evaluate_transformation(calculator::DiabaticCalculator, r)
    Calculators.evaluate_potential!(calculator, r)
    Calculators.eigen!(calculator)
    return calculator.eigen.vectors
end

function evaluate_transformation(calculator::RingPolymerDiabaticCalculator, r)
    RingPolymers.get_centroid!(calculator.tmp_centroid, r)
    Calculators.evaluate_centroid_potential!(calculator, calculator.tmp_centroid)
    Calculators.centroid_eigen!(calculator)
    return calculator.centroid_eigen.vectors
end


text:
using ConcreteStructs, ForwardDiff, IterativeSolvers, JLD2, Lazy, LoopVectorization
using Memoization, Preconditioners, QuadGK, SparseArrays, Base.Threads
using ThreadsX, ToeplitzMatrices

abstract type AbstractGrid{BC,T} <: AbstractVector{T} end

# AbstractArray interface, the rest are in the specific files below
#Base.size(g::AbstractGrid) = g.N

abstract type AbstractField{BC<:AbstractBC} end

demean!(x) = (x .-= mean(x); x)

include("fields/EquispacedValueGrids.jl")
include("fields/Fourier.jl")
include("fields/FiniteDifference.jl")
include("fields/FiniteElement.jl")

deposit!(f::AbstractField, p) = deposit!(f.chargedensity, p)
antideposit(f::AbstractField, p) = antideposit(f.electricfield, p)
electricfield(f::AbstractField, p) = antideposit(f.electricfield, p)
chargedensity(f::AbstractField, p) = antideposit(f.chargedensity, p)
electricfield(f::AbstractField, x::Number) = f.electricfield(x)
chargedensity(f::AbstractField, x::Number) = f.chargedensity(x)
cellsize(l::AbstractField) = l.chargedensity.L / numberofunknowns(l.chargedensity)
domainsize(l::AbstractField) = l.chargedensity.L
numberofunknowns(f::AbstractField) = numberofunknowns(f.chargedensity)

function energy(f::AbstractField{BC}; rtol=100eps()) where {BC}
  bc = BC(0.0, domainsize(f))
  return quadgk(x->f.electricfield(x)^2, 0, domainsize(f), rtol=rtol)[1] / 2
end
function charge(f::AbstractField{BC}) where {BC}
  bc = BC(domainsize(f))
  return mapreduce(b->integral(b, x->1, bc) * weight(b), +, bases(f.chargedensity))
end
energydensity(f::AbstractField; rtol=100eps()) = energy(f; rtol=rtol) / domainsize(f)
chargedensity(f::AbstractField) = charge(f) / domainsize(f)

function cellcentres(f::AbstractField)
  return ((1:numberofunknowns(f)) .- 0.5) * domainsize(f) / numberofunknowns(f)
end

import Base.copy!
function Base.copy!(a::AbstractField, b::AbstractField)
  a.chargedensity .= b.chargedensity # TODO make this and next line good
  a.electricfield .= b.electricfield # TODO make this and prev line good
  return a
end


text:
rootdir = "e:/xxx号机3月29日有风"

machineID ="SJZ-004"
startTime = ""
endTime = ""

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
@doc raw"""
    GammaExponentialKernel([α=1 [,γ=1]])

The gamma exponential kernel is a generalization of the exponential and squared exponential
kernels:

```math
\kappa(\mathbf{x},\mathbf{y}) = \exp\left(-\alpha ||\mathbf{x} - \mathbf{y}||^{\gamma}
\right) \qquad \alpha > 0, \; 0 < \gamma \leq 1
```
where ``\alpha`` is a scaling parameter and ``\gamma`` is a shape parameter.
"""
struct GammaExponentialKernel{T<:AbstractFloat} <: MercerKernel{T}
    alpha::HyperParameter{T}
    gamma::HyperParameter{T}
    GammaExponentialKernel{T}(α::Real, γ::Real) where {T<:AbstractFloat} = new{T}(
        HyperParameter(convert(T,α), interval(OpenBound(zero(T)), nothing)),
        HyperParameter(convert(T,γ), interval(OpenBound(zero(T)), ClosedBound(one(T))))
    )
end
function GammaExponentialKernel(α::T1=1.0, γ::T2=one(T1)) where {T1<:Real,T2<:Real}
    GammaExponentialKernel{floattype(T1,T2)}(α,γ)
end

@inline gammaexponentialkernel(z::T, α::T, γ::T) where {T<:AbstractFloat} = exp(-α*z^γ)

@inline basefunction(::GammaExponentialKernel) = SquaredEuclidean()
@inline function kappa(κ::GammaExponentialKernel{T}, z::T) where {T}
    gammaexponentialkernel(z, getvalue(κ.alpha), getvalue(κ.gamma))
end


text:
# Run the stuff in this file by going to package mode (press ] in REPL) and saying `test`
using Revise
using Test
using AutomotiveInteraction
using AutomotiveSimulator
using AutomotiveVisualization
# using Reel

f = FilteringEnvironment()

    # Test the replay video making
video_trajdata_replay(f,range=1:100,filename="media/replay.mp4")

    # Generate model driven traj and compare to ground truth visually
id_list = [6,19,28,29,34,37,40,42,43,49,50]
compare2truth(f,id_list=id_list,start_frame=101,filename = "media/compare.mp4")

    # multistep_update
final_p_mat,iterwise_p_mat = multistep_update(f,car_id=6,start_frame=1,last_frame=50,num_p=100);
# Reel.extension(m::MIME"image/svg+xml") = "svg"
# plot_pairwise_particles(iterwise_p_mat,filename="media/particles.gif")

    # imitation trajectory with leaders using replay
egoids = [28,29]
new_models,final_particles,mean_dist_mat = obtain_driver_models(f,veh_id_list=egoids,
    num_p=30,ts=1,te=30)
start_frame = 1
scene_list_1 = imitation_with_replay(f,new_models,egoids=egoids,replay_ids=[6,8,13],start_frame=start_frame)
scene_list_2 = f.traj[start_frame:start_frame+length(scene_list_1)-1]
video_overlay_scenelists(scene_list_1,scene_list_2,roadway=f.roadway,filename="replay_imit.mp4",id_list=[6,8,13,28,29])
rmse_pos_dict,rmse_vel_dict = compute_rmse(scene_list_1,scene_list_2,id_list=egoids);

    # Sample vehicles that drive together
timestamped_trajscenes = extract_timespan(f.traj)
sample_simultaneous_vehs(10,50,timestamped_trajscenes,egoid=6)

max_stars_repo_name: sisl/AutomotiveInteraction

max_issues_repo_name: sisl/AutomotiveInteraction

max_forks_repo_name: sisl/AutomotiveInteraction


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMSparse.jl/blob/master/LICENSE

type SparseVectorDOK{Tv,Ti<:Integer} <: AbstractSparseArray{Tv,Ti<:Integer,1}
    data :: Dict{Ti,Tv}
end

function SparseVectorDOK()
    return SparseVectorDOK(Dict{Int64, Float64}())
end

function SparseVectorDOK{Tv,Ti<:Integer}(b::SparseVector{Tv,Ti})
    I, V = findnz(b)
    c = SparseVectorDOK()
    add!(c, I, V)
    return c
end

function SparseVectorDOK{T}(b::Vector{T})
    return SparseVectorDOK(sparsevec(b))
end

function add!{Tv,Ti<:Integer}(b::SparseVectorDOK{Tv,Ti}, i::Ti, v::Tv)
    b.data[i] = get(b, i) + v
    return nothing
end

function add!{Tv,Ti<:Integer}(b::SparseVectorDOK{Tv,Ti}, dofs::Vector{Ti}, data::Vector{Tv})
    @assert length(dofs) == length(data)
    z = Tv(0)
    for i=1:length(dofs)
        @inbounds b.data[dofs[i]] = Base.get(b.data, i, z) + data[i]
    end
    return nothing
end

function get{Tv,Ti<:Integer}(b::SparseVectorDOK{Tv,Ti}, i::Ti)
    return Base.get(b.data, i, Tv(0))
end

function get!{Tv,Ti<:Integer}(b::SparseVectorDOK{Tv,Ti}, dofs::Vector{Ti}, data::Vector{Tv})
    for (i,j) in enumerate(dofs)
        data[i] = get(b, i)
    end
    return nothing
end


text:
include("core.jl")
#creates the beta_v and beta_T from the vector specified
function gerg_betamatrix_from_vector(v,symmetric_op = (a)->a)
    N1 = length(v)
    N = Int64((1+sqrt(1+8*N1))/2) #the real size of the matrix, this will fail if x is not a vector with the values of
    #the diagonal
    A = Array{eltype(v),2}(undef,N,N)
    count=1
   @boundscheck checkbounds(A,N,N)
   @inbounds begin
        for i = 1 : N
            A[i,i] = 1.0
                for j = i+1 : N
                A[i,j] = v[count]
                A[j,i] = symmetric_op(v[count])
                count+=1
            end
        end
   end
    return A
end

struct GERG2008 <: AbstractHelmholtzModel
    N::Int64
    acentric_factor::Vector{Float64}
    molecularWeight::Vector{Float64}
    criticalDensity::Vector{Float64}
    criticalVolume::Vector{Float64}
    criticalTemperature::Vector{Float64}
    criticalPressure::Vector{Float64}
    ideal_iters::Vector{Vector{Int64}}
    nr::Matrix{Float64}
    zeta::Matrix{Float64}
    n0ik::Vector{Vector{Float64}}
    t0ik::Vector{Vector{Float64}}
    d0ik::Vector{Vector{Float64}}
    c0ik::Vector{Vector{Float64}}
    k_pol_ik::Vector{Int64}
    k_exp_ik::Vector{Int64}
    gamma_v::Matrix{Float64}
    gamma_T::Matrix{Float64}
    beta_v::Matrix{Float64}
    beta_T::Matrix{Float64}
    Aij_indices::Vector{Tuple{Int64,Int64,Int64}}
    fij::Vector{Float64}
    dijk::Vector{Vector{Float64}}
    tijk::Vector{Vector{Float64}}
    nijk::Vector{Vector{Float64}}
    etaijk::Vector{Vector{Float64}}
    epsijk::Vector{Vector{Float64}}
    betaijk::Vector{Vector{Float64}}
    gammaijk::Vector{Vector{Float64}}
    k_pol_ijk::Vector{Int64}
    k_exp_ijk::Vector{Int64}
end
function GERG2008(xsel = collect(1:21))

    N=length(xsel)
    acentric_factor = [0.0115478,0.0377215,0.223621,0.099493,0.152291,
                    0.200164,0.183521,0.251506,0.227875,0.301261,
                    0.349469,0.399552,0.44346,0.492328,-0.215993,
                    0.0221798,0.0481621,0.344861,0.0941677,-0.390032,
                    0.0]
    molecularWeight = [16.04246,28.0134,44.0095,30.06904,44.09562,58.1222,58.1222,72.14878,72.14878,86.17536,100.20194,
    114.22852,128.2551,142.28168,2.01588,31.9988,28.0101,18.01528,34.08088,4.002602,39.948]
    criticalTemperature = [190.564,126.192,304.1282,305.322,369.825,407.817,425.125,460.35,469.7,507.82,540.13,
    569.32,594.55,617.7,33.19,154.595,132.86,647.096,373.1,5.1953,150.687]
    criticalDensity =[10.139342719,11.1839,10.624978698,6.87085454,5.000043088,3.86014294,
    3.920016792,3.271,3.215577588,2.705877875,2.315324434,2.056404127,1.81,
    1.64,14.94,13.63,10.85,17.87371609,10.19,17.399,13.407429659] #molar density, in mol/dm3
    criticalPressure=[4599000,3400000,7383000,4872000,4248000,3796000,3640000,
    3370000,3380000,3025000,2740000,2490000,2290000,2110000,1313000,5043000,
    3499000,22064000,8962910,117000,4898000]#critical pressure, Pa
    #iterations for the ideal part, for now i can´t design a general function
    #with alternating indices like in this form.
    iter_ideal = [4, 3, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 4, 2, 2, 3, 2, 0, 0]
    #original order
    iter4 = (findall(z->z==4,iter_ideal[xsel]))
    iter3 = (findall(z->z==3,iter_ideal[xsel]))
    iter2 = (findall(z->z==2,iter_ideal[xsel]))
    iter0 = (findall(z->z==0,iter_ideal[xsel]))
    ideal_iters = [iter4,iter3,iter2,iter0]

    nr = hcat([19.597508817,-83.959667892,3.00088,0.76315,0.00460,8.74432,-4.46921], #methane
    [11.083407489,-22.202102428,2.50031,0.13732,-0.14660,0.90066,0.0] ,     #nitrogen
    [11.925152758,-16.118762264,2.50002,2.04452,-1.06044,2.03366,0.01393],  #carbon dioxide
    [24.675437527,-77.425313760,3.00263,4.33939,1.23722,13.19740,-6.01989], #ethane
    [31.602908195,-84.463284382,3.02939,6.60569,3.19700,19.19210,-8.37267], #propane
    [20.884143364,-91.638478026,3.33944,9.44893,6.89406,24.46180,14.78240], #n-butane
    [20.413726078,-94.467620036,3.06714,8.97575,5.25156,25.14230,16.13880], #iso-butane
    [14.536611217,-89.919548319,3.00000,8.95043,21.83600,33.40320,0.0],     #n-pentane
    [15.449907693,-101.298172792,3.00000,11.76180,20.11010,33.16880,0.0],   #isopentane
    [14.345969349,-96.165722367,3.00000,11.69770,26.81420,38.61640,0.0],    #n-hexane
    [15.063786601,-97.345252349,3.00000,13.72660,30.47070,43.55610,0.0],    #n-heptane
    [15.864687161,-97.370667555,3.00000,15.68650,33.80290,48.17310,0.0],   #n-octane
    [16.313913248,-102.160247463,3.00000,18.02410,38.12350,53.34150,0.0],  #n-nonane
    [15.870791919,-108.858547525,3.00000,21.00690,43.49310,58.36570,0.0],       #n-decane
    [13.796443393,-175.864487294,1.47906,0.95806,0.45444,1.56039,-1.37560], #hydrogen
    [10.001843586,-14.996095135,2.50146,1.07558,1.01334,0.0,0.0],           #oxygen
    [10.813340744,-19.834733959,2.50055,1.02865,0.00493,0.0,0.0],           #carbon monoxide
    [8.203520690,-11.996306443,3.00392,0.01059,0.98763,3.06904,0.0],            #water
    [9.336197742,-16.266508995,3.00000,3.11942,1.00243,0.0,0.0],            #hydrogen sulfide
    [13.628409737,-143.470759602,1.50000,0.0,0.0,0.0,0.0],                  #helium
    [8.316631500,-4.946502600,1.50000,0.0,0.0,0.0,0.0]                      #argon
    )

    zeta = hcat([4.306474465,0.936220902,5.577233895,5.722644361],                  #methane
    [5.251822620,-5.393067706,13.788988208,0],                        #nitrogen
    [3.022758166,-2.844425476,1.589964364,1.121596090],                 #carbon dioxide
    [1.831882406,0.731306621,3.378007481,3.508721939],                  #ethane
    [1.297521801,0.543210978,2.583146083,2.777773271],                  #propane
    [1.101487798,0.431957660,4.502440459,2.124516319],                  #n-butane
    [1.074673199,0.485556021,4.671261865,2.191583480],                  #isobutane
    [0.380391739,1.789520971,3.777411113,0],                          #n-pentane
    [0.635392636,1.977271641,4.169371131,0],                          #isopentane
    [0.359036667,1.691951873,3.596924107,0],                          #n-hexane
    [0.314348398,1.548136560,3.259326458,0],                          #n-heptane
    [0.279143540,1.431644769,2.973845992,0],                          #n-octane
    [0.263819696,1.370586158,2.848860483,0],                          #n-nonane
    [0.267034159,1.353835195,2.833479035,0],                          #n-decane
    [6.891654113,9.847634830,49.765290750,50.367279301],                #hydrogen
    [14.461722565,7.223325463,0,0],                                 #oxygen
    [11.669802800,5.302762306,0,0],                                 #carbon monoxide
    [0.415386589,1.763895929,3.874803739,0],                          #water
    [4.914580541,2.270653980,0,0],                                  #hydrogen sulfide
    [0.0,0,0.0,0.0],                                                  #helium
    [0.0,0,0.0,0.0]                                                   #argon
    )

    n0ik=[[0.57335704239162,-0.16760687523730e1,0.23405291834916,-0.21947376343441,0.16369201404128e-1,
    0.15004406389280e-1,0.98990489492918e-1,0.58382770929055,-0.74786867560390,0.30033302857974,0.20985543806568,
    -0.18590151133061e-1,-0.15782558339049,0.12716735220791,-0.32019743894346e-1,-0.68049729364536e-1,0.24291412853736e-1,
    0.51440451639444e-2,-0.19084949733532e-1,0.55229677241291e-2,-0.44197392976085e-2,0.40061416708429e-1,-0.33752085907575e-1,
    -0.25127658213357e-2],#methane
    [0.59889711801201,-0.16941557480731e1,0.24579736191718,-0.23722456755175,0.17954918715141e-1,0.14592875720215e-1,0.10008065936206,
    0.73157115385532,-0.88372272336366,0.31887660246708,0.20766491728799,-0.19379315454158e-1,-0.16936641554983,0.13546846041701,
    -0.33066712095307e-1,-0.60690817018557e-1,0.12797548292871e-1,0.58743664107299e-2,-0.18451951971969e-1,0.47226622042472e-2,
    -0.52024079680599e-2,0.43563505956635e-1,-0.36251690750939e-1,-0.28974026866543e-2],#nitrogen
    [0.52646564804653,-0.14995725042592e1,0.27329786733782,0.12949500022786,0.15404088341841,-0.58186950946814,-0.18022494838296,
    -0.95389904072812e-1,-0.80486819317679e-2,-0.35547751273090e-1,-0.28079014882405,-0.82435890081677e-1,0.10832427979006e-1,
    -0.67073993161097e-2,-0.46827907600524e-2,-0.28359911832177e-1,0.19500174744098e-1,-0.21609137507166,0.43772794926972,
    -0.22130790113593,0.15190189957331e-1,-0.15380948953300e-1],#carbon_dioxide
    [0.63596780450714,-0.17377981785459e1,0.28914060926272,-0.33714276845694,0.22405964699561e-1,0.15715424886913e-1,0.11450634253745,
    0.10612049379745e1,-0.12855224439423e1,0.39414630777652,0.31390924682041,-0.21592277117247e-1,-0.21723666564905,-0.28999574439489,
    0.42321173025732,0.46434100259260e-1,-0.13138398329741,0.11492850364368e-1,-0.33387688429909e-1,0.15183171583644e-1,
    -0.47610805647657e-2,0.46917166277885e-1,-0.39401755804649e-1,-0.32569956247611e-2],#ethane
    [.10403973107358e1,-0.28318404081403e1,0.84393809606294,-0.76559591850023e-1,0.94697373057280e-1,0.24796475497006e-3,0.27743760422870,
    -0.43846000648377e-1,-0.26991064784350,-0.69313413089860e-1,-0.29632145981653e-1,0.14040126751380e-1],#propane
    [0.10626277411455e1,-0.28620951828350e1,0.88738233403777,-0.12570581155345,0.10286308708106,0.25358040602654e-3,0.32325200233982,
    -0.37950761057432e-1,-0.32534802014452,-0.79050969051011e-1,-0.20636720547775e-1,0.57053809334750e-2],#n-butane
    [0.10429331589100e1,-0.28184272548892e1,0.86176232397850,-0.10613619452487,0.98615749302134e-1,0.23948208682322e-3,0.30330004856950,
    -0.41598156135099e-1,-0.29991937470058,-0.80369342764109e-1,-0.29761373251151e-1,0.13059630303140e-1],#isobutane
    [0.10968643098001e1,-0.29988888298061e1,0.99516886799212,-0.16170708558539,0.11334460072775,
    0.26760595150748e-3,0.40979881986931,-0.40876423083075e-1,-0.38169482469447,-0.10931956843993,
    -0.32073223327990e-1,0.16877016216975e-1],#n-pentane
    [0.10963e1,-0.30402e1,0.10317e1,-0.15410,0.11535,0.29809e-3,0.39571,-0.45881e-1,-0.35804,-0.10107,
    -0.35484e-1,0.18156e-1],#isopentane
    [0.10553238013661e1,-0.26120615890629e1,0.76613882967260,-0.29770320622459,0.11879907733358,0.27922861062617e-3,
    0.46347589844105,0.11433196980297e-1,-0.48256968738131,-0.93750558924659e-1,-0.67273247155994e-2,-0.51141583585428e-2],#n-hexane
    [0.10543747645262e1,-0.26500681506144e1,0.81730047827543,-0.30451391253428,0.12253868710800,0.27266472743928e-3,0.49865825681670,
    -0.71432815084176e-3,-0.54236895525450,-0.13801821610756,-0.61595287380011e-2,0.48602510393022e-3],#n-heptane
    [0.10722544875633e1,-0.24632951172003e1,0.65386674054928,-0.36324974085628,0.12713269626764,0.30713572777930e-3,0.52656856987540,
    0.19362862857653e-1,-0.58939426849155,-0.14069963991934,-0.78966330500036e-2,0.33036597968109e-2],#n-octane
    [0.11151e1,-0.27020e1,0.83416,-0.38828,0.13760,0.28185e-3,0.62037,0.15847e-1,-0.61726,-0.15043,-0.12982e-1,0.44325e-2],#n-nonane
    [0.10461e1,-0.24807e1,0.74372,-0.52579,0.15315,0.32865e-3,0.84178,0.55424e-1,-0.73555,-0.18507,-0.20775e-1,0.12335e-1],#n-decane
    [0.53579928451252e1,-0.62050252530595e1,0.13830241327086,-0.71397954896129e-1,0.15474053959733e-1,-0.14976806405771,
    -0.26368723988451e-1,0.56681303156066e-1,-0.60063958030436e-1,-0.45043942027132,0.42478840244500,-0.21997640827139e-1,
    -0.10499521374530e-1,-0.28955902866816e-2],#hydrogen
    [0.88878286369701,-0.24879433312148e1,0.59750190775886,0.96501817061881e-2,0.71970428712770e-1,0.22337443000195e-3,
    0.18558686391474,-0.38129368035760e-1,-0.15352245383006,-0.26726814910919e-1,-0.25675298677127e-1,0.95714302123668e-2],#oxygen
    [0.90554,-0.24515e1,0.53149,0.24173e-1,0.72156e-1,0.18818e-3,0.19405,-0.43268e-1,-0.12778,-0.27896e-1,-0.34154e-1,0.16329e-1],#carbon_monoxide
    [0.82728408749586,-0.18602220416584e1,-0.11199009613744e1,0.15635753976056,0.87375844859025,-0.36674403715731,
    0.53987893432436e-1,0.10957690214499e1,0.53213037828563e-1,0.13050533930825e-1,-0.41079520434476,0.14637443344120,
    -0.55726838623719e-1,-0.11201774143800e-1,-0.66062758068099e-2,0.46918522004538e-2],#water
    [0.87641,-0.20367e1,0.21634,-0.50199e-1,0.66994e-1,0.19076e-3,0.20227,-0.45348e-2,-0.22230,-0.34714e-1,-0.14885e-1,0.74154e-2],#hydrogen_sulfide
    [-0.45579024006737,0.12516390754925e1,-0.15438231650621e1,0.20467489707221e-1,-0.34476212380781,-0.20858459512787e-1,
    0.16227414711778e-1,-0.57471818200892e-1,0.19462416430715e-1,-0.33295680123020e-1,-0.10863577372367e-1,-0.22173365245954e-1],#helium
    [0.85095714803969,-0.24003222943480e1,0.54127841476466,0.16919770692538e-1,0.68825965019035e-1,0.21428032815338e-3,0.17429895321992,
    -0.33654495604194e-1,-0.13526799857691,-0.16387350791552e-1,-0.24987666851475e-1,0.88769204815709e-2]#argon
    ]

    #most_compounds
    c0ik1=[1,1,2,2,3,3]
    d0ik1=[1,1,1,2,3,7,2,5,1,4,3,4]
    t0ik1=[0.250,1.125,1.500,1.375,0.250,0.875,0.625,1.750,3.625,3.625,14.500,12.000]
    k_exp_ik1 = length(c0ik1)
    k_pol_ik1 = length(d0ik1) - k_exp_ik1

    #methane,nitrogen,ethane
    c0ik2=[1,1,1,1,1,1,2,2,2,2,2,3,3,3,6,6,6,6]
    d0ik2=[1,1,2,2,4,4,1,1,1,2,3,6,2,3,3,4,4,2,3,4,5,6,6,7]
    t0ik2=[0.125,1.125,0.375,1.125,0.625,1.500,0.625,2.625,2.750,2.125,2.000,1.750,4.500,4.750,5.000,4.000,4.500,7.500,
    14.000,11.500,26.000,28.000,30.000,16.000]
    k_exp_ik2 = length(c0ik2)
    k_pol_ik2 = length(d0ik2) - k_exp_ik2

    #carbon_Dioxide
    c0ik3=[1,1,1,1,1,1,2,2,3,3,3,3,3,5,5,5,6,6]
    d0ik3=[1,1,2,3,3,3,4,5,6,6,1,4,1,1,3,3,4,5,5,5,5,5]
    t0ik3=[0.000,1.250,1.625,0.375,0.375,1.375,1.125,1.375,0.125,1.625,3.750,3.500,7.500,8.000,6.000,
    16.000,11.000,24.000,26.000,28.000,24.000,26.000]
    k_exp_ik3 = length(c0ik3)
    k_pol_ik3 = length(d0ik3) - k_exp_ik3
    #hydrogen
    c0ik4=[1,1,1,1,2,2,3,3,5]
    d0ik4=[1,1,2,2,4,1,5,5,5,1,1,2,5,1]
    t0ik4=[0.500,0.625,0.375,0.625,1.125,2.625,0.000,0.250,1.375,
    4.000,4.250,5.000,8.000,8.000]
    k_exp_ik4 = length(c0ik4)
    k_pol_ik4 = length(d0ik4) - k_exp_ik4

    #water
    c0ik5=[1,1,1,2,2,2,3,5,5]
    d0ik5=[1,1,1,2,2,3,4,1,5,5,1,2,4,4,1,1]
    t0ik5=[0.500,1.250,1.875,0.125,1.500,1.000,0.750,1.500,0.625,2.625,5.000,4.000,4.500,3.000,4.000,6.000]
    k_exp_ik5 = length(c0ik5)
    k_pol_ik5 = length(d0ik5) - k_exp_ik5

    #helium
    c0ik6=[1,1,1,1,1,2,3,3]
    d0ik6=[1,1,1,4,1,3,5,5,5,2,1,2]
    t0ik6=[0.000,0.125,0.750,1.000,0.750,2.625,0.125,1.250,2.000,1.000,4.500,5.000]
    k_exp_ik6 = length(c0ik6)
    k_pol_ik6 = length(d0ik6) - k_exp_ik6

    c0ik = [c0ik2,c0ik2,c0ik3,c0ik2,c0ik1,
            c0ik1,c0ik1,c0ik1,c0ik1,c0ik1,
            c0ik1,c0ik1,c0ik1,c0ik1,c0ik4,
            c0ik1,c0ik1,c0ik5,c0ik1,c0ik6,
            c0ik1]

    d0ik = [d0ik2,d0ik2,d0ik3,d0ik2,d0ik1,
            d0ik1,d0ik1,d0ik1,d0ik1,d0ik1,
            d0ik1,d0ik1,d0ik1,d0ik1,d0ik4,
            d0ik1,d0ik1,d0ik5,d0ik1,d0ik6,
            d0ik1]

    t0ik = [t0ik2,t0ik2,t0ik3,t0ik2,t0ik1,
            t0ik1,t0ik1,t0ik1,t0ik1,t0ik1,
            t0ik1,t0ik1,t0ik1,t0ik1,t0ik4,
            t0ik1,t0ik1,t0ik5,t0ik1,t0ik6,
            t0ik1]

    k_pol_ik = vcat(k_pol_ik2,k_pol_ik2,k_pol_ik3,k_pol_ik2,k_pol_ik1,
                k_pol_ik1,k_pol_ik1,k_pol_ik1,k_pol_ik1,k_pol_ik1,
                k_pol_ik1,k_pol_ik1,k_pol_ik1,k_pol_ik1,k_pol_ik4,
                k_pol_ik1,k_pol_ik1,k_pol_ik5,k_pol_ik1,k_pol_ik6,
                k_pol_ik1)

    k_exp_ik = vcat(k_exp_ik2,k_exp_ik2,k_exp_ik3,k_exp_ik2,k_exp_ik1,
                k_exp_ik1,k_exp_ik1,k_exp_ik1,k_exp_ik1,k_exp_ik1,
                k_exp_ik1,k_exp_ik1,k_exp_ik1,k_exp_ik1,k_exp_ik4,
                k_exp_ik1,k_exp_ik1,k_exp_ik5,k_exp_ik1,k_exp_ik6,
                k_exp_ik1)

    fgammat = (a,b) -> 0.5*(a+b)/sqrt(a*b)
    fgammav = (a,b) -> 4*(a+b)/(a^(1/3)+b^(1/3))^3

    gamma_T = mixing_matrix(fgammat,criticalTemperature)

    gamma_v = mixing_matrix(fgammav,1 ./ criticalDensity)

    vector_beta_v=[0.9987213770,0.9995180720,0.9975478660,1.0048270700,0.9791059720,
    1.0112403880,0.9483301200,1.0000000000,0.9580152940,0.9620508310,0.9947406030,
    1.0028522870,1.0330862920,1.0000000000,1.0000000000,0.9973407720,1.0127831690,
    1.0125990870,1.0000000000,1.0346302590,0.9777946340,0.9788801680,0.9744246810,
    0.9960826100,0.9864158300,1.0000000000,1.0000000000,1.0000000000,1.0000000000,
    1.0000000000,1.0000000000,1.0000000000,0.9725320650,1.0000000000,1.0000000000,
    0.9995217700,0.9103942490,0.9695010550,1.0041664120,1.0025257180,0.9968980040,
    1.1747609230,1.0765518820,1.0243114980,1.0607931040,1.0000000000,1.2054699760,
    1.0261693730,1.0000000000,1.0001511320,0.9041421590,1.0000000000,1.0000000000,
    0.9490559590,0.9066305640,0.8466475610,1.0083924280,0.9976072770,0.9991572050,
    1.0000000000,0.9938510090,1.0000000000,1.0000000000,1.0000000000,1.0074697260,
    1.0000000000,0.9956762580,0.9253671710,1.0000000000,1.0000000000,1.0000000000,
    1.0108179090,1.0000000000,1.0000000000,0.9997958680,0.9992431460,1.0449194310,
    1.0404592890,1.0000000000,1.0000000000,1.0000000000,1.0000000000,0.9841042270,
    1.0000000000,1.0000000000,1.0000000000,1.0000000000,0.9368112190,1.0000000000,
    1.0000000000,1.0008804640,1.0000000000,1.0000000000,1.0000000000,1.0000000000,
    1.0000000000,1.0000000000,0.9769519680,1.0000000000,1.0000000000,1.0000000000,
    1.0000000000,0.9081131630,1.0000000000,1.0000000000,1.0000000000,1.0000000000,
    1.0000000000,1.0000000000,1.0000000000,1.0000000000,1.0000000000,1.0000000000,
    1.0000000000,1.0000000000,1.0000000000,1.0129944310,1.0000000000,1.0000000000,
    1.0000000000,1.0000000000,1.0000000000,1.0000000000,1.0000000000,1.0000000000,
    1.0000000000,1.0000000000,1.0000000000,1.0000000000,0.9846132030,1.0000000000,
    1.0000000000,1.0000000000,1.0000000000,1.0000000000,1.0000000000,1.0000000000,
    1.0000000000,1.0000000000,1.0000000000,1.0000000000,1.0000000000,1.0000000000,
    1.0000000000,1.0000000000,1.0000000000,1.0000000000,1.0015163710,1.0000000000,
    1.0000000000,1.0000000000,1.0000000000,0.7544739580,1.0000000000,1.0000000000,
    1.0000000000,1.0000000000,1.0000000000,1.0000000000,1.0000000000,1.0000000000,
    1.0000000000,0.8289671640,1.0000000000,1.0000000000,1.0000000000,1.0000000000,
    1.0000000000,1.0000000000,1.0000000000,1.0000000000,1.0000000000,1.0000000000,
    1.0000000000,1.0000000000,1.0000000000,1.0000000000,1.0000000000,1.0000000000,
    1.0000000000,1.0000000000,1.0000000000,1.6953583820,1.0000000000,1.0000000000,
    1.0000000000,0.9751877660,1.0000000000,1.0000000000,1.0000000000,1.0000000000,
    1.0000000000,1.0000000000,1.0000000000,1.0000000000,1.0000000000,1.0000000000,
    1.0000000000,1.0000000000,0.9997468470,1.0000000000,0.7956603920,1.0000000000,
    1.0000000000,1.0000000000,1.0000000000,1.0000000000,1.0000000000,1.0000000000,
    1.0000000000]

    beta_v=gerg_betamatrix_from_vector(vector_beta_v,x->1/x)

    vector_beta_T = [0.998098830,1.022624490,0.996336508,0.989680305,0.994174910,
    0.980315756,0.992127525,1.0,0.981844797,0.977431529,0.957473785,0.947716769,
    0.937777823,1.0,1.0,0.987411732,1.063333913,1.011090031,1.0,0.990954281,
    1.005894529,1.007671428,1.002677329,0.994515234,0.992868130,1.0,1.0,1.0,1.0,
    1.0,0.956379450,0.957934447,0.946134337,1.0,1.0,0.997190589,1.004692366,
    0.692868765,0.999069843,1.013871147,1.033620538,1.018171004,1.023339824,
    1.027000795,1.019180957,1.0,1.011806317,1.029690780,1.007688620,1.020028790,
    0.942320195,1.0,1.0,0.997372205,1.016034583,0.768377630,0.996512863,
    0.996199694,0.999130554,1.0,0.998688946,1.0,1.0,1.0,0.984068272,1.0,
    0.970918061,0.932969831,1.0,1.0,1.0,0.990197354,1.0,1.0,1.000310289,
    0.998012298,0.996484021,0.994364425,1.0,1.0,1.0,1.0,0.985331233,1.0,1.0,1.0,
    1.0,0.992573556,1.0,1.0,1.000077547,1.0,1.0,1.0,1.0,1.0,1.0,0.993688386,1.0,
    1.0,1.0,1.0,0.985962886,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
    0.974550548,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,0.962006651,1.0,
    1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,0.997641010,
    1.0,1.0,1.0,1.0,0.985891113,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,0.988937417,
    1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
    1.064818089,1.0,1.049594632,0.897162268,0.973091413,1.0,1.0,1.0,1.0,1.0,1.0,
    1.0,1.0,1.0,1.0,1.0,1.0,1.000023103,1.0,1.025536736,1.0,1.0,1.0,1.0,1.0,1.0,
    1.0,1.0]

    beta_T=gerg_betamatrix_from_vector(vector_beta_T,x->1/x)

    indice1 = [1,1,1,1,1,1,1,2,2,4,4,4,5,5,6]
    indice2 = [2,3,4,5,6,7,15,3,4,5,6,7,6,7,7]

    value_eq = 1:length(indice1)
    Aij_indices_mat = zeros(Int64,21,21)
    for i = 1:length(value_eq)
        Aij_indices_mat[indice1[i],indice2[i]] = value_eq[i]
     end

    Aij_indices_mat= Aij_indices_mat[xsel,xsel]
    Aij_indices2 = findall(!iszero, Aij_indices_mat)
    Aij_indices = [tuple(Aij_indices2[i][1],Aij_indices2[i][2],Aij_indices_mat[Aij_indices2[i]]) for i = 1:length(Aij_indices2)]

    fij=[1.0,1.0,1.0,1.0,1.0,0.771035405688,1.0,1.0,1.0,0.130424765150,0.281570073085,0.260632376098,
    0.312572600489e-1,-0.551609771024e-1,-0.551240293009e-1]

    dijk=[[1,4,1,2,2,2,2,2,3],
    [1,2,3,1,2,3],
    [3,4,1,2,2,2,2,2,2,3,3,3],
    [3,3,4,4,4,1,1,1,2],
    [2,3,1,1,1,2],
    [2,2,3,1,2,2],
    [1,3,3,4],
    [1,1,1,2,2,3,3,4,4,4],
    [1,1,1,2,2,3,3,4,4,4],
    [1,1,1,2,2,3,3,4,4,4],
    [1,1,1,2,2,3,3,4,4,4],
    [1,1,1,2,2,3,3,4,4,4],
    [1,1,1,2,2,3,3,4,4,4],
    [1,1,1,2,2,3,3,4,4,4],
    [1,1,1,2,2,3,3,4,4,4]]

    tijk=[[0.000,1.850,7.850,5.400,0.000,0.750,2.800,4.450,4.250],
    [2.600,1.950,0.000,3.950,7.950,8.000],
    [0.650,1.550,3.100,5.900,7.050,3.350,1.200,5.800,2.700,0.450,0.550,1.950],
    [1.850,3.950,0.000,1.850,3.850,5.250,3.850,0.200,6.500],
    [1.850,1.400,3.200,2.500,8.000,3.750],
    [0.000,0.050,0.000,3.650,4.900,4.450],
    [2.0,-1.0,1.75,1.4],
    [1.000,1.550,1.700,0.250,1.350,0.000,1.250,0.000,0.700,5.400],
    [1.000,1.550,1.700,0.250,1.350,0.000,1.250,0.000,0.700,5.400],
    [1.000,1.550,1.700,0.250,1.350,0.000,1.250,0.000,0.700,5.400],
    [1.000,1.550,1.700,0.250,1.350,0.000,1.250,0.000,0.700,5.400],
    [1.000,1.550,1.700,0.250,1.350,0.000,1.250,0.000,0.700,5.400],
    [1.000,1.550,1.700,0.250,1.350,0.000,1.250,0.000,0.700,5.400],
    [1.000,1.550,1.700,0.250,1.350,0.000,1.250,0.000,0.700,5.400],
    [1.000,1.550,1.700,0.250,1.350,0.000,1.250,0.000,0.700,5.400],
    ]

    nijk=[[-0.98038985517335e-2,0.42487270143005e-3,-0.34800214576142e-1,-0.13333813013896,
    -0.11993694974627e-1,0.69243379775168e-1,-0.31022508148249,0.24495491753226,0.22369816716981],
    [-0.10859387354942,0.80228576727389e-1,-0.93303985115717e-2,0.40989274005848e-1,
    -0.24338019772494,0.23855347281124],
    [-0.80926050298746e-3,-0.75381925080059e-3,-0.41618768891219e-1,-0.23452173681569,0.14003840584586,
    0.63281744807738e-1,-0.34660425848809e-1,-0.23918747334251,0.19855255066891e-2,0.61777746171555e1,
    -0.69575358271105e1,0.10630185306388e1],
    [0.13746429958576e-1,-0.74425012129552e-2,-0.45516600213685e-2,-0.54546603350237e-2,0.23682016824471e-2,
    0.18007763721438,-0.44773942932486,0.19327374888200e-1,-0.30632197804624],
    [0.28661625028399,-0.10919833861247,-0.11374032082270e1,0.76580544237358,
    0.42638000926819e-2,0.17673538204534],
    [-0.47376518126608,0.48961193461001,-0.57011062090535e-2,-0.19966820041320,-0.69411103101723,0.69226192739021],
    [-0.25157134971934,-0.62203841111983e-2,0.88850315184396e-1,-0.35592212573239e-1],
    [0.25574776844118e1,-0.79846357136353e1,0.47859131465806e1,-0.73265392369587,0.13805471345312e1,
    0.28349603476365,-0.49087385940425,-0.10291888921447,0.11836314681968,0.55527385721943e-4],
    [0.25574776844118e1,-0.79846357136353e1,0.47859131465806e1,-0.73265392369587,0.13805471345312e1,
    0.28349603476365,-0.49087385940425,-0.10291888921447,0.11836314681968,0.55527385721943e-4],
    [0.25574776844118e1,-0.79846357136353e1,0.47859131465806e1,-0.73265392369587,0.13805471345312e1,
    0.28349603476365,-0.49087385940425,-0.10291888921447,0.11836314681968,0.55527385721943e-4],
    [0.25574776844118e1,-0.79846357136353e1,0.47859131465806e1,-0.73265392369587,0.13805471345312e1,
    0.28349603476365,-0.49087385940425,-0.10291888921447,0.11836314681968,0.55527385721943e-4],
    [0.25574776844118e1,-0.79846357136353e1,0.47859131465806e1,-0.73265392369587,0.13805471345312e1,
    0.28349603476365,-0.49087385940425,-0.10291888921447,0.11836314681968,0.55527385721943e-4],
    [0.25574776844118e1,-0.79846357136353e1,0.47859131465806e1,-0.73265392369587,0.13805471345312e1,
    0.28349603476365,-0.49087385940425,-0.10291888921447,0.11836314681968,0.55527385721943e-4],
    [0.25574776844118e1,-0.79846357136353e1,0.47859131465806e1,-0.73265392369587,0.13805471345312e1,
    0.28349603476365,-0.49087385940425,-0.10291888921447,0.11836314681968,0.55527385721943e-4],
    [0.25574776844118e1,-0.79846357136353e1,0.47859131465806e1,-0.73265392369587,0.13805471345312e1,
    0.28349603476365,-0.49087385940425,-0.10291888921447,0.11836314681968,0.55527385721943e-4]
    ]

    etaijk=[
    [1.000,1.000,0.250,0.000,0.000,0.000,0.000],
    [1.000,0.500,0.000],
    [1.000,1.000,1.000,0.875,0.750,0.500,0.000,0.000,0.000,0.000],
    [0.250,0.250,0.000,0.000],
    [0.250,0.250,0.000,0.000],
    [1.000,1.000,0.875],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[]
    ]

    epsijk=[[0.5,0.5,0.5,0.5,0.5,0.5,0.5],
    [0.5,0.5,0.5],
    [0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5],
    [0.5,0.5,0.5,0.5],
    [0.5,0.5,0.5,0.5],
    [0.5,0.5,0.5],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[]
    ]
    betaijk=[[1.000,1.000,2.500,3.000,3.000,3.000,3.000],
    [1.000,2.000,3.000],
    [1.000,1.000,1.000,1.250,1.500,2.000,3.000,3.000,3.000,3.000],
    [0.750,1.000,2.000,3.000],
    [0.750,1.000,2.000,3.000],
    [1.000,1.000,1.250],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[]
    ]

    gammaijk=[[0.5,0.5,0.5,0.5,0.5,0.5,0.5],
    [0.5,0.5,0.5],
    [0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5],
    [0.5,0.5,0.5,0.5],
    [0.5,0.5,0.5,0.5],
    [0.5,0.5,0.5],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[]
    ]

    length_k_ijk = length(dijk)
    k_pol_ijk = Vector{Int64}(undef,length_k_ijk)
    k_exp_ijk = Vector{Int64}(undef,length_k_ijk)
    for i = 1:length_k_ijk
        k_exp_ijk[i] = length(etaijk[i])
        k_pol_ijk[i] = length(dijk[i]) - k_exp_ijk[i]
    end
    return GERG2008(N,acentric_factor[xsel],molecularWeight[xsel],criticalDensity[xsel],1 ./ criticalDensity[xsel],criticalTemperature[xsel],criticalPressure[xsel],
    ideal_iters,nr[:,xsel],zeta[:,xsel],n0ik[xsel],t0ik[xsel],d0ik[xsel],c0ik[xsel],
    k_pol_ik[xsel],k_exp_ik[xsel],gamma_v[xsel,xsel],gamma_T[xsel,xsel],beta_v[xsel,xsel],beta_T[xsel,xsel],
    Aij_indices,fij,dijk,tijk,nijk,etaijk,epsijk,betaijk,gammaijk, k_pol_ijk,k_exp_ijk)
end

function GERG2008(symsel::Vararg{Symbol})
    symvalues = [:CH4, :N2, :CO2, :C2, :C3, :nC4, :iC4, :nC5, :iC5, :C6, :C7, :C8, :C9, :C10,
    :H2, :O2, :CO, :H2O, :H2S, :He, :Ar]
    xsel = Vector{Int}()
    for i in symsel
        i in symvalues || throw(ArgumentError("invalid rank: $r"))
        push!(xsel,findall(isequal(i),symvalues)[1])

    end
    return GERG2008(xsel)
end

function _f0(model::GERG2008,rho,T,x)

    RR = 8.314472/8.314510
    common_type = promote_type(typeof(rho),typeof(T),eltype(x))

    res = zero(common_type)
    x0 = zero(common_type)  #for comparison
    ao1 =zero(common_type)
    ao2 =zero(common_type)
    ao_zero =zero(common_type)

    for i in model.ideal_iters[1]
        x[i] !=x0 && begin
        ao2 =ao_zero
        ao3 =ao_zero
        delta = rho/model.criticalDensity[i]
        tau = model.criticalTemperature[i]/T
        ao1 =  model.nr[1,i]+ model.nr[2,i]*tau + model.nr[3,i]*log(tau)
        ao2 = model.nr[4,i]*log(abs(sinh(model.zeta[1,i]*tau))) - model.nr[5,i]*log(cosh(model.zeta[2,i]*tau))+
        model.nr[6,i]*log(abs(sinh(model.zeta[3,i]*tau))) - model.nr[7,i]*log(cosh(model.zeta[4,i]*tau))

        ao3 = log(delta)
        a0 = RR*(ao1+ao2)+ao3
        res +=x[i]*(a0+log(x[i]))
    end
    end

    for i =model.ideal_iters[2]
        x[i] !=x0 && begin
        ao2 =ao_zero
        ao3 =ao_zero
        delta = rho/model.criticalDensity[i]
        tau = model.criticalTemperature[i]/T
        ao1 =  model.nr[1,i]+ model.nr[2,i]*tau + model.nr[3,i]*log(tau)
        ao2 = model.nr[4,i]*log(abs(sinh(model.zeta[1,i]*tau))) - model.nr[5,i]*log(cosh(model.zeta[2,i]*tau))+
        model.nr[6,i]*log(abs(sinh(model.zeta[3,i]*tau)))
        ao3 = log(delta)
        a0 = RR*(ao1+ao2)+ao3
        res +=x[i]*(a0+log(x[i]))
        end
    end

    for i in model.ideal_iters[3]
      x[i] !=x0 && begin
      ao2 =ao_zero
      ao3 =ao_zero
      delta = rho/model.criticalDensity[i]
      tau = model.criticalTemperature[i]/T
      ao1 =  model.nr[1,i]+ model.nr[2,i]*tau + model.nr[3,i]*log(tau)
      ao2 = model.nr[4,i]*log(abs(sinh(model.zeta[1,i]*tau))) - model.nr[5,i]*log(cosh(model.zeta[2,i]*tau))
      ao3 = log(delta)
      a0 = RR*(ao1+ao2)+ao3
      res +=x[i]*(a0+log(x[i]))
      end
    end

    for i in model.ideal_iters[4]
            x[i] !=x0 && begin
            ao2 =ao_zero
            ao3 =ao_zero
            delta = rho/model.criticalDensity[i]
            tau = model.criticalTemperature[i]/T
            ao1 =  model.nr[1,i]+ model.nr[2,i]*tau + model.nr[3,i]*log(tau)
            ao3 = log(delta)
            a0 = RR*(ao1)+ao3
            res +=x[i]*(a0+log(x[i]))
            end
    end

        return res
end

_gerg_asymetric_mix_rule(xi,xj,b)= b*(xi+xj)/(xi*b^2+xj)

function _delta(model::GERG2008,rho,T,x)
    rhor = inv(mixing_rule_asymetric(cubic_mean_rule,_gerg_asymetric_mix_rule,
    x, model.criticalVolume ,model.gamma_v,model.beta_v))
    return rho/rhor
end

function _tau(model::GERG2008,rho,T,x)
    Tr = mixing_rule_asymetric(geometric_mean_rule,_gerg_asymetric_mix_rule,
    x,model.criticalTemperature,model.gamma_T,model.beta_T)
    return Tr/T
end

function _fr1(model::GERG2008,delta,tau,x)
common_type = promote_type(typeof(delta),typeof(tau),eltype(x))
res =zero(common_type)
res0 =zero(common_type)
res1 =zero(common_type)
x0 = zero(common_type)

for i = 1:model.N
    x[i]!=x0 && begin
        res1=res0
        for k in 1:model.k_pol_ik[i]
            res1+=model.n0ik[i][k]*
                (delta^model.d0ik[i][k])*
                (tau^model.t0ik[i][k])
        end

        for k in (model.k_pol_ik[i]+1):(model.k_exp_ik[i]+model.k_pol_ik[i])

            res1+=model.n0ik[i][k]*
                (delta^model.d0ik[i][k])*
                (tau^model.t0ik[i][k])*
                exp(-delta^model.c0ik[i][k-model.k_pol_ik[i]])
        end
        res+= x[i]*res1

    end
end

return res
end

function _fr2(model::GERG2008,delta,tau,x)
    common_type = promote_type(typeof(delta),typeof(tau),eltype(x))
    res =zero(common_type)
    res0 =zero(common_type)
    x0 = zero(common_type)

    for kk in 1:length(model.Aij_indices)

        i1,i2,i0 = model.Aij_indices[kk]

            x[i1] != x0 && x[i2] !=x0 && begin
            res1 = res0
            for j=1:model.k_pol_ijk[i0]
                res1 +=model.nijk[i0][j]*
                (delta^model.dijk[i0][j])*
                (tau^model.tijk[i0][j])
            end

            for j in (model.k_pol_ijk[i0]+1):(model.k_pol_ijk[i0]+model.k_exp_ijk[i0])

                idx = j-model.k_pol_ijk[i0]

                res1 +=model.nijk[i0][j]*
                (delta^model.dijk[i0][j])*
                (tau^model.tijk[i0][j])*
            exp(-model.etaijk[i0][idx]*(delta-model.epsijk[i0][idx])^2 -
                model.betaijk[i0][idx]*(delta-model.gammaijk[i0][idx]))
            end
            res+=res1*x[i1]*x[i2]*model.fij[i0]

    end
end
return res
end

function core_ideal_helmholtz(model::GERG2008,v,T,x)
    rho = 1.0e-3/v
    R= Unitful.ustrip(Unitful.R)
    return R*T*_f0(model,rho,T,x)
end

function core_residual_helmholtz(model::GERG2008,v,T,x)
    rho = 1.0e-3/v
    R= Unitful.ustrip(Unitful.R)
    delta = _delta(model,rho,T,x)
    tau = _tau(model,rho,T,x)
    return R*T*(_fr1(model,delta,tau,x)+_fr2(model,delta,tau,x))
end

#an interface to stract properties stored in the model, the minimum are:
#core_helmholtz (the potential energy)
#molecular weight (for the transformations mass/molar)
#compounds_number

function core_helmholtz(model::GERG2008,v,T,x)
    rho = 1.0e-3/v
    R= Unitful.ustrip(Unitful.R)
    delta = _delta(model,rho,T,x)
    tau = _tau(model,rho,T,x)
    return R*T*(_f0(model,rho,T,x)+_fr1(model,delta,tau,x)+_fr2(model,delta,tau,x))
end

compounds_number(model::GERG2008)=model.N
molecular_weight(model::GERG2008) = model.molecularWeight

#functions for a critical model, maybe adding this as a trait?
critical_density(model::GERG2008)=1000.0 .* model.criticalDensity
critical_pressure(model::GERG2008) = model.criticalPressure
critical_temperature(model::GERG2008) = model.criticalTemperature
critical_volume(model::GERG2008) = 1 ./criticalDensity(model)
acentric_factor(model::GERG2008) = model.acentric_factor
covolumes(model::GERG2008) = 0.0778*model.criticalTemperature*Unitful.ustrip(Unitful.R) ./ model.criticalPressure


text:
function write_file(srange,
                    n::Integer,
                    sw::SlidingWindow=SlidingWindow(18),
                    file_name_suffix::String="slidingwindow$(sw.width)")
    file = "data/n$(n)$(file_name_suffix).txt"
    open(file,"a") do io
        write_io(io,srange,n,sw)
    end
end

function write_file(srange,
                    n::Integer,
                    ::MabeyBest,
                    file_name_suffix::String="mabeybest")
    file = "data/n$(n)$(file_name_suffix).txt"
    open(file,"a") do io
        write_io(io,srange,n,MabeyBest())
    end
end

function write_io(io::IO,srange,n::Integer,sw::SlidingWindow=SlidingWindow(18))
    for s in srange
        a = Vector{Int}()
        a,e = search(s,n,a,sw)
        sol = Solution(s,n,a)
        @assert e == err(sol)
        println(io,string(sol),",e=",string(e))
        flush(io)
    end
end

function write_io(io::IO,srange,n::Integer,::MabeyBest)
    for s in srange
        a,e = search(MabeyBest(),s,n)
        sol = Solution(s,n,a,false)
        @assert e == err(sol)
        println(io,string(sol),",e=",string(e))
        flush(io)
    end
end


text:
# Center and/or whiten data matrix
function center{T<:Real}(data::Matrix{T};
                         standardize::Bool=false,
                         bias::Int=0)
    num_samples, num_signals = size(data)
    @inbounds begin
        avgs = mean(data, 1)
        cntr = data .- avgs

        # Standardized moments: divide by the per-signal standard deviation
        if standardize
            stddev = (bias == 1) ? std(data, 1) :
                std(data, vec(avgs), num_samples, num_signals)'
            all(stddev .!= 0) ||
                error("Cannot standardize: variance = 0")
            cntr ./= stddev
        end
    end
    (cntr, num_samples, num_signals)
end

# Centering/whitening for weighted data
function center{T<:Real}(data::Matrix{T},
                         w::Vector{T};
                         standardize::Bool=false,
                         bias::Int=0)
    num_samples, num_signals = size(data)
    length(w) == num_samples || throw(DimensionMismatch("Inconsistent array lengths."))

    # Weighted detrending
    @inbounds begin
        avgs = mean(data, weights(w), 1)
        cntr = data .- avgs
        if standardize
            if bias == 1
                cntr ./= std(data, 1)
            else
                cntr ./= std(data, vec(avgs), num_samples, num_signals)'
            end
        end
    end
    (cntr, num_samples, num_signals)
end


text:
using LinearAlgebra
using SparseArrays
using Plots
using ForwardDiff

"This routine solves Poisson's equation using piecewise linear finite element methods.
The matrix is assembled globally. "

m = 50 # number of points
α = 1.0
β = 2.0

# define spatial grid
x = LinRange(-1,1,m+2) # x_0, x_1, ..., x_m, x_{m+1} = x_0
xint = x[1:end-1]
h = x[2]-x[1]

# construct FEM matrix
A = spdiagm(0=>2*ones(m), 1=>-ones(m-1), -1=>-ones(m-1))
A = (1/h)*A

# assume f(x) = 1.0 for now
b = h*ones(m)
b[1] += α/h
b[m] += β/h
u = A\b

plot(x,[α;u;β])


text:
import JuliaDB: IndexedTable
import ReadOnlyArrays: ReadOnlyArray

struct FilePaths
    corpus::String
    tanzil::String
end
FilePaths() = FilePaths(
    joinpath(@__DIR__, "..", "data", "quranic-corpus-morphology-0.4.txt"),
    joinpath(@__DIR__, "..", "data", "quran-uthmani-final.txt")
)

struct SimpleEncoder
    encode::Dict{Symbol,Symbol}
end
SimpleEncoder() = SimpleEncoder(SIMPLE_ENCODING)

struct ChapterLabel
    english::Array{String,1}
    arabic::Array{String,1}
end
ChapterLabel() = ChapterLabel(EN_CHAPTER_LABELS, AR_CHAPTER_LABELS)

struct QuranData
   path::FilePaths
end
QuranData() = QuranData(FilePaths())

struct CorpusRaw
    data::ReadOnlyArray{String,1,Array{String,1}}
end

struct TanzilRaw
    data::ReadOnlyArray{String,1,Array{String,1}}
end

mutable struct MetaData
    title::String
    author::String
    description::String
    website::String
    language::String
    year::String
    license::String
    version::String
end
MetaData(::Type{CorpusRaw}) = MetaData(
    "Quranic Arabic Corpus (morphology)",
    "Kais Dukes",
    "The Quranic Arabic Corpus includes syntactic and morphological\nannotation of the Quran, and builds on the verified Arabic text\ndistributed by the Tanzil project.",
    "http://corpus.quran.com/",
    "English",
    "2011",
    "GNU General Public License",
    "0.4"
)
MetaData(::Type{TanzilRaw}) = MetaData(
    "Tanzil Quran Text (Uthmani)",
    "Tanzil.net",
    "This copy of quran text is carefully produced, highly\nverified and continuously monitored by a group of specialists\nat Tanzil project.",
    "http://tanzil.net",
    "Arabic",
    "2008-2010",
    "Creative Commons Attribution 3.0",
    "1.0.2"
)

abstract type AbstractQuran end

struct CorpusData <: AbstractQuran
    data::IndexedTable
    meta::MetaData
end
CorpusData(data) = CorpusData(data, MetaData(CorpusRaw))

struct TanzilData <: AbstractQuran
    data::IndexedTable
    meta::MetaData
end
TanzilData(data) = TanzilData(data, MetaData(CorpusRaw))

struct Quran
    corpus::CorpusData
    tanzil::TanzilData
end

struct Chapter <: AbstractQuran
    data::IndexedTable
    numbers::Union{Int64,Array{Int64,1},UnitRange{Int64}}
    verses::Array{Int64,1}
    istanzil::Bool
end
Chapter(data, number, verses) = Chapter(data, number, verses, false)

struct Verse <: AbstractQuran
    data::IndexedTable
    numbers::Union{Int64,Array{Int64,1},UnitRange{Int64}}
    chapters::Union{Int64,Array{Int64,1},UnitRange{Int64}}
    words::Union{Int64,Array{Int64,1},UnitRange{Int64},Missing}
    istanzil::Bool
end

struct Word <: AbstractQuran
    data::IndexedTable
    numbers::Union{Int64,Array{Int64,1},UnitRange{Int64}}
    chapters::Union{Int64,Array{Int64,1},UnitRange{Int64}}
    verses::Union{Int64,Array{Int64,1},UnitRange{Int64}}
    parts::Union{Int64,Array{Int64,1},UnitRange{Int64}}
end

struct Part <: AbstractQuran
    data::IndexedTable
    numbers::Union{Int64,Array{Int64,1},UnitRange{Int64}}
    chapters::Union{Int64,Array{Int64,1},UnitRange{Int64}}
    verses::Union{Int64,Array{Int64,1},UnitRange{Int64}}
    words::Union{Int64,Array{Int64,1},UnitRange{Int64}}
end


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

makedocs(
	sitename="PrincipalMomentAnalysis.jl",
	modules = [PrincipalMomentAnalysis],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages = [
    	"Home" => "index.md",
    	"Tutorial" => "tutorial.md",
    	"Reference" => "reference.md",
    ]
)

deploydocs(
    repo = "github.com/PrincipalMomentAnalysis/PrincipalMomentAnalysis.jl.git",
)


text:
module AircraftGeometry

using Base.Math
using Base.Iterators
using DelimitedFiles
using StaticArrays
using CoordinateTransformations
using Rotations
using LinearAlgebra

using ..AeroMDAO: uniform_spacing, linear_spacing, sine_spacing, cosine_spacing, cosine_interp, splitat, adj3, slope, columns, fwdsum, fwddiv, fwddiff, weighted_vector, vectarray, Point2D, Panel2D, Panel3D, extend_yz, transform, panel_area, panel_normal, wetted_area

abstract type Aircraft end

export Aircraft

## Foil geometry
#==========================================================================================#

include("foil.jl")

# export Foil, kulfan_CST, naca4, camber_CST, paneller, read_foil, split_foil, foil_camthick, camthick_foil, camber_thickness, cosine_foil, camthick_to_CST, coords_to_CST, max_thickness_to_chord_ratio_location

## Fuselage geometry

include("fuselage.jl")

# export Fuselage, projected_area, length, cosine_distribution

## Wing geometry
#==========================================================================================#

aspect_ratio(span, area) = span^2 / area
mean_geometric_chord(span, area) = area / span
taper_ratio(root_chord, tip_chord) = tip_chord / root_chord
area(span, chord) = span * chord
mean_aerodynamic_chord(root_chord, taper_ratio) = (2/3) * root_chord * (1 + taper_ratio + taper_ratio^2)/(1 + taper_ratio)
y_mac(y, b, λ) = y + b / 2 * (1 + 2λ) / 3(1 + λ)
quarter_chord(chord) = 0.25 * chord

include("mesh_tools.jl")
include("halfwing.jl")
include("wing.jl")
include("mesh_wing.jl")

aspect_ratio(wing) = aspect_ratio(span(wing), projected_area(wing))

info(wing :: Union{Wing, HalfWing}) = [ span(wing), projected_area(wing), mean_aerodynamic_chord(wing), aspect_ratio(wing) ]

# export HalfWing, HalfWingSection, Wing, WingSection, mean_aerodynamic_chord, span, aspect_ratio, projected_area, taper_ratio, info, max_tbyc_sweeps, leading_edge, trailing_edge, wing_bounds, chop_leading_edge, chop_trailing_edge, chop_wing, paneller, panel_wing, mesh_horseshoes, mesh_wing, mesh_cambers, make_panels, vlmesh_wing, mean_aerodynamic_center, wetted_area, number_of_spanwise_panels, spanwise_spacing, coordinates

end

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
using Dash
using DashHtmlComponents
using DashCoreComponents

app = dash()

app.layout = html_div() do
    dcc_input(id="first", value = 1),
    dcc_input(id="second", value = 1),
    dcc_input(id="third", value = 1)
end

callback!(
    ClientsideFunction("clientside", "add1_prevent_at_11"),
    app, Output("second","value"), Input("first","value"), State("third","value")
    )

callback!(
    ClientsideFunction("clientside", "add1_prevent_at_11"),
    app, Output("third","value"), Input("second","value"), State("third","value")
    )

run_server(app)


text:
pyQt = pyimport("PyQt5.QtWidgets")
@pydef  mutable struct sd_QTButton <: pyQt.QPushButton
    function  __init__(self,name,parent=nothing)
        name1 = string( name ,"Julia")
        pyQt.QPushButton.__init__(self, name1 ,parent)
    end

end


text:
export get_latest_news, get_news_feeds, get_news_categories, get_news_feeds_and_categories
"""
    get_latest_news(; lang::String=EN, feeds::Array{String,1}, categories::Array{String,1}, ITs::Int64, sortOrder::String)

Gets the latest/popular news from CryptoCompare API.

Input parameters:
1. feeds: Specific news feeds to retrieve news from [ Min length - 1] [ Max length - 1000] [ Default - ALL_NEWS_FEEDS]
2. categories: Category of news articles to return [ Min length - 3] [ Max length - 1000] [ Default - ALL_NEWS_CATEGORIES]
3. excludeCategories: News article categories to exclude from results [ Min length - 3] [ Max length - 1000] [ Default - NO_EXCLUDED_NEWS_CATEGORIES]
4. lTs: Returns news before that timestamp [ Min - 0] [ Default - 0]
5. lang: Preferred language - English (EN) or Portuguese (PT) [ Min length - 1] [ Max length - 4] [ Default - EN]
6. sortOrder: The order to return news articles - latest or popular [ Min length - 1] [ Max length - 8] [ Default - latest]
7. extraParams: The name of your application (we recommend you send it) [ Min length - 1] [ Max length - 2000] [ Default - NotAvailable]
8. sign: If set to true, the server will sign the requests (by default we don't sign them), this is useful for usage in smart contracts [ Default - false]
"""
function get_latest_news(; lang::String="EN", feeds::Array{String,1}=Array{String,1}([]), categories::Array{String,1}=Array{String,1}([]), lTs::Int64=0, sortOrder::String="latest")
    url = build_news_url("latest_news_articles", lang=lang, feeds=feeds, categories=categories, lTs=lTs, sortOrder=sortOrder)

    try
        str = make_API_call(url)
        jobj = JSON.Parser.parse(str)
        datadict = jobj["Data"]

        df = DataFrame(imageurl=String[], id=Int[], guid=String[], published_on=DateTime[], body=String[], downvotes=Int[],
                       url=String[], upvotes=Int[], source=String[], title=String[], tags=String[], categories=String[], lang=String[],
                       source_info=Dict{String,Any}[])

        nrows = size(datadict)[1]

        for i in 1:nrows
            push!(df, datadict[i], cols=:union)
        end

        return df

    catch e
        return e
    end
end

"""
    get_news_feeds(; extraParams="", sign::Bool=false)

Gets the list of feeds from CryptoCompare API.

Input parameters:
1. extraParams: The name of your application (we recommend you send it) [ Min length - 1] [ Max length - 2000] [ Default - NotAvailable]
2. sign: If set to true, the server will sign the requests (by default we don't sign them), this is useful for usage in smart contracts [ Default - false]
"""
function get_news_feeds(; extraParams="", sign::Bool=false)

    url = build_news_url("feeds"; extraParams=extraParams, sign=sign)

    try
        str = make_API_call(url)
        jobj = JSON.Parser.parse(str)
        datadict = jobj
        df = DataFrame(key=String[], name=String[], lang=String[], img=String[])
        nrows = size(datadict)[1]

        for i in 1:nrows
            push!(df.key, datadict[i]["key"])
            push!(df.name, datadict[i]["name"])
            push!(df.lang, datadict[i]["lang"])
            push!(df.img, datadict[i]["img"])
        end

        return df

    catch e
        return e
    end
end

"""
    get_news_feeds_and_categories(; extraParams::String="", sign::Bool=false)

Gets the list of feeds and categories from CryptoCompare API.

Input parameters:
1. extraParams: The name of your application (we recommend you send it) [ Min length - 1] [ Max length - 2000] [ Default - NotAvailable]
2. sign: If set to true, the server will sign the requests (by default we don't sign them), this is useful for usage in smart contracts [ Default - false]
"""
function get_news_feeds_and_categories(; extraParams::String="", sign::Bool=false)
    url = build_news_url("feedsandcategories"; extraParams=extraParams, sign=sign)

    try
        str = make_API_call(url)
        jobj = JSON.Parser.parse(str)
        datadict = jobj["Data"]
        datadict_feeds = datadict["Feeds"]
        datadict_categories = datadict["Categories"]
        df_feeds = DataFrame(key=String[], name=String[], lang=String[], img=String[])
        df_categories = DataFrame(categoryName=String[], excludedPhrases=Array{String,1}[], includedPhrases=Array{String,1}[], wordsAssociatedWithCategory=Array{String,1}[])

        nrows = size(datadict_categories)[1]
        for i in 1:nrows
            push!(df_categories, datadict_categories[i], cols=:union)
        end

        nrows = size(datadict_feeds)[1]
        for i in 1:nrows
            push!(df_feeds, datadict_feeds[i], cols=:union)
        end

        return df_feeds, df_categories

    catch e
        return e
    end
end

"""
    get_news_categories(; extraParams::String="", sign::Bool=false)

Gets the list of news categories from CryptoCompare API.

Input parameters:
1. extraParams: The name of your application (we recommend you send it) [ Min length - 1] [ Max length - 2000] [ Default - NotAvailable]
2. sign: If set to true, the server will sign the requests (by default we don't sign them), this is useful for usage in smart contracts [ Default - false]
"""
function get_news_categories(; extraParams="", sign::Bool=false)
    url = build_news_url("categories"; extraParams=extraParams, sign=sign)

    try
        str = make_API_call(url)
        jobj = JSON.Parser.parse(str)
        datadict = jobj
        df = DataFrame(categoryName=String[], excludedPhrases=Array{String,1}[], includedPhrases=Array{String,1}[], wordsAssociatedWithCategory=Array{String,1}[])
        nrows = size(datadict)[1]

        for i in 1:nrows
            push!(df, datadict[i], cols=:union)
        end

        return df

    catch e
        return e
    end
end


text:
"""
    Viewer(args...; kwargs...)

    Call napari.Viewer(args...; kwargs...)
"""
Viewer(args...; kwargs...) = napari.Viewer(args...; kwargs...)

"""
    Viewer(f::Function, args...; kwargs...)

    Call f(Napari.Viewer(args...; kwargs...). Enables the following
    code pattern:

    ```julia
    Napari.Viewer() do v
        @add_image( v, zeros(1024,1024) )
        @add_iamge( v, Napari.astronaut() )
    end
    ```
"""
Viewer(f::Function, args...; kwargs...) = f(Viewer(args...; kwargs...))


text:
function calculate_num_error_prop(chain)
    num_error = get(chain, :numerical_error)
    return 100 * sum(num_error.numerical_error.data) / length(num_error.numerical_error.data)
end

#function calculate_experiment_results(@nospecialize(chain), @nospecialize(logger), include_vars)
function calculate_experiment_results(chain, logger, include_vars)
    logdir = logger.logdir
    has_num_error =  (:numerical_error in keys(chain))

    # Store the chain in several formats.
    # NOTE: For now, we save in JLS to start with. We comment out the next two lines for now.
    # JLSO.save(joinpath(logdir, "chain.jlso"), :chain => chain) # As JLSO, most robust
    # CSV.write(joinpath(logdir, "chain.csv"), DataFrame(chain)) # As CSV, in case everything else fails
    serialize(joinpath(logdir, "chain.jls"), chain)            # Fast but flimsy, no version control

    # summary statistics
    sum_stats = describe(chain[include_vars])
    param_names = sum_stats[1][:, 1]
    param_mean = sum_stats[1][:, 2]
    param_sd = sum_stats[1][:, 3]
    param_ess = sum_stats[1][:, 6]
    param_rhat = sum_stats[1][:, 7]
    param_esspersec = sum_stats[1][:, 8]

    CSV.write(
        joinpath(logdir, "sumstats.csv"),
        DataFrame(
            Parameter = param_names,
            Mean = param_mean,
            StdDev = param_sd,
            ESS = param_ess,
            Rhat = param_rhat,
            ESSpersec = param_esspersec,
            Num_error = has_num_error ? calculate_num_error_prop(chain) * ones(length(include_vars)) : missing # TODO: verify the "ones" logic here?
        )
    )
    # Save the preconditioner if one is available
    if :samplerstate in keys(chain.info)
        state = chain.info.samplerstate

        # Check if there's a Hamiltonian in the state
        if :hamiltonian in propertynames(state)
            h = state.hamiltonian

            # Get the metric
            mass = h.metric.M⁻¹
            mass_shape = size(mass)

            # If it's diagonal, make it dense
            m = if length(mass_shape) == 1
                diagm(mass)
            else
                mass
            end

            # Write it to disk
            writedlm(joinpath(logdir,"preconditioner.csv"), m, ',')
            end

        # cumulative averages
        for var in include_vars
            values = chain[var]
            cum_average = MCMCChains.cummean(values)
            writedlm(
                joinpath(logdir, "cumaverage_$(var).csv"),
                cum_average,
                ','
            )
        end
    end

    # Log the ESS/sec and rhat.  Nice to show as summary results from tensorboard
    for (i, name) = enumerate(param_names)
        TensorBoardLogger.log_value(
            logger,
            "$(name)_ess_per_sec",
            param_ess[i],
        )
        TensorBoardLogger.log_value(
            logger,
            "$(name)_rhat",
            param_rhat[i],
        )
    end

    # Make a succinct results JSON
    parameters = JSON.parsefile(joinpath(logdir, "parameters.json"))

    for (i, pname) in enumerate(param_names)
        parameters["$(pname)_mean"] = param_mean[i]
        parameters["$(pname)_std"] = param_sd[i]
        parameters["$(pname)_ess"] = param_ess[i]
        parameters["$(pname)_rhat"] = param_rhat[i]
        parameters["$(pname)_esspersec"] = param_esspersec[i]
    end

    # Add time difference if there is such information
    if :start_time in keys(chain.info)
        parameters["time_elapsed"] = chain.info.stop_time - chain.info.start_time
    end

    result_path = joinpath(logdir, "result.json")
    @info "Storing results at $result_path"
    open(result_path, "w") do io
        JSON.print(io, parameters)
    end
end

max_stars_repo_name: HighDimensionalEconLab/HMCExamples

max_issues_repo_name: HighDimensionalEconLab/HMCExamples

max_forks_repo_name: HighDimensionalEconLab/HMCExamples


text:
function fill_string(data, value=""; method="mostfrequent")
    nonempty = [x for x in data if x != "?"]

    if method == "sampling"
        data = [x == "?" ? sample(nonempty) : x for x in data]
    elseif method == "constant"
        data = [x == "?" ? value : x for x in data]
    elseif method == "mostfrequent"
        mostfrequent = frequentvalue(nonempty)
        data = [x == "?" ? mostfrequent : x for x in data]
    else
        throw("Undefined filling strategy")
    end
end

function fill_int(data, value=0; method="mean")
    nonmissing = [x for x in data if x != "?"]
    ints = [typeof(x) == String ? Int64(parse(Float64, x)) : Int64(x)
            for x in nonmissing]

    if method == "sampling"
        data = [x == "?" ? sample(ints) : x for x in data]
    elseif method == "mean"
        meanvalue = Int64(round(mean(ints)))
        data = [x == "?" ? meanvalue : x for x in data]
    elseif method == "constant"
        data = [x == "?" ? value : x for x in data]
    elseif method == "mostfrequent"
        mostfrequent = frequentvalue(nonempty)
        data = [x == "?" ? mostfrequent : x for x in data]
    else
        throw("Undefined filling strategy")
    end

    data = [typeof(x) == String ? Int64(parse(Float64, x)) : Int64(x)
            for x in data]
end

function fill_bin_category(data, value=""; method="mostfrequent")
    nonempty = [x for x in data if x != "?"]

    if method == "sampling"
        data = [x == "?" ? sample(nonempty) : x for x in data]
    elseif method == "constant"
        data = [x == "?" ? value : x for x in data]
    elseif method == "mostfrequent"
        mostfrequent = frequentvalue(nonempty)
        data = [x == "?" ? mostfrequent : x for x in data]
    else
        throw("Undefined filling strategy")
    end
end

function fill_category(data, value=""; method="mostfrequent")
    nonempty = [x for x in data if x != "?"]

    if method == "sampling"
        data = [x == "?" ? sample(nonempty) : x for x in data]
    elseif method == "constant"
        data = [x == "?" ? value : x for x in data]
    elseif method == "mostfrequent"
        mostfrequent = frequentvalue(nonempty)
        data = [x == "?" ? mostfrequent : x for x in data]
    else
        throw("Undefined filling strategy")
    end
end

function fill_date(data, value=""; method="mostfrequent")
    dates = [Date(x) for x in data if x != "?"]

    if method == "sampling"
       data = [x == "?" ? sample(dates) : x for x in data]
    elseif method == "constant"
        data = [x == "?" ? value : x for x in data]
    elseif method == "mostfrequent"
        mostfrequent = frequentvalue(dates)
        data = [x == "?" ? mostfrequent : x for x in data]
    else
        throw("Undefined filling strategy")
    end
end

function fill_timestamp(data, value=""; method="mostfrequent")
    timestamps = [DateTime(x) for x in data if x != "?"]

    if method == "sampling"
       data = [x == "?" ? sample(timestamps) : x for x in data]
    elseif method == "constant"
        data = [x == "?" ? value : x for x in data]
    elseif method == "mostfrequent"
        mostfrequent = frequentvalue(timestamps)
        data = [x == "?" ? mostfrequent : x for x in data]
    else
        throw("Undefined filling strategy")
    end
end

function fill_text(data, fname, value=""; method="constant")
    ids = [i for i in 1:length(data[fname]) if data[fname][i] != "?"]

    if method == "drop"
        data = Dict(fname => value[ids] for (fname, value) in data)
    elseif method == "constant"
        data[fname] = replace(data[fname], "?" => "")

        return data
    else
        throw("Undefined filling strategy")
    end
end

function fill_image(data, fname, value=""; method="drop")
    ids = [i for i in 1:length(data[fname]) if data[fname][i] != "?"]

    if method == "drop"
        data = Dict(fname => value[ids] for (fname, value) in data)
    elseif method == "constant"
        data[fname] = [x == "?" ? value : x for x in data[fname]]

        return data
    else
        throw("Undefined filling strategy")
    end
end

function remove_invalid(data, fname)
    data[fname] = [isvalid(docs) ? docs : string(c for c in docs if isvalid(c))
                   for docs in data[fname]]
    data
end

max_stars_repo_name: alperkesen/AutoML

max_issues_repo_name: alperkesen/AutoML


text:
using BlockArrays
using ControlSystems
using LinearAlgebra
using MatrixEquations
using PredictiveControl
using Test
using SafeTestsets

# Some utilities for testing
include( "../src/utilities.jl" )
include( "testUtils.jl" )

@testset "PredictiveControl" begin
    @testset "Problem Types" begin
        @safetestset "LTI CLQR" begin include( "typeTests.jl" ) end
    end

    @testset "Fully condensed problem" begin
        @safetestset "Initial Propagation Matrix" begin include( "./matrices/FullCondensed_InitialPropagation.jl" ) end
        @safetestset "Prediction Matrix"          begin include( "./matrices/FullCondensed_Prediction.jl" ) end
        @safetestset "Constraint Matrix"          begin include( "./matrices/FullCondensed_Constraints.jl" ) end
    end

    @testset "Analysis" begin
        @safetestset "Condition number bound" begin include( "analysis/conditioning_bound.jl" ) end
    end

    @testset "Solvers" begin
        @testset "Iteration Utilities" begin
            @safetestset "apply" begin include( "solvers/iterationUtils/apply.jl" )    end
            @safetestset "apply" begin include( "solvers/iterationUtils/stopping.jl" ) end
        end

        @testset "Fast Gradient Method" begin
            @safetestset "FGM Algorithm"        begin include( "solvers/fastGradientMethod/fgm_algorithm.jl" ) end
            @safetestset "FGM Δ - Cold Start"   begin include( "solvers/fastGradientMethod/fgm_colddelta.jl" ) end
            @safetestset "FGM upper iter bound" begin include( "solvers/fastGradientMethod/fgm_iterationbound.jl" ) end
        end
    end
end


text:
using Base: isexported

"""
    stdlibobject(name::Symbol)

If `name` is an export of Base or Core, return `Nullable{Any}(x)` where `x` is
the object itself. Otherwise, return `Nullable{Any}()`.
"""
function stdlibobject(ex::Symbol)
    if isexported(Base, ex) && isdefined(Base, ex)
        Nullable{Any}(getfield(Base, ex))
    elseif isexported(Core, ex) && isdefined(Core, ex)
        Nullable{Any}(getfield(Core, ex))
    else
        Nullable{Any}()
    end
end

"""
    parsetype(ctx::LintContext, ex::Expr) :: Type

Obtain a supertype of the type represented by `ex`.
"""
function parsetype(ctx::LintContext, ex)
    obj = abstract_eval(ctx, ex)
    if !isnull(obj) && isa(get(obj), Type)
        get(obj)
    elseif isexpr(ex, :curly)
        obj = abstract_eval(ctx, ex.args[1])
        if !isnull(obj) && isa(get(obj), Type) && get(obj) !== Union
            get(obj)
        else
            Any
        end
    else
        Any
    end
end

function guesstype(ex::Symbol, ctx::LintContext)
    result = lookup(ctx, ex)
    if isnull(result)
        Any  # conservative guess
    else
        get(result).typeactual
    end
end

function guesstype(ex::Expr, ctx::LintContext)::Type
    ex = ExpressionUtils.expand_trivial_calls(ex)

    if isexpr(ex, :tuple)
        ts = Type[]
        for a in ex.args
            push!(ts, guesstype(a, ctx))
        end
        return Tuple{ts...}
    end

    if isexpr(ex, :(::)) && length(ex.args) == 2
        return parsetype(ctx, ex.args[2])
    end

    if isexpr(ex, :block)
        return isempty(ex.args) ? Void : guesstype(ex.args[end], ctx)
    end

    if isexpr(ex, :return)
        tmp = guesstype(ex.args[1], ctx)
        return tmp
    end

    if isexpr(ex, :call)
        fn = ex.args[1]
        if any(x -> isexpr(x, :kw) || isexpr(x, :(...)), ex.args[2:end])
            # TODO: smarter way to deal with kw/vararg
            return Any
        end
        argtypes = map(x -> guesstype(x, ctx), ex.args[2:end])

        # infer return types of Base functions
        obj = abstract_eval(ctx, fn)
        type_argtypes = [isa(t, Type) ? t : Any for t in argtypes]
        if !isnull(obj)
            inferred = StaticTypeAnalysis.infertype(get(obj), type_argtypes)
            if inferred ≠ Any
                return inferred
            end
        end
    end

    if isexpr(ex, :macrocall)
        if ex.args[1] == Symbol("@sprintf")
            return String
        elseif ex.args[1] == Symbol("@compat")
            return guesstype(ex.args[2], ctx)
        end
    end

    if isexpr(ex, :curly)
        return Type
    end

    if isexpr(ex, :ref) # it could be a ref a[b] or an array Int[1,2,3], Vector{Int}[]
        if isexpr(ex.args[1], :curly) # must be a datatype, right?
            elt = abstract_eval(ctx, ex.args[1])
            if !isnull(elt) && isa(get(elt), Type)
                return Vector{get(elt)}
            else
                return Vector
            end
        end

        if isa(ex.args[1], Symbol)
            what = registersymboluse(ex.args[1], ctx)
            if what <: Type
                elt = abstract_eval(ctx, ex.args[1])
                if !isnull(elt) && isa(get(elt), Type)
                    return Vector{get(elt)}
                else
                    return Vector
                end
            end
        end
        # not symbol, or symbol but it refers to a variable
        partyp = guesstype(ex.args[1], ctx)
        if isa(partyp, Symbol)
            # we are in a context of a constructor of a new type, so it's
            # difficult to figure out the content
            return Any
        elseif partyp <: AbstractArray && !(partyp <: Range)
            eletyp = StaticTypeAnalysis.eltype(partyp)
            try
                nd = ndims(partyp) # This may throw if we couldn't infer the dimension
                tmpdim = nd - (length(ex.args)-1)
                if tmpdim < 0
                    if nd == 0 && ex.args[2] == 1 # ok to do A[1] for a 0-dimensional array
                        return eletyp
                    else
                        msg(ctx, :E436, ex, "more indices than dimensions")
                        return Any
                    end
                end

                for i in 2:length(ex.args)
                    if ex.args[i] == :(:) || isexpr(ex.args[i], :call) &&
                            ex.args[i].args[1] == :Colon
                        tmpdim += 1
                    end
                end
                if tmpdim != 0
                    return Array{eletyp, tmpdim} # is this strictly right?
                else
                    return eletyp
                end
            end
            return Any
        else
            argtypes = [guesstype(x, ctx) for x in ex.args]
            type_argtypes = [isa(t, Type) ? t : Any for t in argtypes]
            inferred = StaticTypeAnalysis.infertype(getindex, type_argtypes)
            if ctx.versionreachable(VERSION) && inferred == Union{}
                indtypes = if length(type_argtypes) == 1
                    "no indices"
                else
                    string("index types ", join(type_argtypes[2:end], ", "))
                end
                msg(ctx, :E522, ex,
                    string("indexing $(type_argtypes[1]) with ",
                           indtypes,
                           " is not supported"))
            end
            return inferred
        end
    end

    if isexpr(ex, :comparison)
        return Bool
    end

    # simple if statement e.g. test ? 0 : 1
    if isexpr(ex, :if) && 2 ≤ length(ex.args) ≤ 3
        tt = guesstype(ex.args[2], ctx)
        ft = if length(ex.args) == 3
            guesstype(ex.args[3], ctx)
        else
            Void
        end
        if tt == ft
            # we need this case because tt and ft might be symbols
            return tt
        elseif isa(tt, Type) && isa(ft, Type)
            return Union{tt, ft}
        else
            return Any
        end
    end

    if isexpr(ex, :(->))
        return Function
    end

    return Any
end

guesstype(ex, ctx::LintContext) = lexicaltypeof(ex)


text:
mutable struct GetIndex <: DiffableFunction
    ind
    grad_field :: GradField
    GetIndex(ind) = new(ind, GradField())
end

function forward(f::GetIndex, x)
    return x[f.ind...]
end

function backward(f::GetIndex, gy)
    x = f.grad_field.inputs[1]
    return GetIndexGrad(f.ind, size(x))(gy)
end

mutable struct GetIndexGrad <: DiffableFunction
    ind
    in_shape
    grad_field :: GradField
    GetIndexGrad(ind, in_shape) = new(ind, in_shape, GradField())
end

function add_at(arr::Vector, ind, val)
    arr[ind...] += val
    return arr
end

function add_at(arr::AbstractArray, ind, val)
    arr[ind...] .+= val
    return arr
end

function forward(f::GetIndexGrad, gy)
    gx = zeros(f.in_shape)
    return add_at(gx, f.ind, gy)
end

function backward(::GetIndexGrad, ggx)
    return GetIndex(ind, size(ggx))(ggx)
end

"""
    Base.getindex(x::Variable, ind...)
return `x.values[ind...]` as `Variable`.

# Examples
```julia-repl
julia> x = Variable(rand(2, 2))
name: nothing
values: [0.7050007249265509 0.5075375401538957; 0.9953109600473362 0.8447135817368259]
creator: User-Defined(nothing)

julia> x[1, 2]
name: nothing
values: 0.5075375401538957
creator: JITrench.GetIndex

julia> x[1, :]
name: nothing
values: [0.7050007249265509, 0.5075375401538957]
creator: JITrench.GetIndex
```
"""
Base.getindex(x::Variable, ind...) = GetIndex(ind)(x)

is_support(::typeof(Base.getindex)) = true
get_jt_struct(::typeof(Base.getindex)) = GetIndex


text:
function part1(input)
    sues = reduce(hcat, parse_line.(eachline(input)))
    true_sue = Dict("children" => 3,
                    "cats" => 7,
                    "samoyeds" => 2,
                    "pomeranians" => 3,
                    "akitas" => 0,
                    "vizslas" => 0,
                    "goldfish" => 5,
                    "trees" => 3,
                    "cars" => 2,
                    "perfumes" => 1)
    for (i, sue) in enumerate(sues)
        all(true_sue[key] == value for (key, value) in sue) && return i
    end
end

function part2(input)
    sues = reduce(hcat, parse_line.(eachline(input)))
    true_sue = Dict("children" => ==(3),
                    "cats" => >(7),
                    "samoyeds" => ==(2),
                    "pomeranians" => <(3),
                    "akitas" => ==(0),
                    "vizslas" => ==(0),
                    "goldfish" => <(5),
                    "trees" => >(3),
                    "cars" => ==(2),
                    "perfumes" => ==(1))
    for (i, sue) in enumerate(sues)
        all(true_sue[key](value) for (key, value) in sue) && return i
    end
end

function parse_line(line)
    facts = Dict{String, Int}()
    for part in split(line, ", ")
        pieces = split(part, ": ")
        facts[pieces[end - 1]] = parse(Int, pieces[end])
    end
    return facts
end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Author: Maur좩cio Carvalho Maciel
#* $Id$
#*--------------------------------------------------------------------
type sepComp
	sepComp()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of chemical components"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet stream",
				:PosX=>0,
				:PosY=>0.5001,
				:Symbol=>"_{Inlet}"
			)),
			streamPH(Dict{Symbol,Any}(
				:Brief=>"Overhead_Outlet stream",
				:PosX=>0.5,
				:PosY=>0,
				:Symbol=>"_{Overhead}"
			)),
			streamPH(Dict{Symbol,Any}(
				:Brief=>"Bottom_Outlet stream",
				:PosX=>0.5,
				:PosY=>1,
				:Symbol=>"_{Bottom}"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Rate of heat supply",
				:PosX=>1,
				:PosY=>0.7559,
				:Symbol=>"_{in}"
			)),
			fill(fraction(Dict{Symbol,Any}(
				:Brief=>"Fraction to Overhead_Outlet",
				:Symbol=>"\\phi_{overhead}"
			)),(NComp)),
			fill(fraction(Dict{Symbol,Any}(
				:Brief=>"Fraction to Bottom_Outlet",
				:Symbol=>"\\phi_{bottom}"
			)),(NComp)),
			[
				:(Overhead_Outlet.F*Overhead_Outlet.z = Inlet.F * Inlet.z*Overhead_Splits),
				:(sum(Bottom_Outlet.z) = 1),
				:(Overhead_Splits+Bottom_Splits = 1),
				:(sum(Overhead_Outlet.z) = 1),
				:(Overhead_Outlet.F*Overhead_Outlet.z + Bottom_Outlet.F*Bottom_Outlet.z = Inlet.F*Inlet.z),
				:(Inlet.F*Inlet.h = Overhead_Outlet.F*Overhead_Outlet.h + Bottom_Outlet.F*Bottom_Outlet.h - InletQ),
			],
			[
				"Composition Overhead","Bottom Composition Constraints","Sum Of Splits","Overhead Composition Constraints","Global Composition","Energy Balance",
			],
			[:PP,:NComp,],
			[:Inlet,:Overhead_Outlet,:Bottom_Outlet,:InletQ,:Overhead_Splits,:Bottom_Splits,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	Inlet::stream
	Overhead_Outlet::streamPH
	Bottom_Outlet::streamPH
	InletQ::power
	Overhead_Splits::Array{fraction}
	Bottom_Splits::Array{fraction}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export sepComp
function setEquationFlow(in::sepComp)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	#"Overhead Pressure"
	#	Overhead_Outlet.P = Inlet.P;
	#"Bottom Pressure"
	#	Bottom_Outlet.P = Inlet.P;
	#"Overhead Temperature"
	#	Overhead_Outlet.T = Inlet.T;
	#"Bottom Temperature"
	#	Bottom_Outlet.T = Inlet.T;

end
function atributes(in::sepComp,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/SepComp"
	fields[:Brief]="Model of a separator of components"
	fields[:Info]="== Assumptions ==
* thermodynamics equilibrium

== Specify ==
* the inlet stream
* (NComp) Overhead_Splits or (NComp) Bottom_Splits
* the Pressure and Temperature of the outlet streams
"
	drive!(fields,_)
	return fields
end
sepComp(_::Dict{Symbol,Any})=begin
	newModel=sepComp()
	newModel.attributes=atributes(newModel,_)
	newModel
end