text:
type BehaviorParticleBelief <: BehaviorBelief
    gen::UniformIDMMOBIL
    physical::MLPhysicalState
    particles::Vector{Vector{BehaviorModel}} # First index is the position in physical.cars
    weights::Vector{Vector{Float64}}   # Second index is the particle number
end

function rand(rng::AbstractRNG,
                b::BehaviorParticleBelief,
                s::MLState=MLState(b.physical, Array(CarState, length(b.physical.cars))))
    s.crashed = b.physical.crashed
    resize!(s.cars, length(b.physical.cars))
    for i in 1:length(s.cars)
        particle = sample(rng, b.particles[i], WeightVec(b.weights[i]))
        s.cars[i] = CarState(b.physical.cars[i], particle)
    end
    return s
end

function rand(rng::AbstractRNG,
                b::BehaviorParticleBelief,
                sample_noises::Vector,
                s::MLState=MLState(b.physical, Array(CarState, length(b.physical.cars))))

    s.crashed = b.physical.crashed
    resize!(s.cars, length(b.physical.cars))
    for i in 1:length(s.cars)
        particle = sample(rng, b.particles[i], WeightVec(b.weights[i]))
        nudged = clip(particle+sample_noises[i].*randn(rng,9), b.gen)
        s.cars[i] = CarState(b.physical.cars[i], nudged)
    end
    return s
end

# action(p::Policy, b::BehaviorParticleBelief) = action(p, b.physical)

function most_likely_state(b::BehaviorParticleBelief)
    s = MLState(b.physical, Array(CarState, length(b.physical.cars)))
    for i in 1:length(s.cars)
        ml_ind = indmax(b.weights[i])
        behavior = b.particles[i][ml_ind]
        s.cars[i] = CarState(b.physical.cars[i], behavior)
    end
    return s
end

param_means(b::BehaviorParticleBelief) = [sum(b.weights[i].*b.particles[i])/sum(b.weights[i]) for i in 1:length(b.particles)]
function param_stds(b::BehaviorParticleBelief)
    means = param_means(b)
    stds = Array(BehaviorModel, length(b.physical.cars))
    for i in 1:length(b.physical.cars)
        stds[i] = sqrt(sum(b.weights[i].*(b.particles[i].-means[i]).^2)/sum(b.weights[i]))
    end
    return stds
end

function weights_from_particles!(b::BehaviorParticleBelief,
                                 problem::NoCrashProblem,
                                 o::MLPhysicalState,
                                 particles,
                                 p::WeightUpdateParams)

    b.physical = o
    resize!(b.weights, length(o.cars))
    resize!(b.particles, length(o.cars))
    for i in 1:length(o.cars)
        # make sure we're not going to be allocating a bunch of memory in the loop
        if isdefined(b.particles, i)
            sizehint!(b.particles[i], length(particles))
            resize!(b.particles[i], 0)
        else
            b.particles[i] = Array(BehaviorModel, length(particles))
            resize!(b.particles[i], 0)
        end
        if isdefined(b.weights, i)
            sizehint!(b.weights[i], length(particles))
            resize!(b.weights[i], 0)
        else
            b.weights[i] = Array(Float64, length(particles))
            resize!(b.weights[i], 0)
        end
    end
    for (j, sp) in enumerate(particles)
        isp = 1
        io = 1
        while io <= length(o.cars) && isp <= length(sp.cars)
            co = o.cars[io]
            csp = sp.cars[isp]
            if co.id == csp.id
                if abs(co.x-csp.x) < 0.2*problem.dmodel.phys_param.lane_length
                    @assert isa(csp.behavior, IDMMOBILBehavior)
                    a = csp.behavior.p_idm.a
                    dt = problem.dmodel.phys_param.dt
                    veld = TriangularDist(csp.vel-a*dt/2.0, csp.vel+a*dt/2.0, csp.vel)
                    proportional_likelihood = Distributions.pdf(veld, co.vel)
                    if proportional_likelihood > 0.0
                        if co.y == csp.y
                            push!(b.particles[io], csp.behavior)
                            push!(b.weights[io], proportional_likelihood)
                        elseif abs(co.y - csp.y) < 1.0
                            push!(b.particles[io], csp.behavior)
                            push!(b.weights[io], p.wrong_lane_factor*proportional_likelihood)
                        end # if greater than one lane apart, do nothing
                    end
                end
                io += 1
                isp += 1
            elseif co.id < csp.id
                io += 1
            else
                @assert co.id > csp.id
                isp += 1
            end
        end
    end

    return b
end

type BehaviorParticleUpdater <: Updater{BehaviorParticleBelief}
    problem::Nullable{NoCrashProblem}
    nb_sims::Int
    p_resample_noise::Float64
    resample_noise_factor::Float64
    params::WeightUpdateParams
    rng::AbstractRNG
end
function set_problem!(u::BehaviorParticleUpdater, p::Union{POMDP,MDP})
    u.problem = Nullable{NoCrashProblem}(p)
end
function set_rng!(u::BehaviorParticleUpdater, rng::AbstractRNG)
    u.rng = rng
end

function update(up::BehaviorParticleUpdater,
                b_old::BehaviorParticleBelief,
                a::MLAction,
                o::MLPhysicalState,
                b_new::BehaviorParticleBelief=BehaviorParticleBelief(get(up.problem).dmodel.behaviors, o,
                                                    Array(Vector{BehaviorModel}, length(o.cars)),
                                                    Array(Vector{Float64}, length(o.cars))))

    gen = get(up.problem).dmodel.behaviors
    particles = Array(MLState, up.nb_sims)
    min_std = 0.001*IDMMOBILBehavior(gen.max_idm-gen.min_idm, gen.max_mobil-gen.min_mobil, 0)
    stds = param_stds(b_old)
    for i in 1:length(stds)
        stds[i] = max(stds[i], min_std)
    end
    for i in 1:up.nb_sims
        if rand(up.rng) < up.p_resample_noise
            s = rand(up.rng, b_old, up.resample_noise_factor.*stds)
        else
            s = rand(up.rng, b_old)
        end
        particles[i] = generate_s(get(up.problem), s, a, up.rng)
    end

    weights_from_particles!(b_new, get(up.problem), o, particles, up.params)

    for i in 1:length(o.cars)
        if isempty(b_new.weights[i])
            b_new.particles[i] = [rand(up.rng, b_new.gen) for i in 1:up.nb_sims]
            b_new.weights[i] = ones(up.nb_sims)
        end
    end

    return b_new
end

function initialize_belief(up::BehaviorParticleUpdater, distribution)
    gen = get(up.problem).dmodel.behaviors
    states = [rand(up.rng, distribution) for i in 1:up.nb_sims]
    particles = Array(Array{BehaviorModel}, length(first(states).cars))
    weights = Array(Array{Float64}, length(first(states).cars))
    for i in 1:length(first(states).cars)
        particles[i] = Array(BehaviorModel, length(states))
        weights[i] = ones(length(states))
        for (j,s) in enumerate(states)
            particles[i][j] = s.cars[i].behavior
        end
    end
    return BehaviorParticleBelief(gen, MLPhysicalState(first(states)), particles, weights)
end

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
export
    make

function make(env_id::String, env_params::Dict)
    try
        if env_id == "DeterministicSingleStepDebugEnv"
            return DeterministicSingleStepDebugEnv(env_params)
        elseif env_id == "NGSIMEnv"
            return NGSIMEnv(env_params)
        elseif env_id == "VectorizedNGSIMEnv"
            return VectorizedNGSIMEnv(env_params)
        elseif env_id == "MultiagentNGSIMEnv"
            return MultiagentNGSIMEnv(env_params)
	elseif env_id == "MultiagentNGSIMEnvVideoMaker"
	    	#println("RAUNAK make.jl says video maker")
		return MultiagentNGSIMEnvVideoMaker(env_params)
        else
            throw(ArgumentError("Invalid env_id: $(env_id)"))
        end
    catch e
        println("exception raised while making environment")
        println(backtrace(e))
        println(e)
        rethrow(e)
    end
end

ext: jl


text:
###
### LongSequence specific specializations of src/biosequence/indexing.jl
###
### This file is a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

# assumes `i` is positive and `bitsof(A)` is a power of 2

@inline function bitindex(seq::LongSequence, i::Integer)
    return bitindex(BitsPerSymbol(seq), encoded_data_eltype(seq), i + first(seq.part) - 1)
end

@inline function Base.getindex(seq::LongSequence, part::UnitRange)
    return LongSequence(seq, part)
end
@inline function Base.view(seq::LongSequence, part::UnitRange)
    return getindex(seq, part)
end

# Set a single sequence position to a single symbol value.
function Base.setindex!(seq::LongSequence, x, i::Integer)
    @boundscheck checkbounds(seq, i)
    orphan!(seq)
    return unsafe_setindex!(seq, x, i)
end

# Set multiple sequence positions to a single symbol value.
function Base.setindex!(seq::LongSequence{A}, x, locs::AbstractVector{<:Integer}) where {A}
    @boundscheck checkbounds(seq, locs)
    orphan!(seq)
    return unsafe_setindex!(seq, x, locs)
end

function Base.setindex!(seq::LongSequence{A}, x, locs::AbstractVector{Bool}) where {A}
    @boundscheck checkbounds(seq, locs)
    orphan!(seq)
    return unsafe_setindex!(seq, x, locs)
end

function Base.setindex!(seq::LongSequence{A},
                        other::LongSequence{A},
                        locs::AbstractVector{<:Integer}) where {A}
    @boundscheck checkbounds(seq, locs)
    checkdimension(other, locs)
    orphan!(seq)
    return unsafe_setindex!(seq, other, locs)
end

function Base.setindex!(seq::LongSequence{A},
                        other::LongSequence{A},
                        locs::AbstractVector{Bool}) where {A}
    @boundscheck checkbounds(seq, locs)
    checkdimension(other, locs)
    orphan!(seq)
    return unsafe_setindex!(seq, other, locs)
end

function Base.setindex!(seq::LongSequence{A},
                        other::LongSequence{A},
                        locs::UnitRange{<:Integer}) where {A}
    @boundscheck checkbounds(seq, locs)
    checkdimension(other, locs)
    return copyto!(seq, locs.start, other, 1, length(locs))
end

function Base.setindex!(seq::LongSequence{A},
			            other::LongSequence{A}, ::Colon) where {A}
    return setindex!(seq, other, 1:lastindex(seq))
end

function Base.setindex!(seq::LongSequence{A}, x, ::Colon) where {A}
    return setindex!(seq, x, 1:lastindex(seq))
end

# These are "unsafe" because of no bounds check and no orphan! call
@inline function unsafe_setindex!(seq::LongSequence{A}, x, i::Integer) where {A}
    bin = enc64(seq, x)
    return encoded_setindex!(seq, bin, i)
end

@inline function unsafe_setindex!(seq::LongSequence{A}, x, locs::AbstractVector{<:Integer}) where {A}
    bin = enc64(seq, x)
    for i in locs
        encoded_setindex!(seq, bin, i)
    end
    return seq
end

function unsafe_setindex!(seq::LongSequence{A}, x, locs::AbstractVector{Bool}) where {A}
    bin = enc64(seq, x)
    i = j = 0
    while true
        i = findnext(locs, i + 1)
        if i === nothing
            break
        end
        encoded_setindex!(seq, bin, i)
    end
    return seq
end

function unsafe_setindex!(seq::LongSequence{A}, other::LongSequence{A}, locs::AbstractVector{<:Integer}) where{A}
    for (i, x) in zip(locs, other)
        unsafe_setindex!(seq, x, i)
    end
    return seq
end

function unsafe_setindex!(seq::LongSequence{A},
                          other::LongSequence{A},
                          locs::AbstractVector{Bool}) where {A}
    i = j = 0
    while true
        i = findnext(locs, i + 1)
        if i === nothing
            break
        end
        unsafe_setindex!(seq, other[j += 1], i)
    end
    return seq
end

function enc64(::LongSequence{A}, x) where {A}
    #TODO: Resolve these two use cases of A().
    return UInt64(encode(A(), convert(eltype(A()), x)))
end

@inline function encoded_setindex!(seq::LongSequence{A},
				   bin::UInt64, i::Integer) where {A}
    j, r = bitindex(seq, i)
    data = encoded_data(seq)
    @inbounds data[j] = (bin << r) | (data[j] & ~(bindata_mask(seq) << r))
    return seq
end

#=
function Base.iterate(seq::LongSequence{A}, off::Int=first(seq.part)-1) where {A <: Alphabet}
    off == last(seq.part) && return nothing
    bps = bits_per_symbol(A())
    @inbounds chunk = seq.data[(off >>> index_shift(BitsPerSymbol(A()))) + 1]
    shift = (unsigned(off) * bps) & 63
    encoding = (chunk >>> shift) & bitmask(A())
    return decode(A(), encoding), off + 1
end

=#


text:
# Copyright 2015, 2016 Martin Holters
# See accompanying license file.

using ACME

function diodeclipper(::Type{Circuit})
    j_in = voltagesource()
    r1 = resistor(1e3)
    c1 = capacitor(47e-9)
    d1 = diode(is=1e-15)
    d2 = diode(is=1.8e-15)
    j_out = voltageprobe()

    circ = Circuit()

    connect!(circ, j_in[:+], r1[1])
    connect!(circ, j_in[:-], :gnd)
    connect!(circ, r1[2], c1[1], d1[:+], d2[:-], j_out[:+])
    connect!(circ, :gnd, c1[2], d1[:-], d2[:+], j_out[:-])

    return circ
end

diodeclipper{T<:DiscreteModel}(::Type{T}=DiscreteModel; fs=44100) =
    T(diodeclipper(Circuit), 1/fs)


text:
"""
    ModelFrame(formula, data; model=StatisticalModel, contrasts=Dict())

Wrapper that encapsulates a `FormulaTerm`, schema, data table, and model type.

This wrapper encapsulates all the information that's required to transform data
of the same structure as the wrapped data frame into a model matrix (the
`FormulaTerm`), as well as the information about how that formula term was
instantiated (the schema and model type)

Creating a model frame involves first extracting the [`schema`](@ref) for the
data (using any contrasts provided as hints), and then applying that schema with
[`apply_schema`](@ref) to the formula in the context of the provided model type.

# Constructors

```julia
ModelFrame(f::FormulaTerm, data; model::Type{M} = StatisticalModel, contrasts::Dict = Dict())
```

# Fields

* `f::FormulaTerm`: Formula whose left hand side is the *response* and right hand
  side are the *predictors*.
* `schema::Any`: The schema that was applied to generate `f`.
* `data::D`: The data table being modeled.  The only restriction is that `data`
  is a table (`Tables.istable(data) == true`)
* `model::Type{M}`: The type of the model that will be fit from this model frame.

# Examples

```julia
julia> df = (x = 1:4, y = 5:8)
julia> mf = ModelFrame(@formula(y ~ 1 + x), df)
```
"""
mutable struct ModelFrame{D,M}
    f::FormulaTerm
    schema
    data::D
    model::Type{M}
end

## copied from DataFrames:
function _nonmissing!(res, col)
    # workaround until JuliaLang/julia#21256 is fixed
    eltype(col) >: Missing || return
    res .&= .!ismissing.(col)
end

function missing_omit(d::T) where T<:ColumnTable
    nonmissings = trues(length(first(d)))
    for col in d
        _nonmissing!(nonmissings, col)
    end

    rows = findall(nonmissings)
    d_nonmissing =
        NamedTuple{Tables.names(T)}(tuple((copyto!(similar(col,
                                                           Base.nonmissingtype(eltype(col)),
                                                           length(rows)),
                                                   view(col, rows)) for col in d)...))
    d_nonmissing, nonmissings
end

missing_omit(data::T, formula::AbstractTerm) where T<:ColumnTable =
    missing_omit(NamedTuple{tuple(termvars(formula)...)}(data))

function ModelFrame(f::FormulaTerm, data::ColumnTable;
                    model::Type{M}=StatisticalModel, contrasts=Dict{Symbol,Any}()) where M
    data, _ = missing_omit(data, f)

    sch = schema(f, data, contrasts)
    f = apply_schema(f, sch, M)

    ModelFrame(f, sch, data, model)
end

ModelFrame(f::FormulaTerm, data; model=StatisticalModel, contrasts=Dict{Symbol,Any}()) =
    ModelFrame(f, columntable(data); model=model, contrasts=contrasts)

StatsBase.modelmatrix(f::FormulaTerm, data; kwargs...) = modelmatrix(f.rhs, data; kwargs...)

function StatsBase.modelmatrix(t::Union{AbstractTerm, TupleTerm}, data;
                               hints=Dict{Symbol,Any}(), mod::Type{M}=StatisticalModel) where M
    Tables.istable(data) ||
        throw(ArgumentError("expected data in a Table, got $(typeof(data))"))
    t = has_schema(t) ? t : apply_schema(t, schema(t, data, hints), M)
    modelcols(collect_matrix_terms(t), columntable(data))
end
function StatsBase.response(f::FormulaTerm, data;
                            hints=Dict{Symbol,Any}(),
                            mod::Type{M}=StatisticalModel) where M
    Tables.istable(data) ||
        throw(ArgumentError("expected data in a Table, got $(typeof(data))"))
    f = has_schema(f) ? f : apply_schema(f, schema(f, data, hints), M)
    modelcols(f.lhs, columntable(data))
end

StatsBase.modelmatrix(mf::ModelFrame; data=mf.data) = modelcols(mf.f.rhs, data)
StatsBase.response(mf::ModelFrame; data=mf.data) = modelcols(mf.f.lhs, data)

StatsBase.coefnames(mf::ModelFrame) = vectorize(coefnames(mf.f.rhs))

"""
    setcontrasts!(mf::ModelFrame; kwargs...)
    setcontrasts!(mf::ModelFrame, contrasts::Dict{Symbol})

Update the contrasts used for coding categorical variables in
[`ModelFrame`](@ref) in place.  This is accomplished by computing a new schema
based on the provided contrasts and the `ModelFrame`'s data, and applying it to
the `ModelFrame`'s `FormulaTerm`.

Note that only the `ModelFrame` itself is mutated: because `AbstractTerm`s are
immutable, any changes will produce a copy.

"""
setcontrasts!(mf::ModelFrame; kwargs...) = setcontrasts!(mf, Dict(kwargs))
function setcontrasts!(mf::ModelFrame, contrasts::Dict{Symbol})
    # TODO: don't consume the whole table again if it's not needed
    new_schema = schema([term(k) for k in keys(contrasts)], mf.data, contrasts)

    # warn of no-op for keys that dont correspond to known terms from old schema
    unknown_keys = [k for k in keys(new_schema) if !haskey(mf.schema, k)]
    if !isempty(unknown_keys)
        unknown_keys_str = join(unknown_keys, ", ", " and ")
        @warn "setcontrasts! for terms " * unknown_keys_str *
            " has no effect since they are not found in original schema"
    end

    # apply only the re-mapped terms
    mf.f = apply_schema(mf.f, new_schema, mf.model)
    merge!(mf.schema, new_schema)
    mf
end

"""
    ModelMatrix(mf::ModelFrame)

Convert a `ModelFrame` into a numeric matrix suitable for modeling

# Fields

* `m::AbstractMatrix{<:AbstractFloat}`: the generated numeric matrix
* `assign::Vector{Int}` the index of the term corresponding to each column of `m`.

# Constructors

```julia
ModelMatrix(mf::ModelFrame)
# Specify the type of the resulting matrix (default Matrix{Float64})
ModelMatrix{T <: AbstractMatrix{<:AbstractFloat}}(mf::ModelFrame)
```
"""
mutable struct ModelMatrix{T <: AbstractMatrix{<:AbstractFloat}}
    m::T
    assign::Vector{Int}
end

Base.size(mm::ModelMatrix, dim...) = size(mm.m, dim...)

asgn(f::FormulaTerm) = asgn(f.rhs)
asgn(mt::MatrixTerm) = asgn(mt.terms)
asgn(t) = mapreduce(((i,t), ) -> i*ones(width(t)),
                    append!,
                    enumerate(vectorize(t)),
                    init=Int[])

function ModelMatrix{T}(mf::ModelFrame) where T<:AbstractMatrix{<:AbstractFloat}
    mat = modelmatrix(mf)
    ModelMatrix(convert(T, mat), asgn(mf.f))
end

ModelMatrix(mf::ModelFrame) = ModelMatrix{Matrix{Float64}}(mf)


text:
using Documenter, MiniFB
using Literate

out_path = joinpath(@__DIR__,"src/example")
rm("gzexamples"; force=true, recursive=true)
rm(out_path; force=true, recursive=true )

config = Dict{Any, Any}("documenter"=>false, "execute"=>false, "credit"=>false)
Literate.markdown(joinpath(@__DIR__, "../example/Noise.jl"), out_path; config=config)
Literate.markdown(joinpath(@__DIR__,"../example/Plasma.jl"), out_path; config=config)
Literate.markdown(joinpath(@__DIR__,"../example/Image.jl"), out_path; config=config)

makedocs(;
    modules=[MiniFB],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Examples" => Any[
            "Noise" => "example/Noise.md",
            "Plasma" => "example/Plasma.md",
            "Image" => "example/Image.md",
        ],
        "API" => "api.md",
        "Generated Code" => "gen.md"
    ],
    sitename="MiniFB",
    authors="Avik Sengupta"
)


text:
"""
    CR_gammaReg_inv(topredict::String, varnames::Vector{Symbol}, train::DataFrame, test::DataFrame)

Tests and ranks (RMSD) all possible set of variables fitting a gamma regression (inverse link).

# Arguments

- `topredict::String`: The name of the column of the dataframe of interest that is to be predicted.
- `varnames::Vector{Symbol}`: The explanatory variables to include in the comprehensive research.
- `train::DataFrame`: The dataset used to train the model.
- `test::DataFrame`: The dataset used to test the model.
"""
function CR_gammaReg_inv(topredict::String, varnames::Vector{Symbol}, train::DataFrame, test::DataFrame)
    diagnostic = DataFrame(Variables = Array{Symbol}[], Deviance = Float64[], RMSD = Float64[])
    model_var = collect(combinations(varnames))

    for i in 1:length(model_var)
        vars = model_var[i]
        model = glm(Term(Symbol(topredict)) ~ sum(Term.(vars)), train, Gamma(), InverseLink())
        predictions = convert(Vector{Float64}, GLM.predict(model, test))
        RMSD = rmsd(predictions, test[:, Symbol(topredict)])
        push!(diagnostic, [model_var[i], deviance(model), RMSD])
    end

    sort!(diagnostic, [:RMSD, :Deviance])
end

"""
    stepRMSD_gammaReg_inv(topredict::String, varnames::Vector{Symbol}, initial_list::Vector{Symbol}, train::DataFrame, test::DataFrame, verbose::Bool=true)

Selects the best set of explanatory variables after performing a forward feature selection with gamma regression (inverse link) selecting on the RMSD.

# Arguments

- `topredict::String`: The name of the column of the dataframe of interest that is to be predicted.
- `varnames::Vector{Symbol}`: The explanatory variables to include in the regression.
- `initial_list::T where T<:Vector{Any}`: list of features to start with.
- `train::DataFrame`: The dataset used to train the model.
- `test::DataFrame`: The dataset used to test the model.
- `verbose::Bool=true`: whether to print the sequence of inclusions and exclusions. Default = true.
"""

function stepRMSD_gammaReg_inv(topredict::String, varnames::Vector{Symbol}, initial_list::Union{Vector{Symbol}, Vector{Any}}, train::DataFrame, test::DataFrame, verbose::Bool=true)
    included = initial_list
    while true
        changed=false

        # forward step
        excluded = setdiff(varnames, included)
        new_RMSD = DataFrame(var = [excluded[i] for i in 1:length(excluded)])
        vec_RMSD = [] #Will store the RMSD scores to then add them to the new_RMSD df
        for i=1:length(excluded)
            concat = vcat(included,excluded[i]) #concatenates included and one variable from excluded
            model = glm(Term(Symbol(topredict)) ~ sum(Term.(concat)), train, Gamma(), InverseLink())
            predictions = convert(Vector{Float64},GLM.predict(model, test))
            RMSD = rmsd(predictions, test[:, Symbol(topredict)])
            push!(vec_RMSD, RMSD)
        end
        new_RMSD.RMSD = vec_RMSD

        best_RMSD = minimum(new_RMSD.RMSD)

        if included != []
            model = glm(Term(Symbol(topredict)) ~ sum(Term.(included)), train, Gamma(), InverseLink())
            predictions = convert(Vector{Float64},GLM.predict(model, test))
            RMSD = rmsd(predictions, test[:, Symbol(topredict)])
            ex_RMSD = RMSD
        else
            ex_RMSD = 10^10
        end
        if  best_RMSD < ex_RMSD
            best_feature = new_RMSD[new_RMSD.RMSD .== best_RMSD,:var][1] #[1] to extract a Symbol and not a Vector
            push!(included, best_feature)
            changed=true
            if verbose
                println("Add $(best_feature) with RMSD $(best_RMSD)")
            end
        end

        if changed == false
            println("RMSD: $(ex_RMSD)")
            return convert(Vector{Symbol}, included)
            break
        end
    end
    return convert(Vector{Symbol}, included)
end

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
function hh2d(A::HMat)
    if A.blockType == DENSE
        return A.DMat;
    elseif A.blockType == LOWRANK
        return A.UMat*A.VMat';
    elseif A.blockType == HMAT
        D = zeros(A.height,A.width);
        hoffset = 0;
        for i = 1:4
            woffset = 0;
            for j = 1:4
                C = A.childHMat[i,j];
                D[hoffset+(1:C.height),woffset+(1:C.width)] = hh2d(C);
                woffset += C.width;
                if j == 4
                    hoffset += C.height;
                end
            end
        end
        return D;
    end
end

max_stars_repo_name: YingzhouLi/HMat

max_issues_repo_name: YingzhouLi/HMat


text:
export Limit, MeasurementDependency, CoilDescription, minimalHeader

function GeneralParameters(str::String)
  return GeneralParameters( parse_string(str) )
end

function addToDict!(params::Dict, e::XMLElement, paramName::String, ::Type{T},
                    dictParamName=paramName) where T
  u = get_elements_by_tagname(e, paramName)
  if !isempty(u)
    if length(u) == 1
      params[dictParamName] = customParse(T,u[1])
    else
      tmp = T[]
      for l=1:length(u)
        push!(tmp, customParse(T,u[l]))
      end
      params[dictParamName] = tmp
    end
  end
end

struct Limit
  minimum::Int
  maximum::Int
  center::Int
end

struct MeasurementDependency
  dependencyType::String
  measurementID::String
end

struct CoilDescription
  number::Int
  name::String
end

customParse(::Type{T}, node) where T = parse(T,content(node))
customParse(::Type{String}, node) where T = content(node)
function customParse(::Type{Vector{T}}, node) where T
  x = content(get_elements_by_tagname(node, "x")[1])
  y = content(get_elements_by_tagname(node, "y")[1])
  z = content(get_elements_by_tagname(node, "z")[1])
  return [parse(T,x),parse(T,y),parse(T,z)]
end
function customParse(::Type{MeasurementDependency}, node)
  a = content(get_elements_by_tagname(node, "dependencyType")[1])
  b = content(get_elements_by_tagname(node, "measurementID")[1])
  return MeasurementDependency(a, b)
end
function customParse(::Type{CoilDescription}, node)
  a = content(get_elements_by_tagname(node, "coilNumber")[1])
  b = content(get_elements_by_tagname(node, "coilName")[1])
  return CoilDescription(parse(Int,a), b)
end
function customParse(::Type{Limit}, node)
  a = content(get_elements_by_tagname(node, "minimum")[1])
  b = content(get_elements_by_tagname(node, "maximum")[1])
  c = content(get_elements_by_tagname(node, "center")[1])
  return Limit(parse(Int,a), parse(Int,b), parse(Int,c))
end

function GeneralParameters(xdoc::XMLDocument)
    params = Dict{String,Any}()

    # SubjectInformation
    e = get_elements_by_tagname(LightXML.root(xdoc),"subjectInformation")
    if !isempty(e)
      addToDict!(params, e[1], "patientName", String)
      addToDict!(params, e[1], "patientWeight_kg", Float64)
      addToDict!(params, e[1], "patientID", String)
      addToDict!(params, e[1], "patientBirthdate", String)
      addToDict!(params, e[1], "patientGender", String)
    end

    # StudyInformation
    e = get_elements_by_tagname(LightXML.root(xdoc),"studyInformation")
    if !isempty(e)
      addToDict!(params, e[1], "studyDate", String)
      addToDict!(params, e[1], "studyDate", String)
      addToDict!(params, e[1], "studyID", String)
      addToDict!(params, e[1], "accessionNumber", Int)
      addToDict!(params, e[1], "referringPhysicianName", String)
      addToDict!(params, e[1], "studyDescription", String)
      addToDict!(params, e[1], "studyInstanceUID", String)
    end

    # MeasurementInformation
    e = get_elements_by_tagname(LightXML.root(xdoc),"measurementInformation")
    if !isempty(e)
      addToDict!(params, e[1], "measurementID", String)
      addToDict!(params, e[1], "seriesDate", String)
      addToDict!(params, e[1], "seriesTime", String)
      addToDict!(params, e[1], "patientPosition", String) # non-optional?
      addToDict!(params, e[1], "initialSeriesNumber", Int)
      addToDict!(params, e[1], "protocolName", String)
      addToDict!(params, e[1], "seriesDescription", String)
      addToDict!(params, e[1], "seriesInstanceUIDRoot", String)
      addToDict!(params, e[1], "frameOfReferenceUID", String)

      addToDict!(params, e[1], "measurementDependency", MeasurementDependency)
      addToDict!(params, e[1], "referencedImageSequence", String)
    end

    # AcquisitionSystemInformation
    e = get_elements_by_tagname(LightXML.root(xdoc),"acquisitionSystemInformation")
    if !isempty(e)
      addToDict!(params, e[1], "systemVendor", String)
      addToDict!(params, e[1], "systemModel", String)
      addToDict!(params, e[1], "systemFieldStrength_T", Float64)
      addToDict!(params, e[1], "relativeReceiverNoiseBandwidth", Float64)
      addToDict!(params, e[1], "receiverChannels", Int)
      addToDict!(params, e[1], "institutionName", String)
      addToDict!(params, e[1], "stationName", String)
      addToDict!(params, e[1], "coilLabel", CoilDescription)
    end

    # ExperimentalConditions
    e = get_elements_by_tagname(LightXML.root(xdoc),"experimentalConditions")
    if !isempty(e)
      addToDict!(params, e[1], "H1resonanceFrequency_Hz", Int)
    end

    # Encoding
    e = get_elements_by_tagname(LightXML.root(xdoc),"encoding")
    if !isempty(e)
      if !isempty(e[1]["encodedSpace"])
        params["encodedSize"] = customParse(Vector{Int}, e[1]["encodedSpace"][1]["matrixSize"][1])
        params["encodedFOV"] = customParse(Vector{Float64}, e[1]["encodedSpace"][1]["fieldOfView_mm"][1])
      end
      if !isempty(e[1]["reconSpace"])
        params["reconSize"] = customParse(Vector{Int}, e[1]["reconSpace"][1]["matrixSize"][1])
        params["reconFOV"] = customParse(Vector{Float64}, e[1]["reconSpace"][1]["fieldOfView_mm"][1])
      end

      d = e[1]["encodingLimits"]
      if !isempty(d)
        addToDict!(params, d[1], "kspace_encoding_step_0", Limit, "enc_lim_kspace_encoding_step_0")
        addToDict!(params, d[1], "kspace_encoding_step_1", Limit, "enc_lim_kspace_encoding_step_1")
        addToDict!(params, d[1], "kspace_encoding_step_2", Limit, "enc_lim_kspace_encoding_step_2")
        addToDict!(params, d[1], "average", Limit, "enc_lim_average")
        addToDict!(params, d[1], "slice", Limit, "enc_lim_slice")
        addToDict!(params, d[1], "contrast", Limit, "enc_lim_contrast")
        addToDict!(params, d[1], "phase", Limit, "enc_lim_phase")
        addToDict!(params, d[1], "repetition", Limit, "enc_lim_repetition")
        addToDict!(params, d[1], "set", Limit, "enc_lim_set")
        addToDict!(params, d[1], "segment", Limit, "enc_lim_segment")
      end

      addToDict!(params, e[1], "trajectory", String)

      d = e[1]["trajectoryDescription"]
      if !isempty(d)
        y = Dict{String,Any}()
        y["identifier"] = content(d[1]["identifier"][1])
        for q in d[1]["userParameterLong"]
          y[content(q["name"][1])] = parse(Int,content(q["value"][1]))
        end
        for q in d[1]["userParameterDouble"]
          y[content(q["name"][1])] = parse(Float64,content(q["value"][1]))
        end
        for q in d[1]["userParameterString"]
          y[content(q["name"][1])] = content(q["value"][1])
        end
        if !isempty(d[1]["comment"])
          y["comment"] = content(d[1]["comment"][1])
        end
        params["trajectoryDescription"] = y
      end

      d = e[1]["parallelImaging"]
      if !isempty(d)
        a = parse(Int,content(d[1]["accelerationFactor"][1]["kspace_encoding_step_1"][1]))
        b = parse(Int,content(d[1]["accelerationFactor"][1]["kspace_encoding_step_2"][1]))
        params["accelerationFactor"] = (a,b)

        addToDict!(params, d[1], "calibrationMode", String)
        addToDict!(params, d[1], "interleavingDimension", String)
      end
      addToDict!(params, e[1], "echoTrainLength", Int)
    end

    # Support wrongly created files (which are downloaded in the tests):
    e = get_elements_by_tagname(LightXML.root(xdoc),"parallelImaging")
    if !isempty(e)
      a = parse(Int,content(e[1]["accelerationFactor"][1]["kspace_encoding_step_1"][1]))
      b = parse(Int,content(e[1]["accelerationFactor"][1]["kspace_encoding_step_2"][1]))
      params["accelerationFactor"] = (a,b)

      addToDict!(params, e[1], "calibrationMode", String)
      addToDict!(params, e[1], "interleavingDimension", String)
    end

    # SequenceParameters
    e = get_elements_by_tagname(LightXML.root(xdoc),"sequenceParameters")
    if !isempty(e)
      addToDict!(params, e[1], "TR", Float64)
      addToDict!(params, e[1], "TE", Float64)
      addToDict!(params, e[1], "TI", Float64)
      addToDict!(params, e[1], "flipAngle_deg", Float64)
      addToDict!(params, e[1], "sequence_type", String)
      addToDict!(params, e[1], "echo_spacing", Float64)
    end

    # waveformInformation
    e = get_elements_by_tagname(LightXML.root(xdoc),"waveformInformation")
    if !isempty(e)
      addToDict!(params, e[1], "waveformName", Float64)
      addToDict!(params, e[1], "waveformType", Float64)

      if !isempty(e[1]["userParameters"])
          d = e[1]["userParameters"]
          y = Dict{String,Any}()
          for q in d["userParameterLong"]
            y[content(q["name"][1])] = parse(Int,content(q["value"][1]))
          end
          for q in d["userParameterDouble"]
            y[content(q["name"][1])] = parse(Float64,content(q["value"][1]))
          end
          for q in d["userParameterString"]
            y[content(q["name"][1])] = content(q["value"][1])
          end
          params["waveformUserParameters"] = y
      end
    end

    # UserParameters
    e = get_elements_by_tagname(LightXML.root(xdoc),"userParameters")
    if !isempty(e)
        d = e[1]
        y = Dict{String,Any}()
        for q in d["userParameterLong"]
          y[content(q["name"][1])] = parse(Int,content(q["value"][1]))
        end
        for q in d["userParameterDouble"]
          y[content(q["name"][1])] = parse(Float64,content(q["value"][1]))
        end
        for q in d["userParameterString"]
          y[content(q["name"][1])] = content(q["value"][1])
        end
        params["userParameters"] = y
    end

    return params
end

function generateGroup(params, paramVec, node, groupName, prefix="")
  if any( [haskey(params,prefix*n) for n in paramVec] )
    xs = new_child(node, groupName)
    for p in paramVec
      p_ = prefix*p
      if haskey(params,p_)
        insertNode(xs, p, params[p_])
      end
    end
  end
end

function insertNode(xs, paramName::String, param::T) where T
  xsp = new_child(xs, paramName)
  add_text(xsp, string(param))
end

function insertNode(xs, paramName::String, param::Vector{T}) where T <: Number

  if paramName in ["TE", "TR"]
    for p in param
      insertNode(xs, paramName, p)
    end
  else
    xsp = new_child(xs, paramName)
    insertNode(xsp, "x", param[1])
    insertNode(xsp, "y", param[2])
    insertNode(xsp, "z", param[3])
  end
end

function insertNode(xs, paramName::String, param::Vector{T}) where T
  for p in param
    insertNode(xs, paramName, p)
  end
end

function insertNode(xs, paramName::String, param::Tuple{T,T}) where T <: Integer
  xsp = new_child(xs, "accelerationFactor")
  insertNode(xsp, "kspace_encoding_step_1", param[1])
  insertNode(xsp, "kspace_encoding_step_2", param[2])
end

function insertNode(xs, paramName::String, param::MeasurementDependency)
  xsp = new_child(xs, paramName)
  insertNode(xsp, "dependencyType", param.dependencyType)
  insertNode(xsp, "measurementID", param.measurementID)
end

function insertNode(xs, paramName::String, param::CoilDescription)
  xsp = new_child(xs, paramName)
  insertNode(xsp, "coilNumber", param.number)
  insertNode(xsp, "coilName", param.name)
end

function insertNode(xs, paramName::String, param::Limit)
  xsp = new_child(xs, paramName)
  insertNode(xsp, "minimum", param.minimum)
  insertNode(xsp, "maximum", param.maximum)
  insertNode(xsp, "center", param.center)
end

function insertNode(xs, paramName::String, param_::Dict{String,Any})
  xsp = new_child(xs, paramName)
  param = deepcopy(param_)
  if haskey(param, "comment")
    insertNode(xsp, "comment", param["comment"])
    delete!(param,"comment")
  end
  if haskey(param, "identifier")
    insertNode(xsp, "identifier", param["identifier"])
    delete!(param,"identifier")
  end
  for (key,val) in param
    if typeof(val) <: Integer
      xspp = new_child(xsp, "userParameterLong")
      insertNode(xspp, "name", key)
      insertNode(xspp, "value", val)
    elseif typeof(val) <: AbstractFloat
      xspp = new_child(xsp, "userParameterDouble")
      insertNode(xspp, "name", key)
      insertNode(xspp, "value", val)
    elseif typeof(val) <: AbstractString
      xspp = new_child(xsp, "userParameterString")
      insertNode(xspp, "name", key)
      insertNode(xspp, "value", val)
    end
  end
end

function GeneralParametersToXML(params::Dict{String,Any})

  xdoc = XMLDocument()
  xroot = create_root(xdoc, "ismrmrdHeader")
  set_attributes(xroot, Dict{String,String}("xmlns"=>"http://www.ismrm.org/ISMRMRD",
                         "xmlns:xsi"=>"http://www.w3.org/2001/XMLSchema-instance",
                         "xmlns:xs"=>"http://www.w3.org/2001/XMLSchema",
                         "xsi:schemaLocation"=>"http://www.ismrm.org/ISMRMRD ismrmrd.xsd"))

  p = ["patientName","patientWeight_kg", "patientID", "patientBirthdate", "patientGender"]
  generateGroup(params, p, xroot, "subjectInformation")

  p = ["studyDate", "studyDate", "studyID", "accessionNumber","referringPhysicianName",
       "studyDescription", "studyInstanceUID"]
  generateGroup(params, p, xroot, "studyInformation")

  p = ["measurementID", "seriesDate", "seriesTime", "patientPosition", "initialSeriesNumber",
       "protocolName", "seriesDescription", "seriesInstanceUIDRoot", "frameOfReferenceUID",
       "measurementDependency", "referencedImageSequence"]
  generateGroup(params, p, xroot, "measurementInformation")

  p = ["systemVendor", "systemModel", "systemFieldStrength_T",
       "relativeReceiverNoiseBandwidth", "receiverChannels",
       "institutionName", "stationName", "coilLabel"]
  generateGroup(params, p, xroot, "acquisitionSystemInformation")

  p = ["H1resonanceFrequency_Hz"]
  generateGroup(params, p, xroot, "experimentalConditions")

  pEnc = ["encodedSize", "encodedFOV"]
  pRecon = ["reconSize", "reconFOV"]
  pEncLim = ["kspace_encoding_step_0", "kspace_encoding_step_1",
             "kspace_encoding_step_2", "average", "slice",
             "contrast", "phase", "repetition", "set",
             "segment"]

  pElse = ["trajectory", "trajectoryDescription", "echoTrainLength"]

  pPI = ["accelerationFactor", "calibrationMode", "interleavingDimension"]

  p = union(pEnc, pRecon, pEncLim, pElse, pPI)

  if any( [haskey(params,n) for n in p] )
    xs = new_child(xroot, "encoding")

    if haskey(params, "encodedSize") || haskey(params, "encodedFOV")
      xsp = new_child(xs, "encodedSpace")
      insertNode(xsp, "matrixSize", params["encodedSize"])
      insertNode(xsp, "fieldOfView_mm", params["encodedFOV"])
    end

    if haskey(params, "reconSize") || haskey(params, "reconFOV")
      xsp = new_child(xs, "reconSpace")
      insertNode(xsp, "matrixSize", params["reconSize"])
      insertNode(xsp, "fieldOfView_mm", params["reconFOV"])
    end

    generateGroup(params, pEncLim, xs, "encodingLimits", "enc_lim_")

    for n in pElse
      if haskey(params, n)
        insertNode(xs, n, params[n])
      end
    end

    generateGroup(params, pPI, xs, "parallelImaging")
  end

  p = ["TR", "TE", "TI", "flipAngle_deg", "sequence_type", "echo_spacing"]
  generateGroup(params, p, xroot, "sequenceParameters")

  p = ["waveformName", "waveformType", "waveformUserParameters"]
  generateGroup(params, p, xroot, "waveformInformation")

  p = ["userParameters"]
  generateGroup(params, p, xroot, "userParameters")

  return string(xdoc)
end

function minimalHeader(encodingSize::NTuple{3,Int},fov::NTuple{3,AbstractFloat};f_res::Integer=1,tr_name::AbstractString="cartesian",numChannels::Int=1)
  params = Dict{String,Any}()
  params["H1resonanceFrequency_Hz"] = f_res
  params["encodedSize"] = collect(encodingSize)
  params["encodedFOV"] = collect(fov)
  params["trajectory"] = tr_name
  params["receiverChannels"] = numChannels

  return params
end


text:
module DictDotDot
export dotset!
get_path_segments(path::String) = split(path, ".")

function dotset!(root::AbstractDict, dottedkeys::String, value)
    dict = root
    paths = get_path_segments(dottedkeys)
    parentdict = dict
    for i = 1:length(paths)
        p = paths[i]
        if i == length(paths)
            dict[p] = value
        else
            if haskey(dict, p)
            else
                if length(keys(dict)) == 0
                    dict[p] = Dict()
                else
                    newchild = Dict(String(p) => Dict())
                    dict = merge(Dict{Any,Any}(dict), newchild)
                    if (i > 1)
                        parentdict[paths[i-1]] = dict
                    end
                end
            end
            parentdict = dict
            dict = dict[p]
        end
    end
    root
end
end # module

max_stars_repo_name: nampdn/DictDotDot

max_issues_repo_name: nampdn/DictDotDot

max_forks_repo_name: nampdn/DictDotDot


text:
function Sigmoid()
    return Sigmoid(())
end

function Sigmoid(arg0::jdouble, arg1::jdouble)
    return Sigmoid((jdouble, jdouble), arg0, arg1)
end

function value(obj::Sigmoid, arg0::DerivativeStructure)
    return jcall(obj, "value", DerivativeStructure, (DerivativeStructure,), arg0)
end

function value(obj::Sigmoid, arg0::jdouble)
    return jcall(obj, "value", jdouble, (jdouble,), arg0)
end


text:
using UAParser, YAML, Base.Test

#Tests look a little funky, because YAML reads in file as Dict
#UAParser uses composite type; tests build Dict for comparison
#Test files where tests are commented out due to not implementing jsParseBits from ua-parser Python library
#These user-agents look to be older, may not still be out in wild

#Test 1: Validation of parsedevice
test_device = YAML.load(open(joinpath(dirname(@__FILE__), "data", "test_device.yaml")));

for test_case in test_device["test_cases"]
  @test test_case["family"] == parsedevice(test_case["user_agent_string"]).family
end

#Test 2: Validation of parseos
test_os = YAML.load(open(joinpath(dirname(@__FILE__), "data", "test_user_agent_parser_os.yaml")));

for value in test_os["test_cases"]
  @test Dict{Any, Any}("major" => value["major"],
      "minor" => value["minor"],
      "patch" => value["patch"],
      "patch_minor" => value["patch_minor"],
      "family" => value["family"]) ==
      Dict{Any, Any}("major" => parseos(value["user_agent_string"]).major,
       "minor" => parseos(value["user_agent_string"]).minor,
       "patch" => parseos(value["user_agent_string"]).patch,
       "patch_minor" => parseos(value["user_agent_string"]).patch_minor,
       "family" => parseos(value["user_agent_string"]).family)
end

#Test 3: Additional validation of parseos
test_os_2 = YAML.load(open(joinpath(dirname(@__FILE__), "data", "additional_os_tests.yaml")));

for value in test_os_2["test_cases"]
@test Dict{Any, Any}("major" => value["major"], "minor" => value["minor"], "patch" => value["patch"], "patch_minor" => value["patch_minor"], "family" => value["family"]) ==
     Dict{Any, Any}("major" => parseos(value["user_agent_string"]).major,
       "minor" => parseos(value["user_agent_string"]).minor,
       "patch" => parseos(value["user_agent_string"]).patch,
       "patch_minor" => parseos(value["user_agent_string"]).patch_minor,
       "family" => parseos(value["user_agent_string"]).family)
end

#Test 4: Validation of parseuseragent
test_ua = YAML.load(open(joinpath(dirname(@__FILE__), "data", "test_user_agent_parser.yaml")));

for value in test_ua["test_cases"]
@test Dict{Any, Any}("major" => value["major"], "minor" => value["minor"], "patch" => value["patch"], "family" => value["family"]) ==
    Dict{Any, Any}("major" => parseuseragent(value["user_agent_string"]).major,
     "minor" => parseuseragent(value["user_agent_string"]).minor,
     "patch" => parseuseragent(value["user_agent_string"]).patch,
     "family" => parseuseragent(value["user_agent_string"]).family)
end

#Test 5: Additional validation of parseuseragent
test_ua_2 = YAML.load(open(joinpath(dirname(@__FILE__), "data", "firefox_user_agent_strings.yaml")));

for value in test_ua_2["test_cases"]
@test Dict{Any, Any}("major" => value["major"], "minor" => value["minor"], "patch" => value["patch"], "family" => value["family"]) ==
     Dict{Any, Any}("major" => parseuseragent(value["user_agent_string"]).major,
     "minor" => parseuseragent(value["user_agent_string"]).minor,
     "patch" => parseuseragent(value["user_agent_string"]).patch,
     "family" => parseuseragent(value["user_agent_string"]).family)

end


text:
using Test
using Dates
using LinearAlgebra
using MPI
using Printf
using StaticArrays
using UnPack

using ClimateMachine
using ClimateMachine.ConfigTypes
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.MPIStateArrays
using ClimateMachine.ODESolvers
using ClimateMachine.GenericCallbacks
using ClimateMachine.Atmos
using ClimateMachine.BalanceLaws
using ClimateMachine.Orientations
using ClimateMachine.VariableTemplates
using Thermodynamics
using ClimateMachine.TurbulenceClosures
using ClimateMachine.VTK

import ClimateMachine.BalanceLaws: source, prognostic_vars

using CLIMAParameters
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

import CLIMAParameters
# Assume zero reference temperature
CLIMAParameters.Planet.T_0(::EarthParameterSet) = 0

if !@isdefined integration_testing
    const integration_testing = parse(
        Bool,
        lowercase(get(ENV, "JULIA_CLIMA_INTEGRATION_TESTING", "false")),
    )
end

include("mms_solution_generated.jl")

import Thermodynamics: total_specific_enthalpy
using ClimateMachine.Atmos

total_specific_enthalpy(ts::PhaseDry{FT}, e_tot::FT) where {FT <: Real} =
    zero(FT)

function mms2_init_state!(problem, bl, state::Vars, aux::Vars, localgeo, t)
    (x1, x2, x3) = localgeo.coord
    state.ρ = ρ_g(t, x1, x2, x3, Val(2))
    state.ρu = SVector(
        U_g(t, x1, x2, x3, Val(2)),
        V_g(t, x1, x2, x3, Val(2)),
        W_g(t, x1, x2, x3, Val(2)),
    )
    state.energy.ρe = E_g(t, x1, x2, x3, Val(2))
end

struct MMSSource{N} <: TendencyDef{Source} end

prognostic_vars(::MMSSource{N}) where {N} = (Mass(), Momentum(), Energy())

function source(::Mass, s::MMSSource{N}, m, args) where {N}
    @unpack aux, t = args
    x1, x2, x3 = aux.coord
    return Sρ_g(t, x1, x2, x3, Val(N))
end
function source(::Momentum, s::MMSSource{N}, m, args) where {N}
    @unpack aux, t = args
    x1, x2, x3 = aux.coord
    return SVector(
        SU_g(t, x1, x2, x3, Val(N)),
        SV_g(t, x1, x2, x3, Val(N)),
        SW_g(t, x1, x2, x3, Val(N)),
    )
end
function source(::Energy, s::MMSSource{N}, m, args) where {N}
    @unpack aux, t = args
    x1, x2, x3 = aux.coord
    return SE_g(t, x1, x2, x3, Val(N))
end

function mms3_init_state!(problem, bl, state::Vars, aux::Vars, localgeo, t)
    (x1, x2, x3) = localgeo.coord
    state.ρ = ρ_g(t, x1, x2, x3, Val(3))
    state.ρu = SVector(
        U_g(t, x1, x2, x3, Val(3)),
        V_g(t, x1, x2, x3, Val(3)),
        W_g(t, x1, x2, x3, Val(3)),
    )
    state.energy.ρe = E_g(t, x1, x2, x3, Val(3))
end

# initial condition

function test_run(mpicomm, ArrayType, dim, topl, warpfun, N, timeend, FT, dt)

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
        meshwarp = warpfun,
    )

    physics = AtmosPhysics{FT}(
        param_set;
        ref_state = NoReferenceState(),
        turbulence = ConstantDynamicViscosity(FT(μ_exact), WithDivergence()),
        moisture = DryModel(),
    )
    if dim == 2
        problem = AtmosProblem(
            boundaryconditions = (InitStateBC(),),
            init_state_prognostic = mms2_init_state!,
        )
        model = AtmosModel{FT}(
            AtmosLESConfigType,
            physics;
            problem = problem,
            orientation = NoOrientation(),
            source = (MMSSource{2}(),),
        )
    else
        problem = AtmosProblem(
            boundaryconditions = (InitStateBC(),),
            init_state_prognostic = mms3_init_state!,
        )
        model = AtmosModel{FT}(
            AtmosLESConfigType,
            physics;
            problem = problem,
            orientation = NoOrientation(),
            source = (MMSSource{3}(),),
        )
    end
    show_tendencies(model)

    dg = DGModel(
        model,
        grid,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )

    Q = init_ode_state(dg, FT(0))
    Qcpu = init_ode_state(dg, FT(0); init_on_cpu = true)
    @test euclidean_distance(Q, Qcpu) < sqrt(eps(FT))

    lsrk = LSRK54CarpenterKennedy(dg, Q; dt = dt, t0 = 0)

    eng0 = norm(Q)
    @info @sprintf """Starting
    norm(Q₀) = %.16e""" eng0

    # Set up the information callback
    starttime = Ref(now())
    cbinfo = GenericCallbacks.EveryXWallTimeSeconds(60, mpicomm) do (s = false)
        if s
            starttime[] = now()
        else
            energy = norm(Q)
            @info @sprintf(
                """Update
                simtime = %.16e
                runtime = %s
                norm(Q) = %.16e""",
                ODESolvers.gettime(lsrk),
                Dates.format(
                    convert(Dates.DateTime, Dates.now() - starttime[]),
                    Dates.dateformat"HH:MM:SS",
                ),
                energy
            )
        end
    end

    solve!(Q, lsrk; timeend = timeend, callbacks = (cbinfo,))
    # solve!(Q, lsrk; timeend=timeend, callbacks=(cbinfo, cbvtk))

    # Print some end of the simulation information
    engf = norm(Q)
    Qe = init_ode_state(dg, FT(timeend))

    engfe = norm(Qe)
    errf = euclidean_distance(Q, Qe)
    @info @sprintf """Finished
    norm(Q)                 = %.16e
    norm(Q) / norm(Q₀)      = %.16e
    norm(Q) - norm(Q₀)      = %.16e
    norm(Q - Qe)            = %.16e
    norm(Q - Qe) / norm(Qe) = %.16e
    """ engf engf / eng0 engf - eng0 errf errf / engfe
    errf
end

let
    ClimateMachine.init()
    ArrayType = ClimateMachine.array_type()

    mpicomm = MPI.COMM_WORLD

    polynomialorder = 4
    base_num_elem = 4

    expected_result = [
        1.6931876910307017e-01 5.4603193051929394e-03 2.3307776694542282e-04
        3.3983777728925593e-02 1.7808380837573065e-03 9.176181458773599e-5
    ]
    lvls = integration_testing ? size(expected_result, 2) : 1

    @testset "mms_bc_atmos" begin
        for FT in (Float64,) #Float32)
            result = zeros(FT, lvls)
            for dim in 2:3
                for l in 1:lvls
                    if dim == 2
                        Ne = (
                            2^(l - 1) * base_num_elem,
                            2^(l - 1) * base_num_elem,
                        )
                        brickrange = (
                            range(FT(0); length = Ne[1] + 1, stop = 1),
                            range(FT(0); length = Ne[2] + 1, stop = 1),
                        )
                        topl = BrickTopology(
                            mpicomm,
                            brickrange,
                            periodicity = (false, false),
                        )
                        dt = 1e-2 / Ne[1]
                        warpfun =
                            (x1, x2, _) -> begin
                                (x1 + sin(x1 * x2), x2 + sin(2 * x1 * x2), 0)
                            end

                    elseif dim == 3
                        Ne = (
                            2^(l - 1) * base_num_elem,
                            2^(l - 1) * base_num_elem,
                        )
                        brickrange = (
                            range(FT(0); length = Ne[1] + 1, stop = 1),
                            range(FT(0); length = Ne[2] + 1, stop = 1),
                            range(FT(0); length = Ne[2] + 1, stop = 1),
                        )
                        topl = BrickTopology(
                            mpicomm,
                            brickrange,
                            periodicity = (false, false, false),
                        )
                        dt = 5e-3 / Ne[1]
                        warpfun =
                            (x1, x2, x3) -> begin
                                (
                                    x1 +
                                    (x1 - 1 / 2) * cos(2 * π * x2 * x3) / 4,
                                    x2 + exp(sin(2π * (x1 * x2 + x3))) / 20,
                                    x3 + x1 / 4 + x2^2 / 2 + sin(x1 * x2 * x3),
                                )
                            end
                    end
                    timeend = 1
                    nsteps = ceil(Int64, timeend / dt)
                    dt = timeend / nsteps

                    @info (ArrayType, FT, dim, nsteps, dt)
                    result[l] = test_run(
                        mpicomm,
                        ArrayType,
                        dim,
                        topl,
                        warpfun,
                        polynomialorder,
                        timeend,
                        FT,
                        dt,
                    )
                    @test result[l] ≈ FT(expected_result[dim - 1, l])
                end
                if integration_testing
                    @info begin
                        msg = ""
                        for l in 1:(lvls - 1)
                            rate = log2(result[l]) - log2(result[l + 1])
                            msg *= @sprintf(
                                "\n  rate for level %d = %e\n",
                                l,
                                rate
                            )
                        end
                        msg
                    end
                end
            end
        end
    end
end


text:
using Plots, DelimitedFiles
using ProgressMeter
using Random, Distributions, StatsBase
Random.seed!(21);

mutable struct ParticleGroup
    locations;
    weights;
end

struct NoiseDists
    actual;
    estimated;
    perturb;
end

function elevAboveSeaLevel(x)
    sin_term = sin(0.01 * x + 3);
    cos_term = cos((0.01 * x - 2) * (0.005 * x - 5));
    prefix = (100000 * x) / (1 + (x + 500) ^ 2);
    elev = prefix * (cos_term - sin_term);
    return elev;
end

function get_measurement(x, distribution)
    msmt = elevAboveSeaLevel(x);
    msmt += rand(distribution);
    return msmt;
end

function get_likelihood(x, msmt, distribution)
    ytrue = elevAboveSeaLevel(x);
    residual = msmt - ytrue;
    prob = pdf(distribution, residual);
    return prob;
end

function init_particles(x0, x1, Np)
    scale = rand(Uniform(), Np);
    locations = x0 .+ (x1 - x0) .* scale;
    return ParticleGroup(locations, ones(Np));
end

function init_histories(Nt, Np)
    xtrue = zeros(Nt);
    x = zeros(Nt, Np);
    return (xtrue, x);
end

function update_weights!(particles, xtrue, noise)
    msmt = get_measurement(xtrue, noise.actual);
    particles.weights = get_likelihood.(
        particles.locations,
        msmt,
        noise.estimated
    );
end

function update_locations!(particles)
    wsample!(
        particles.locations,
        particles.weights,
        particles.locations,
        replace=true
    );
end

function move!(particles, noise, dt, velocity)
    particles.locations .+= dt * velocity;
    particles.locations .+= rand.(noise.perturb);
end

function forward!(particles, xtrue, noise, dt, velocity)
    update_weights!(particles, xtrue, noise);
    update_locations!(particles);
    move!(particles, noise, dt, velocity);
end

function run_filter(x0, x1, xtrue, Np, Nt, dt, velocity, σ_true, σ_est, σ_pert)
    particles = init_particles(x0, x1, Np);
    xtrue_history, x_history = init_histories(Nt, Np);
    noise = NoiseDists(map((σ) -> Normal(0, σ), [σ_true, σ_est, σ_pert])...);
    @showprogress 1 "Running pfilter..." for t = 1 : Nt
        forward!(particles, xtrue, noise, dt, velocity);
        xtrue += dt * velocity;
        xtrue_history[t] = xtrue;
        x_history[t, :] = particles.locations;
    end
    return (xtrue_history, x_history);
end

function make_marker(fpath, scale=5)
    verts = scale .* readdlm(fpath);
    x, y = verts[:, 1], verts[:, 2];
    return [Shape(x, y)];
end

function make_mountain(x0, x1)
    x = collect(x0 : 1 : x1);
    y = elevAboveSeaLevel.(x);
    return (x, y);
end

function make_animation(x, y, x0, y0, x1, y1, x_history, xtrue_history, aeroplane, Nt)
    p = Progress(Nt, 1, "Generating animation...")
    anim = @animate for t = 1 : Nt
        xest_t = x_history[t, :];
        yest_t = elevAboveSeaLevel.(xest_t);
        plot(x, y, xlims=(x0, x1), ylims=(y0, y1), legend=false);
        scatter!(xest_t, yest_t, markersize=1, color="darkblue");
        scatter!([xtrue_history[t] - 50], [y1 - 50], shape=aeroplane, ms=20);
        next!(p)
    end
    return anim;
end

function main()
    Nt, Np = 200, 5000;
    dt, velocity = 0.01, 200;
    x0, x1, xtrue = 0, 1000, 100;
    σ_true, σ_est, σ_pert  = 50, 40, 0.1 * dt * velocity;

    args = (x0, x1, xtrue, Np, Nt, dt, velocity, σ_true, σ_est, σ_pert)
    xtrue_history, x_history = run_filter(args...);

    x, y = make_mountain(x0, x1);
    y0, y1 = quantile(y[:], [0, 1]) .+ [-10, +75];
    aeroplane = make_marker("../data/airplane_verts.txt")

    anim = make_animation(x, y, x0, y0, x1, y1, x_history, xtrue_history, aeroplane, Nt);
    gif(anim, "../data/demo.gif", fps=20);
end

main();

max_stars_repo_name: jhultman/pfilter

max_issues_repo_name: jhultman/pfilter

max_forks_repo_name: jhultman/pfilter


text:
"""
    BernoulliCodeDistribution

Distribution where bits are each independent Bernoulli trials (coin flips)
"""
mutable struct BernoulliCodeDistribution <: AbstractBinaryVectorDistribution
    p::Vector{Float64}
    # I::Vector{Int} # moving this to metadata
    metadata::Dict{Any,Any}
    cache::Dict{Any,Any}

    function BernoulliCodeDistribution(P::Vector{Float64}; kwargs...)
        new(P, Dict(kwargs), Dict(:pdf=>spzeros(2^length(P))))
    end
end

function show(io::IO, P::BernoulliCodeDistribution)
    println(io,"Bernoulli Code Distribution")
    println(io,"N_neurons: $(n_bits(P))")
    println(io,"p_i's:     $(P.p)")
    show_metadata(io, P)
end

==(P1::BernoulliCodeDistribution, P2::BernoulliCodeDistribution) = P1.p == P2.p

n_bits(B::BernoulliCodeDistribution) = length(B.p)

function pdf(Pr::BernoulliCodeDistribution, x)
    if length(x) != n_bits(Pr)
        error("BernoulliCodeDistribution pdf: out of domain error")
    else
        # minor inefficiency: this won't differentiate between a 0.0 value that
        # is stored vs. not stored in the sparse vector. But, this would only
        # come up if one of the p_i = 0.
        # also, might be worth splitting this back into two cases (Vector{Bool}
        # vs. BitVector)
        # idx = 1 + dot([2^i for i = 0:(n_bits(Pr) - 1)], x)
        idx = 1 + _binary_to_int(x)
        if Pr.cache[:pdf][idx] == 0.0
            Pr.cache[:pdf][idx] = prod([x[i] ? Pr.p[i] : (1-Pr.p[i]) for i = 1:length(Pr.p)])
        end
        return Pr.cache[:pdf][idx]
        # if Pr.cache[:pdf][idx] > 0.0
        #     return Pr.cache[:pdf][idx]
        # else
        #     Pr.cache[:pdf][idx] = prod([x[i] ? Pr.p[i] : (1-Pr.p[i]) for i = 1:length(Pr.p)])
        #     return Pr.cache[:pdf][idx]
        # end
    end
    # length(x) == length(Pr.p) ? prod([x[i] ? Pr.p[i] : (1-Pr.p[i]) for i = 1:length(Pr.p)]) : error("BernoulliCodeDistribution: out of domain error")
end

function expectation_matrix(B::BernoulliCodeDistribution)
    em = B.p * B.p' # bits are independent coin flips
    em[1:(n_bits(B)+1):end] = B.p # set the diagonal to the probabilities instead of squares
    return em
end

# I'm going to be lazy here and allow this value to be computed each time. The
# parentheses around the second assignment expression are just for human
# readability.
entropy(P::BernoulliCodeDistribution) = (P.metadata[:entropy] = -sum([(p_i * log(p_i) + (1-p_i) * log(1-p_i)) for p_i in filter(x -> 0 < x < 1, P.p)]))
entropy2(P::BernoulliCodeDistribution) = (P.metadata[:entropy2] = -sum([(p_i * log2(p_i) + (1-p_i) * log2(1-p_i)) for p_i in filter(x -> 0 < x < 1, P.p)]))

random(B::BernoulliCodeDistribution, n_samples=1) = rand(n_bits(B), n_samples) .< B.p

max_stars_repo_name: sekunder/BinaryVectorProbability

max_issues_repo_name: sekunder/BinaryVectorProbability

max_forks_repo_name: sekunder/BinaryVectorProbability


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

function array_expr_inliner(node, cbdata, index, top_level, read)
    if isa(node, Expr)
        if node.head == :ref
            arr = cbdata[node.args[1]]
            idxs = node.args[2:end]
            map!((i) -> AstWalk(i, array_expr_inliner, cbdata), idxs)
            if length(idxs) == ndims(arr)
                flattened_idx = :($(idxs[end]) - 1)
                for i in length(idxs)-1:-1:1
                    flattened_idx = :(($(idxs[i]) - 1) + $(size(arr, i)) * ($flattened_idx))
                end
                node.args = [node.args[1], :(1 + $flattened_idx)]
            else
                # 1-d indexing
                @assert length(idxs) == 1 node
            end
        elseif node.head == :call && node.args[1] == :size
            arr = cbdata[node.args[2]]
            dim = node.args[3]
            return size(arr, dim)
        elseif node.head == :call && node.args[1] == :length
            if haskey(cbdata, node.args[2])
                arr = cbdata[node.args[2]]
                return length(arr)
            end
        elseif node.head in [:loophead, :parallel_loophead, :loopend, :parallel_loopend]
            return node
        end
    end
    ASTWALK_RECURSE
end

function inline_array_expressions(expr, env)
    AstWalk(expr, array_expr_inliner, env)
end


text:
module MyModule

export MyType, foo

struct MyType
    x
end

bar(x) = 2x

foo(a::MyType) = bar(a.x) + 1

show(io::IO, a::MyType) = print(io, "MyType $(a.x)")

end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


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: Estefane Horn, N�bia do Carmo Ferreira
#*$Id$
#*-------------------------------------------------------------------
#-------------------------------------------------------------------
#* Model of a valve (simplified)
#*--------------------------------------------------------------------
#*
#* Author: Paula B. Staudt
#*--------------------------------------------------------------------
type valve
	valve()=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",
				:Lower=>1
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Valid=>["linear", "parabolic", "equal", "quick", "hyperbolic"],
				:Default=>"linear"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Valid Phases for Flash Calculation",
				:Valid=>["Vapour-Only", "Liquid-Only"],
				:Default=>"Liquid-Only"
			)),
			dens_mass(Dict{Symbol,Any}(
				:Brief=>"Water Mass Density at 60 F",
				:Hidden=>true
			)),
			flow_mass(Dict{Symbol,Any}(
				:DisplayUnit=>"kg/s"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Pressure Drop",
				:DisplayUnit=>"kPa",
				:Symbol=>"\\Delta P"
			)),
			flow_vol(Dict{Symbol,Any}(
				:Brief=>"Volumetric Flow"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Opening Function",
				:Hidden=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Valve Flow Coefficient",
				:Unit=>"gal/min/psi^0.5"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Valve Flow Coefficient",
				:Unit=>"m^2"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Valve Gas Flow Coefficient",
				:Unit=>"ft^3/h/psi"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Liquid-gas Coefficient Ratio",
				:Unit=>"(ft^3/gal)*(min/h)/(psi^.5)"
			)),
			fraction(Dict{Symbol,Any}(
				:Brief=>"Actual valve stem position"
			)),
			DanaReal(),
			vol_mol(Dict{Symbol,Any}(
				:Brief=>"Mixture Molar Volume"
			)),
			dens_mass(Dict{Symbol,Any}(
				:Brief=>"Mixture Mass Density"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet stream",
				:PosX=>0,
				:PosY=>0.7365,
				:Symbol=>"_{in}"
			)),
			streamPH(Dict{Symbol,Any}(
				:Brief=>"Outlet stream",
				:PosX=>1,
				:PosY=>0.7365,
				:Symbol=>"_{out}"
			)),
			[
				:(Outlet.P = Inlet.P - Pdrop),
				:(Outlet.h = Inlet.h),
				:(Outlet.F = Inlet.F),
				:(Outlet.z = Inlet.z),
				:(Cv1=Cv*"1/(gal/min/psi^0.5)"*2.3837e-5*"m^2"),
				:(Fvol = fc*(Cv/sqrt(1/rho60F))*sqrt(Pdrop/rho)),
				:(C*Cv=Cg),
				:(a=1/(1.6764e-2*C*"1/((ft^3/gal)*(min/h)/(psi^.5))")*sqrt(Pdrop/Inlet.P)),
				:(Fvol = fc*(Cv/sqrt(1/rho60F))*sqrt(Pdrop/rho)),
				:(C*Cv=Cg),
				:(a=1/(1.6764e-2*C*"1/((ft^3/gal)*(min/h)/(psi^.5))")*sqrt(Pdrop/Inlet.P)),
				:(rho = PP.LiquidDensity(Inlet.T,Inlet.P,Inlet.z)),
				:(vm = PP.LiquidVolume(Inlet.T,Inlet.P,Inlet.z)),
				:(C*Cv=Cg),
				:(a=1/(1.6764e-2*C*"1/((ft^3/gal)*(min/h)/(psi^.5))")*sqrt(Pdrop/Inlet.P)),
				:(Fvol = fc*0.13446*"psi^.5"*Cg*sqrt(Inlet.P/1000*rho60F/rho)*sin(a*"rad")),
				:(Fvol = fc*0.13446*Cv*sqrt(Inlet.P*rho60F/rho)),
				:(Fvol = fc*(Cv/sqrt(1/rho60F))*sqrt(Pdrop/rho)),
				:(C*Cv=Cg),
				:(a=1/(1.6764e-2*C*"1/((ft^3/gal)*(min/h)/(psi^.5))")*sqrt(Pdrop/Inlet.P)),
				:(rho = PP.VapourDensity(Inlet.T,Inlet.P,Inlet.z)),
				:(vm = PP.VapourVolume(Inlet.T,Inlet.P,Inlet.z)),
				:(Fvol = Inlet.F*vm),
				:(W=Fvol*rho),
				:(fc = StemPosition),
				:(fc = StemPosition^2),
				:(fc = StemPosition^2/(2-StemPosition^4)^(1/2)),
				:(fc = 10*StemPosition/sqrt(1+99*StemPosition^2)),
				:(fc = 0.1*StemPosition/sqrt(1-0.99*StemPosition^2)),
			],
			[
				"Pressure Drop","Enthalpy Balance","Molar Balance","Outlet Composition","","Valve Equation - Liquid Flow","Liquid-gas Coefficient Ratio","","Valve Equation - Liquid Flow","Liquid-gas Coefficient Ratio","","Liquid Mass Density","Liquid Molar Volume","Liquid-gas Coefficient Ratio","","Valve Equation - Vapour Flow","Valve Equation - Vapour Flow","Valve Equation - Vapour Flow","Liquid-gas Coefficient Ratio","","Vapour Mass Density","Vapour Molar Volume","Calculate Mass Flow","","Opening Equation","Opening Equation","Opening Equation","Opening Equation","Opening Equation",
			],
			[:PP,:NComp,:ValveType,:ValidPhases,:rho60F,],
			[:W,:Pdrop,:Fvol,:fc,:Cv,:Cv1,:Cg,:C,:StemPosition,:a,:vm,:rho,:Inlet,:Outlet,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	ValveType::DanaSwitcher
	ValidPhases::DanaSwitcher
	rho60F::dens_mass
	W::flow_mass
	Pdrop::press_delta
	Fvol::flow_vol
	fc::positive
	Cv::DanaReal
	Cv1::DanaReal
	Cg::DanaReal
	C::DanaReal
	StemPosition::fraction
	a::DanaReal
	vm::vol_mol
	rho::dens_mass
	Inlet::stream
	Outlet::streamPH
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export valve
function set(in::valve)
	rho60F = 984.252 * "kg/m^3"

end
function setEquationFlow(in::valve)
	#"First order valve dynamics"
	#	Tau*diff(StemPosition) = vsp-StemPosition;
	#"Flow Signal"
	#	vsp = vsignal;
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	let switch=ValidPhases
		#############################################################
		if switch=="Liquid-Only"
			if Pdrop > 0 * "atm"
				addEquation(6)
				addEquation(7)
				addEquation(8)
			else
				addEquation(9)
				addEquation(10)
				addEquation(11)
			end
			addEquation(12)
			addEquation(13)
			############################################################

		elseif switch=="Vapour-Only"
			if Pdrop > 0 * "atm"
				#Update for gas flow !!!!
				addEquation(14)
				addEquation(15)
				if 1.5708 > a
					#Fvol = fc*Cg*sqrt(Inlet.P/1000*rho60F/rho);####rho60f/rho ok!!!!
					#Fvol = fc*Cv*sqrt(Pdrop/1000*rho60F/rho);
					#W = fc*Cv1*sqrt(Pdrop/1000*rho);
					addEquation(16)
				else
					addEquation(17)
				end
			else
				addEquation(18)
				addEquation(19)
				addEquation(20)
			end
			addEquation(21)
			#rho=3.708741*'kg/m^3';
			addEquation(22)
		end
	end
	######################################################
	addEquation(23)
	addEquation(24)
	let switch=ValveType
		#Update the valve Type
		if switch=="linear"
			addEquation(25)
		elseif switch=="parabolic"
			addEquation(26)
		elseif switch=="equal"
			addEquation(27)
		elseif switch=="quick"
			addEquation(28)
		elseif switch=="hyperbolic"
			addEquation(29)
		end
	end
end
function atributes(in::valve,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/Valve"
	fields[:Brief]="Model of a valve."
	fields[:Info]="== Model of valves ==
* Linear;
* Parabolic;
* Equal;
* Quick;
* Hyperbolic.

== Assumptions ==
* First Order Dynamic;
* Only Liquid or Only Vapour;
* Isentalpic.

== Specify ==
* the valve type;
* the Valve Coefficient (Cv);
* the valve time constant (Tau).
"
	drive!(fields,_)
	return fields
end
valve(_::Dict{Symbol,Any})=begin
	newModel=valve()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
mutable struct RecordWithoutSubrecords{T,C} <: Record
   io       :: IO    # underlying I/O stream
   reclen   :: T     # length of this record
   nleft    :: T     # bytes left in this record
   writable :: Bool  # whether this record is writable
   convert  :: C     # convert method
end

function Record( f::FortranFile{SequentialAccess{WithoutSubrecords{T}},C} ) where {T,C}
## constructor for readable records
   conv = f.convert
   reclen = conv.onread( read(f.io, T) ) # read leading record marker
   RecordWithoutSubrecords{T,C}(f.io, reclen, reclen, false, conv)
end

function Record( f::FortranFile{SequentialAccess{WithoutSubrecords{T}},C}, towrite::Integer ) where {T,C}
## constructor for writable records
   conv = f.convert
   write(f.io, conv.onwrite( convert(T, towrite) )) # write leading record marker
   RecordWithoutSubrecords{T,C}(f.io, towrite, towrite, true, conv)
end

function Base.unsafe_read( rec::RecordWithoutSubrecords, p::Ptr{UInt8}, n::UInt )
   if (n > rec.nleft)
      fthrow("attempting to read beyond record end")
   end
   unsafe_read( rec.io, p, n )
   rec.nleft -= n
   nothing
end

function Base.unsafe_write( rec::RecordWithoutSubrecords, p::Ptr{UInt8}, n::UInt )
   if (n > rec.nleft); fthrow("attempting to write beyond record end"); end
   nwritten = unsafe_write( rec.io, p, n )
   rec.nleft -= n
   return nwritten
end

function Base.close( rec::RecordWithoutSubrecords{T} ) where {T}
   if rec.writable
      @assert rec.nleft == 0
      write(rec.io, rec.convert.onwrite( convert(T, rec.reclen)) ) # write trailing record marker
   else
      skip(rec.io, rec.nleft)
      reclen = rec.convert.onread( read(rec.io, T) ) # read trailing record marker
      if reclen != rec.reclen
         fthrow("trailing record marker doesn't match")
      end
   end
   nothing
end


text:
import Pkg
Pkg.activate(".")
using ProjectMoana
using ArgParse
import ColorTypes:RGB
import Images:save

function parse_commandline()
    s = ArgParseSettings(description = "This program make an averege image from 4 images. Try me!",
                               usage = "usage: [--help] [--file_in FILE_IN]",
                              epilog = "Let's try again!")

    @add_arg_table s begin
        "--file_in"
            help = "input PFM file name"
            required = true
    end

    return parse_args(s)
end

function main()
    params = parse_commandline()

    fileinput=params["file_in"]

    img0 = read_pfm_image(string(fileinput, "000.pfm"))
    img1 = read_pfm_image(string(fileinput, "001.pfm"))
    img2 = read_pfm_image(string(fileinput, "002.pfm"))
    img3 = read_pfm_image(string(fileinput, "003.pfm"))

    out=HdrImage(img0.width,img0.height)

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

            color_r = (img0.pixels[ get_pixel(img0, x, y)].r + img1.pixels[ get_pixel(img1, x, y)].r + img2.pixels[ get_pixel(img2, x, y)].r + img3.pixels[ get_pixel(img3, x, y)].r) / 4.
            color_g = (img0.pixels[ get_pixel(img0, x, y)].g + img1.pixels[ get_pixel(img1, x, y)].g + img2.pixels[ get_pixel(img2, x, y)].g + img3.pixels[ get_pixel(img3, x, y)].g) / 4.
            color_b = (img0.pixels[ get_pixel(img0, x, y)].b + img1.pixels[ get_pixel(img1, x, y)].b + img2.pixels[ get_pixel(img2, x, y)].b + img3.pixels[ get_pixel(img3, x, y)].b) / 4.
            color = RGB(color_r,color_g,color_b)
            set_pixel(out, x, y, color)

        end
    end

    write(string(fileinput, ".pfm"), out)
    println("$(string(fileinput, ".pfm")) has been written to disk.")
end

main()


text:
tests = [
    "types",
	"utils",
	"sample",
	"fit",
	"discrete",
	"univariate",
	"truncate",
	"multinomial",
	"dirichlet",
	"mvnormal",
	"conjugates",
	"kolmogorov",
	"edgeworth",
	"matrix"]

println("Running tests:")

for t in tests
	test_fn = joinpath("test", "$t.jl")
    println(" * $test_fn")
    include(test_fn)
end


text:
struct RiskPriors{T} <: AbstractRisk{T}
  sparks::Union{Nothing, Vector{UnivariateDistribution}}
  susceptibility::Union{Nothing, Vector{UnivariateDistribution}}
  infectivity::Union{Nothing, Vector{UnivariateDistribution}}
  transmissibility::Union{Nothing, Vector{UnivariateDistribution}}
  latency::Union{Nothing, Vector{UnivariateDistribution}}
  removal::Union{Nothing, Vector{UnivariateDistribution}}

  RiskPriors{SEIR}(pϵ, ps, pκ, pt, pl, pr) = new{SEIR}(pϵ, ps, pκ, pt, pl, pr)
  RiskPriors{SEI}(pϵ, ps, pκ, pt, pl) = new{SEI}(pϵ, ps, pκ, pt, pl, nothing)
  RiskPriors{SIR}(pϵ, ps, pκ, pt, pr) = new{SIR}(pϵ, ps, pκ, pt, nothing, pr)
  RiskPriors{SI}(pϵ, ps, pκ, pt) = new{SI}(pϵ, ps, pκ, pt, nothing, nothing)
end

function Base.show(io::IO, x::RiskPriors{T}) where T <: DiseaseStateSequence
  return print(io, "$T model risk function priors")
end